allow systems to accept named parameters

Hamiltonian value functions may now have signatures that depend on
different parameters. On finalization, Builders now inspect all of their
value functions and store the names of the parameters on which the value
function depends. Then, when evaluating the Hamiltonian, a particular
value function is only passed the parameters on which it depends.

'System.hamiltonian' has been updated to accept an extra keyword-only
parameter, 'params', which is a dictionary mapping parameter names to
values. This is mutually exclusive with the existing '*args'. All
top-level API that takes an 'args' parameter now also takes a
keyword-only parameter, 'params'

kwant/_common.py

@@ -10,6 +10,7 @@ import subprocess
 import os
 import numpy as np
 import numbers
+import inspect
 
 __all__ = ['version', 'KwantDeprecationWarning', 'UserCodeError']
 
@@ -130,3 +131,26 @@ def ensure_rng(rng=None):
         return rng
     raise ValueError("Expecting a seed or an object that offers the "
                      "numpy.random.RandomState interface.")
+
+
+def get_parameters(func):
+    """Get the names of the parameters to 'func' and whether it takes kwargs.
+
+    Returns
+    -------
+    names : list
+        Positional, keyword and keyword only parameter names in the order
+        that they appear in the signature of 'func'.
+    takes_kwargs : bool
+        True if 'func' takes '**kwargs'.
+    """
+    sig = inspect.signature(func)
+    pars = sig.parameters
+
+    # Signature.parameters is an *ordered mapping*
+    names = [k for (k, v) in pars.items()
+             if v.kind in (inspect.Parameter.POSITIONAL_OR_KEYWORD,
+                           inspect.Parameter.KEYWORD_ONLY)]
+    takes_kwargs = any(i.kind is inspect.Parameter.VAR_KEYWORD
+                       for i in pars.values())
+    return names, takes_kwargs

kwant/_system.pyx

@@ -21,7 +21,7 @@ msg = ('Hopping from site {0} to site {1} does not match the '
        'dimensions of onsite Hamiltonians of these sites.')
 
 @cython.boundscheck(False)
