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.
9__all__ = ['SparseSolver', 'SMatrix', 'GreensFunction']
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
22LinearSys = namedtuple('LinearSys', ['lhs', 'rhs', 'indices', 'num_orb'])
25class SparseSolver(metaclass=abc.ABCMeta):
26 """Solver class for computing physical quantities based on solving
27 a liner system of equations.
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
33 - `_factorize` and `_solve_linear_sys`, that solve a linear system of
34 equations
36 and the following properties:
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'}.
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 """
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`.
55 Parameters
56 ----------
57 a : a scipy.sparse.coo_matrix sparse matrix.
59 Returns
60 -------
61 factorized_a : object
62 factorized lhs to be used with `_solve_linear_sys`.
63 """
64 pass
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.
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
81 Returns
82 -------
83 output : NumPy matrix
84 Solution to the system of equations.
85 """
86 pass
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.
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'.
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.
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.
134 Notes
135 -----
136 All the leads should implement a method `modes` if `realspace=False`
137 and a method `selfenergy`.
139 The system of equations that is created will be described in detail
140 elsewhere.
141 """
143 syst = sys # ensure consistent naming across function bodies
144 ensure_isinstance(syst, system.System)
146 splhsmat = getattr(sp, self.lhsformat + '_matrix')
147 sprhsmat = getattr(sp, self.rhsformat + '_matrix')
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]
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.')
167 offsets = np.empty(norb.shape[0] + 1, int)
168 offsets[0] = 0
169 offsets[1 :] = np.cumsum(norb)
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
211 if len(u) == 0:
212 rhs.append(None)
213 continue
215 indices.append(np.arange(lhs.shape[0], lhs.shape[0] + nprop))
217 u_out, ulinv_out = u[:, nprop:], ulinv[:, nprop:]
218 u_in, ulinv_in = u[:, :nprop], ulinv[:, :nprop]
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)]
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))
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]))
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
241 lhs = sp.bmat([[lhs, vdaguout_sp], [v_sp, lead_mat]],
242 format=self.lhsformat)
244 if leadnum in in_leads and nprop > 0:
245 vdaguin_sp = transf.T * sp.csc_matrix(
246 -np.dot(svd_v, u_in))
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)]
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))
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,))
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])
288 zero_mat = sprhsmat((zero_rows, mats[0].shape[1]))
289 bmat = [[mats[0]], [mats[1]], [zero_mat]]
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')
298 return LinearSys(lhs, rhs, indices, num_orb), lead_info
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.
307 An alias exists for this common name: ``kwant.smatrix``.
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'.
334 Returns
335 -------
336 output : `~kwant.solvers.common.SMatrix`
337 See the notes below and `~kwant.solvers.common.SMatrix`
338 documentation.
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`.
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.
350 Both `in_leads` and `out_leads` must be sorted and may only contain
351 unique entries.
352 """
354 syst = sys # ensure consistent naming across function bodies
355 ensure_isinstance(syst, system.System)
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.")
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 )
381 linsys, lead_info = self._make_linear_sys(syst, in_leads, energy, args,
382 check_hermiticity, False,
383 params=params)
385 kept_vars = np.concatenate([coords for i, coords in
386 enumerate(linsys.indices) if i in
387 out_leads])
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)
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)
400 return SMatrix(data, lead_info, out_leads, in_leads, check_hermiticity)
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.
409 An alias exists for this common name: ``kwant.greens_function``.
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'.
436 Returns
437 -------
438 output : `~kwant.solvers.common.GreensFunction`
439 See the notes below and `~kwant.solvers.common.GreensFunction`
440 documentation.
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.
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.
455 Both `in_leads` and `out_leads` must be sorted and may only contain
456 unique entries.
457 """
459 syst = sys # ensure consistent naming across function bodies
460 ensure_isinstance(syst, system.System)
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.")
477 linsys, lead_info = self._make_linear_sys(syst, in_leads, energy, args,
478 check_hermiticity, True,
479 params=params)
481 kept_vars = np.concatenate([coords for i, coords in
482 enumerate(linsys.indices) if i in
483 out_leads])
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)
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)
496 return GreensFunction(data, lead_info, out_leads, in_leads,
497 check_hermiticity)
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.
505 An alias exists for this common name: ``kwant.ldos``.
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'.
524 Returns
525 -------
526 ldos : a NumPy array
527 Local density of states at each orbital of the system.
528 """
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.")
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.")
542 linsys = self._make_linear_sys(syst, range(len(syst.leads)), energy,
543 args, check_hermiticity,
544 params=params)[0]
546 ldos = np.zeros(linsys.num_orb, float)
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
552 factored = self._factorized(linsys.lhs)
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)
561 return ldos * (0.5 / np.pi)
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.
570 An alias exists for this common name: ``kwant.wave_function``.
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'.
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.
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::
599 left lead SR right lead
600 /---------\ /---\ /---------\
601 ...-3-2-1-0-X-X-X-0-1-2-3-...
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.
607 Examples
608 --------
609 >>> wf = kwant.solvers.default.wave_function(some_syst, some_energy)
610 >>> wfs_of_lead_2 = wf(2)
612 """
613 return WaveFunction(self, sys, energy, args, check_hermiticity, params)
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)
621 modes_leads = [
622 i for i, l in enumerate(syst.leads)
623 if not system.is_selfenergy_lead(l)
624 ]
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
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)
641 result = self.solve(self.factorized_h, self.rhs[lead],
642 slice(self.num_orb))
643 return result.transpose()
646class BlockResult(metaclass=abc.ABCMeta):
647 """
648 ABC for a linear system solution with variable grouping.
649 """
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])
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)
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])
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])
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)]
689 @abc.abstractmethod
690 def num_propagating(self, lead):
691 """Return the number of propagating modes in the lead."""
692 pass
694 @abc.abstractmethod
695 def _transmission(self, lead_out, lead_in):
696 """Return transmission from lead_in to lead_out.
698 The calculation is expected to be done directly from the corresponding
699 matrix elements.
700 """
701 pass
703 def transmission(self, lead_out, lead_in):
704 """Return transmission from lead_in to lead_out.
706 If the option ``current_conserving`` has been enabled for this object,
707 this method will deduce missing transmission values whenever possible.
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).
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)
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)))
739 raise ValueError("Insufficient matrix elements to compute "
740 "transmission({0}, {1}).".format(*chosen))
742 def conductance_matrix(self):
743 """Return the conductance matrix.
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.
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
761class SMatrix(BlockResult):
762 """A scattering matrix.
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.)
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`.
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.
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 """
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 )
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]
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])
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])
881 def _transmission(self, lead_out, lead_in):
882 return np.linalg.norm(self.submatrix(lead_out, lead_in)) ** 2
884 def transmission(self, lead_out, lead_in):
885 """Return transmission from lead_in to lead_out.
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)
898 def num_propagating(self, lead):
899 """Return the number of propagating modes in the lead."""
900 return self.sizes[lead]
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))
907class GreensFunction(BlockResult):
908 """
909 Retarded Green's function.
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).
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`.
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 """
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)
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)
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))
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)
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)
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)
985 result += 2 * np.trace(np.dot(gamma, gf)).imag + N
987 return result
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))