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-2014 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__all__ = ['SparseSolver', 'SMatrix', 'GreensFunction'] 

10 

11from collections import namedtuple 

12from itertools import product 

13import abc 

14from numbers import Integral 

15import numpy as np 

16import scipy.sparse as sp 

17from .._common import ensure_isinstance, deprecate_args 

18from .. import system, physics 

19from functools import reduce 

20 

21 

22LinearSys = namedtuple('LinearSys', ['lhs', 'rhs', 'indices', 'num_orb']) 

23 

24 

25class SparseSolver(metaclass=abc.ABCMeta): 

26 """Solver class for computing physical quantities based on solving 

27 a liner system of equations. 

28 

29 `SparseSolver` is an abstract base class. It cannot be instantiated as it 

30 does not specify the actual linear solve step. In order to be directly 

31 usable, a derived class must implement the methods 

32 

33 - `_factorize` and `_solve_linear_sys`, that solve a linear system of 

34 equations 

35 

36 and the following properties: 

37 

38 - `lhsformat`, `rhsformat`: Sparse matrix format of left and right hand 

39 sides of the linear system, respectively. Must be one of {'coo', 'csc', 

40 'csr'}. 

41 

42 - `nrhs`: Number of right hand side vectors that should be solved at once 

43 in a call to `solve_linear_sys`, when the full solution is computed (i.e. 

44 kept_vars covers all entries in the solution). This should be not too big 

45 too avoid excessive memory usage, but for some solvers not too small for 

46 performance reasons. 

47 """ 

48 

49 @abc.abstractmethod 

50 def _factorized(self, a): 

51 """ 

52 Return a preprocessed version of a matrix for the use with 

53 `solve_linear_sys`. 

54 

55 Parameters 

56 ---------- 

57 a : a scipy.sparse.coo_matrix sparse matrix. 

58 

59 Returns 

60 ------- 

61 factorized_a : object 

62 factorized lhs to be used with `_solve_linear_sys`. 

63 """ 

64 pass 

65 

66 @abc.abstractmethod 

67 def _solve_linear_sys(self, factorized_a, b, kept_vars): 

68 """ 

69 Solve the linar system `a x = b`, returning the part of 

70 the result indicated in kept_vars. 

71 

72 Parameters 

73 ---------- 

74 factorized : object 

75 The result of calling `_factorized` for the matrix a. 

76 b : sparse matrix 

77 The right hand side. Its format must match `rhsformat`. 

78 kept_vars : slice object or sequence of integers 

79 A sequence of numbers of variables to keep in the solution 

80 

81 Returns 

82 ------- 

83 output : NumPy matrix 

84 Solution to the system of equations. 

85 """ 

86 pass 

87 

88 def _make_linear_sys(self, sys, in_leads, energy=0, args=(), 

89 check_hermiticity=True, realspace=False, 

90 *, params=None): 

91 """Make a sparse linear system of equations defining a scattering 

92 problem. 

93 

94 Parameters 

95 ---------- 

96 sys : kwant.system.FiniteSystem 

97 Low level system, containing the leads and the Hamiltonian of a 

98 scattering region. 

99 in_leads : sequence of integers 

100 Numbers of leads in which current or wave function is injected. 

101 energy : number 

102 Excitation energy at which to solve the scattering problem. 

103 args : tuple, defaults to empty 

104 Positional arguments to pass to the ``hamiltonian`` method. 

105 Deprecated in favor of 'params' (and mutually exclusive with it). 

106 check_hermiticity : bool 

107 Check if Hamiltonian matrices are in fact Hermitian. 

108 Enables deduction of missing transmission coefficients. 

109 realspace : bool 

110 Calculate Green's function between the outermost lead 

111 sites, instead of lead modes. This is almost always 

112 more computationally expensive and less stable. 

113 params : dict, optional 

114 Dictionary of parameter names and their values. Mutually exclusive 

115 with 'args'. 

116 

117 Returns 

118 ------- 

119 (lhs, rhs, indices, num_orb) : LinearSys 

120 `lhs` is a scipy.sparse.csc_matrix, containing the left hand side 

121 of the system of equations. `rhs` is a list of matrices with the 

122 right hand side, with each matrix corresponding to one lead 

123 mentioned in `in_leads`. `indices` is a list of arrays of variables 

124 in the system of equations corresponding to the the outgoing modes 

125 in each lead, or the indices of variables, on which a lead defined 

126 via self-energy adds the self-energy. Finally `num_orb` is the 

127 total number of degrees of freedom in the scattering region. 

128 

129 lead_info : list of objects 

130 Contains one entry for each lead. If `realspace=False`, this is an 

131 instance of `~kwant.physics.PropagatingModes` for all leads that 

132 have modes, otherwise the lead self-energy matrix. 

133 

134 Notes 

135 ----- 

136 All the leads should implement a method `modes` if `realspace=False` 

137 and a method `selfenergy`. 

138 

139 The system of equations that is created will be described in detail 

140 elsewhere. 

141 """ 

142 

143 syst = sys # ensure consistent naming across function bodies 

144 ensure_isinstance(syst, system.System) 

145 

146 splhsmat = getattr(sp, self.lhsformat + '_matrix') 

147 sprhsmat = getattr(sp, self.rhsformat + '_matrix') 

148 

149 if not syst.lead_interfaces: 149 ↛ 150line 149 didn't jump to line 150, because the condition on line 149 was never true

150 raise ValueError('System contains no leads.') 