-def make_sparse(ham, args, CGraph gr, diag,
+def make_sparse(ham, args, params, CGraph gr, diag,
                 gint [:] from_sites, n_by_to_site,
                 gint [:] to_norb, gint [:] to_off,
                 gint [:] from_norb, gint [:] from_off):
@@ -73,7 +73,7 @@ def make_sparse(ham, args, CGraph gr, diag,
             if ts not in n_by_to_site:
                 continue
             n_ts = n_by_to_site[ts]
-            h = matrix(ham(ts, fs, *args), complex)
+            h = matrix(ham(ts, fs, *args, params=params), complex)
             if h.shape[0] != to_norb[n_ts] or h.shape[1] != from_norb[n_fs]:
                 raise ValueError(msg.format(fs, ts))
             for i in range(h.shape[0]):
@@ -98,7 +98,7 @@ def make_sparse(ham, args, CGraph gr, diag,
 
 
 @cython.boundscheck(False)
-def make_sparse_full(ham, args, CGraph gr, diag,
+def make_sparse_full(ham, args, params, CGraph gr, diag,
                      gint [:] to_norb, gint [:] to_off,
                      gint [:] from_norb, gint [:] from_off):
     """For internal use by hamiltonian_submatrix."""
@@ -143,7 +143,7 @@ def make_sparse_full(ham, args, CGraph gr, diag,
         for ts in nbors.data[:nbors.size]:
             if ts < fs:
                 continue
-            h = matrix(ham(ts, fs, *args), complex)
+            h = matrix(ham(ts, fs, *args, params=params), complex)
             if h.shape[0] != to_norb[ts] or h.shape[1] != from_norb[fs]:
                 raise ValueError(msg.format(fs, ts))
             for i in range(h.shape[0]):
@@ -168,7 +168,7 @@ def make_sparse_full(ham, args, CGraph gr, diag,
 
 
 @cython.boundscheck(False)
-def make_dense(ham, args, CGraph gr, diag,
+def make_dense(ham, args, params, CGraph gr, diag,
                gint [:] from_sites, n_by_to_site,
                gint [:] to_norb, gint [:] to_off,
                gint [:] from_norb, gint [:] from_off):
@@ -197,7 +197,7 @@ def make_dense(ham, args, CGraph gr, diag,
             if ts not in n_by_to_site:
                 continue
             n_ts = n_by_to_site[ts]
-            h = matrix(ham(ts, fs, *args), complex)
+            h = matrix(ham(ts, fs, *args, params=params), complex)
             if h.shape[0] != to_norb[n_ts] or h.shape[1] != from_norb[n_fs]:
                 raise ValueError(msg.format(fs, ts))
             h_sub_view[to_off[n_ts] : to_off[n_ts + 1],
@@ -206,7 +206,7 @@ def make_dense(ham, args, CGraph gr, diag,
 
 
 @cython.boundscheck(False)
-def make_dense_full(ham, args, CGraph gr, diag,
+def make_dense_full(ham, args, params, CGraph gr, diag,
                     gint [:] to_norb, gint [:] to_off,
                     gint [:] from_norb, gint [:] from_off):
     """For internal use by hamiltonian_submatrix."""
@@ -230,7 +230,7 @@ def make_dense_full(ham, args, CGraph gr, diag,
         for ts in nbors.data[:nbors.size]:
             if ts < fs:
                 continue
-            h = mat = matrix(ham(ts, fs, *args), complex)
+            h = mat = matrix(ham(ts, fs, *args, params=params), complex)
             h_herm = mat.transpose().conjugate()
             if h.shape[0] != to_norb[ts] or h.shape[1] != from_norb[fs]:
                 raise ValueError(msg.format(fs, ts))
@@ -243,19 +243,23 @@ def make_dense_full(ham, args, CGraph gr, diag,
 
 @cython.embedsignature(True)
 def hamiltonian_submatrix(self, args=(), to_sites=None, from_sites=None,
-                          sparse=False, return_norb=False):
+                          sparse=False, return_norb=False, *, params=None):
     """Return a submatrix of the system Hamiltonian.
 
     Parameters
     ----------
     args : tuple, defaults to empty
-        Positional arguments to pass to the ``hamiltonian`` method.
+        Positional arguments to pass to the ``hamiltonian`` method. Mutually
+        exclusive with 'params'.
     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``.
+    params : dict, optional
+        Dictionary of parameter names and their values. Mutually exclusive
+        with 'args'.
 
     Returns
     -------
@@ -286,7 +290,8 @@ def hamiltonian_submatrix(self, args=(), to_sites=None, from_sites=None,
         diag = n * [None]
         from_norb = np.empty(n, gint_dtype)
         for site in range(n):
-            diag[site] = h = matrix(ham(site, site, *args), complex)
+            diag[site] = h = matrix(ham(site, site, *args, params=params),
+                                    complex)
             from_norb[site] = h.shape[0]
     else:
         diag = len(from_sites) * [None]
@@ -294,7 +299,7 @@ def hamiltonian_submatrix(self, args=(), to_sites=None, from_sites=None,
         for n_site, site in enumerate(from_sites):
             if site < 0 or site >= n:
                 raise IndexError('Site number out of range.')
-            diag[n_site] = h = matrix(ham(site, site, *args),
+            diag[n_site] = h = matrix(ham(site, site, *args, params=params),
                                       complex)
             from_norb[n_site] = h.shape[0]
     from_off = np.empty(from_norb.shape[0] + 1, gint_dtype)
@@ -308,14 +313,14 @@ def hamiltonian_submatrix(self, args=(), to_sites=None, from_sites=None,
         if to_sites is None:
             to_norb = np.empty(n, gint_dtype)
             for site in range(n):
-                h = matrix(ham(site, site, *args), complex)
+                h = matrix(ham(site, site, *args, params=params), complex)
                 to_norb[site] = h.shape[0]
         else:
             to_norb = np.empty(len(to_sites), gint_dtype)
             for n_site, site in enumerate(to_sites):
                 if site < 0 or site >= n:
                     raise IndexError('Site number out of range.')
-                h = matrix(ham(site, site, *args), complex)
+                h = matrix(ham(site, site, *args, params=params), complex)
                 to_norb[n_site] = h.shape[0]
         to_off = np.empty(to_norb.shape[0] + 1, gint_dtype)
         to_off[0] = 0
@@ -324,7 +329,7 @@ def hamiltonian_submatrix(self, args=(), to_sites=None, from_sites=None,
 
     if to_sites is from_sites is None:
         func = make_sparse_full if sparse else make_dense_full
-        mat = func(ham, args, self.graph, diag, to_norb, to_off,
+        mat = func(ham, args, params, self.graph, diag, to_norb, to_off,
                    from_norb, from_off)
     else:
         if to_sites is None:
@@ -340,7 +345,7 @@ def hamiltonian_submatrix(self, args=(), to_sites=None, from_sites=None,
             from_sites = np.asarray(from_sites, gint_dtype)
 
         func = make_sparse if sparse else make_dense
-        mat = func(ham, args, self.graph, diag, from_sites,
+        mat = func(ham, args, params, self.graph, diag, from_sites,
                    n_by_to_site, to_norb, to_off, from_norb, from_off)
     return (mat, to_norb, from_norb) if return_norb else mat
 

kwant/builder.py

@@ -13,7 +13,7 @@ import abc
 import warnings
 import operator
 import collections
-from functools import total_ordering
+from functools import total_ordering, wraps
 from itertools import islice, chain
 import tinyarray as ta
 import numpy as np
@@ -21,7 +21,7 @@ from scipy import sparse
 from . import system, graph, KwantDeprecationWarning, UserCodeError
 from .operator import Density
 from .physics import DiscreteSymmetry
-from ._common import ensure_isinstance
+from ._common import ensure_isinstance, get_parameters
 
 
 
@@ -490,8 +490,8 @@ class HermConjOfFunc:
     def __init__(self, function):
         self.function = function
 
-    def __call__(self, i, j, *args):
-        return herm_conj(self.function(j, i, *args))
+    def __call__(self, i, j, *args, **kwargs):
+        return herm_conj(self.function(j, i, *args, **kwargs))
 
 
 ################ Leads
@@ -582,8 +582,8 @@ class SelfEnergyLead(Lead):
     Parameters
     ----------
     selfenergy_func : function
-        Function which returns the self energy matrix for the interface sites
-        given the energy and optionally a sequence of extra arguments.
+        Has the same signature as `selfenergy` (without the ``self``
+        parameter) and returns the self energy matrix for the interface sites.
     interface : sequence of `Site` instances
     """
     def __init__(self, selfenergy_func, interface):
@@ -594,8 +594,8 @@ class SelfEnergyLead(Lead):
         """Trivial finalization: the object is returned itself."""
         return self
 
-    def selfenergy(self, energy, args=()):
-        return self.selfenergy_func(energy, args)
+    def selfenergy(self, energy, args=(), *, params=None):
+        return self.selfenergy_func(energy, args, params=params)
 
 
 class ModesLead(Lead):
@@ -604,9 +604,9 @@ class ModesLead(Lead):
     Parameters
     ----------
     modes_func : function
-        Function which returns the modes of the lead as a tuple of
-        `~kwant.physics.PropagatingModes` and `~kwant.physics.StabilizedModes`
-        given the energy and optionally a sequence of extra arguments.
+        Has the same signature as `modes` (without the ``self`` parameter)
+        and returns the modes of the lead as a tuple of
+        `~kwant.physics.PropagatingModes` and `~kwant.physics.StabilizedModes`.
     interface :
         sequence of `Site` instances
 
@@ -619,11 +619,11 @@ class ModesLead(Lead):
         """Trivial finalization: the object is returned itself."""
         return self
 
-    def modes(self, energy, args=()):
-        return self.modes_func(energy, args)
+    def modes(self, energy, args=(), *, params=None):
+        return self.modes_func(energy, args, params=params)
 
-    def selfenergy(self, energy, args=()):
-        stabilized = self.modes(energy, args)[1]
+    def selfenergy(self, energy, args=(), *, params=None):
+        stabilized = self.modes(energy, args, params=params)[1]
         return stabilized.selfenergy()
 
 
@@ -1458,6 +1458,22 @@ class Builder:
 
             lead_interfaces.append(np.array(interface))
 
+        #### Find parameters taken by all value functions
+        hoppings = [self._get_edge(sites[tail], sites[head])
+                           for tail, head in g]
+        onsite_hamiltonians = [self.H[site][1] for site in sites]
+
+        _ham_param_map = {}
+        for hams, skip in [(onsite_hamiltonians, 1), (hoppings, 2)]:
+            for ham in hams:
+                if (not callable(ham) or ham is Other or
+                    ham in _ham_param_map):
+                    continue
+                # parameters come in the same order as in the function signature
+                params, takes_kwargs = get_parameters(ham)
+                params = params[skip:]  # remove site argument(s)
+                _ham_param_map[ham] = (params, takes_kwargs)
+
         #### Assemble and return result.
         result = FiniteSystem()
         result.graph = g
@@ -1465,9 +1481,9 @@ class Builder:
         result.site_ranges = _site_ranges(sites)
         result.id_by_site = id_by_site
         result.leads = finalized_leads
-        result.hoppings = [self._get_edge(sites[tail], sites[head])
-                           for tail, head in g]
-        result.onsite_hamiltonians = [self.H[site][1] for site in sites]
+        result.hoppings = hoppings
+        result.onsite_hamiltonians = onsite_hamiltonians
+        result._ham_param_map = _ham_param_map
         result.lead_interfaces = lead_interfaces
         result.symmetry = self.symmetry
         return result
@@ -1592,6 +1608,18 @@ class Builder:
                 tail, head = sym.to_fd(tail, head)
             hoppings.append(self._get_edge(tail, head))
 
+        #### Find parameters taken by all value functions
+        _ham_param_map = {}
+        for hams, skip in [(onsite_hamiltonians, 1), (hoppings, 2)]:
+            for ham in hams:
+                if (not callable(ham) or ham is Other or
+                    ham in _ham_param_map):
+                    continue
+                # parameters come in the same order as in the function signature
+                params, takes_kwargs = get_parameters(ham)
+                params = params[skip:]  # remove site argument(s)
+                _ham_param_map[ham] = (params, takes_kwargs)
+
         #### Assemble and return result.
         result = InfiniteSystem()
         result.cell_size = cell_size
@@ -1601,6 +1629,7 @@ class Builder:
         result.graph = g
         result.hoppings = hoppings
         result.onsite_hamiltonians = onsite_hamiltonians
+        result._ham_param_map = _ham_param_map
         result.symmetry = self.symmetry
         return result
 
@@ -1613,24 +1642,24 @@ def _raise_user_error(exc, func):
     raise UserCodeError(msg.format(func.__name__)) from exc
 
 
-def discrete_symmetry(self, args=()):
+def discrete_symmetry(self, args=(), *, params=None):
     if self._cons_law is not None:
         eigvals, eigvecs = self._cons_law
-        eigvals = eigvals.tocoo(args)
+        eigvals = eigvals.tocoo(args, params=params)
         if not np.allclose(eigvals.data, np.round(eigvals.data)):
             raise ValueError("Conservation law must have integer eigenvalues.")
         eigvals = np.round(eigvals).tocsr()
         # Avoid appearance of zero eigenvalues
         eigvals = eigvals + 0.5 * sparse.identity(eigvals.shape[0])
         eigvals.eliminate_zeros()
-        eigvecs = eigvecs.tocoo(args)
+        eigvecs = eigvecs.tocoo(args, params=params)
         projectors = [eigvecs.dot(eigvals == val)
                       for val in sorted(np.unique(eigvals.data))]
     else:
         projectors = None
 
     def evaluate(op):
-        return None if op is None else op.tocoo(args)
+        return None if op is None else op.tocoo(args, params=params)
 
     return DiscreteSymmetry(projectors, *(evaluate(symm) for symm in
                                           self._symmetries))
@@ -1651,15 +1680,17 @@ def _transfer_symmetry(syst, builder):
         syst._cons_law = None
 
     elif callable(cons_law):
-        def vals(site, *args):
+        @wraps(cons_law)
+        def vals(site, *args, **kwargs):
             if site.family.norbs == 1:
-                return cons_law(site, *args)
-            return np.diag(np.linalg.eigvalsh(cons_law(site, *args)))
+                return cons_law(site, *args, **kwargs)
+            return np.diag(np.linalg.eigvalsh(cons_law(site, *args, **kwargs)))
 
-        def vecs(site, *args):
+        @wraps(cons_law)
+        def vecs(site, *args, **kwargs):
             if site.family.norbs == 1:
                 return 1
-            return np.linalg.eigh(cons_law(site, *args))[1]
+            return np.linalg.eigh(cons_law(site, *args, **kwargs))[1]
 
         syst._cons_law = operator(vals), operator(vecs)
 
@@ -1687,6 +1718,7 @@ def _transfer_symmetry(syst, builder):
                                                     builder.particle_hole,
                                                     builder.chiral)]
 
+
 class FiniteSystem(system.FiniteSystem):
     """Finalized `Builder` with leads.
 
@@ -1702,14 +1734,25 @@ class FiniteSystem(system.FiniteSystem):
         The inverse of ``sites``; maps from ``i`` to ``sites[i]``.
     """
 
-    def hamiltonian(self, i, j, *args):
+    def hamiltonian(self, i, j, *args, params=None):
+        if args and params:
+            raise TypeError("'args' and 'params' are mutually exclusive.")
         if i == j:
             value = self.onsite_hamiltonians[i]
             if callable(value):
-                try:
-                    value = value(self.sites[i], *args)
-                except Exception as exc:
-                    _raise_user_error(exc, value)
+                if params:
+                    param_names, takes_kwargs = self._ham_param_map[value]
+                    if not takes_kwargs:
+                        params = {pn: params[pn] for pn in param_names}
+                    try:
+                            value = value(self.sites[i], **params)
+                    except Exception as exc:
+                        _raise_user_error(exc, value)
+                else:
+                    try:
+                        value = value(self.sites[i], *args)
+                    except Exception as exc:
+                        _raise_user_error(exc, value)
         else:
             edge_id = self.graph.first_edge_id(i, j)
             value = self.hoppings[edge_id]
@@ -1720,10 +1763,19 @@ class FiniteSystem(system.FiniteSystem):
                 value = self.hoppings[edge_id]
             if callable(value):
                 sites = self.sites
-                try:
-                    value = value(sites[i], sites[j], *args)
-                except Exception as exc:
-                    _raise_user_error(exc, value)
+                if params:
+                    param_names, takes_kwargs = self._ham_param_map[value]
+                    if not takes_kwargs:
+                        params = {pn: params[pn] for pn in param_names}
+                    try:
+                        value = value(sites[i], sites[j], **params)
+                    except Exception as exc:
+                        _raise_user_error(exc, value)
+                else:
+                    try:
+                        value = value(sites[i], sites[j], *args)
+                    except Exception as exc:
+                        _raise_user_error(exc, value)
             if conj:
                 value = herm_conj(value)
         return value
@@ -1759,16 +1811,28 @@ class InfiniteSystem(system.InfiniteSystem):
     hoppings. Each of these three subsequences is individually sorted.
     """
 
-    def hamiltonian(self, i, j, *args):
+    def hamiltonian(self, i, j, *args, params=None):
+        if args and params:
+            raise TypeError("'args' and 'params' are mutually exclusive.")
         if i == j:
             if i >= self.cell_size:
                 i -= self.cell_size
             value = self.onsite_hamiltonians[i]
             if callable(value):
-                try:
-                    value = value(self.symmetry.to_fd(self.sites[i]), *args)
-                except Exception as exc:
-                    _raise_user_error(exc, value)
+                site = self.symmetry.to_fd(self.sites[i])
+                if params:
+                    param_names, takes_kwargs = self._ham_param_map[value]
+                    if not takes_kwargs:
+                        params = {pn: params[pn] for pn in param_names}
+                    try:
+                        value = value(site, **params)
+                    except Exception as exc:
+                        _raise_user_error(exc, value)
+                else:
+                    try:
+                        value = value(site, *args)
+                    except Exception as exc:
+                        _raise_user_error(exc, value)
         else:
             edge_id = self.graph.first_edge_id(i, j)
             value = self.hoppings[edge_id]
@@ -1779,13 +1843,20 @@ class InfiniteSystem(system.InfiniteSystem):
                 value = self.hoppings[edge_id]
             if callable(value):
                 sites = self.sites
-                site_i = sites[i]
-                site_j = sites[j]
-                site_i, site_j = self.symmetry.to_fd(site_i, site_j)
-                try:
-                    value = value(site_i, site_j, *args)
-                except Exception as exc:
-                    _raise_user_error(exc, value)
+                site_i, site_j = self.symmetry.to_fd(sites[i], sites[j])
+                if params:
+                    param_names, takes_kwargs = self._ham_param_map[value]
+                    if not takes_kwargs:
+                        params = {pn: params[pn] for pn in param_names}
+                    try:
+                        value = value(site_i, site_j, **params)
+                    except Exception as exc:
+                        _raise_user_error(exc, value)
+                else:
+                    try:
+                        value = value(site_i, site_j, *args)
+                    except Exception as exc:
+                        _raise_user_error(exc, value)
             if conj:
                 value = herm_conj(value)
         return value

kwant/operator.pyx

@@ -24,7 +24,7 @@ from .graph.core cimport EdgeIterator
 from .graph.defs cimport gint
 from .graph.defs import gint_dtype
 from .system import InfiniteSystem
-from ._common import UserCodeError
+from ._common import UserCodeError, get_parameters
 
 
 ################ Generic Utility functions
@@ -92,7 +92,7 @@ cdef int _check_onsite(complex[:, :] M, gint norbs,
     return 0
 
 
-cdef int _check_ham(complex[:, :] H, ham, args,
+cdef int _check_ham(complex[:, :] H, ham, args, params,
                     gint a, gint a_norbs, gint b, gint b_norbs,
                     int check_hermiticity) except -1:
     "Check Hamiltonian matrix for correct shape and hermiticity."
@@ -100,7 +100,8 @@ cdef int _check_ham(complex[:, :] H, ham, args,
         raise UserCodeError(_shape_msg.format('Hamiltonian'))
     if check_hermiticity:
         # call the "partner" element if we are not on the diagonal
-        H_conj = H if a == b else ta.matrix(ham(b, a, *args), complex)
+        H_conj = H if a == b else ta.matrix(ham(b, a, *args, params=params),
+                                                complex)
         if not _is_herm_conj(H_conj, H):
             raise ValueError(_herm_msg.format('Hamiltonian'))
     return 0
@@ -237,11 +238,18 @@ def _normalize_onsite(syst, onsite, check_hermiticity):
     If `onsite` is a function or a mapping (dictionary) then a function
     is returned.
     """
+    parameter_info =  ((), False)
+
     if callable(onsite):
+        # make 'onsite' compatible with hamiltonian value functions
+        parameters, takes_kwargs = get_parameters(onsite)
+        parameters = parameters[1:]  # skip 'site' parameter
+        parameter_info = (tuple(parameters), takes_kwargs)
         try:
             _sites = syst.sites
-            def _onsite(site_id, *args):
-                return onsite(_sites[site_id], *args)
+            @ft.wraps(onsite)
+            def _onsite(site_id, *args, **kwargs):
+                return onsite(_sites[site_id], *args, **kwargs)
         except AttributeError:
             _onsite = onsite
     elif isinstance(onsite, collections.Mapping):
@@ -280,7 +288,7 @@ def _normalize_onsite(syst, onsite, check_hermiticity):
                    'different numbers of orbitals on different sites')
             raise ValueError(msg)
 
-    return _onsite
+    return _onsite, parameter_info
 
 
 cdef class BlockSparseMatrix:
@@ -396,7 +404,7 @@ cdef class _LocalOperator:
     """
 
     cdef public int check_hermiticity, sum
-    cdef public object syst, onsite
+    cdef public object syst, onsite, _onsite_params_info
     cdef public gint[:, :]  where, _site_ranges
     cdef public BlockSparseMatrix _bound_onsite, _bound_hamiltonian
 
@@ -410,7 +418,8 @@ cdef class _LocalOperator:
                              'the site families (lattices).')
 
         self.syst = syst
-        self.onsite = _normalize_onsite(syst, onsite, check_hermiticity)
+        self.onsite, self._onsite_params_info = \
+            _normalize_onsite(syst, onsite, check_hermiticity)
         self.check_hermiticity = check_hermiticity
         self.sum = sum
         self._site_ranges = np.asarray(syst.site_ranges, dtype=gint_dtype)
@@ -419,8 +428,8 @@ cdef class _LocalOperator:
         self._bound_hamiltonian = None
 
     @cython.embedsignature
-    def __call__(self, bra, ket=None, args=()):
-        r"""Return matrix elements of the operator.
+    def __call__(self, bra, ket=None, args=(), *, params=None):
+        r"""Return the matrix elements of the operator.
 
         An operator ``A`` can be called like
 
@@ -449,6 +458,10 @@ cdef class _LocalOperator:
         args : tuple, optional
             The arguments to pass to the system. Used to evaluate
             the ``onsite`` elements and, possibly, the system Hamiltonian.
+            Mutually exclusive with 'params'.
+        params : dict, optional
+            Dictionary of parameter names and their values. Mutually exclusive
+            with 'args'.
 
         Returns
         -------
@@ -456,10 +469,12 @@ cdef class _LocalOperator:
         otherwise `complex`. If this operator was created with ``sum=True``,
         then a single value is returned, otherwise an array is returned.
         """
-        if (self._bound_onsite or self._bound_hamiltonian) and args:
-            raise ValueError('Extra arguments are already bound to this '
-                             'operator. You should call this operator '
-                             'without providing `args`.')
+        if (self._bound_onsite or self._bound_hamiltonian) and (args or params):
+            raise ValueError("Extra arguments are already bound to this "
+                             "operator. You should call this operator "
+                             "providing neither 'args' nor 'params'.")
+        if args and params:
+            raise TypeError("'args' and 'params' are mutually exclusive.")
         if bra is None:
             raise TypeError('bra must be an array')
         bra = np.asarray(bra, dtype=complex)
@@ -479,14 +494,15 @@ cdef class _LocalOperator:
             where = where.reshape(-1)
 
         result = np.zeros((self.where.shape[0],), dtype=complex)
-        self._operate(out_data=result, bra=bra, ket=ket, args=args, op=MAT_ELS)
+        self._operate(out_data=result, bra=bra, ket=ket, args=args,
+                      params=params, op=MAT_ELS)
         # if everything is Hermitian then result is real if bra == ket
         if self.check_hermiticity and bra is ket:
             result = result.real
         return np.sum(result) if self.sum else result
 
     @cython.embedsignature
-    def act(self, ket, args=()):
+    def act(self, ket, args=(), *, params=None):
         """Act with the operator on a wavefunction.
 
         For an operator :math:`Q_{iαβ}` and ``ket`` :math:`ψ_β`
@@ -498,16 +514,21 @@ cdef class _LocalOperator:
             Wavefunctions defined over all the orbitals of the system.
         args : tuple
             The extra arguments to the Hamiltonian value functions and
-            the operator ``onsite`` function.
+            the operator ``onsite`` function. Mutually exclusive with 'params'.
+        params : dict, optional
+            Dictionary of parameter names and their values. Mutually exclusive
+            with 'args'.
 
         Returns
         -------
         Array of `complex`.
         """
-        if (self._bound_onsite or self._bound_hamiltonian) and args:
-            raise ValueError('Extra arguments are already bound to this '
-                             'operator. You should call this operator '
-                             'without providing `args`.')
+        if (self._bound_onsite or self._bound_hamiltonian) and (args or params):
+            raise ValueError("Extra arguments are already bound to this "
+                             "operator. You should call this operator "
+                             "providing neither 'args' nor 'params'.")
+        if args and params:
+            raise TypeError("'args' and 'params' are mutually exclusive.")
 
         if ket is None:
             raise TypeError('ket must be an array')
@@ -516,32 +537,36 @@ cdef class _LocalOperator:
         if ket.shape != (tot_norbs,):
             raise ValueError('ket vector is incorrect shape')
         result = np.zeros((tot_norbs,), dtype=np.complex)
-        self._operate(out_data=result, bra=None, ket=ket, args=args, op=ACT)
+        self._operate(out_data=result, bra=None, ket=ket, args=args,
+                      params=params, op=ACT)
         return result
 
     @cython.embedsignature
-    def bind(self, args=()):
+    def bind(self, args=(), *, params=None):
         """Bind the given arguments to this operator.
 
         Returns a copy of this operator that does not need to be passed extra
         arguments when subsequently called or when using the ``act`` method.
         """
+        if args and params:
+            raise TypeError("'args' and 'params' are mutually exclusive.")
         # generic creation of new instance
         cls = self.__class__
         q = cls.__new__(cls)
         q.syst = self.syst
         q.onsite = self.onsite
+        q._onsite_params_info = self._onsite_params_info
         q.where = self.where
         q.sum = self.sum
         q._site_ranges = self._site_ranges
         q.check_hermiticity = self.check_hermiticity
         if callable(self.onsite):
-            q._bound_onsite = self._eval_onsites(args)
+            q._bound_onsite = self._eval_onsites(args, params)
         # NOTE: subclasses should populate `bound_hamiltonian` if needed
         return q
 
     def _operate(self, complex[:] out_data, complex[:] bra, complex[:] ket,
-                 args, operation op):
+                 args, operation op, *, params=None):
         """Do an operation with the operator.
 
         Parameters
@@ -555,38 +580,47 @@ cdef class _LocalOperator:
             If `op` is `ACT` then `bra` is None.
         args : tuple
             The extra arguments to the Hamiltonian value functions and
-            the operator ``onsite`` function.
+            the operator ``onsite`` function. Mutually exclusive with 'params'.
         op : operation
             The operation to perform.
             `MAT_ELS`: calculate matrix elements between `bra` and `ket`
             `ACT`: act on `ket` with the operator
+        params : dict, optional
+            Dictionary of parameter names and their values. Mutually exclusive
+            with 'args'.
         """
         raise NotImplementedError()
 
-    cdef BlockSparseMatrix _eval_onsites(self, args):
+    cdef BlockSparseMatrix _eval_onsites(self, args, params):
         """Evaluate the onsite matrices on all elements of `where`"""
         assert callable(self.onsite)
+        assert not (args and params)
+        params = params or {}
         matrix = ta.matrix
         onsite = self.onsite
         check_hermiticity = self.check_hermiticity
 
+        param_names, takes_kwargs = self._onsite_params_info
+        if params and not takes_kwargs:
+            params = {pn: params[pn] for pn in param_names}
+
         def get_onsite(a, a_norbs, b, b_norbs):
-            mat = matrix(onsite(a, *args), complex)
+            mat = matrix(onsite(a, *args, **params), complex)
             _check_onsite(mat, a_norbs, check_hermiticity)
             return mat
 
         offsets, norbs = _get_all_orbs(self.where, self._site_ranges)
         return  BlockSparseMatrix(self.where, offsets, norbs, get_onsite)
 
-    cdef BlockSparseMatrix _eval_hamiltonian(self, args):
+    cdef BlockSparseMatrix _eval_hamiltonian(self, args, params):
         """Evaluate the Hamiltonian on all elements of `where`."""
         matrix = ta.matrix
         hamiltonian = self.syst.hamiltonian
         check_hermiticity = self.check_hermiticity
 
         def get_ham(a, a_norbs, b, b_norbs):
-            mat = matrix(hamiltonian(a, b, *args), complex)
-            _check_ham(mat, hamiltonian, args,
+            mat = matrix(hamiltonian(a, b, *args, params=params), complex)
+            _check_ham(mat, hamiltonian, args, params,
                        a, a_norbs, b, b_norbs, check_hermiticity)
             return mat
 
@@ -641,7 +675,7 @@ cdef class Density(_LocalOperator):
     @cython.boundscheck(False)
     @cython.wraparound(False)
     def _operate(self, complex[:] out_data, complex[:] bra, complex[:] ket,
-                 args, operation op):
+                 args, operation op, *, params=None):
         matrix = ta.matrix
         cdef int unique_onsite = not callable(self.onsite)
         # prepare onsite matrices
@@ -655,7 +689,7 @@ cdef class Density(_LocalOperator):
         elif self._bound_onsite:
             M_a_blocks = self._bound_onsite
         else:
-            M_a_blocks = self._eval_onsites(args)
+            M_a_blocks = self._eval_onsites(args, params)
 
         # loop-local variables
         cdef gint a, a_s, a_norbs
@@ -688,15 +722,17 @@ cdef class Density(_LocalOperator):
     @cython.wraparound(False)
     @cython.cdivision(True)
     @cython.embedsignature
-    def tocoo(self, args=()):
+    def tocoo(self, args=(), *, params=None):
         """Convert the operator to coordinate format sparse matrix."""
         cdef int blk, blk_size, n_blocks, n, k = 0
         cdef int [:, :] offsets, shapes
         cdef int [:] row, col
-        if self._bound_onsite and args:
+        if self._bound_onsite and (args or params):
            raise ValueError("Extra arguments are already bound to this "
                             "operator. You should call this operator "
-                            "without providing 'args'.")
+                            "providing neither 'args' nor 'params'.")
+        if args and params:
+            raise TypeError("'args' and 'params' are mutually exclusive.")
 
         if not callable(self.onsite):
             offsets = _get_all_orbs(self.where, self._site_ranges)[0]
@@ -709,7 +745,7 @@ cdef class Density(_LocalOperator):
             if self._bound_onsite is not None:
                 onsite_matrix = self._bound_onsite
             else:
-                onsite_matrix = self._eval_onsites(args)
+                onsite_matrix = self._eval_onsites(args, params)
             data = onsite_matrix.data
             offsets = np.asarray(onsite_matrix.block_offsets)
             shapes = np.asarray(onsite_matrix.block_shapes)
@@ -781,20 +817,20 @@ cdef class Current(_LocalOperator):
                          check_hermiticity=check_hermiticity, sum=sum)
 
     @cython.embedsignature
-    def bind(self, args=()):
+    def bind(self, args=(), *, params=None):
         """Bind the given arguments to this operator.
 
         Returns a copy of this operator that does not need to be passed extra
         arguments when subsequently called or when using the ``act`` method.
         """
-        q = super().bind(args)
-        q._bound_hamiltonian = self._eval_hamiltonian(args)
+        q = super().bind(args, params=params)
+        q._bound_hamiltonian = self._eval_hamiltonian(args, params)
         return q
 
     @cython.boundscheck(False)
     @cython.wraparound(False)
     def _operate(self, complex[:] out_data, complex[:] bra, complex[:] ket,
-                 args, operation op):
+                 args, operation op, *, params=None):
         # prepare onsite matrices and hamiltonians
         cdef int unique_onsite = not callable(self.onsite)
         cdef complex[:, :] _tmp_mat
@@ -808,12 +844,12 @@ cdef class Current(_LocalOperator):
         elif self._bound_onsite:
             M_a_blocks = self._bound_onsite
         else:
-            M_a_blocks = self._eval_onsites(args)
+            M_a_blocks = self._eval_onsites(args, params)
 
         if self._bound_hamiltonian:
             H_ab_blocks = self._bound_hamiltonian
         else:
-            H_ab_blocks = self._eval_hamiltonian(args)
+            H_ab_blocks = self._eval_hamiltonian(args, params)
 
         # main loop
         cdef gint a, a_s, a_norbs, b, b_s, b_norbs
@@ -904,20 +940,20 @@ cdef class Source(_LocalOperator):
                          check_hermiticity=check_hermiticity, sum=sum)
 
     @cython.embedsignature
-    def bind(self, args=()):
+    def bind(self, args=(), *, params=None):
         """Bind the given arguments to this operator.
 
         Returns a copy of this operator that does not need to be passed extra
         arguments when subsequently called or when using the ``act`` method.
         """
-        q = super().bind(args)
-        q._bound_hamiltonian = self._eval_hamiltonian(args)
+        q = super().bind(args, params=params)
+        q._bound_hamiltonian = self._eval_hamiltonian(args, params)
         return q
 
     @cython.boundscheck(False)
     @cython.wraparound(False)
     def _operate(self, complex[:] out_data, complex[:] bra, complex[:] ket,
-                 args, operation op):
+                 args, operation op, *, params=None):
         # prepare onsite matrices and hamiltonians
         cdef int unique_onsite = not callable(self.onsite)
         cdef complex[:, :] _tmp_mat
@@ -931,12 +967,12 @@ cdef class Source(_LocalOperator):
         elif self._bound_onsite:
             M_a_blocks = self._bound_onsite
         else:
-            M_a_blocks = self._eval_onsites(args)
+            M_a_blocks = self._eval_onsites(args, params)
 
         if self._bound_hamiltonian:
             H_aa_blocks = self._bound_hamiltonian
         else:
-            H_aa_blocks = self._eval_hamiltonian(args)
+            H_aa_blocks = self._eval_hamiltonian(args, params)
 
         # main loop
         cdef gint a, a_s, a_norbs

kwant/physics/dispersion.py

@@ -25,6 +25,10 @@ class Bands:
         calculated.
     args : tuple, defaults to empty
         Positional arguments to pass to the ``hamiltonian`` method.
+        Mutually exclusive with 'params'.
+    params : dict, optional
+        Dictionary of parameter names and their values. Mutually exclusive
+        with 'args'.
 
     Notes
     -----
@@ -42,13 +46,13 @@ class Bands:
     >>> pyplot.show()
     """
 
-    def __init__(self, sys, args=()):
+    def __init__(self, sys, args=(), *, params=None):
         syst = sys
         ensure_isinstance(syst, system.InfiniteSystem)
-        self.ham = syst.cell_hamiltonian(args)
+        self.ham = syst.cell_hamiltonian(args, params=params)
         if not np.allclose(self.ham, self.ham.T.conj()):
             raise ValueError('The cell Hamiltonian is not Hermitian.')
-        hop = syst.inter_cell_hopping(args)
+        hop = syst.inter_cell_hopping(args, params=params)
         self.hop = np.empty(self.ham.shape, dtype=complex)
         self.hop[:, : hop.shape[1]] = hop
         self.hop[:, hop.shape[1]:] = 0

kwant/plotter.py

@@ -1592,7 +1592,7 @@ def map(sys, value, colorbar=True, cmap=None, vmin=None, vmax=None, a=None,
 
 
 def bands(sys, args=(), momenta=65, file=None, show=True, dpi=None,
-          fig_size=None, ax=None):
+          fig_size=None, ax=None, *, params=None):
     """Plot band structure of a translationally invariant 1D system.
 
     Parameters
@@ -1601,6 +1601,7 @@ def bands(sys, args=(), momenta=65, file=None, show=True, dpi=None,
         A system bands of which are to be plotted.
     args : tuple, defaults to empty
         Positional arguments to pass to the ``hamiltonian`` method.
+        Mutally exclusive with 'params'.
     momenta : int or 1D array-like
         Either a number of sampling points on the interval [-pi, pi], or an
         array of points at which the band structure has to be evaluated.
@@ -1619,6 +1620,9 @@ def bands(sys, args=(), momenta=65, file=None, show=True, dpi=None,
         If `ax` is not `None`, no new figure is created, but the plot is done
         within the existing Axes `ax`. in this case, `file`, `show`, `dpi`
         and `fig_size` are ignored.
+    params : dict, optional
+        Dictionary of parameter names and their values. Mutually exclusive
+        with 'args'.
 
     Returns
     -------
@@ -1639,7 +1643,7 @@ def bands(sys, args=(), momenta=65, file=None, show=True, dpi=None,
     if momenta.ndim != 1:
         momenta = np.linspace(-np.pi, np.pi, momenta)
 
-    bands = physics.Bands(syst, args)
+    bands = physics.Bands(syst, args, params=params)
     energies = [bands(k) for k in momenta]
 
     if ax is None:

kwant/solvers/common.py

@@ -97,7 +97,8 @@ class SparseSolver(metaclass=abc.ABCMeta):
         pass
 
     def _make_linear_sys(self, sys, in_leads, energy=0, args=(),
-                         check_hermiticity=True, realspace=False):
+                         check_hermiticity=True, realspace=False,
+                         *, params=None):
         """Make a sparse linear system of equations defining a scattering
         problem.
 
@@ -112,6 +113,7 @@ class SparseSolver(metaclass=abc.ABCMeta):
             Excitation energy at which to solve the scattering problem.
         args : tuple, defaults to empty
             Positional arguments to pass to the ``hamiltonian`` method.
+            Mutually exclusive with 'params'.
         check_hermiticity : bool
             Check if Hamiltonian matrices are in fact Hermitian.
             Enables deduction of missing transmission coefficients.
@@ -119,6 +121,9 @@ class SparseSolver(metaclass=abc.ABCMeta):
             Calculate Green's function between the outermost lead
             sites, instead of lead modes. This is almost always
             more computationally expensive and less stable.
+        params : dict, optional
+            Dictionary of parameter names and their values. Mutually exclusive
+            with 'args'.
 
         Returns
         -------
@@ -155,7 +160,8 @@ class SparseSolver(metaclass=abc.ABCMeta):
         if not syst.lead_interfaces:
             raise ValueError('System contains no leads.')
         lhs, norb = syst.hamiltonian_submatrix(args, sparse=True,
-                                              return_norb=True)[:2]
+                                               return_norb=True,
+                                               params=params)[:2]
         lhs = getattr(lhs, 'to' + self.lhsformat)()
         lhs = lhs - energy * sp.identity(lhs.shape[0], format=self.lhsformat)
         num_orb = lhs.shape[0]
@@ -182,7 +188,7 @@ class SparseSolver(metaclass=abc.ABCMeta):
         for leadnum, interface in enumerate(syst.lead_interfaces):
             lead = syst.leads[leadnum]
             if not realspace:
-                prop, stab = lead.modes(energy, args)
+                prop, stab = lead.modes(energy, args, params=params)
                 lead_info.append(prop)
                 u = stab.vecs
                 ulinv = stab.vecslmbdainv
@@ -241,7 +247,7 @@ class SparseSolver(metaclass=abc.ABCMeta):
                 else:
                     rhs.append(None)
             else:
-                sigma = lead.selfenergy(energy, args)
+                sigma = lead.selfenergy(energy, args, params=params)
                 lead_info.append(sigma)
                 coords = np.r_[tuple(slice(offsets[i], offsets[i + 1])
                                       for i in interface)]
@@ -291,7 +297,8 @@ class SparseSolver(metaclass=abc.ABCMeta):
         return LinearSys(lhs, rhs, indices, num_orb), lead_info
 
     def smatrix(self, sys, energy=0, args=(),
-                out_leads=None, in_leads=None, check_hermiticity=True):
+                out_leads=None, in_leads=None, check_hermiticity=True,
+                *, params=None):
         """
         Compute the scattering matrix of a system.
 
@@ -304,6 +311,7 @@ class SparseSolver(metaclass=abc.ABCMeta):
             Excitation energy at which to solve the scattering problem.
         args : tuple, defaults to empty
             Positional arguments to pass to the ``hamiltonian`` method.
+            Mutually exclusive with 'params'.
         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
@@ -315,6 +323,9 @@ class SparseSolver(metaclass=abc.ABCMeta):
         check_hermiticity : ``bool``
             Check if the Hamiltonian matrices are Hermitian.
             Enables deduction of missing transmission coefficients.
+        params : dict, optional
+            Dictionary of parameter names and their values. Mutually exclusive
+            with 'args'.
 
         Returns
         -------
@@ -355,7 +366,8 @@ class SparseSolver(metaclass=abc.ABCMeta):
             raise ValueError("No output is requested.")
 
         linsys, lead_info = self._make_linear_sys(syst, in_leads, energy, args,
-                                                  check_hermiticity, False)
+                                                  check_hermiticity, False,
+                                                  params=params)
 
         kept_vars = np.concatenate([coords for i, coords in
                                     enumerate(linsys.indices) if i in
@@ -377,7 +389,8 @@ class SparseSolver(metaclass=abc.ABCMeta):
         return SMatrix(data, lead_info, out_leads, in_leads, check_hermiticity)
 
     def greens_function(self, sys, energy=0, args=(),
-                        out_leads=None, in_leads=None, check_hermiticity=True):
+                        out_leads=None, in_leads=None, check_hermiticity=True,
+                        *, params=None):
         """
         Compute the retarded Green's function of the system between its leads.
 
@@ -390,6 +403,7 @@ class SparseSolver(metaclass=abc.ABCMeta):
             Excitation energy at which to solve the scattering problem.
         args : tuple, defaults to empty
             Positional arguments to pass to the ``hamiltonian`` method.
+            Mutually exclusive with 'params'.
         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
@@ -401,6 +415,9 @@ class SparseSolver(metaclass=abc.ABCMeta):
         check_hermiticity : ``bool``
             Check if the Hamiltonian matrices are Hermitian.
             Enables deduction of missing transmission coefficients.
+        params : dict, optional
+            Dictionary of parameter names and their values. Mutually exclusive
+            with 'args'.
 
         Returns
         -------
@@ -444,7 +461,8 @@ class SparseSolver(metaclass=abc.ABCMeta):
             raise ValueError("No output is requested.")
 
         linsys, lead_info = self._make_linear_sys(syst, in_leads, energy, args,
-                                                  check_hermiticity, True)
+                                                  check_hermiticity, True,
+                                                  params=params)
 
         kept_vars = np.concatenate([coords for i, coords in
                                     enumerate(linsys.indices) if i in
@@ -466,7 +484,8 @@ class SparseSolver(metaclass=abc.ABCMeta):
         return GreensFunction(data, lead_info, out_leads, in_leads,
                               check_hermiticity)
 
-    def ldos(self, sys, energy=0, args=(), check_hermiticity=True):
+    def ldos(self, sys, energy=0, args=(), check_hermiticity=True,
+             *, params=None):
         """
         Calculate the local density of states of a system at a given energy.
 
@@ -479,9 +498,13 @@ class SparseSolver(metaclass=abc.ABCMeta):
             Excitation energy at which to solve the scattering problem.
         args : tuple of arguments, or empty tuple
             Positional arguments to pass to the function(s) which
-            evaluate the hamiltonian matrix elements
+            evaluate the hamiltonian matrix elements.  Mutually exclusive
+            with 'params'.
         check_hermiticity : ``bool``
             Check if the Hamiltonian matrices are Hermitian.
+        params : dict, optional
+            Dictionary of parameter names and their values. Mutually exclusive
+            with 'args'.
 
         Returns
         -------
@@ -502,7 +525,8 @@ class SparseSolver(metaclass=abc.ABCMeta):
                                           "self-energy is not implemented yet.")
 
         linsys = self._make_linear_sys(syst, range(len(syst.leads)), energy,
-                                       args, check_hermiticity)[0]
+                                       args, check_hermiticity,
+                                       params=params)[0]
 
         ldos = np.zeros(linsys.num_orb, float)
 
@@ -523,7 +547,8 @@ class SparseSolver(metaclass=abc.ABCMeta):
 
         return ldos * (0.5 / np.pi)
 
-    def wave_function(self, sys, energy=0, args=(), check_hermiticity=True):
+    def wave_function(self, sys, energy=0, args=(), check_hermiticity=True,
+                      *, params=None):
         """
         Return a callable object for the computation of the wave function
         inside the scattering region.
@@ -535,9 +560,13 @@ class SparseSolver(metaclass=abc.ABCMeta):
             calculated.
         args : tuple of arguments, or empty tuple
             Positional arguments to pass to the function(s) which
-            evaluate the hamiltonian matrix elements
+            evaluate the hamiltonian matrix elements. Mutually exclusive
+            with 'params'.
         check_hermiticity : ``bool``
             Check if the Hamiltonian matrices are Hermitian.
+        params : dict, optional
+            Dictionary of parameter names and their values. Mutually exclusive
+            with 'args'.
 
         Notes
         -----
@@ -555,11 +584,11 @@ class SparseSolver(metaclass=abc.ABCMeta):
         >>> wfs_of_lead_2 = wf(2)
 
         """
-        return WaveFunction(self, sys, energy, args, check_hermiticity)
+        return WaveFunction(self, sys, energy, args, check_hermiticity, params)
 
 
 class WaveFunction:
-    def __init__(self, solver, sys, energy, args, check_hermiticity):
+    def __init__(self, solver, sys, energy, args, check_hermiticity, params):
         syst = sys  # ensure consistent naming across function bodies
         ensure_isinstance(syst, system.System)
         for lead in syst.leads:
@@ -569,7 +598,8 @@ class WaveFunction:
                        ' are not available yet.')
                 raise NotImplementedError(msg)
         linsys = solver._make_linear_sys(syst, range(len(syst.leads)), energy,
-                                         args, check_hermiticity)[0]
+                                         args, check_hermiticity,
+                                         params=params)[0]
         self.solve = solver._solve_linear_sys
         self.rhs = linsys.rhs
         self.factorized_h = solver._factorized(linsys.lhs)

kwant/solvers/tests/_test_sparse.py

@@ -23,8 +23,9 @@ class LeadWithOnlySelfEnergy:
     def __init__(self, lead):
         self.lead = lead
 
-    def selfenergy(self, energy, args=()):
+    def selfenergy(self, energy, args=(), *, params=None):
         assert args == ()
+        assert params == None
         return self.lead.selfenergy(energy)
 
 
@@ -500,3 +501,41 @@ def test_wavefunc_ldos_consistency(wave_function, ldos):
     raises(ValueError, check, syst.precalculate(what='selfenergy'))
     syst.leads[0] = LeadWithOnlySelfEnergy(syst.leads[0])
     raises(NotImplementedError, check, syst)
+
+
+def test_arg_passing(wave_function, ldos, smatrix):
+
+    def onsite(site, a, b):
+        return site.pos[0] + site.pos[1] + a + b
+
+    def hopping(site1, site2, a, b):
+        return b - a
+
+    W = 3
+    L = 4
+
+    syst = kwant.Builder()
+    syst[(square(i, j) for i in range(L) for j in range(W))] = onsite
+    syst[square.neighbors()] = hopping
+
+    lead = kwant.Builder(kwant.TranslationalSymmetry((-1, 0)))
+    lead[(square(0, j) for j in range(W))] = onsite
+    lead[square.neighbors()] = hopping
+
+    syst.attach_lead(lead)
+    syst.attach_lead(lead.reversed())
+
+    fsyst = syst.finalized()
+
+    # compare results to those when we pass `args` only
+    args = (1, 3)
+    params = dict(a=1, b=3)
+    np.testing.assert_array_equal(
+        wave_function(fsyst, args=args)(0),
+        wave_function(fsyst, params=params)(0))
+    np.testing.assert_array_equal(
+        ldos(fsyst, args=args),
+        ldos(fsyst, params=params))
+    np.testing.assert_array_equal(
+        smatrix(fsyst, args=args).data,
+        smatrix(fsyst, params=params).data)

kwant/solvers/tests/test_mumps.py

@@ -112,3 +112,7 @@ def test_wavefunc_ldos_consistency():
     for opts in opt_list:
         options(**opts)
         _test_sparse.test_wavefunc_ldos_consistency(wave_function, ldos)
+
+
+def test_arg_passing():
+    _test_sparse.test_arg_passing(wave_function, ldos, smatrix)

kwant/solvers/tests/test_sparse.py

@@ -59,3 +59,6 @@ def test_ldos():
 
 def test_wavefunc_ldos_consistency():
     _test_sparse.test_wavefunc_ldos_consistency(wave_function, ldos)
+
+def test_arg_passing():
+    _test_sparse.test_arg_passing(wave_function, ldos, smatrix)

kwant/system.py

@@ -47,7 +47,7 @@ class System(metaclass=abc.ABCMeta):
     numbers of orbitals.
     """
     @abc.abstractmethod
-    def hamiltonian(self, i, j, *args):
+    def hamiltonian(self, i, j, *args, params=None):
         """Return the hamiltonian matrix element for sites ``i`` and ``j``.
 
         If ``i == j``, return the on-site Hamiltonian of site ``i``.
@@ -59,7 +59,7 @@ class System(metaclass=abc.ABCMeta):
         """
         pass
 
-    def discrete_symmetry(self, args):
+    def discrete_symmetry(self, args, *, params=None):
         """Return the discrete symmetry of the system."""
         # Avoid the circular import.
         from .physics import DiscreteSymmetry
@@ -76,9 +76,11 @@ class FiniteSystem(System, metaclass=abc.ABCMeta):
     Attributes
     ----------
     leads : sequence of leads
-        Each lead has to provide a method
-        ``selfenergy(energy, args)``.
-        It may provide ``modes(energy, args)`` as well.
+        Each lead has to provide a method ``selfenergy`` that has
+        the same signature as `InfiniteSystem.selfenergy` (without the
+        ``self`` parameter). It may also provide ``modes`` that has the
+        same signature as `InfiniteSystem.modes` (without the ``self``
+        parameter).
     lead_interfaces : sequence of sequences of integers
         Each sub-sequence contains the indices of the system sites
         to which the lead is connected.
@@ -99,7 +101,7 @@ class FiniteSystem(System, metaclass=abc.ABCMeta):
     """
 
     def precalculate(self, energy=0, args=(), leads=None,
-                     what='modes'):
+                     what='modes', *, params=None):
         """
         Precalculate modes or self-energies in the leads.
 
@@ -113,13 +115,17 @@ class FiniteSystem(System, metaclass=abc.ABCMeta):
             Energy at which the modes or self-energies have to be
             evaluated.
         args : sequence
-            Additional parameters required for calculating the Hamiltionians
+            Additional parameters required for calculating the Hamiltionians.
+            Mutually exclusive with 'params'.
         leads : sequence of integers or None
             Numbers of the leads to be precalculated. If ``None``, all are
             precalculated.
         what : 'modes', 'selfenergy', 'all'
             The quantitity to precompute. 'all' will compute both
             modes and self-energies. Defaults to 'modes'.
+        params : dict, optional
+            Dictionary of parameter names and their values. Mutually exclusive
+            with 'args'.
 
         Returns
         -------
@@ -147,12 +153,12 @@ class FiniteSystem(System, metaclass=abc.ABCMeta):
                 continue
             modes, selfenergy = None, None
             if what in ('modes', 'all'):
-                modes = lead.modes(energy, args)
+                modes = lead.modes(energy, args, params=params)
             if what in ('selfenergy', 'all'):
                 if modes:
                     selfenergy = modes[1].selfenergy()
                 else:
-                    selfenergy = lead.selfenergy(energy, args)
+                    selfenergy = lead.selfenergy(energy, args, params=params)
             new_leads.append(PrecalculatedLead(modes, selfenergy))
         result.leads = new_leads
         return result
@@ -200,29 +206,29 @@ class InfiniteSystem(System, metaclass=abc.ABCMeta):
     exchanged, as well as of site 3 and 4.
 
     """
