diff --git a/doc/source/images/ab_ring.py.diff b/doc/source/images/ab_ring.py.diff index d4f2c825bfc628c74edc75ffe0919a26969ae6a0..ce73ce4a322b60683d521ace7b67a6ba7e9e79ce 100644 --- a/doc/source/images/ab_ring.py.diff +++ b/doc/source/images/ab_ring.py.diff @@ -80,7 +80,7 @@ # compute conductance @@ -87,18 +133,31 @@ - smatrix = kwant.solve(sys, energy, args=[flux]) + smatrix = kwant.smatrix(sys, energy, args=[flux]) data.append(smatrix.transmission(1, 0)) - pyplot.figure() diff --git a/doc/source/images/graphene.py.diff b/doc/source/images/graphene.py.diff index 59f9d43e6bf319913fefd96996c7b50a69124050..741576cda2885471d88c694892569aefe7d4b025 100644 --- a/doc/source/images/graphene.py.diff +++ b/doc/source/images/graphene.py.diff @@ -9,7 +9,7 @@ # Define the graphene lattice @@ -100,22 +101,40 @@ - smatrix = kwant.solve(sys, energy) + smatrix = kwant.smatrix(sys, energy) data.append(smatrix.transmission(0, 1)) - pyplot.figure() diff --git a/doc/source/images/quantum_well.py.diff b/doc/source/images/quantum_well.py.diff index 9262fd0197e4e0816ce55440156d62d2a259e4e0..c1bcd61be6cd89cd8053e7d10c6027f851393933 100644 --- a/doc/source/images/quantum_well.py.diff +++ b/doc/source/images/quantum_well.py.diff @@ -9,7 +9,7 @@ def make_system(a=1, t=1.0, W=10, L=30, L_well=10): @@ -52,19 +53,25 @@ - smatrix = kwant.solve(sys, energy, args=[-welldepth]) + smatrix = kwant.smatrix(sys, energy, args=[-welldepth]) data.append(smatrix.transmission(1, 0)) - pyplot.figure() diff --git a/doc/source/images/spin_orbit.py.diff b/doc/source/images/spin_orbit.py.diff index d16e4a8556610742e45fe840972a40f47edb12e2..5f737a9dd5ad555c8aae2c3baa2ab92bc1a8bf92 100644 --- a/doc/source/images/spin_orbit.py.diff +++ b/doc/source/images/spin_orbit.py.diff @@ -9,7 +9,7 @@ # define Pauli-matrices for convenience sigma_0 = tinyarray.array([[1, 0], [0, 1]]) @@ -67,19 +68,24 @@ - smatrix = kwant.solve(sys, energy) + smatrix = kwant.smatrix(sys, energy) data.append(smatrix.transmission(1, 0)) - pyplot.figure() diff --git a/doc/source/reference/kwant.rst b/doc/source/reference/kwant.rst index 13481d3d6f603e381b627d743b582a363eabb504..9c3b6a92dd08bc1d5b6c7487c35ca8e66c0c1bbe 100644 --- a/doc/source/reference/kwant.rst +++ b/doc/source/reference/kwant.rst @@ -38,4 +38,4 @@ From ``kwant.solvers.default`` ------------------------------ .. autosummary:: - solve + smatrix diff --git a/doc/source/reference/kwant.solvers.rst b/doc/source/reference/kwant.solvers.rst index 586b27e4562189882e9a5023eeddf606b79bda5f..94a3369b758d66000e9b9d05ed477dc998b4ce68 100644 --- a/doc/source/reference/kwant.solvers.rst +++ b/doc/source/reference/kwant.solvers.rst @@ -26,18 +26,26 @@ reasons to use another. The following functions are provided. .. autosummary:: :toctree: generated/ - solve + smatrix + greens_function wave_function ldos -``solve`` returns an object of the following type: +``smatrix`` returns an object of the following type: .. module:: kwant.solvers .. autosummary:: :toctree: generated/ - kwant.solvers.common.BlockResult + kwant.solvers.common.SMatrix + +The analog of ``smatrix``, ``greens_function`` accordingly returns: + +.. autosummary:: + :toctree: generated/ + + kwant.solvers.common.GreensFunction Being just a thin wrapper around other solvers, the default solver selectively imports their functionality. To find out the origin of any function in this diff --git a/doc/source/tutorial/ab_ring.py b/doc/source/tutorial/ab_ring.py index 8b302964ad3b040633852004eb608a200a620212..69aeb9ccc4851b500dcedf9a615168ab1c96a320 100644 --- a/doc/source/tutorial/ab_ring.py +++ b/doc/source/tutorial/ab_ring.py @@ -94,7 +94,7 @@ def plot_conductance(sys, energy, fluxes): normalized_fluxes = [flux / (2 * pi) for flux in fluxes] data = [] for flux in fluxes: - smatrix = kwant.solve(sys, energy, args=[flux]) + smatrix = kwant.smatrix(sys, energy, args=[flux]) data.append(smatrix.transmission(1, 0)) pyplot.figure() diff --git a/doc/source/tutorial/graphene.py b/doc/source/tutorial/graphene.py index abdb352179103005d19ac2a3cfc5af7bde27e6bb..2c652d0bec918dfe0bc9e3bea12071e2d98087fe 100644 --- a/doc/source/tutorial/graphene.py +++ b/doc/source/tutorial/graphene.py @@ -113,7 +113,7 @@ def plot_conductance(sys, energies): # Compute transmission as a function of energy data = [] for energy in energies: - smatrix = kwant.solve(sys, energy) + smatrix = kwant.smatrix(sys, energy) data.append(smatrix.transmission(0, 1)) pyplot.figure() diff --git a/doc/source/tutorial/quantum_well.py b/doc/source/tutorial/quantum_well.py index fb901bb616629399099ba74a0b738d13adc4f66b..50001f70c5cf8af5e7860683cfeec0e4201f9f9b 100644 --- a/doc/source/tutorial/quantum_well.py +++ b/doc/source/tutorial/quantum_well.py @@ -54,7 +54,7 @@ def plot_conductance(sys, energy, welldepths): # Compute conductance data = [] for welldepth in welldepths: - smatrix = kwant.solve(sys, energy, args=[-welldepth]) + smatrix = kwant.smatrix(sys, energy, args=[-welldepth]) data.append(smatrix.transmission(1, 0)) pyplot.figure() diff --git a/doc/source/tutorial/quantum_wire.py b/doc/source/tutorial/quantum_wire.py index b1067a785dabbd9ad8fa15aeb9d3f7b1b9c9e983..35e1ac6a73114618c51dde629a11acdf0e22ca0c 100644 --- a/doc/source/tutorial/quantum_wire.py +++ b/doc/source/tutorial/quantum_wire.py @@ -97,8 +97,8 @@ data = [] for ie in xrange(100): energy = ie * 0.01 - # compute the scattering matrix at energy energy - smatrix = kwant.solve(sys, energy) + # compute the scattering matrix at a given energy + smatrix = kwant.smatrix(sys, energy) # compute the transmission probability from lead 0 to # lead 1 diff --git a/doc/source/tutorial/quantum_wire_revisited.py b/doc/source/tutorial/quantum_wire_revisited.py index 84ac6afcd5e04075309b9deb726bbccfc5c1ec90..c52425686b5446e18204dc5b6f5713b186a7558a 100644 --- a/doc/source/tutorial/quantum_wire_revisited.py +++ b/doc/source/tutorial/quantum_wire_revisited.py @@ -54,7 +54,7 @@ def plot_conductance(sys, energies): # Compute conductance data = [] for energy in energies: - smatrix = kwant.solve(sys, energy) + smatrix = kwant.smatrix(sys, energy) data.append(smatrix.transmission(1, 0)) pyplot.figure() diff --git a/doc/source/tutorial/spin_orbit.py b/doc/source/tutorial/spin_orbit.py index f3b9f2ca98d9f768489ad04d0dd2471e467712e6..8f3a35a681557af5b7d89a09a5bf0f8a91160258 100644 --- a/doc/source/tutorial/spin_orbit.py +++ b/doc/source/tutorial/spin_orbit.py @@ -72,7 +72,7 @@ def plot_conductance(sys, energies): # Compute conductance data = [] for energy in energies: - smatrix = kwant.solve(sys, energy) + smatrix = kwant.smatrix(sys, energy) data.append(smatrix.transmission(1, 0)) pyplot.figure() diff --git a/doc/source/tutorial/superconductor_transport.py b/doc/source/tutorial/superconductor_transport.py index c9b1ad47fbfefe8f862ed2c935f48fc3662838df..3f95387515c3967c25a8869184fb9a0030d2c0ce 100644 --- a/doc/source/tutorial/superconductor_transport.py +++ b/doc/source/tutorial/superconductor_transport.py @@ -87,7 +87,7 @@ def plot_conductance(sys, energies): # Compute conductance data = [] for energy in energies: - smatrix = kwant.solve(sys, energy) + smatrix = kwant.smatrix(sys, energy) # Conductance is N - R_ee + R_he data.append(smatrix.submatrix(0, 0).shape[0] - smatrix.transmission(0, 0) + diff --git a/doc/source/tutorial/tutorial1.rst b/doc/source/tutorial/tutorial1.rst index 86157680ae1a03498608d604b637307302ed3b94..71efbed5c1d79d9eed39189105671c362041884b 100644 --- a/doc/source/tutorial/tutorial1.rst +++ b/doc/source/tutorial/tutorial1.rst @@ -142,11 +142,12 @@ its conductance as a function of energy: :start-after: #HIDDEN_BEGIN_buzn :end-before: #HIDDEN_END_buzn -We use ``kwant.solve`` which is a short name for `kwant.solvers.default.solve` -of the default solver module `kwant.solvers.default`. ``kwant.solve`` computes -the scattering matrix ``smatrix`` solving a sparse linear system. ``smatrix`` -itself allows to directly compute the total transmission probability from lead -0 to lead 1 as ``smatrix.transmission(1, 0)``. +We use ``kwant.smatrix`` which is a short name for +`kwant.solvers.default.smatrix` of the default solver module +`kwant.solvers.default`. ``kwant.smatrix`` computes the scattering matrix +``smatrix`` solving a sparse linear system. ``smatrix`` itself allows to +directly compute the total transmission probability from lead 0 to lead 1 as +``smatrix.transmission(1, 0)``. Finally we can use `matplotlib` to make a plot of the computed data (although writing to file and using an external viewer such as diff --git a/doc/source/tutorial/tutorial2.rst b/doc/source/tutorial/tutorial2.rst index ce779b88818ced22df479c69fc14ccf91e0a7ecd..4750451e60d81202efef92f0c5f03695fb31ff8b 100644 --- a/doc/source/tutorial/tutorial2.rst +++ b/doc/source/tutorial/tutorial2.rst @@ -163,7 +163,7 @@ Kwant now allows us to pass a function as a value to For each lattice point, the corresponding site is then passed as the first argument to the function `onsite`. The values of any additional parameters, which can be used to alter the Hamiltonian matrix elements -at a later stage, are specified later during the call to `solve`. +at a later stage, are specified later during the call to `smatrix`. Note that we had to define `onsite`, as it is not possible to mix values and functions as in ``sys[...] = 4 * t + potential``. @@ -179,7 +179,7 @@ Finally, we compute the transmission probability: :start-after: #HIDDEN_BEGIN_sqvr :end-before: #HIDDEN_END_sqvr -``kwant.solve`` allows us to specify a list, `args`, that will be passed as +``kwant.smatrix`` allows us to specify a list, `args`, that will be passed as additional arguments to the functions that provide the Hamiltonian matrix elements. In this example we are able to solve the system for different depths of the potential well by passing the potential value. We obtain the result: diff --git a/doc/source/tutorial/tutorial5.rst b/doc/source/tutorial/tutorial5.rst index ff078636ef379181a19600971dca182fe9477884..63e891bce6c0f35247dd9c07fda273e71f0e2fb1 100644 --- a/doc/source/tutorial/tutorial5.rst +++ b/doc/source/tutorial/tutorial5.rst @@ -168,9 +168,9 @@ above the gap. At the gap edge, we observe a resonant Andreev reflection. - It is in fact possible to separate electron and hole degrees of freedom in the scattering matrix, even if one uses matrices for these degrees of freedom. In the solve step, - `~kwant.solvers.common.SparseSolver.solve` returns an array containing + `~kwant.solvers.common.SparseSolver.smatrix` returns an array containing the transverse wave functions of the lead modes (in - `BlockResult.lead_info <kwant.solvers.common.BlockResult.lead_info>`. + `SMatrix.lead_info <kwant.solvers.common.SMatrix.lead_info>`. By inspecting the wave functions, electron and hole wave functions can be distinguished (they only have entries in either the electron part *or* the hole part. If you encounter modes diff --git a/doc/source/whatsnew/1.0.rst b/doc/source/whatsnew/1.0.rst index 36a09ebb54e09ec2d52d792d9a475d7dac355bf2..c9a3bc69e6e6aabdfea195efa13937c9865d0014 100644 --- a/doc/source/whatsnew/1.0.rst +++ b/doc/source/whatsnew/1.0.rst @@ -54,6 +54,13 @@ Some renames * ``wave_func`` has been renamed to `~kwant.solvers.default.wave_function`, * ``MonatomicLattice`` has been renamed to `~kwant.lattice.Monatomic`, * ``PolyatomicLattice`` has been renamed to `~kwant.lattice.Polyatomic`. +* ``solve`` was split into two functions: `~kwant.solvers.default.smatrix`, and + `~kwant.solvers.default.greens_function`. The former calculates the + scattering matrix, the latter the retarded Green's function between the sites + adjacent to the leads. It is temporarily not possible to mix self-energy and + modes leads within the same system. +* The object that contained the results, `BlockResult` was also split into + `~kwant.solvers.common.SMatrix` and `~kwant.solvers.common.GreensFunction`. Band structure plots -------------------- @@ -98,7 +105,7 @@ regular function parameters (e.g. global variables). Now the value functions may accept arbitrary arguments after the `Site` arguments. These extra arguments can be specified when -`~kwant.solvers.default.solve` is called by setting the arguments: +`~kwant.solvers.default.smatrix` is called by setting the arguments: args A tuple of values to be passed as the positional arguments to the @@ -114,10 +121,10 @@ the following prototype:: ... then the values of `t`, `B` and `pot` for which to solve the system can be -passed to `solve` like this:: +passed to `~kwant.solvers.default.smatrix` like this:: - kwant.solve(sys, energy, - args=(2., 3., 4.)) + kwant.smatrix(sys, energy, + args=(2., 3., 4.)) With many parameters it can be less error-prone to collect all of them into a single object and pass this object as the single argument. Such a parameter @@ -137,7 +144,7 @@ collection could be a dictionary, or a class instance, for example:: params = SimpleNamespace(t=1, mu=2) for params.B in B_values: - kwant.solve(sys, energy, args=[params]) + kwant.smatrix(sys, energy, args=[params]) Arguments can be passed in an equivalent way to `~kwant.solvers.default.wave_function`, @@ -168,7 +175,7 @@ space, as well as their velocities and momenta. All these quantities were previously not directly available. The second object contains the propagating and evanescent modes in the compressed format expected by the sparse solver (previously this was the sole output of `~kwant.physics.modes`). Accordingly, -the `lead_info` attribute of `~kwant.solvers.default.BlockResult` contains the +the `lead_info` attribute of `~kwant.solvers.common.SMatrix` contains the real space information about the modes in the leads (a list of `~kwant.physics.PropagatingModes` objects). diff --git a/examples/advanced_construction.py b/examples/advanced_construction.py index f47d4100f5f3379d313b51d3f375a739996aab6e..2165b4047c7dea89c6d6253c050b88d934316f56 100644 --- a/examples/advanced_construction.py +++ b/examples/advanced_construction.py @@ -48,7 +48,7 @@ def make_system(R): def main(): sys = make_system(100) - print kwant.solve(sys, 1.1, args=[0.1]).transmission(0, 1) + print kwant.smatrix(sys, 1.1, args=[0.1]).transmission(0, 1) if __name__ == '__main__': diff --git a/examples/square.py b/examples/square.py index 0d31155ba0dcb9ada42db3ae4fffef2a30992071..c7df5d54200b247cb7584db9a7e0a9699ebf2f58 100644 --- a/examples/square.py +++ b/examples/square.py @@ -99,7 +99,7 @@ class System(kwant.system.FiniteSystem): def main(): sys = System((10, 5), 1) energies = [0.04 * i for i in xrange(100)] - data = [kwant.solve(sys, energy).transmission(1, 0) + data = [kwant.smatrix(sys, energy).transmission(1, 0) for energy in energies] pyplot.plot(energies, data) diff --git a/kwant/__init__.py b/kwant/__init__.py index 5fdc46367ed608c6826130ce21085890824fdac5..0f11f86aa02e5235dd3000707f361be0f57c52d1 100644 --- a/kwant/__init__.py +++ b/kwant/__init__.py @@ -19,7 +19,7 @@ from .lattice import TranslationalSymmetry __all__.append('TranslationalSymmetry') # Make kwant.solvers.default.solve available as kwant.solve. -solve = solvers.default.solve +smatrix = solvers.default.smatrix __all__.append('solve') # Importing plotter might not work, but this does not have to be a problem -- diff --git a/kwant/physics/noise.py b/kwant/physics/noise.py index 92b9a9cbbbd0135ca168f4777a917508335665d5..847e2910d2e1265492d19b809cf29f7cf087e4db 100644 --- a/kwant/physics/noise.py +++ b/kwant/physics/noise.py @@ -16,7 +16,7 @@ def two_terminal_shotnoise(smatrix): Parameters ---------- - smatrix : `~kwant.solvers.common.BlockResult` instance + smatrix : `~kwant.solvers.common.SMatrix` instance A two terminal scattering matrix. Returns @@ -25,6 +25,10 @@ def two_terminal_shotnoise(smatrix): Shot noise measured in noise quanta `2 e^3 |V| / pi hbar`. """ + if hasattr(smatrix.lead_info[0], 'shape'): + raise NotImplementedError("Noise expressions in terms of Green's " + "functions are not implemented.") + if len(smatrix.lead_info) != 2: raise ValueError("Only works for two-terminal systems!") diff --git a/kwant/physics/tests/test_noise.py b/kwant/physics/tests/test_noise.py index adc71d78e4c809e15623a2742b1ad2b2cda9b877..18b2c3aedf99b169fb824a3d87c01fc6d9ab1ef5 100644 --- a/kwant/physics/tests/test_noise.py +++ b/kwant/physics/tests/test_noise.py @@ -36,7 +36,7 @@ def test_multiterminal_input(): sys = twoterminal_system() sys.attach_lead(sys.leads[0].builder) - sol = kwant.solve(sys.finalized(), out_leads=[0], in_leads=[0]) + sol = kwant.smatrix(sys.finalized(), out_leads=[0], in_leads=[0]) assert_raises(ValueError, two_terminal_shotnoise, sol) @@ -45,7 +45,7 @@ def test_twoterminal(): fsys = twoterminal_system().finalized() - sol = kwant.solve(fsys) + sol = kwant.smatrix(fsys) t = sol.submatrix(1, 0) Tn = np.linalg.eigvalsh(np.dot(t, t.conj().T)) noise_should_be = np.sum(Tn * (1 - Tn)) diff --git a/kwant/solvers/common.py b/kwant/solvers/common.py index bd839057c302bb9cd4a9657656fc516811606dd8..c8b95ba5521fd16e07eb05d648f0d84323327967 100644 --- a/kwant/solvers/common.py +++ b/kwant/solvers/common.py @@ -6,7 +6,7 @@ # the AUTHORS file at the top-level directory of this distribution and at # http://kwant-project.org/authors. -__all__ = ['SparseSolver', 'BlockResult'] +__all__ = ['SparseSolver', 'SMatrix', 'GreensFunction'] from collections import namedtuple import abc @@ -90,7 +90,7 @@ class SparseSolver(object): """ pass - def _make_linear_sys(self, sys, in_leads, energy=0, force_realspace=False, + def _make_linear_sys(self, sys, in_leads, energy=0, realspace=False, check_hermiticity=True, args=()): """Make a sparse linear system of equations defining a scattering problem. @@ -104,7 +104,7 @@ class SparseSolver(object): Numbers of leads in which current or wave function is injected. energy : number Excitation energy at which to solve the scattering problem. - force_realspace : bool + realspace : bool Calculate Green's function between the outermost lead sites, instead of lead modes. This is almost always more computationally expensive and less stable. @@ -126,20 +126,18 @@ class SparseSolver(object): total number of degrees of freedom in the scattering region. lead_info : list of objects - Contains one entry for each lead. For a lead defined as a - tight-binding system, this is an instance of - `~kwant.physics.PropagatingModes` with a corresponding format, - otherwise the lead self-energy matrix. + Contains one entry for each lead. If `realspace=False`, this is an + instance of `~kwant.physics.PropagatingModes` with a corresponding + format, otherwise the lead self-energy matrix. Notes ----- - Both incoming and outgoing leads can be defined via either self-energy, - or a low-level translationally invariant system. - The system of equations that is created is described in - kwant/doc/other/linear_system.pdf + All the leads should implement a method `modes` if `realspace=False` + and a method `selfenergy`. + The system of equations that is created will be described in detail + elsewhere. """ - splhsmat = getattr(sp, self.lhsformat + '_matrix') sprhsmat = getattr(sp, self.rhsformat + '_matrix') @@ -170,7 +168,7 @@ class SparseSolver(object): lead_info = [] for leadnum, interface in enumerate(sys.lead_interfaces): lead = sys.leads[leadnum] - if hasattr(lead, 'modes') and not force_realspace: + if not realspace: modes, (u, ulinv, nprop, svd_v) = lead.modes(energy, args=args) lead_info.append(modes) @@ -275,11 +273,10 @@ class SparseSolver(object): return LinearSys(lhs, rhs, indices, num_orb), lead_info - def solve(self, sys, energy=0, out_leads=None, in_leads=None, - force_realspace=False, check_hermiticity=True, - args=()): + def smatrix(self, sys, energy=0, out_leads=None, in_leads=None, + check_hermiticity=True, args=()): """ - Compute the scattering matrix or Green's function between leads. + Compute the scattering matrix of a system. Parameters ---------- @@ -296,10 +293,6 @@ class SparseSolver(object): Numbers of leads in which current or wave function is injected. None is interpreted as all leads. Default is ``None`` and means "all leads". - force_realspace : ``bool`` - Calculate Green's function between the outermost lead - sites, instead of lead modes. This is almost always - more computationally expensive and less stable. check_hermiticity : ``bool`` Check if the Hamiltonian matrices are Hermitian. args : tuple, defaults to empty @@ -307,34 +300,19 @@ class SparseSolver(object): Returns ------- - output : `~kwant.solvers.common.BlockResult` - See the notes below and `~kwant.solvers.common.BlockResult` + output : `~kwant.solvers.common.SMatrix` + See the notes below and `~kwant.solvers.common.SMatrix` documentation. Notes ----- This function can be used to calculate the conductance and other transport properties of a system. See the documentation for its output - type, `~kwant.solvers.common.BlockResult`. + type, `~kwant.solvers.common.SMatrix`. - It returns an object encapsulating the Green's function elements - between the desired leads. For leads defined as a self-energy, the - result is just the real-space retarded Green's function between the - system sites interfacing the leads in `in_leads` and those interfacing - the leads in `out_leads`. If, as is usually the case, the leads are - defined as tight-binding systems, then the Green's function from - incoming to outgoing modes is returned (more commonly known as the - scattering matrix). If some leads are defined via a self-energy, and - some as tight-binding systems, the result has Green's function's - elements between modes and sites. The returned object also contains - information about the modes or self-energies of the leads. - - If `force_realspace` is set to ``True`` all leads will be treated as if - they would be defined in terms of self energies. The returned Green's - function will be thus the retarded Green's function between sites in - real space, just like if the tradidional RGF algorithm would have been - used. Enabling this option is more computationally expensive and can - be less stable. + The returned object contains the scattering matrix elements from the + `in_leads` to the `out_leads` as well as information about the lead + modes. Both `in_leads` and `out_leads` must be sorted and may only contain unique entries. @@ -353,8 +331,7 @@ class SparseSolver(object): if len(in_leads) == 0 or len(out_leads) == 0: raise ValueError("No output is requested.") - linsys, lead_info = self._make_linear_sys(sys, in_leads, energy, - force_realspace, + linsys, lead_info = self._make_linear_sys(sys, in_leads, energy, False, check_hermiticity, args) kept_vars = np.concatenate([vars for i, vars in @@ -365,7 +342,7 @@ class SparseSolver(object): len_rhs = sum(i.shape[1] for i in linsys.rhs) len_kv = len(kept_vars) if not(len_rhs and len_kv): - return BlockResult(np.zeros((len_kv, len_rhs)), + return SMatrix(np.zeros((len_kv, len_rhs)), lead_info, out_leads, in_leads) # See comment about zero-shaped sparse matrices at the top of common.py. @@ -374,8 +351,90 @@ class SparseSolver(object): flhs = self._factorized(linsys.lhs) data = self._solve_linear_sys(flhs, rhs, kept_vars) - return BlockResult(data, lead_info, out_leads, in_leads) + return SMatrix(data, lead_info, out_leads, in_leads) + + def greens_function(self, sys, energy=0, out_leads=None, in_leads=None, + check_hermiticity=True, args=()): + """ + Compute the retarded Green's function of the system between its leads. + + Parameters + ---------- + sys : `kwant.system.FiniteSystem` + Low level system, containing the leads and the Hamiltonian of a + scattering region. + energy : number + Excitation energy at which to solve the scattering problem. + out_leads : sequence of integers or ``None`` + Numbers of leads where current or wave function is extracted. None + is interpreted as all leads. Default is ``None`` and means "all + leads". + in_leads : sequence of integers or ``None`` + Numbers of leads in which current or wave function is injected. + None is interpreted as all leads. Default is ``None`` and means + "all leads". + check_hermiticity : ``bool`` + Check if the Hamiltonian matrices are Hermitian. + args : tuple, defaults to empty + Positional arguments to pass to the ``hamiltonian`` method. + + Returns + ------- + output : `~kwant.solvers.common.GreensFunction` + See the notes below and `~kwant.solvers.common.GreensFunction` + documentation. + + Notes + ----- + This function can be used to calculate the conductance and other + transport properties of a system. It is often slower and less stable + than the scattering matrix-based calculation executed by + `~kwant.smatrix`, and is currently provided mostly for testing + purposes and compatibility with RGF code. + + It returns an object encapsulating the Green's function elements + between the system sites interfacing the leads in `in_leads` and those + interfacing the leads in `out_leads`. The returned object also + contains a list with self-energies of the leads. + + Both `in_leads` and `out_leads` must be sorted and may only contain + unique entries. + """ + + n = len(sys.lead_interfaces) + if in_leads is None: + in_leads = range(n) + if out_leads is None: + out_leads = range(n) + if sorted(in_leads) != in_leads or sorted(out_leads) != out_leads or \ + len(set(in_leads)) != len(in_leads) or \ + len(set(out_leads)) != len(out_leads): + raise ValueError("Lead lists must be sorted and " + "with unique entries.") + if len(in_leads) == 0 or len(out_leads) == 0: + raise ValueError("No output is requested.") + + linsys, lead_info = self._make_linear_sys(sys, in_leads, energy, True, + check_hermiticity, args) + + kept_vars = np.concatenate([vars for i, vars in + enumerate(linsys.indices) if i in + out_leads]) + + # Do not perform factorization if no calculation is to be done. + len_rhs = sum(i.shape[1] for i in linsys.rhs) + len_kv = len(kept_vars) + if not(len_rhs and len_kv): + return GreensFunction(np.zeros((len_kv, len_rhs)), + lead_info, out_leads, in_leads) + + # See comment about zero-shaped sparse matrices at the top of common.py. + rhs = sp.bmat([[i for i in linsys.rhs if i.shape[1]]], + format=self.rhsformat) + flhs = self._factorized(linsys.lhs) + data = self._solve_linear_sys(flhs, rhs, kept_vars) + return GreensFunction(data, lead_info, out_leads, in_leads) def ldos(self, fsys, energy=0, args=()): """ @@ -481,52 +540,16 @@ class WaveFunction(object): class BlockResult(object): """ - Solution of a transport problem, subblock of retarded Green's function - or scattering matrix. + Container for a linear system solution with variable grouping. - Transport properties can be easily accessed using the - `~BlockResult.transmission` method (don't be fooled by the name, - it can also compute reflection, which is just transmission from one - lead back into the same lead.) - - `BlockResult` however also allows for a more direct access to the result: - The data stored in `BlockResult` is either a real-space Green's - function (e.g. if ``force_realspace=True`` in - `~kwant.solvers.default.solve`) or a scattering matrix with respect to - lead modes. The details of this data can be directly accessed through - the instance variables `data` and `lead_info`. Subblocks of data - corresponding to particular leads are conveniently obtained by - `~BlockResult.submatrix`. - - Instance Variables - ------------------ - data : NumPy matrix - a matrix containing all the requested matrix elements of Green's - function. - lead_info : list of data - a list with output of `kwant.physics.modes` for each lead - defined as a builder, and self-energy for each lead defined as - self-energy term. - out_leads, in_leads : list of integers - indices of the leads where current is extracted (out) or injected - (in). Only those are listed for which BlockResult contains the - calculated result. + This class is not intended to be used directly. """ - - def __init__(self, data, lead_info, out_leads, in_leads): + def __init__(self, data, lead_info, out_leads, in_leads, sizes): self.data = data self.lead_info = lead_info self.out_leads = out_leads self.in_leads = in_leads - - sizes = [] - for i in self.lead_info: - if isinstance(i, physics.PropagatingModes): - sizes.append(len(i.momenta) // 2) - else: - sizes.append(i.shape[0]) self._sizes = np.array(sizes) - self._in_offsets = np.zeros(len(self.in_leads) + 1, int) self._in_offsets[1 :] = np.cumsum(self._sizes[self.in_leads]) self._out_offsets = np.zeros(len(self.out_leads) + 1, int) @@ -539,16 +562,15 @@ class BlockResult(object): return self.out_block_coords(lead_out), self.in_block_coords(lead_in) def out_block_coords(self, lead_out): - """Return a slice corresponding to the rows in the block corresponding - to lead_out + """Return a slice with the rows in the block corresponding to lead_out. """ lead_out = self.out_leads.index(lead_out) return slice(self._out_offsets[lead_out], self._out_offsets[lead_out + 1]) def in_block_coords(self, lead_in): - """Return a slice corresponding to the columns in the block - corresponding to lead_in + """ + Return a slice with the columns in the block corresponding to lead_in. """ lead_in = self.in_leads.index(lead_in) return slice(self._in_offsets[lead_in], @@ -558,44 +580,124 @@ class BlockResult(object): """Return the matrix elements from lead_in to lead_out.""" return self.data[self.block_coords(lead_out, lead_in)] + +class SMatrix(BlockResult): + """A scattering matrix. + + Transport properties can be easily accessed using the + `~SMatrix.transmission` method (don't be fooled by the name, + it can also compute reflection, which is just transmission from one + lead back into the same lead.) + + `SMatrix` however also allows for a more direct access to the result: The + data stored in `SMatrix` is a scattering matrix with respect to lead modes + and these modes themselves. The details of this data can be directly + accessed through the instance variables `data` and `lead_info`. Subblocks + of data corresponding to particular leads are conveniently obtained by + `~SMatrix.submatrix`. + + Instance Variables + ------------------ + data : NumPy array + a matrix containing all the requested matrix elements of the scattering + matrix. + lead_info : list of data + a list containing `kwant.physics.PropagatingModes` for each lead. + out_leads, in_leads : list of integers + indices of the leads where current is extracted (out) or injected + (in). Only those are listed for which SMatrix contains the + calculated result. + """ + + def __init__(self, data, lead_info, out_leads, in_leads): + sizes = [len(i.momenta) // 2 for i in lead_info] + super(SMatrix, self).__init__(data, lead_info, out_leads, in_leads, + sizes) + + def transmission(self, lead_out, lead_in): + """Return transmission from lead_in to lead_out.""" + return np.linalg.norm(self.submatrix(lead_out, lead_in)) ** 2 + + def __repr__(self): + return "SMatrix(data=%r, lead_info=%r, " \ + "out_leads=%r, in_leads=%r)" % (self.data, self.lead_info, + self.out_leads, self.in_leads) + + +class GreensFunction(BlockResult): + """ + Retarded Green's function. + + Transport properties can be easily accessed using the + `~GreensFunction.transmission` method (don't be fooled by the name, it can + also compute reflection, which is just transmission from one lead back into + the same lead). + + `GreensFunction` however also allows for a more direct access to the + result: The data stored in `GreensFunction` is the real-space Green's + function. The details of this data can be directly accessed through the + instance variables `data` and `lead_info`. Subblocks of data corresponding + to particular leads are conveniently obtained by + `~GreensFunction.submatrix`. + + Instance Variables + ------------------ + data : NumPy array + a matrix containing all the requested matrix elements of Green's + function. + lead_info : list of matrices + a list with self-energies of each lead. + out_leads, in_leads : list of integers + indices of the leads where current is extracted (out) or injected + (in). Only those are listed for which SMatrix contains the + calculated result. + """ + + def __init__(self, data, lead_info, out_leads, in_leads): + sizes = [i.shape[0] for i in lead_info] + super(GreensFunction, self).__init__(data, lead_info, out_leads, + in_leads, sizes) + def _a_ttdagger_a_inv(self, lead_out, lead_in): + """Return t * t^dagger in a certain basis.""" gf = self.submatrix(lead_out, lead_in) factors = [] for lead, gf2 in ((lead_out, gf), (lead_in, gf.conj().T)): possible_se = self.lead_info[lead] - if not isinstance(possible_se, physics.PropagatingModes): - # Lead is a "self energy lead": multiply gf2 with a gamma - # matrix. - factors.append(1j * (possible_se - possible_se.conj().T)) + factors.append(1j * (possible_se - possible_se.conj().T)) factors.append(gf2) return reduce(np.dot, factors) def transmission(self, lead_out, lead_in): """Return transmission from lead_in to lead_out.""" - if isinstance(self.lead_info[lead_out], physics.PropagatingModes) and \ - isinstance(self.lead_info[lead_in], physics.PropagatingModes): - return np.linalg.norm(self.submatrix(lead_out, lead_in)) ** 2 - else: - result = np.trace(self._a_ttdagger_a_inv(lead_out, lead_in)).real - if lead_out == lead_in: - # For reflection we have to be more careful - gamma = 1j * (self.lead_info[lead_in] - - self.lead_info[lead_in].conj().T) - gf = self.submatrix(lead_out, lead_in) - - # The number of channels is given by the number of - # nonzero eigenvalues of Gamma - # rationale behind the threshold from - # Golub; van Loan, chapter 5.5.8 - eps = np.finfo(gamma.dtype).eps * 1000 - N = np.sum(np.linalg.eigvalsh(gamma) > - eps * np.linalg.norm(gamma, np.inf)) - - result += 2 * np.trace(np.dot(gamma, gf)).imag + N - - return result + gf = self.submatrix(lead_out, lead_in) + factors = [] + for lead, gf2 in ((lead_out, gf), (lead_in, gf.conj().T)): + self_en = self.lead_info[lead] + factors.append(1j * (self_en - self_en.conj().T)) + factors.append(gf2) + attdagainv = reduce(np.dot, factors) + + result = np.trace(attdagainv).real + if lead_out == lead_in: + # For reflection we have to be more careful + gamma = 1j * (self.lead_info[lead_in] - + self.lead_info[lead_in].conj().T) + gf = self.submatrix(lead_out, lead_in) + + # The number of channels is given by the number of + # nonzero eigenvalues of Gamma + # rationale behind the threshold from + # Golub; van Loan, chapter 5.5.8 + eps = np.finfo(gamma.dtype).eps * 1000 + N = np.sum(np.linalg.eigvalsh(gamma) > + eps * np.linalg.norm(gamma, np.inf)) + + result += 2 * np.trace(np.dot(gamma, gf)).imag + N + + return result def __repr__(self): - return "BlockResult(data=%r, lead_info=%r, " \ + return "GreensFunction(data=%r, lead_info=%r, " \ "out_leads=%r, in_leads=%r)" % (self.data, self.lead_info, self.out_leads, self.in_leads) diff --git a/kwant/solvers/default.py b/kwant/solvers/default.py index e92fdbcaeaaef376a0cddddec9c461839b2b2bf3..56260d720da54fddb0cccb1f5c4b0059c27cd69f 100644 --- a/kwant/solvers/default.py +++ b/kwant/solvers/default.py @@ -6,7 +6,7 @@ # the AUTHORS file at the top-level directory of this distribution and at # http://kwant-project.org/authors. -__all__ = ['solve', 'ldos', 'wave_function'] +__all__ = ['smatrx', 'ldos', 'wave_function', 'greens_function'] # MUMPS usually works best. Use SciPy as fallback. try: @@ -16,6 +16,7 @@ except ImportError: hidden_instance = smodule.Solver() -solve = hidden_instance.solve +smatrix = hidden_instance.smatrix ldos = hidden_instance.ldos wave_function = hidden_instance.wave_function +greens_function = hidden_instance.greens_function diff --git a/kwant/solvers/mumps.py b/kwant/solvers/mumps.py index d99e615db58d39f1c82258ebf9add10c038fc7bc..ab491e4e5ca27e1eaf9cedaf072fcdd396f0dfa6 100644 --- a/kwant/solvers/mumps.py +++ b/kwant/solvers/mumps.py @@ -6,7 +6,8 @@ # the AUTHORS file at the top-level directory of this distribution and at # http://kwant-project.org/authors. -__all__ = ['solve', 'ldos', 'wave_function', 'options', 'Solver'] +__all__ = ['smatrix', 'ldos', 'wave_function', 'greens_function', 'options', + 'Solver'] import numpy as np from . import common @@ -122,7 +123,8 @@ class Solver(common.SparseSolver): default_solver = Solver() -solve = default_solver.solve +smatrix = default_solver.smatrix +greens_function = default_solver.greens_function ldos = default_solver.ldos wave_function = default_solver.wave_function options = default_solver.options diff --git a/kwant/solvers/sparse.py b/kwant/solvers/sparse.py index 3c7f157da48fe7fd38846cb20eb72b3d437b8c0c..b2712f38c17a2fe930e9025f805d1902a1557ffa 100644 --- a/kwant/solvers/sparse.py +++ b/kwant/solvers/sparse.py @@ -6,7 +6,7 @@ # the AUTHORS file at the top-level directory of this distribution and at # http://kwant-project.org/authors. -__all__ = ['solve', 'ldos', 'wave_function', 'Solver'] +__all__ = ['smatrix', 'greens_function', 'ldos', 'wave_function', 'Solver'] import warnings import numpy as np @@ -115,6 +115,7 @@ class Solver(common.SparseSolver): default_solver = Solver() -solve = default_solver.solve +smatrix = default_solver.smatrix +greens_function = default_solver.greens_function ldos = default_solver.ldos wave_function = default_solver.wave_function diff --git a/kwant/solvers/tests/_test_sparse.py b/kwant/solvers/tests/_test_sparse.py index 9f8b2d839699132da3d36049c8ae6ddecf55ded2..1a3e49af250d22539a98fb2cd6c1df4042d550a3 100644 --- a/kwant/solvers/tests/_test_sparse.py +++ b/kwant/solvers/tests/_test_sparse.py @@ -281,7 +281,7 @@ def test_tricky_singular_hopping(solve): # Test equivalence between self-energy and scattering matrix representations. # Also check that transmission works. -def test_selfenergy(solve): +def test_selfenergy(greens_function, smatrix): np.random.seed(4) system = kwant.Builder() left_lead = kwant.Builder(kwant.TranslationalSymmetry((-1,))) @@ -299,12 +299,12 @@ def test_selfenergy(solve): system.attach_lead(right_lead) fsys = system.finalized() - t = solve(fsys, 0, [1], [0]).data + t = smatrix(fsys, 0, [1], [0]).data eig_should_be = np.linalg.eigvals(t * t.conjugate().transpose()) n_eig = len(eig_should_be) fsys.leads[1] = LeadWithOnlySelfEnergy(fsys.leads[1]) - sol = solve(fsys, 0, [1], [0]) + sol = greens_function(fsys, 0, [1], [0]) ttdagnew = sol._a_ttdagger_a_inv(1, 0) eig_are = np.linalg.eigvals(ttdagnew) t_should_be = np.sum(eig_are) @@ -314,7 +314,7 @@ def test_selfenergy(solve): assert_almost_equal(t_should_be, sol.transmission(1, 0)) fsys.leads[0] = LeadWithOnlySelfEnergy(fsys.leads[0]) - sol = solve(fsys, 0, [1], [0]) + sol = greens_function(fsys, 0, [1], [0]) ttdagnew = sol._a_ttdagger_a_inv(1, 0) eig_are = np.linalg.eigvals(ttdagnew) t_should_be = np.sum(eig_are) diff --git a/kwant/solvers/tests/test_mumps.py b/kwant/solvers/tests/test_mumps.py index 3c78eb7cd86c26a562c83a180f420877e1ca79e3..9802e1fee996c4e1d781dbe882370eba71cd36c2 100644 --- a/kwant/solvers/tests/test_mumps.py +++ b/kwant/solvers/tests/test_mumps.py @@ -10,7 +10,7 @@ from numpy.testing.decorators import skipif try: from kwant.solvers.mumps import \ - solve, ldos, wave_function, options, reset_options + smatrix, greens_function, ldos, wave_function, options, reset_options import _test_sparse _no_mumps = False except ImportError: @@ -30,7 +30,7 @@ def test_output(): for opts in opt_list: reset_options() options(**opts) - _test_sparse.test_output(solve) + _test_sparse.test_output(smatrix) @skipif(_no_mumps) @@ -38,7 +38,7 @@ def test_one_lead(): for opts in opt_list: reset_options() options(**opts) - _test_sparse.test_one_lead(solve) + _test_sparse.test_one_lead(smatrix) @skipif(_no_mumps) @@ -46,7 +46,7 @@ def test_smatrix_shape(): for opts in opt_list: reset_options() options(**opts) - _test_sparse.test_smatrix_shape(solve) + _test_sparse.test_smatrix_shape(smatrix) @skipif(_no_mumps) @@ -54,14 +54,14 @@ def test_two_equal_leads(): for opts in opt_list: reset_options() options(**opts) - _test_sparse.test_two_equal_leads(solve) + _test_sparse.test_two_equal_leads(smatrix) @skipif(_no_mumps) def test_graph_system(): for opts in opt_list: reset_options() options(**opts) - _test_sparse.test_graph_system(solve) + _test_sparse.test_graph_system(smatrix) @skipif(_no_mumps) @@ -69,7 +69,7 @@ def test_singular_graph_system(): for opts in opt_list: reset_options() options(**opts) - _test_sparse.test_singular_graph_system(solve) + _test_sparse.test_singular_graph_system(smatrix) @skipif(_no_mumps) @@ -77,7 +77,7 @@ def test_tricky_singular_hopping(): for opts in opt_list: reset_options() options(**opts) - _test_sparse.test_tricky_singular_hopping(solve) + _test_sparse.test_tricky_singular_hopping(smatrix) @skipif(_no_mumps) @@ -85,7 +85,7 @@ def test_selfenergy(): for opts in opt_list: reset_options() options(**opts) - _test_sparse.test_selfenergy(solve) + _test_sparse.test_selfenergy(greens_function, smatrix) @skipif(_no_mumps) @@ -93,7 +93,7 @@ def test_selfenergy_reflection(): for opts in opt_list: reset_options() options(**opts) - _test_sparse.test_selfenergy_reflection(solve) + _test_sparse.test_selfenergy_reflection(greens_function) @skipif(_no_mumps) @@ -101,7 +101,7 @@ def test_very_singular_leads(): for opts in opt_list: reset_options() options(**opts) - _test_sparse.test_very_singular_leads(solve) + _test_sparse.test_very_singular_leads(smatrix) @skipif(_no_mumps) diff --git a/kwant/solvers/tests/test_sparse.py b/kwant/solvers/tests/test_sparse.py index 8b0a2375e36218843779ba285f891067879be561..399a1a9cd334406b029a899a36660bf21348add4 100644 --- a/kwant/solvers/tests/test_sparse.py +++ b/kwant/solvers/tests/test_sparse.py @@ -7,48 +7,48 @@ # http://kwant-project.org/authors. from nose.plugins.skip import Skip, SkipTest -from kwant.solvers.sparse import solve, ldos, wave_function +from kwant.solvers.sparse import smatrix, greens_function, ldos, wave_function import kwant.solvers.sparse import _test_sparse def test_output(): - _test_sparse.test_output(solve) + _test_sparse.test_output(smatrix) def test_one_lead(): - _test_sparse.test_one_lead(solve) + _test_sparse.test_one_lead(smatrix) def test_smatrix_shape(): - _test_sparse.test_smatrix_shape(solve) + _test_sparse.test_smatrix_shape(smatrix) def test_two_equal_leads(): - _test_sparse.test_two_equal_leads(solve) + _test_sparse.test_two_equal_leads(smatrix) def test_graph_system(): - _test_sparse.test_graph_system(solve) + _test_sparse.test_graph_system(smatrix) def test_singular_graph_system(): - _test_sparse.test_singular_graph_system(solve) + _test_sparse.test_singular_graph_system(smatrix) def test_tricky_singular_hopping(): - _test_sparse.test_tricky_singular_hopping(solve) + _test_sparse.test_tricky_singular_hopping(smatrix) def test_selfenergy(): - _test_sparse.test_selfenergy(solve) + _test_sparse.test_selfenergy(greens_function, smatrix) def test_selfenergy_reflection(): - _test_sparse.test_selfenergy_reflection(solve) + _test_sparse.test_selfenergy_reflection(greens_function) def test_very_singular_leads(): - _test_sparse.test_very_singular_leads(solve) + _test_sparse.test_very_singular_leads(smatrix) def test_umfpack_del(): diff --git a/kwant/tests/test_comprehensive.py b/kwant/tests/test_comprehensive.py index 0bb59a05b14b05c64de540dbcacbc4b1549d0bd8..4b0086f75ead185004c8a25113253cc53f1d011c 100644 --- a/kwant/tests/test_comprehensive.py +++ b/kwant/tests/test_comprehensive.py @@ -41,7 +41,7 @@ def test_qhe(W=16, L=8): # reciprocal_phis = numpy.linspace(0.1, 7, 200) # conductances = [] # for phi in 1 / reciprocal_phis: - # smatrix = kwant.solve(sys, 1.0, args=[phi, ""]) + # smatrix = kwant.smatrix(sys, 1.0, args=[phi, ""]) # conductances.append(smatrix.transmission(1, 0)) # pyplot.plot(reciprocal_phis, conductances) # pyplot.show() @@ -51,5 +51,5 @@ def test_qhe(W=16, L=8): ((5.2, 5.5), 3, 1e-1)]: for r_phi in r_phis: args = (1.0 / r_phi, "") - T_actual = kwant.solve(sys, 1.0, args=args).transmission(1, 0) + T_actual = kwant.smatrix(sys, 1.0, args=args).transmission(1, 0) assert abs(T_nominal - T_actual) < max_err diff --git a/kwant/tests/test_system.py b/kwant/tests/test_system.py index 898cdf69aa34893fa25bfe5f4e2c55764389b192..8359b9ab7f9b00d78d8c0820125e77f0c65d1e7d 100644 --- a/kwant/tests/test_system.py +++ b/kwant/tests/test_system.py @@ -63,11 +63,12 @@ def test_hamiltonian_submatrix(): lead[gr(0), gr(1)] = np.random.randn(2, 2) sys.attach_lead(lead) sys2 = sys.finalized() - smatrix = kwant.solve(sys2, .1).data + smatrix = kwant.smatrix(sys2, .1).data sys3 = sys2.precalculate(.1, calculate_selfenergy=False) - smatrix2 = kwant.solve(sys3, .1).data + smatrix2 = kwant.smatrix(sys3, .1).data np.testing.assert_almost_equal(smatrix, smatrix2) - assert_raises(ValueError, kwant.solve, sys3, 0.2, None, None, True) + assert_raises(ValueError, kwant.solvers.default.greens_function, sys3, 0.2, + None, None) # Test for shape errors. sys[gr(0), gr(2)] = np.array([[1, 2]])