151 lhs, norb = syst.hamiltonian_submatrix(args, sparse=True, 

152 return_norb=True, 

153 params=params)[:2] 

154 lhs = getattr(lhs, 'to' + self.lhsformat)() 

155 lhs = lhs - energy * sp.identity(lhs.shape[0], format=self.lhsformat) 

156 num_orb = lhs.shape[0] 

157 

158 if check_hermiticity and len(lhs.data): 158 ↛ 167line 158 didn't jump to line 167, because the condition on line 158 was never false

159 rtol = 1e-13 

160 atol = 1e-300 

161 tol = rtol * np.max(np.abs(lhs.data)) + atol 

162 if np.any(np.abs((lhs - lhs.T.conj()).data) > tol): 162 ↛ 163line 162 didn't jump to line 163, because the condition on line 162 was never true

163 raise ValueError('System Hamiltonian is not Hermitian. ' 

164 'Use option `check_hermiticity=False` ' 

165 'if this is intentional.') 

166 

167 offsets = np.empty(norb.shape[0] + 1, int) 

168 offsets[0] = 0 

169 offsets[1 :] = np.cumsum(norb) 

170 

171 # Process the leads, generate the eigenvector matrices and lambda 

172 # vectors. Then create blocks of the linear system and add them 

173 # step by step. 

174 indices = [] 

175 rhs = [] 

176 lead_info = [] 

177 for leadnum, interface in enumerate(syst.lead_interfaces): 

178 lead = syst.leads[leadnum] 

179 if not realspace: 

180 if system.is_selfenergy_lead(lead): 

181 # We make equivalent checks to this in 'smatrix', 

182 # 'greens_function' and 'wave_function' (where we 

183 # can give more informative error messages). 

184 # Putting this check here again is just a sanity check, 

185 assert leadnum not in in_leads 

186 sigma = np.asarray(lead.selfenergy(energy, args, params=params)) 

187 rhs.append(None) 

188 lead_info.append(sigma) 

189 indices.append(None) 

190 # add selfenergy to the LHS 

191 coords = np.r_[tuple(slice(offsets[i], offsets[i + 1]) 

192 for i in interface)] 

193 if sigma.shape != 2 * coords.shape: 193 ↛ 194line 193 didn't jump to line 194, because the condition on line 193 was never true

194 raise ValueError( 

195 f"Self-energy dimension for lead {leadnum} " 

196 "does not match the total number of orbitals " 

197 "of the sites for which it is defined." 

198 ) 

199 y, x = np.meshgrid(coords, coords) 

200 sig_sparse = splhsmat((sigma.flat, [x.flat, y.flat]), 

201 lhs.shape) 

202 lhs = lhs + sig_sparse 

203 else: 

204 prop, stab = lead.modes(energy, args, params=params) 

205 lead_info.append(prop) 

206 u = stab.vecs 

207 ulinv = stab.vecslmbdainv 

208 nprop = stab.nmodes 

209 svd_v = stab.sqrt_hop 

210 

211 if len(u) == 0: 

212 rhs.append(None) 

213 continue 

214 

215 indices.append(np.arange(lhs.shape[0], lhs.shape[0] + nprop)) 

216 

217 u_out, ulinv_out = u[:, nprop:], ulinv[:, nprop:] 

218 u_in, ulinv_in = u[:, :nprop], ulinv[:, :nprop] 

219 

220 # Construct a matrix of 1's that translates the 

221 # inter-cell hopping to a proper hopping 

222 # from the system to the lead. 

223 iface_orbs = np.r_[tuple(slice(offsets[i], offsets[i + 1]) 

224 for i in interface)] 

225 

226 n_lead_orbs = svd_v.shape[0] 

227 if n_lead_orbs != len(iface_orbs): 227 ↛ 228line 227 didn't jump to line 228, because the condition on line 227 was never true

228 msg = ('Lead {0} has hopping with dimensions ' 

229 'incompatible with its interface dimension.') 

230 raise ValueError(msg.format(leadnum)) 

231 

232 coords = np.r_[[np.arange(len(iface_orbs))], [iface_orbs]] 

233 transf = sp.csc_matrix((np.ones(len(iface_orbs)), coords), 

234 shape=(iface_orbs.size, lhs.shape[0])) 

235 

236 v_sp = sp.csc_matrix(svd_v.T.conj()) * transf 

237 vdaguout_sp = (transf.T * 

238 sp.csc_matrix(np.dot(svd_v, u_out))) 

239 lead_mat = - ulinv_out 

240 

241 lhs = sp.bmat([[lhs, vdaguout_sp], [v_sp, lead_mat]], 

242 format=self.lhsformat) 

243 

244 if leadnum in in_leads and nprop > 0: 

245 vdaguin_sp = transf.T * sp.csc_matrix( 

246 -np.dot(svd_v, u_in)) 

247 

248 # defer formation of the real matrix until the proper 

249 # system size is known 

250 rhs.append((vdaguin_sp, ulinv_in)) 

251 else: 

252 rhs.append(None) 

253 else: 

254 sigma = np.asarray(lead.selfenergy(energy, args, params=params)) 

255 lead_info.append(sigma) 

256 coords = np.r_[tuple(slice(offsets[i], offsets[i + 1]) 

257 for i in interface)] 

258 

259 if sigma.shape != 2 * coords.shape: 259 ↛ 260line 259 didn't jump to line 260, because the condition on line 259 was never true