-    def cell_hamiltonian(self, args=(), sparse=False):
+    def cell_hamiltonian(self, args=(), sparse=False, *, params=None):
         """Hamiltonian of a single cell of the infinite system."""
         cell_sites = range(self.cell_size)
         return self.hamiltonian_submatrix(args, cell_sites, cell_sites,
-                                          sparse=sparse)
+                                          sparse=sparse, params=params)
 
-    def inter_cell_hopping(self, args=(), sparse=False):
+    def inter_cell_hopping(self, args=(), sparse=False, *, params=None):
         """Hopping Hamiltonian between two cells of the infinite system."""
         cell_sites = range(self.cell_size)
         interface_sites = range(self.cell_size, self.graph.num_nodes)
         return self.hamiltonian_submatrix(args, cell_sites, interface_sites,
-                                          sparse=sparse)
+                                          sparse=sparse, params=params)
 
-    def modes(self, energy=0, args=()):
+    def modes(self, energy=0, args=(), *, params=None):
         """Return mode decomposition of the lead
 
         See documentation of `~kwant.physics.PropagatingModes` and
         `~kwant.physics.StabilizedModes` for the return format details.
         """
         from . import physics   # Putting this here avoids a circular import.
-        ham = self.cell_hamiltonian(args)
-        hop = self.inter_cell_hopping(args)
-        symmetries = self.discrete_symmetry(args)
+        ham = self.cell_hamiltonian(args, params=params)
+        hop = self.inter_cell_hopping(args, params=params)
+        symmetries = self.discrete_symmetry(args, params=params)
         broken = symmetries.validate(ham)
         if broken is not None:
             raise ValueError("Cell Hamiltonian breaks " + broken.lower())
