diff --git a/doc/source/tutorial/closed_system.py b/doc/source/tutorial/closed_system.py index 8aebe582ae70fc2c02a0b0a243dc95a35b8ce590..abe9114bca37624cb6893e5359dfb6159dad843e 100644 --- a/doc/source/tutorial/closed_system.py +++ b/doc/source/tutorial/closed_system.py @@ -63,7 +63,7 @@ def plot_spectrum(sys, Bfields): energies = [] for B in Bfields: # Obtain the Hamiltonian as a dense matrix - ham_mat = sys.hamiltonian_submatrix(sparse=True, args=[B]) + ham_mat = sys.hamiltonian_submatrix(args=[B], sparse=True) # we only calculate the 15 lowest eigenvalues ev = sla.eigsh(ham_mat, k=15, which='SM', return_eigenvectors=False) diff --git a/examples/advanced_construction.py b/examples/advanced_construction.py index 2165b4047c7dea89c6d6253c050b88d934316f56..dcfcfa92e6291fc577cd529732310b26d7b79c00 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.smatrix(sys, 1.1, args=[0.1]).transmission(0, 1) + print kwant.smatrix(sys, 1.1, [0.1]).transmission(0, 1) if __name__ == '__main__': diff --git a/kwant/_system.pyx b/kwant/_system.pyx index 1975cb038baf97076a80b5b40134a99b476df812..ac82fd8e32705d0840d0f9baaed73be4308fb77d 100644 --- a/kwant/_system.pyx +++ b/kwant/_system.pyx @@ -226,21 +226,20 @@ def make_dense_full(ham, args, CGraph gr, diag, @cython.embedsignature(True) -def hamiltonian_submatrix(self, to_sites=None, from_sites=None, - sparse=False, return_norb=False, - args=()): +def hamiltonian_submatrix(self, args=(), to_sites=None, from_sites=None, + sparse=False, return_norb=False): """Return a submatrix of the system Hamiltonian. Parameters ---------- + args : tuple, defaults to empty + Positional arguments to pass to the ``hamiltonian`` method. to_sites : sequence of sites or None (default) from_sites : sequence of sites or None (default) sparse : bool Whether to return a sparse or a dense matrix. Defaults to `False`. return_norb : bool Whether to return arrays of numbers of orbitals. Defaults to `False`. - args : tuple, defaults to empty - Positional arguments to pass to the ``hamiltonian`` method. Returns ------- diff --git a/kwant/physics/dispersion.py b/kwant/physics/dispersion.py index 6ce27149dbdcd17a139f527cfdbd42259062d8be..1c8bb2042dafc62616a5e47a4f1b64842d58042c 100644 --- a/kwant/physics/dispersion.py +++ b/kwant/physics/dispersion.py @@ -40,10 +40,10 @@ class Bands(object): """ def __init__(self, sys, args=()): - self.ham = sys.cell_hamiltonian(args=args) + self.ham = sys.cell_hamiltonian(args) if not np.allclose(self.ham, self.ham.T.conj()): raise ValueError('The cell Hamiltonian is not Hermitian.') - hop = sys.inter_cell_hopping(args=args) + hop = sys.inter_cell_hopping(args) self.hop = np.empty(self.ham.shape, dtype=complex) self.hop[:, : hop.shape[1]] = hop self.hop[:, hop.shape[1]:] = 0 diff --git a/kwant/plotter.py b/kwant/plotter.py index ff88c64e2028847f307fac96f456a0cc690b3653..036a0c6ed9f1c28b630238049faa2ecf4ae24e29 100644 --- a/kwant/plotter.py +++ b/kwant/plotter.py @@ -1533,7 +1533,7 @@ def map(sys, value, colorbar=True, cmap=None, vmin=None, vmax=None, a=None, return output_fig(fig, file=file, show=show) -def bands(sys, momenta=65, args=(), file=None, show=True, dpi=None, +def bands(sys, args=(), momenta=65, file=None, show=True, dpi=None, fig_size=None, ax=None): """Plot band structure of a translationally invariant 1D system. @@ -1580,7 +1580,7 @@ def bands(sys, momenta=65, args=(), file=None, show=True, dpi=None, if momenta.ndim != 1: momenta = np.linspace(-np.pi, np.pi, momenta) - bands = physics.Bands(sys, args=args) + bands = physics.Bands(sys, args) energies = [bands(k) for k in momenta] if ax is None: diff --git a/kwant/solvers/common.py b/kwant/solvers/common.py index 0ab39f87b04e279c3ba1eb867f014d8f318bc1a5..b628b5af2d7689100951d222d46153ae8b58e53e 100644 --- a/kwant/solvers/common.py +++ b/kwant/solvers/common.py @@ -90,8 +90,8 @@ class SparseSolver(object): """ pass - def _make_linear_sys(self, sys, in_leads, energy=0, realspace=False, - check_hermiticity=True, args=()): + def _make_linear_sys(self, sys, in_leads, energy=0, args=(), + realspace=False, check_hermiticity=True): """Make a sparse linear system of equations defining a scattering problem. @@ -143,9 +143,8 @@ class SparseSolver(object): if not sys.lead_interfaces: raise ValueError('System contains no leads.') - lhs, norb = sys.hamiltonian_submatrix(sparse=True, - return_norb=True, - args=args)[:2] + lhs, norb = sys.hamiltonian_submatrix(args, sparse=True, + return_norb=True)[:2] lhs = getattr(lhs, 'to' + self.lhsformat)() lhs = lhs - energy * sp.identity(lhs.shape[0], format=self.lhsformat) num_orb = lhs.shape[0] @@ -169,7 +168,7 @@ class SparseSolver(object): for leadnum, interface in enumerate(sys.lead_interfaces): lead = sys.leads[leadnum] if not realspace: - prop, stab = lead.modes(energy, args=args) + prop, stab = lead.modes(energy, args) lead_info.append(prop) u, ulinv, nprop, svd_v = \ stab.vecs, stab.vecslmbdainv, stab.nmodes, stab.sqrt_hop @@ -275,8 +274,8 @@ class SparseSolver(object): return LinearSys(lhs, rhs, indices, num_orb), lead_info - def smatrix(self, sys, energy=0, out_leads=None, in_leads=None, - check_hermiticity=True, args=()): + def smatrix(self, sys, energy=0, args=(), + out_leads=None, in_leads=None, check_hermiticity=True): """ Compute the scattering matrix of a system. @@ -287,6 +286,8 @@ class SparseSolver(object): scattering region. energy : number Excitation energy at which to solve the scattering problem. + args : tuple, defaults to empty + Positional arguments to pass to the ``hamiltonian`` method. 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 @@ -297,8 +298,6 @@ class SparseSolver(object): "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 ------- @@ -333,8 +332,8 @@ 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, False, - check_hermiticity, args) + linsys, lead_info = self._make_linear_sys(sys, in_leads, energy, args, + False, check_hermiticity) kept_vars = np.concatenate([vars for i, vars in enumerate(linsys.indices) if i in @@ -355,8 +354,8 @@ class SparseSolver(object): 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=()): + def greens_function(self, sys, energy=0, args=(), + out_leads=None, in_leads=None, check_hermiticity=True): """ Compute the retarded Green's function of the system between its leads. @@ -367,6 +366,8 @@ class SparseSolver(object): scattering region. energy : number Excitation energy at which to solve the scattering problem. + args : tuple, defaults to empty + Positional arguments to pass to the ``hamiltonian`` method. 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 @@ -377,8 +378,6 @@ class SparseSolver(object): "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 ------- @@ -416,8 +415,8 @@ 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, True, - check_hermiticity, args) + linsys, lead_info = self._make_linear_sys(sys, in_leads, energy, args, + True, check_hermiticity) kept_vars = np.concatenate([vars for i, vars in enumerate(linsys.indices) if i in @@ -465,8 +464,7 @@ class SparseSolver(object): "tight binding systems.") linsys, lead_info = \ - self._make_linear_sys(fsys, xrange(len(fsys.leads)), energy, - args=args) + self._make_linear_sys(fsys, xrange(len(fsys.leads)), energy, args) ldos = np.zeros(linsys.num_orb, float) factored = None @@ -527,8 +525,7 @@ class WaveFunction(object): ' are not available yet.' raise NotImplementedError(msg) linsys, lead_info = \ - solver._make_linear_sys(sys, xrange(len(sys.leads)), energy, - args=args) + solver._make_linear_sys(sys, xrange(len(sys.leads)), energy, args) self.solve = solver._solve_linear_sys self.rhs = linsys.rhs self.factorized_h = solver._factorized(linsys.lhs) diff --git a/kwant/solvers/tests/_test_sparse.py b/kwant/solvers/tests/_test_sparse.py index 1a3e49af250d22539a98fb2cd6c1df4042d550a3..1d14883f33c8e1e2a1af3c7699b0e12895b92726 100644 --- a/kwant/solvers/tests/_test_sparse.py +++ b/kwant/solvers/tests/_test_sparse.py @@ -37,7 +37,7 @@ def assert_modes_equal(modes1, modes2): # Test output sanity: that an error is raised if no output is requested, # and that solving for a subblock of a scattering matrix is the same as taking # a subblock of the full scattering matrix. -def test_output(solve): +def test_output(smatrix): np.random.seed(3) system = kwant.Builder() left_lead = kwant.Builder(kwant.TranslationalSymmetry((-1,))) @@ -55,19 +55,19 @@ def test_output(solve): system.attach_lead(right_lead) fsys = system.finalized() - result1 = solve(fsys) + result1 = smatrix(fsys) s, modes1 = result1.data, result1.lead_info assert s.shape == 2 * (sum(len(i.momenta) for i in modes1) // 2,) s1 = result1.submatrix(1, 0) - result2 = solve(fsys, 0, [1], [0]) + result2 = smatrix(fsys, 0, (), [1], [0]) s2, modes2 = result2.data, result2.lead_info assert s2.shape == (len(modes2[1].momenta) // 2, len(modes2[0].momenta) // 2) assert_almost_equal(s1, s2) assert_almost_equal(np.dot(s.T.conj(), s), np.identity(s.shape[0])) - assert_raises(ValueError, solve, fsys, 0, []) - modes = solve(fsys).lead_info + assert_raises(ValueError, smatrix, fsys, out_leads=[]) + modes = smatrix(fsys).lead_info h = fsys.leads[0].cell_hamiltonian() t = fsys.leads[0].inter_cell_hopping() modes1 = kwant.physics.modes(h, t)[0] @@ -79,7 +79,7 @@ def test_output(solve): # Test that a system with one lead has unitary scattering matrix. -def test_one_lead(solve): +def test_one_lead(smatrix): np.random.seed(3) system = kwant.Builder() lead = kwant.Builder(kwant.TranslationalSymmetry((-1,))) @@ -95,14 +95,14 @@ def test_one_lead(solve): system.attach_lead(lead) fsys = system.finalized() - s = solve(fsys).data + s = smatrix(fsys).data assert_almost_equal(np.dot(s.conjugate().transpose(), s), np.identity(s.shape[0])) # Test that a system with one lead with no propagating modes has a # 0x0 S-matrix. -def test_smatrix_shape(solve): +def test_smatrix_shape(smatrix): system = kwant.Builder() lead0 = kwant.Builder(kwant.TranslationalSymmetry((-1,))) lead1 = kwant.Builder(kwant.TranslationalSymmetry((1,))) @@ -123,30 +123,30 @@ def test_smatrix_shape(solve): lead0_val = 4 lead1_val = 4 - s = solve(fsys, energy=1.0, out_leads=[1], in_leads=[0]).data + s = smatrix(fsys, 1.0, (), [1], [0]).data assert s.shape == (0, 0) lead0_val = 2 lead1_val = 2 - s = solve(fsys, energy=1.0, out_leads=[1], in_leads=[0]).data + s = smatrix(fsys, 1.0, (), [1], [0]).data assert s.shape == (1, 1) lead0_val = 4 lead1_val = 2 - s = solve(fsys, energy=1.0, out_leads=[1], in_leads=[0]).data + s = smatrix(fsys, 1.0, (), [1], [0]).data assert s.shape == (1, 0) lead0_val = 2 lead1_val = 4 - s = solve(fsys, energy=1.0, out_leads=[1], in_leads=[0]).data + s = smatrix(fsys, 1.0, (), [1], [0]).data assert s.shape == (0, 1) # Test that a translationally invariant system with two leads has only # transmission and that transmission does not mix modes. -def test_two_equal_leads(solve): +def test_two_equal_leads(smatrix): def check_fsys(): - sol = solve(fsys) + sol = smatrix(fsys) s, leads = sol.data, sol.lead_info assert_almost_equal(np.dot(s.conjugate().transpose(), s), np.identity(s.shape[0])) @@ -183,7 +183,7 @@ def test_two_equal_leads(solve): # Test a more complicated graph with non-singular hopping. -def test_graph_system(solve): +def test_graph_system(smatrix): np.random.seed(11) system = kwant.Builder() lead = kwant.Builder(kwant.TranslationalSymmetry((-1, 0))) @@ -202,7 +202,7 @@ def test_graph_system(solve): system.attach_lead(lead.reversed()) fsys = system.finalized() - result = solve(fsys) + result = smatrix(fsys) s, leads = result.data, result.lead_info assert_almost_equal(np.dot(s.conjugate().transpose(), s), np.identity(s.shape[0])) @@ -216,7 +216,7 @@ def test_graph_system(solve): # Test a system with singular hopping. -def test_singular_graph_system(solve): +def test_singular_graph_system(smatrix): np.random.seed(11) system = kwant.Builder() @@ -235,7 +235,7 @@ def test_singular_graph_system(solve): system.attach_lead(lead.reversed()) fsys = system.finalized() - result = solve(fsys) + result = smatrix(fsys) s, leads = result.data, result.lead_info assert_almost_equal(np.dot(s.conjugate().transpose(), s), np.identity(s.shape[0])) @@ -251,7 +251,7 @@ def test_singular_graph_system(solve): # This test features inside the cell Hamiltonian a hopping matrix with more # zero eigenvalues than the lead hopping matrix. Older version of the # sparse solver failed here. -def test_tricky_singular_hopping(solve): +def test_tricky_singular_hopping(smatrix): system = kwant.Builder() lead = kwant.Builder(kwant.TranslationalSymmetry((4, 0))) @@ -274,7 +274,7 @@ def test_tricky_singular_hopping(solve): system.leads.append(kwant.builder.BuilderLead(lead, interface)) fsys = system.finalized() - s = solve(fsys, -1.3).data + s = smatrix(fsys, -1.3).data assert_almost_equal(np.dot(s.conjugate().transpose(), s), np.identity(s.shape[0])) @@ -299,12 +299,12 @@ def test_selfenergy(greens_function, smatrix): system.attach_lead(right_lead) fsys = system.finalized() - t = smatrix(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 = greens_function(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(greens_function, smatrix): assert_almost_equal(t_should_be, sol.transmission(1, 0)) fsys.leads[0] = LeadWithOnlySelfEnergy(fsys.leads[0]) - sol = greens_function(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) @@ -324,7 +324,7 @@ def test_selfenergy(greens_function, smatrix): assert_almost_equal(t_should_be, sol.transmission(1, 0)) -def test_selfenergy_reflection(solve): +def test_selfenergy_reflection(smatrix): np.random.seed(4) system = kwant.Builder() left_lead = kwant.Builder(kwant.TranslationalSymmetry((-1,))) @@ -339,15 +339,15 @@ def test_selfenergy_reflection(solve): system.attach_lead(left_lead) fsys = system.finalized() - t = solve(fsys, 0, [0], [0]) + t = smatrix(fsys, 0, (), [0], [0]) fsys.leads[0] = LeadWithOnlySelfEnergy(fsys.leads[0]) - sol = solve(fsys, 0, [0], [0]) + sol = smatrix(fsys, 0, (), [0], [0]) assert_almost_equal(sol.transmission(0,0), t.transmission(0,0)) -def test_very_singular_leads(solve): +def test_very_singular_leads(smatrix): sys = kwant.Builder() gr = kwant.lattice.chain() left_lead = kwant.Builder(kwant.TranslationalSymmetry((-1,))) @@ -358,7 +358,7 @@ def test_very_singular_leads(solve): sys.attach_lead(left_lead) sys.attach_lead(right_lead) fsys = sys.finalized() - leads = solve(fsys).lead_info + leads = smatrix(fsys).lead_info assert [len(i.momenta) for i in leads] == [0, 4] diff --git a/kwant/system.py b/kwant/system.py index 2a69650a9439cd56802cb93a7772b3ce585145c9..545739d926a700fdcf9b19cc0cd9244812d1d969 100644 --- a/kwant/system.py +++ b/kwant/system.py @@ -180,18 +180,18 @@ class InfiniteSystem(System): """ __metaclass__ = abc.ABCMeta - def cell_hamiltonian(self, sparse=False, args=()): + def cell_hamiltonian(self, args=(), sparse=False): """Hamiltonian of a single cell of the infinite system.""" cell_sites = xrange(self.cell_size) - return self.hamiltonian_submatrix(cell_sites, cell_sites, - sparse=sparse, args=args) + return self.hamiltonian_submatrix(args, cell_sites, cell_sites, + sparse=sparse) - def inter_cell_hopping(self, sparse=False, args=()): + def inter_cell_hopping(self, args=(), sparse=False): """Hopping Hamiltonian between two cells of the infinite system.""" cell_sites = xrange(self.cell_size) neighbor_sites = xrange(self.cell_size, self.graph.num_nodes) - return self.hamiltonian_submatrix(cell_sites, neighbor_sites, - sparse=sparse, args=args) + return self.hamiltonian_submatrix(args, cell_sites, neighbor_sites, + sparse=sparse) def modes(self, energy=0, args=()): """Return mode decomposition of the lead @@ -199,13 +199,13 @@ class InfiniteSystem(System): See documentation of `~kwant.physics.PropagatingModes` and `~kwant.physics.StabilizedModes` for the return format details. """ - ham = self.cell_hamiltonian(args=args) + ham = self.cell_hamiltonian(args) shape = ham.shape assert len(shape) == 2 assert shape[0] == shape[1] # Subtract energy from the diagonal. ham.flat[::ham.shape[0] + 1] -= energy - return physics.modes(ham, self.inter_cell_hopping(args=args)) + return physics.modes(ham, self.inter_cell_hopping(args)) def selfenergy(self, energy=0, args=()): """Return self-energy of a lead. @@ -214,13 +214,13 @@ class InfiniteSystem(System): ``sum(len(self.hamiltonian(i, i)) for i in range(self.graph.num_nodes - self.cell_size))``. """ - ham = self.cell_hamiltonian(args=args) + ham = self.cell_hamiltonian(args) shape = ham.shape assert len(shape) == 2 assert shape[0] == shape[1] # Subtract energy from the diagonal. ham.flat[::ham.shape[0] + 1] -= energy - return physics.selfenergy(ham, self.inter_cell_hopping(args=args)) + return physics.selfenergy(ham, self.inter_cell_hopping(args)) class PrecalculatedLead(object): diff --git a/kwant/tests/test_comprehensive.py b/kwant/tests/test_comprehensive.py index 4b0086f75ead185004c8a25113253cc53f1d011c..158fffe0935a4641976738e5705acde3812e81ab 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.smatrix(sys, 1.0, args=[phi, ""]) + # smatrix = kwant.smatrix(sys, 1.0, [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.smatrix(sys, 1.0, args=args).transmission(1, 0) + T_actual = kwant.smatrix(sys, 1.0, 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 c541957b2fef739de453bb37b06003ed336c953c..2f3d6bb5f5d9bc78d8473007123d6200a68f8fff 100644 --- a/kwant/tests/test_system.py +++ b/kwant/tests/test_system.py @@ -37,10 +37,10 @@ def test_hamiltonian_submatrix(): mat = mat[:, perm] np.testing.assert_array_equal(mat, mat_should_be) - mat = sys2.hamiltonian_submatrix(perm[[0, 1]], perm[[2]]) + mat = sys2.hamiltonian_submatrix((), perm[[0, 1]], perm[[2]]) np.testing.assert_array_equal(mat, mat_should_be[:2, 2:3]) - mat = sys2.hamiltonian_submatrix(perm[[0, 1]], perm[[2]], sparse=True) + mat = sys2.hamiltonian_submatrix((), perm[[0, 1]], perm[[2]], sparse=True) mat = mat.todense() np.testing.assert_array_equal(mat, mat_should_be[:2, 2:3]) @@ -67,8 +67,7 @@ def test_hamiltonian_submatrix(): sys3 = sys2.precalculate(.1, calculate_selfenergy=False) smatrix2 = kwant.smatrix(sys3, .1).data np.testing.assert_almost_equal(smatrix, smatrix2) - assert_raises(ValueError, kwant.solvers.default.greens_function, sys3, 0.2, - None, None) + assert_raises(ValueError, kwant.solvers.default.greens_function, sys3, 0.2) # Test for shape errors. sys[gr(0), gr(2)] = np.array([[1, 2]]) @@ -91,7 +90,7 @@ def test_hamiltonian_submatrix(): sys[(gr(i) for i in xrange(3))] = onsite sys[((gr(i), gr(i + 1)) for i in xrange(2))] = hopping sys2 = sys.finalized() - mat = sys2.hamiltonian_submatrix(args=(2, 1)) + mat = sys2.hamiltonian_submatrix((2, 1)) mat_should_be = [[5, 1, 0], [1, 4, 1.], [0, 1, 3]] # Sorting is required due to unknown compression order of builder.