260 msg = ('Self-energy dimension for lead {0} does not ' 

261 'match the total number of orbitals of the ' 

262 'sites for which it is defined.') 

263 raise ValueError(msg.format(leadnum)) 

264 

265 y, x = np.meshgrid(coords, coords) 

266 sig_sparse = splhsmat((sigma.flat, [x.flat, y.flat]), 

267 lhs.shape) 

268 lhs = lhs + sig_sparse # __iadd__ is not implemented in v0.7 

269 indices.append(coords) 

270 if leadnum in in_leads: 

271 # defer formation of true rhs until the proper system 

272 # size is known 

273 rhs.append((coords,)) 

274 

275 # Resize the right-hand sides to be compatible with the full lhs 

276 for i, mats in enumerate(rhs): 

277 if isinstance(mats, tuple): 

278 if len(mats) == 1: 

279 # self-energy lead 

280 l = mats[0].shape[0] 

281 rhs[i] = sprhsmat((-np.ones(l), [mats[0], np.arange(l)]), 

282 shape=(lhs.shape[0], l)) 

283 else: 

284 # lead with modes 

285 zero_rows = (lhs.shape[0] - mats[0].shape[0] - 

286 mats[1].shape[0]) 

287 

288 zero_mat = sprhsmat((zero_rows, mats[0].shape[1])) 

289 bmat = [[mats[0]], [mats[1]], [zero_mat]] 

290 

291 rhs[i] = sp.bmat(bmat, format=self.rhsformat) 

292 elif mats is None: 292 ↛ 296line 292 didn't jump to line 296, because the condition on line 292 was never false

293 # A lead with no rhs. 

294 rhs[i] = np.zeros((lhs.shape[0], 0)) 

295 else: 

296 raise RuntimeError('Unknown right-hand side format') 

297 

298 return LinearSys(lhs, rhs, indices, num_orb), lead_info 

299 

300 @deprecate_args 

301 def smatrix(self, sys, energy=0, args=(), 

302 out_leads=None, in_leads=None, check_hermiticity=True, 

303 *, params=None): 

304 """ 

305 Compute the scattering matrix of a system. 

306 

307 An alias exists for this common name: ``kwant.smatrix``. 

308 

309 Parameters 

310 ---------- 

311 sys : `kwant.system.FiniteSystem` 

312 Low level system, containing the leads and the Hamiltonian of a 

313 scattering region. 

314 energy : number 

315 Excitation energy at which to solve the scattering problem. 

316 args : tuple, defaults to empty 

317 Positional arguments to pass to the ``hamiltonian`` method. 

318 Deprecated in favor of 'params' (and mutually exclusive with it). 

319 out_leads : sequence of integers or ``None`` 

320 Numbers of leads where current or wave function is extracted. None 

321 is interpreted as all leads. Default is ``None`` and means "all 

322 leads". 

323 in_leads : sequence of integers or ``None`` 

324 Numbers of leads in which current or wave function is injected. 

325 None is interpreted as all leads. Default is ``None`` and means 

326 "all leads". 

327 check_hermiticity : ``bool`` 

328 Check if the Hamiltonian matrices are Hermitian. 

329 Enables deduction of missing transmission coefficients. 

330 params : dict, optional 

331 Dictionary of parameter names and their values. Mutually exclusive 

332 with 'args'. 

333 

334 Returns 

335 ------- 

336 output : `~kwant.solvers.common.SMatrix` 

337 See the notes below and `~kwant.solvers.common.SMatrix` 

338 documentation. 

339 

340 Notes 

341 ----- 

342 This function can be used to calculate the conductance and other 

343 transport properties of a system. See the documentation for its output 

344 type, `~kwant.solvers.common.SMatrix`. 

345 

346 The returned object contains the scattering matrix elements from the 

347 `in_leads` to the `out_leads` as well as information about the lead 

348 modes. 

349 

350 Both `in_leads` and `out_leads` must be sorted and may only contain 

351 unique entries. 

352 """ 

353 

354 syst = sys # ensure consistent naming across function bodies 

355 ensure_isinstance(syst, system.System) 

356 

357 n = len(syst.lead_interfaces) 

358 if in_leads is None: 

359 in_leads = list(range(n)) 

360 else: 

361 in_leads = list(in_leads) 

362 if out_leads is None: 

363 out_leads = list(range(n)) 

364 else: 

365 out_leads = list(out_leads) 

366 if (np.any(np.diff(in_leads) <= 0) or np.any(np.diff(out_leads) <= 0)): 366 ↛ 367line 366 didn't jump to line 367, because the condition on line 366 was never true

367 raise ValueError("Lead lists must be sorted and " 

368 "with unique entries.") 

369 if len(in_leads) == 0 or len(out_leads) == 0: 

370 raise ValueError("No output is requested.") 

371 

372 for direction, leads in [("in", in_leads), ("out", out_leads)]: 

373 for lead in leads: 

374 if system.is_selfenergy_lead(syst.leads[lead]): 

375 raise ValueError( 

376 f"lead {lead} is only defined by a self-energy, " 

377 "so we cannot calculate scattering matrix elements for it. " 

378 f"Specify '{direction}_leads' without '{lead}'." 

379 ) 

380 

381 linsys, lead_info = self._make_linear_sys(syst, in_leads, energy, args, 

382 check_hermiticity, False, 

383 params=params) 

384 

385 kept_vars = np.concatenate([coords for i, coords in 

386 enumerate(linsys.indices) if i in 

387 out_leads]) 