@@ -240,7 +246,7 @@ class InfiniteSystem(System, metaclass=abc.ABCMeta):
             symmetries.particle_hole = symmetries.chiral = None
         return physics.modes(ham, hop, discrete_symmetry=symmetries)
 
-    def selfenergy(self, energy=0, args=()):
+    def selfenergy(self, energy=0, args=(), *, params=None):
         """Return self-energy of a lead.
 
         The returned matrix has the shape (s, s), where s is
@@ -248,13 +254,14 @@ class InfiniteSystem(System, metaclass=abc.ABCMeta):
         self.cell_size))``.
         """
         from . import physics   # Putting this here avoids a circular import.
-        ham = self.cell_hamiltonian(args)
+        ham = self.cell_hamiltonian(args, params=params)
         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))
+        return physics.selfenergy(ham,
+                                  self.inter_cell_hopping(args, params=params))
 
 
 class PrecalculatedLead:
@@ -277,7 +284,7 @@ class PrecalculatedLead:
         self._modes = modes
         self._selfenergy = selfenergy
 
-    def modes(self, energy=0, args=()):
+    def modes(self, energy=0, args=(), *, params=None):
         if self._modes is not None:
             return self._modes
         else:
@@ -285,7 +292,7 @@ class PrecalculatedLead:
                              "Consider using precalculate() with "
                              "what='modes' or what='all'")
 