388 

389 # Do not perform factorization if no calculation is to be done. 

390 len_rhs = sum(i.shape[1] for i in linsys.rhs) 

391 len_kv = len(kept_vars) 

392 if not(len_rhs and len_kv): 

393 return SMatrix(np.zeros((len_kv, len_rhs)), lead_info, 

394 out_leads, in_leads, check_hermiticity) 

395 

396 rhs = sp.bmat([linsys.rhs], format=self.rhsformat) 

397 flhs = self._factorized(linsys.lhs) 

398 data = self._solve_linear_sys(flhs, rhs, kept_vars) 

399 

400 return SMatrix(data, lead_info, out_leads, in_leads, check_hermiticity) 

401 

402 @deprecate_args 

403 def greens_function(self, sys, energy=0, args=(), 

404 out_leads=None, in_leads=None, check_hermiticity=True, 

405 *, params=None): 

406 """ 

407 Compute the retarded Green's function of the system between its leads. 

408 

409 An alias exists for this common name: ``kwant.greens_function``. 

410 

411 Parameters 

412 ---------- 

413 sys : `kwant.system.FiniteSystem` 

414 Low level system, containing the leads and the Hamiltonian of a 

415 scattering region. 

416 energy : number 

417 Excitation energy at which to solve the scattering problem. 

418 args : tuple, defaults to empty 

419 Positional arguments to pass to the ``hamiltonian`` method. 

420 Deprecated in favor of 'params' (and mutually exclusive with it). 

421 out_leads : sequence of integers or ``None`` 

422 Numbers of leads where current or wave function is extracted. None 

423 is interpreted as all leads. Default is ``None`` and means "all 

424 leads". 

425 in_leads : sequence of integers or ``None`` 

426 Numbers of leads in which current or wave function is injected. 

427 None is interpreted as all leads. Default is ``None`` and means 

428 "all leads". 

429 check_hermiticity : ``bool`` 

430 Check if the Hamiltonian matrices are Hermitian. 

431 Enables deduction of missing transmission coefficients. 

432 params : dict, optional 

433 Dictionary of parameter names and their values. Mutually exclusive 

434 with 'args'. 

435 

436 Returns 

437 ------- 

438 output : `~kwant.solvers.common.GreensFunction` 

439 See the notes below and `~kwant.solvers.common.GreensFunction` 

440 documentation. 

441 

442 Notes 

443 ----- 

444 This function can be used to calculate the conductance and other 

445 transport properties of a system. It is often slower and less stable 

446 than the scattering matrix-based calculation executed by 

447 `~kwant.solvers.common.smatrix`, and is currently provided mostly for testing 

448 purposes and compatibility with RGF code. 

449 

450 It returns an object encapsulating the Green's function elements 

451 between the system sites interfacing the leads in `in_leads` and those 

452 interfacing the leads in `out_leads`. The returned object also 

453 contains a list with self-energies of the leads. 

454 

455 Both `in_leads` and `out_leads` must be sorted and may only contain 

456 unique entries. 

457 """ 

458 

459 syst = sys # ensure consistent naming across function bodies 

460 ensure_isinstance(syst, system.System) 

461 

462 n = len(syst.lead_interfaces) 

463 if in_leads is None: 

464 in_leads = list(range(n)) 

465 else: 

466 in_leads = list(in_leads) 

467 if out_leads is None: 

468 out_leads = list(range(n)) 

469 else: 

470 out_leads = list(out_leads) 

471 if (np.any(np.diff(in_leads) <= 0) or np.any(np.diff(out_leads) <= 0)): 471 ↛ 472line 471 didn't jump to line 472, because the condition on line 471 was never true

472 raise ValueError("Lead lists must be sorted and " 

473 "with unique entries.") 

474 if len(in_leads) == 0 or len(out_leads) == 0: 474 ↛ 475line 474 didn't jump to line 475, because the condition on line 474 was never true

475 raise ValueError("No output is requested.") 

476 

477 linsys, lead_info = self._make_linear_sys(syst, in_leads, energy, args, 

478 check_hermiticity, True, 

479 params=params) 

480 

481 kept_vars = np.concatenate([coords for i, coords in 

482 enumerate(linsys.indices) if i in 

483 out_leads]) 

484 

485 # Do not perform factorization if no calculation is to be done. 

486 len_rhs = sum(i.shape[1] for i in linsys.rhs) 

487 len_kv = len(kept_vars) 

488 if not(len_rhs and len_kv): 488 ↛ 489line 488 didn't jump to line 489, because the condition on line 488 was never true

489 return GreensFunction(np.zeros((len_kv, len_rhs)), lead_info, 

490 out_leads, in_leads, check_hermiticity) 

491 

492 rhs = sp.bmat([linsys.rhs], format=self.rhsformat) 

493 flhs = self._factorized(linsys.lhs) 

494 data = self._solve_linear_sys(flhs, rhs, kept_vars) 

495 

496 return GreensFunction(data, lead_info, out_leads, in_leads, 

497 check_hermiticity) 

498 

499 @deprecate_args 

500 def ldos(self, sys, energy=0, args=(), check_hermiticity=True, 

501 *, params=None): 

502 """ 

503 Calculate the local density of states of a system at a given energy. 

504 

505 An alias exists for this common name: ``kwant.ldos``. 

506 

507 Parameters 

508 ---------- 

509 sys : `kwant.system.FiniteSystem` 

510 Low level system, containing the leads and the Hamiltonian of the 

511 scattering region. 

512 energy : number 

513 Excitation energy at which to solve the scattering problem. 

514 args : tuple of arguments, or empty tuple 

515 Positional arguments to pass to the function(s) which 

516 evaluate the hamiltonian matrix elements. 

517 Deprecated in favor of 'params' (and mutually exclusive with it). 

518 check_hermiticity : ``bool`` 

519 Check if the Hamiltonian matrices are Hermitian. 

520 params : dict, optional 

521 Dictionary of parameter names and their values. Mutually exclusive 

522 with 'args'. 

523 

524 Returns 

525 ------- 

526 ldos : a NumPy array 

527 Local density of states at each orbital of the system. 

528 """ 

529 

530 syst = sys # ensure consistent naming across function bodies 

531 ensure_isinstance(syst, system.System) 

532 if not check_hermiticity: 532 ↛ 533line 532 didn't jump to line 533, because the condition on line 532 was never true

533 raise NotImplementedError("ldos for non-Hermitian Hamiltonians " 

534 "is not implemented yet.") 

535 

536 for lead in syst.leads: 

537 if system.is_selfenergy_lead(lead): 

538 # TODO: fix this 

539 raise NotImplementedError("ldos for leads with only " 

540 "self-energy is not implemented yet.") 

541 

542 linsys = self._make_linear_sys(syst, range(len(syst.leads)), energy, 

543 args, check_hermiticity, 

544 params=params)[0] 

545 

546 ldos = np.zeros(linsys.num_orb, float) 

547 

548 # Do not perform factorization if no further calculation is needed. 

549 if not sum(i.shape[1] for i in linsys.rhs): 

550 return ldos 

551 

552 factored = self._factorized(linsys.lhs) 

553 

554 rhs = sp.bmat([linsys.rhs], format=self.rhsformat) 

555 for j in range(0, rhs.shape[1], self.nrhs): 

556 jend = min(j + self.nrhs, rhs.shape[1]) 

557 psi = self._solve_linear_sys(factored, rhs[:, j:jend], 

558 slice(linsys.num_orb)) 

559 ldos += np.sum(np.square(abs(psi)), axis=1) 

560 

561 return ldos * (0.5 / np.pi) 

562 

563 @deprecate_args 

564 def wave_function(self, sys, energy=0, args=(), check_hermiticity=True, 

565 *, params=None): 

566 r""" 

567 Return a callable object for the computation of the wave function 

568 inside the scattering region. 

569 

570 An alias exists for this common name: ``kwant.wave_function``. 

571 

572 Parameters 

573 ---------- 

574 sys : `kwant.system.FiniteSystem` 

575 The low level system for which the wave functions are to be 

576 calculated. 

577 args : tuple of arguments, or empty tuple 

578 Positional arguments to pass to the function(s) which 

579 evaluate the hamiltonian matrix elements. 

580 Deprecated in favor of 'params' (and mutually exclusive with it). 

581 check_hermiticity : ``bool`` 

582 Check if the Hamiltonian matrices are Hermitian. 

583 params : dict, optional 

584 Dictionary of parameter names and their values. Mutually exclusive 

585 with 'args'. 

586 

587 Notes 

588 ----- 

589 The returned object can be itself called like a function. Given a lead 

590 number, it returns a 2d NumPy array that contains the wave function 

591 within the scattering region due to each incoming mode of the given 

592 lead. Index 0 is the mode number, index 1 is the orbital number. 

593 

594 The modes appear in the same order as the negative velocity modes in 

595 `kwant.physics.PropagatingModes`. In Kwant's convention leads are attached 

596 so that their translational symmetry points *away* from the scattering 

597 region:: 

598 

599 left lead SR right lead 

600 /---------\ /---\ /---------\ 

601 ...-3-2-1-0-X-X-X-0-1-2-3-... 

602 

603 This means that incoming modes (coming from infinity towards the 

604 scattering region) have *negative* velocity with respect to the 

605 lead's symmetry direction. 

606 

607 Examples 

608 -------- 

609 >>> wf = kwant.solvers.default.wave_function(some_syst, some_energy) 

610 >>> wfs_of_lead_2 = wf(2) 

611 

612 """ 

613 return WaveFunction(self, sys, energy, args, check_hermiticity, params) 

614 

615 

616class WaveFunction: 

617 def __init__(self, solver, sys, energy, args, check_hermiticity, params): 

618 syst = sys # ensure consistent naming across function bodies 

619 ensure_isinstance(syst, system.System) 

620 

621 modes_leads = [ 

622 i for i, l in enumerate(syst.leads) 

623 if not system.is_selfenergy_lead(l) 

624 ] 

625 

626 linsys = solver._make_linear_sys(syst, modes_leads, energy, 

627 args, check_hermiticity, 

628 params=params)[0] 

629 self.modes_leads = modes_leads 

630 self.solve = solver._solve_linear_sys 

631 self.rhs = linsys.rhs 

632 self.factorized_h = solver._factorized(linsys.lhs) 

633 self.num_orb = linsys.num_orb 

634 

635 def __call__(self, lead): 

636 if lead not in self.modes_leads: 

637 msg = ('Scattering wave functions for leads with only ' 

638 'self-energy are undefined.') 

639 raise ValueError(msg) 

640 

641 result = self.solve(self.factorized_h, self.rhs[lead], 

642 slice(self.num_orb)) 

643 return result.transpose() 

644 