-    def selfenergy(self, energy=0, args=()):
+    def selfenergy(self, energy=0, args=(), *, params=None):
         if self._selfenergy is not None:
             return self._selfenergy
         else:

kwant/tests/test_builder.py

@@ -964,3 +964,79 @@ def test_discrete_symmetries():
     p2 = np.zeros((4, 2))
     p2[1, 0] = p2[3, 1] = 1
     assert np.allclose(sym.projectors[1].todense(), p2)
+
+    # test parameter passing to conservation_law
+    syst = builder.Builder(conservation_law=lambda site, b: b)
+    syst[lat2(1)] = 0
+    sym = syst.finalized().discrete_symmetry(params=dict(a=None, b=1))
+    [proj] = sym.projectors
+    assert np.allclose(proj.todense(), [[1]])
+
+
+def test_argument_passing():
+
+    chain = kwant.lattice.chain()
+
+    # Test for passing parameters to hamiltonian matrix elements
+    def onsite(site, p1, p2=1):
+        return p1 + p2
+
+    def hopping(site1, site2, p1, p2=1):
+        return p1 - p2
+
+    def fill_syst(syst):
+        syst[(chain(i) for i in range(3))] = onsite
+        syst[chain.neighbors()] = hopping
+        return syst.finalized()
+
+    syst = fill_syst(kwant.Builder())
+    inf_syst = fill_syst(kwant.Builder(kwant.TranslationalSymmetry((-3,))))
+
+    args= (2, 1)
+    params = dict(p1=2, p2=1)
+
+    np.testing.assert_array_equal(
+        syst.hamiltonian_submatrix(args=args),
+        syst.hamiltonian_submatrix(params=params))
+    np.testing.assert_array_equal(
+        inf_syst.cell_hamiltonian(args=args),
+        inf_syst.cell_hamiltonian(params=params))
+    np.testing.assert_array_equal(
+        inf_syst.inter_cell_hopping(args=args),
+        inf_syst.inter_cell_hopping(params=params))
+    np.testing.assert_array_equal(
+        inf_syst.selfenergy(args=args),
+        inf_syst.selfenergy(params=params))
+    np.testing.assert_array_equal(
+        inf_syst.modes(args=args)[0].wave_functions,
+        inf_syst.modes(params=params)[0].wave_functions)
+
+    # test that mixing 'args' and 'params' raises TypeError
+    with raises(TypeError):
+        syst.hamiltonian(0, 0, *args, params=params)
+    with raises(TypeError):
+        inf_syst.hamiltonian(0, 0, *args, params=params)
+
+
+    # Some common, some different args for value functions
+    def onsite2(site, a, b):
+        return site.pos + a + b
+
+    def hopping2(site1, site2, a, c, b):
+        return a + b + c
+
+    syst = kwant.Builder()
+    syst[(chain(i) for i in range(3))] = onsite2
+    syst[((chain(i), chain(i + 1)) for i in range(2))] = hopping2
+    fsyst = syst.finalized()
+
+    def expected_hamiltonian(a, b, c):
+        return [[a + b, a + b + c, 0],
+                [a + b + c, 1 + a + b, a + b + c],
+                [0, a+ b + c, 2 + a + b]]
+
+    params = dict(a=1, b=2, c=3)
+    np.testing.assert_array_equal(
+        fsyst.hamiltonian_submatrix(params=params),
+        expected_hamiltonian(**params)
+    )