645 

646class BlockResult(metaclass=abc.ABCMeta): 

647 """ 

648 ABC for a linear system solution with variable grouping. 

649 """ 

650 

651 def __init__(self, data, lead_info, out_leads, in_leads, sizes, 

652 current_conserving=False): 

653 self.data = data 

654 self.lead_info = lead_info 

655 self.out_leads = out_leads = list(out_leads) 

656 self.in_leads = in_leads = list(in_leads) 

657 self.sizes = sizes = np.array(sizes) 

658 self.current_conserving = current_conserving 

659 self.in_offsets = np.zeros(len(in_leads) + 1, int) 

660 self.in_offsets[1 :] = np.cumsum(self.sizes[in_leads]) 

661 self.out_offsets = np.zeros(len(self.out_leads) + 1, int) 

662 self.out_offsets[1 :] = np.cumsum(self.sizes[out_leads]) 

663 

664 def block_coords(self, lead_out, lead_in): 

665 """ 

666 Return slices corresponding to the block from lead_in to lead_out. 

667 """ 

668 return self.out_block_coords(lead_out), self.in_block_coords(lead_in) 

669 

670 def out_block_coords(self, lead_out): 

671 """Return a slice with the rows in the block corresponding to lead_out. 

672 """ 

673 lead_out = self.out_leads.index(lead_out) 

674 return slice(self.out_offsets[lead_out], 

675 self.out_offsets[lead_out + 1]) 

676 

677 def in_block_coords(self, lead_in): 

678 """ 

679 Return a slice with the columns in the block corresponding to lead_in. 

680 """ 

681 lead_in = self.in_leads.index(lead_in) 

682 return slice(self.in_offsets[lead_in], 

683 self.in_offsets[lead_in + 1]) 

684 

685 def submatrix(self, lead_out, lead_in): 

686 """Return the matrix elements from lead_in to lead_out.""" 

687 return self.data[self.block_coords(lead_out, lead_in)] 

688 

689 @abc.abstractmethod 

690 def num_propagating(self, lead): 

691 """Return the number of propagating modes in the lead.""" 

692 pass 

693 

694 @abc.abstractmethod 

695 def _transmission(self, lead_out, lead_in): 

696 """Return transmission from lead_in to lead_out. 

697 

698 The calculation is expected to be done directly from the corresponding 

699 matrix elements. 

700 """ 

701 pass 

702 

703 def transmission(self, lead_out, lead_in): 

704 """Return transmission from lead_in to lead_out. 

705 

706 If the option ``current_conserving`` has been enabled for this object, 

707 this method will deduce missing transmission values whenever possible. 

708 

709 Current conservation is enabled by default for objects returned by 

710 `~kwant.solvers.default.smatrix` and 

711 `~kwant.solvers.default.greens_function` whenever the Hamiltonian has 

712 been verified to be Hermitian (option ``check_hermiticity``, enabled by 

713 default). 

714 

715 """ 

716 chosen = [lead_out, lead_in] 

717 present = self.out_leads, self.in_leads 

718 # OK means: chosen (outgoing, incoming) lead is present. 

719 ok = [int(c in p) for c, p in zip(chosen, present)] 

720 if all(ok): 

721 return self._transmission(lead_out, lead_in) 

722 

723 if self.current_conserving: 723 ↛ 739line 723 didn't jump to line 739, because the condition on line 723 was never false

724 all_but_one = len(self.lead_info) - 1 

725 s_t = self._transmission 

726 if sum(ok) == 1: 

727 # Calculate the transmission element by summing over a single 

728 # row or column. 

729 sum_dim, ok_dim = ok 

730 if len(present[sum_dim]) == all_but_one: 

731 return (self.num_propagating(chosen[ok_dim]) - sum( 

732 s_t(*chosen) for chosen[sum_dim] in present[sum_dim])) 

733 elif all(len(p) == all_but_one for p in present): 733 ↛ exit,   733 ↛ 7362 missed branches: 1) line 733 didn't finish the generator expression on line 733, 2) line 733 didn't jump to line 736, because the condition on line 733 was never true

734 # Calculate the transmission element at the center of a single 

735 # missing row and a single missing column. 

736 return (sum(s_t(o, i) - self.num_propagating(o) if o == i else 0 

737 for o, i in product(*present))) 

738 

739 raise ValueError("Insufficient matrix elements to compute " 

740 "transmission({0}, {1}).".format(*chosen)) 

741 

742 def conductance_matrix(self): 

743 """Return the conductance matrix. 

744 

745 This is the matrix :math:`C` such that :math:`I = CV` where :math:`I` 

746 and :math:`V` are, respectively, the vectors of currents and voltages 

747 for each lead. 

748 

749 This matrix is useful for calculating non-local resistances. See 

750 Section 2.4 of the book by S. Datta. 

751 """ 

752 n = len(self.lead_info) 

753 rn = range(n) 

754 result = np.array([[-self.transmission(i, j) if i != j else 0 

755 for j in rn] for i in rn]) 

756 # Set the diagonal elements such that the sums of columns are zero. 

757 result.flat[::n+1] = -result.sum(axis=0) 

758 return result 

759 

760 

761class SMatrix(BlockResult): 

762 """A scattering matrix. 

763 

764 Transport properties can be easily accessed using the 

765 `~SMatrix.transmission` method (don't be fooled by the name, 

766 it can also compute reflection, which is just transmission from one 

767 lead back into the same lead.) 

768 

769 `SMatrix` however also allows for a more direct access to the result: The 

770 data stored in `SMatrix` is a scattering matrix with respect to lead modes 

771 and these modes themselves. The details of this data can be directly 

772 accessed through the instance variables `data` and `lead_info`. Subblocks 

773 of data corresponding to particular leads are conveniently obtained by 

774 `~SMatrix.submatrix`. 

775 

776 `SMatrix` also respects the conservation laws present in the lead, such as 

777 spin conservation, if they are declared during system construction. If 

778 queried with length-2 sequence the first number is the number of the lead, 

779 and the second number is the index of the corresponding conserved 

780 quantity. For example ``smatrix.transmission((1, 3), (0, 2))`` is 

781 transmission from block 2 of the conserved quantity in lead 0 to the block 

782 3 of the conserved quantity in lead 1. 

783 

784 Attributes 

785 ---------- 

786 data : NumPy array 

787 a matrix containing all the requested matrix elements of the scattering 

788 matrix. 

789 lead_info : list of data 

790 a list containing `kwant.physics.PropagatingModes` for each lead. 

791 If a lead is a selfenergy lead, then the corresponding entry 

792 in lead_info is a selfenergy. 

793 out_leads, in_leads : sequence of integers 

794 indices of the leads where current is extracted (out) or injected 

795 (in). Only those are listed for which SMatrix contains the 

796 calculated result. 

797 """ 

798 

799 def __init__(self, data, lead_info, out_leads, in_leads, 

800 current_conserving=False): 

801 # An equivalent condition to this is checked in 'Solver.smatrix', 

802 # but we add it here again as a sanity check. If 'lead_info' is not 

803 # a 'PropagatingModes' that means that the corresponding lead is a 

804 # selfenergy lead. Scattering matrix elements to/from a selfenergy lead 

805 # are not defined. 

806 assert not any( 

807 not isinstance(info, physics.PropagatingModes) 

808 and leadnum in (out_leads + in_leads) 

809 for leadnum, info in enumerate(lead_info) 

810 ) 

811 

812 # 'getattr' in order to handle self-energy leads, for which 

813 # 'info' is just a matrix (so no 'momenta'). 

814 sizes = [ 

815 len(getattr(info, "momenta", [])) // 2 for info in lead_info 

816 ] 

817 super().__init__(data, lead_info, out_leads, in_leads, 

818 sizes, current_conserving) 

819 # in_offsets marks beginnings and ends of blocks in the scattering 

820 # matrix corresponding to the in leads. The end of the block 

821 # corresponding to lead i is taken as the beginning of the block of 

822 # i+1. Same for out leads. For each lead block of the scattering 

823 # matrix, we want to pick out the computed blocks of the conservation 

824 # law. The offsets of these symmetry blocks are stored in 

825 # block_offsets, for all leads. List of lists containing the sizes of 

826 # symmetry blocks in all leads. leads_block_sizes[i][j] is the number 

827 # of left or right moving modes in symmetry block j of lead 

828 # i. len(leads_block_sizes[i]) is the number of blocks in lead i. 

829 leads_block_sizes = [] 

830 for info in self.lead_info: 

831 # If a lead does not have block structure, append None. 

832 leads_block_sizes.append(getattr(info, 'block_nmodes', None)) 

833 self.leads_block_sizes = leads_block_sizes 

834 block_offsets = [] 

835 for lead_block_sizes in self.leads_block_sizes: # Cover all leads 

836 if lead_block_sizes is None: 

837 block_offsets.append(lead_block_sizes) 

838 else: 

839 block_offset = np.zeros(len(lead_block_sizes) + 1, int) 

840 block_offset[1:] = np.cumsum(lead_block_sizes) 

841 block_offsets.append(block_offset) 

842 # Symmetry block offsets for all leads - or None if lead does not have 

843 # blocks. 

844 block_offsets = np.array(block_offsets, dtype=object) 

845 # Pick out symmetry block offsets for in and out leads 

846 self._in_block_offsets = block_offsets[list(self.in_leads)] 

847 self._out_block_offsets = block_offsets[list(self.out_leads)] 

848 # Block j of in lead i starts at in_block_offsets[i][j] 

849 

850 def out_block_coords(self, lead_out): 

851 """Return a slice with the rows in the block corresponding to lead_out. 

852 """ 

853 if isinstance(lead_out, Integral): 853 ↛ 858line 853 didn't jump to line 858, because the condition on line 853 was never false

854 lead_out = self.out_leads.index(lead_out) 

855 return slice(self.out_offsets[lead_out], 

856 self.out_offsets[lead_out + 1]) 

857 else: 

858 lead_ind, block_ind = lead_out 

859 lead_ind = self.out_leads.index(lead_ind) 

860 return slice(self.out_offsets[lead_ind] + 

861 self._out_block_offsets[lead_ind][block_ind], 

862 self.out_offsets[lead_ind] + 

863 self._out_block_offsets[lead_ind][block_ind + 1]) 

864 

865 def in_block_coords(self, lead_in): 

866 """ 

867 Return a slice with the columns in the block corresponding to lead_in. 

868 """ 

869 if isinstance(lead_in, Integral): 869 ↛ 874line 869 didn't jump to line 874, because the condition on line 869 was never false

870 lead_in = self.in_leads.index(lead_in) 

871 return slice(self.in_offsets[lead_in], 

872 self.in_offsets[lead_in + 1]) 

873 else: 