kwant/tests/test_operator.py

@@ -404,3 +404,83 @@ def test_tocoo():
                      check_hermiticity=False)
     op = op.bind(args=(1,))
     raises(ValueError, op.tocoo, [1])
+
+
+def test_arg_passing():
+    lat1 = kwant.lattice.chain(norbs=1)
+
+    syst = kwant.Builder()
+    syst[lat1(0)] = syst[lat1(1)] = lambda s0, a, b: s0.pos + a + b
+    syst[lat1.neighbors()] = lambda s0, s1, a, b: a - b
+    fsyst = syst.finalized()
+
+    wf = np.ones(len(fsyst.sites))
+
+    for A in opservables:
+        op = A(fsyst)
+        canonical_args = (1, 2)
+        params = dict(a=1, b=2)
+        call_should_be = op(wf, args=canonical_args)
+        act_should_be = op.act(wf, args=canonical_args)
+        has_tocoo = hasattr(op, 'tocoo')
+        if has_tocoo:
+            tocoo_should_be = op.tocoo(args=canonical_args).todense()
+
+        with raises(TypeError) as exc:
+            op(wf, args=canonical_args, params=params)
+        assert 'mutually exclusive' in str(exc)
+        with raises(TypeError) as exc:
+            op.act(wf, args=canonical_args, params=params)
+        assert 'mutually exclusive' in str(exc)
+        with raises(TypeError) as exc:
+            op.bind(args=canonical_args, params=params)
+        assert 'mutually exclusive' in str(exc)
+        if has_tocoo:
+            with raises(TypeError) as exc:
+                op.tocoo(args=canonical_args, params=params)
+            assert 'mutually exclusive' in str(exc)
+
+        np.testing.assert_array_equal(
+            call_should_be, op(wf, params=params))
+        np.testing.assert_array_equal(
+            act_should_be, op.act(wf, params=params))
+        if has_tocoo:
+            np.testing.assert_array_equal(
+                tocoo_should_be, op.tocoo(params=params).todense())
+        # after binding
+        op2 = op.bind(params=params)
+        np.testing.assert_array_equal(
+            call_should_be, op2(wf))
+        np.testing.assert_array_equal(
+            act_should_be, op2.act(wf))
+        if has_tocoo:
+            np.testing.assert_array_equal(
+                tocoo_should_be, op2.tocoo().todense())
+
+        # system and onsite having different args
+        def onsite(site, flip):
+            return -1 if flip else 1
+
+        op = A(fsyst, onsite=onsite)
+        params['flip'] = True
+        call_should_be = -call_should_be
+        act_should_be = -act_should_be
+        if has_tocoo:
+            tocoo_should_be = -tocoo_should_be
+
+        np.testing.assert_array_equal(
+            call_should_be, op(wf, params=params))
+        np.testing.assert_array_equal(
+            act_should_be, op.act(wf, params=params))
+        if has_tocoo:
+            np.testing.assert_array_equal(
+                tocoo_should_be, op.tocoo(params=params).todense())
+        # after binding
+        op2 = op.bind(params=params)
+        np.testing.assert_array_equal(
+            call_should_be, op2(wf))
+        np.testing.assert_array_equal(
+            act_should_be, op2.act(wf))
+        if has_tocoo:
+            np.testing.assert_array_equal(
+                tocoo_should_be, op2.tocoo().todense())

kwant/tests/test_system.py

@@ -83,27 +83,6 @@ def test_hamiltonian_submatrix():
     raises(ValueError, syst2.hamiltonian_submatrix)
     raises(ValueError, syst2.hamiltonian_submatrix, sparse=True)
 
-    # Test for passing parameters to hamiltonian matrix elements
-    def onsite(site, p1, p2=0):
-        return site.pos + p1 + p2
-
-    def hopping(site1, site2, p1, p2=0):
-        return p1 - p2
-
-    syst = kwant.Builder()
-    syst[(chain(i) for i in range(3))] = onsite
-    syst[((chain(i), chain(i + 1)) for i in range(2))] = hopping
-    syst2 = syst.finalized()
-    mat = syst2.hamiltonian_submatrix((2, 1))
-    mat_should_be = [[3, 1, 0], [1, 4, 1], [0, 1, 5]]
-
-    # Sorting is required due to unknown compression order of builder.
-    onsite_hamiltonians = mat.flat[::4]
-    perm = np.argsort(onsite_hamiltonians)
-    mat = mat[perm, :]
-    mat = mat[:, perm]
-    np.testing.assert_array_equal(mat, mat_should_be)
-
 
 def test_pickling():
     syst = kwant.Builder()