874 lead_ind, block_ind = lead_in 

875 lead_ind = self.in_leads.index(lead_ind) 

876 return slice(self.in_offsets[lead_ind] + 

877 self._in_block_offsets[lead_ind][block_ind], 

878 self.in_offsets[lead_ind] + 

879 self._in_block_offsets[lead_ind][block_ind + 1]) 

880 

881 def _transmission(self, lead_out, lead_in): 

882 return np.linalg.norm(self.submatrix(lead_out, lead_in)) ** 2 

883 

884 def transmission(self, lead_out, lead_in): 

885 """Return transmission from lead_in to lead_out. 

886 

887 If `lead_in` or `lead_out` are a length-2 sequence, the first number is 

888 the number of the lead, and the second number indexes the eigenvalue of 

889 the conserved quantity in that lead (e.g. spin) if it was specified. 

890 """ 

891 # If transmission is between leads and not subblocks, we can use 

892 # unitarity (like BlockResult does). 

893 if isinstance(lead_in, Integral) and isinstance(lead_out, Integral): 893 ↛ 896line 893 didn't jump to line 896, because the condition on line 893 was never false

894 return super().transmission(lead_out, lead_in) 

895 else: 

896 return self._transmission(lead_out, lead_in) 

897 

898 def num_propagating(self, lead): 

899 """Return the number of propagating modes in the lead.""" 

900 return self.sizes[lead] 

901 

902 def __repr__(self): 

903 return ("SMatrix(data=%r, lead_info=%r, out_leads=%r, in_leads=%r)" % 

904 (self.data, self.lead_info, self.out_leads, self.in_leads)) 

905 

906 

907class GreensFunction(BlockResult): 

908 """ 

909 Retarded Green's function. 

910 

911 Transport properties can be easily accessed using the 

912 `~GreensFunction.transmission` method (don't be fooled by the name, it can 

913 also compute reflection, which is just transmission from one lead back into 

914 the same lead). 

915 

916 `GreensFunction` however also allows for a more direct access to the 

917 result: The data stored in `GreensFunction` is the real-space Green's 

918 function. The details of this data can be directly accessed through the 

919 instance variables `data` and `lead_info`. Subblocks of data corresponding 

920 to particular leads are conveniently obtained by 

921 `~GreensFunction.submatrix`. 

922 

923 Attributes 

924 ---------- 

925 data : NumPy array 

926 a matrix containing all the requested matrix elements of Green's 

927 function. 

928 lead_info : list of matrices 

929 a list with self-energies of each lead. 

930 out_leads, in_leads : sequence of integers 

931 indices of the leads where current is extracted (out) or injected 

932 (in). Only those are listed for which SMatrix contains the 

933 calculated result. 

934 """ 

935 

936 def __init__(self, data, lead_info, out_leads, in_leads, 

937 current_conserving=False): 

938 sizes = [i.shape[0] for i in lead_info] 

939 super().__init__( 

940 data, lead_info, out_leads, in_leads, sizes, current_conserving) 

941 

942 def num_propagating(self, lead): 

943 """Return the number of propagating modes in the lead.""" 

944 gamma = 1j * (self.lead_info[lead] - 

945 self.lead_info[lead].conj().T) 

946 

947 # The number of channels is given by the number of 

948 # nonzero eigenvalues of Gamma 

949 # rationale behind the threshold from 

950 # Golub; van Loan, chapter 5.5.8 

951 eps = np.finfo(gamma.dtype).eps * 1000 

952 return np.sum(np.linalg.eigvalsh(gamma) > 

953 eps * np.linalg.norm(gamma, np.inf)) 

954 

955 def _a_ttdagger_a_inv(self, lead_out, lead_in): 

956 """Return t * t^dagger in a certain basis.""" 

957 gf = self.submatrix(lead_out, lead_in) 

958 factors = [] 

959 for lead, gf2 in ((lead_out, gf), (lead_in, gf.conj().T)): 

960 possible_se = self.lead_info[lead] 

961 factors.append(1j * (possible_se - possible_se.conj().T)) 

962 factors.append(gf2) 

963 return reduce(np.dot, factors) 

964 

965 def _transmission(self, lead_out, lead_in): 

966 attdagainv = self._a_ttdagger_a_inv(lead_out, lead_in) 

967 result = np.trace(attdagainv).real 

968 if lead_out == lead_in: 

969 # For reflection we have to be more careful 

970 sigma = self.lead_info[lead_in] 

971 gamma = 1j * (sigma - sigma.conj().T) 

972 gf = self.submatrix(lead_out, lead_in) 

973 

974 # The number of channels is given by the number of 

975 # nonzero eigenvalues of Gamma 

976 # rationale behind the threshold from 

977 # Golub; van Loan, chapter 5.5.8 

978 # We use ‖Σ‖, not ‖Γ‖, for the tolerance as ‖Γ‖~0 when there 

979 # are no open modes. 

980 eps = ( 

981 1e6 * np.finfo(gamma.dtype).eps * np.linalg.norm(sigma, np.inf) 

982 ) 

983 N = np.sum(np.linalg.eigvalsh(gamma) > eps) 

984 

985 result += 2 * np.trace(np.dot(gamma, gf)).imag + N 

986 

987 return result 

988 

989 def __repr__(self): 

990 return ("GreensFunction(data=%r, lead_info=%r, " 

991 "out_leads=%r, in_leads=%r)" % 

992 (self.data, self.lead_info, self.out_leads, self.in_leads))