diff --git a/doc/source/code/figure/ab_ring.py.diff b/doc/source/code/figure/ab_ring.py.diff
index 50a10d1920c488552360b3bfc388d661da86d88f..c4fb178b462033ee9f0cdc906c4cfb5f349e403b 100644
--- a/doc/source/code/figure/ab_ring.py.diff
+++ b/doc/source/code/figure/ab_ring.py.diff
@@ -143,7 +143,7 @@
      normalized_fluxes = [flux / (2 * pi) for flux in fluxes]
      data = []
      for flux in fluxes:
-         smatrix = kwant.smatrix(syst, energy, args=[flux])
+         smatrix = kwant.smatrix(syst, energy, params=dict(phi=flux))
          data.append(smatrix.transmission(1, 0))
  
 -    pyplot.figure()
diff --git a/doc/source/code/figure/closed_system.py.diff b/doc/source/code/figure/closed_system.py.diff
index 5f9cc10122b239ebd33488bbc433868dbba02a55..6ad93be562a7e793ca1f825fa8af2ec5641ef114 100644
--- a/doc/source/code/figure/closed_system.py.diff
+++ b/doc/source/code/figure/closed_system.py.diff
@@ -68,7 +68,7 @@
      energies = []
      for B in Bfields:
          # Obtain the Hamiltonian as a dense matrix
-         ham_mat = syst.hamiltonian_submatrix(args=[B], sparse=True)
+         ham_mat = syst.hamiltonian_submatrix(params=dict(B=B), sparse=True)
  
          # we only calculate the 15 lowest eigenvalues
          ev = sla.eigsh(ham_mat.tocsc(), k=15, sigma=0,
@@ -105,7 +105,7 @@
 +    size = (_defs.figwidth_in, _defs.figwidth_in)
 +
      # Calculate the wave functions in the system.
-     ham_mat = syst.hamiltonian_submatrix(sparse=True, args=[B])
+     ham_mat = syst.hamiltonian_submatrix(sparse=True, params=dict(B=B))
      evals, evecs = sorted_eigs(sla.eigsh(ham_mat.tocsc(), k=20, sigma=0))
  
      # Plot the probability density of the 10th eigenmode.
@@ -124,12 +124,12 @@
 +    size = (_defs.figwidth_in, _defs.figwidth_in)
 +
      # Calculate the wave functions in the system.
-     ham_mat = syst.hamiltonian_submatrix(sparse=True, args=[B])
+     ham_mat = syst.hamiltonian_submatrix(sparse=True, params=dict(B=B))
      evals, evecs = sorted_eigs(sla.eigsh(ham_mat.tocsc(), k=20, sigma=0))
  
      # Calculate and plot the local current of the 10th eigenmode.
      J = kwant.operator.Current(syst)
-     current = J(evecs[:, 9], args=[B])
+     current = J(evecs[:, 9], params=dict(B=B))
 -    kwant.plotter.current(syst, current, colorbar=False)
 +    for extension in ('pdf', 'png'):
 +        kwant.plotter.current(
diff --git a/doc/source/code/figure/quantum_well.py.diff b/doc/source/code/figure/quantum_well.py.diff
index 1c81f1c54c37e24e231b15d6144c3cac7207e54e..8b1be4bb17cd8513816464e1782381842961ab88 100644
--- a/doc/source/code/figure/quantum_well.py.diff
+++ b/doc/source/code/figure/quantum_well.py.diff
@@ -59,7 +59,7 @@
      # Compute conductance
      data = []
      for welldepth in welldepths:
-         smatrix = kwant.smatrix(syst, energy, args=[-welldepth])
+         smatrix = kwant.smatrix(syst, energy, params=dict(pot=-welldepth))
          data.append(smatrix.transmission(1, 0))
  
 -    pyplot.figure()
diff --git a/doc/source/pre/whatsnew/1.4.rst b/doc/source/pre/whatsnew/1.4.rst
index e29c11afcf895e160b60bf18d0d8b9f66da69257..33ec8a92589754f355eca696b135bbc9c6e36500 100644
--- a/doc/source/pre/whatsnew/1.4.rst
+++ b/doc/source/pre/whatsnew/1.4.rst
@@ -170,6 +170,27 @@ parameters that the system (and its leads) expects::
 
 This is a provisional API that may be changed in a future version of Kwant.
 
+Passing system arguments via ``args`` is deprecated in favor of ``params``
+--------------------------------------------------------------------------
+It is now deprecated to pass arguments to systems by providing the
+``args`` parameter (in ``kwant.smatrix`` and elsewhere). This is
+error prone and requires that all value functions take the same
+formal parameters, even if they do not depend on all of them. The
+preferred way of passing parameters to Kwant systems is by passing
+a dictionary using ``params``::
+
+  def onsite(site, magnetic_field, voltage):
+    return magnetic_field * sigma_z + voltage * sigma_0
+
+  syst = make_system(onsite).finalized()
+
+  kwant.smatrix(syst, params=dict(magnetic_field=0.5, voltage=0.2))
+
+  # Compare this to the deprecated 'args'
+  kwant.smatrix(syst, args=(0.5, 0.2))
+
+The ability to provide ``args`` will be removed in a future Kwant version.
+
 Interpolated density plots
 --------------------------
 A new function `~kwant.plotter.density` has been added that can be used to
diff --git a/doc/source/tutorial/discretize.rst b/doc/source/tutorial/discretize.rst
index 1ed5015cd2203e100d1d4e4cf9f21ae283716f04..b5edb0157f8a4ba5c3542b525328a48d852a8d42 100644
--- a/doc/source/tutorial/discretize.rst
+++ b/doc/source/tutorial/discretize.rst
@@ -160,9 +160,9 @@ energy eigenstates:
 
 .. image:: /code/figure/discretizer_gs.*
 
-Note in the above that we provided the function ``V`` to
-``syst.hamiltonian_submatrix`` using ``params=dict(V=potential)``, rather than
-via ``args``.
+Note in the above that we pass the spatially varying potential *function*
+to our system via a parameter called ``V``, because the symbol $V$
+was used in the intial, symbolic, definition of the Hamiltonian.
 
 In addition, the function passed as ``V`` expects two input parameters ``x``
 and ``y``, the same as in the initial continuum Hamiltonian.
diff --git a/doc/source/tutorial/spin_potential_shape.rst b/doc/source/tutorial/spin_potential_shape.rst
index bcf95dacc7627cf3822474009480287802cc6889..02ac1fbad566ad8daa570f26677e86c3d83fd01d 100644
--- a/doc/source/tutorial/spin_potential_shape.rst
+++ b/doc/source/tutorial/spin_potential_shape.rst
@@ -184,10 +184,12 @@ Finally, we compute the transmission probability:
     :start-after: #HIDDEN_BEGIN_sqvr
     :end-before: #HIDDEN_END_sqvr
 
-``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:
+``kwant.smatrix`` allows us to specify a dictionary, `params`, that contains
+the additional arguments required by 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 (remember above
+we defined our `onsite` function that takes a parameter named `pot`).
+We obtain the result:
 
 .. image:: /code/figure/quantum_well_result.*
 
diff --git a/kwant/_common.py b/kwant/_common.py
index 71b4f7826e3d963278a72b199fc9827227fee141..e4623b8d09c0d8a307fcab8dff2b6571fb4112c3 100644
--- a/kwant/_common.py
+++ b/kwant/_common.py
@@ -12,6 +12,7 @@ import numbers
 import inspect
 import warnings
 import importlib
+import functools
 from contextlib import contextmanager
 
 __all__ = ['KwantDeprecationWarning', 'UserCodeError']
@@ -41,6 +42,46 @@ class UserCodeError(Exception):
     pass
 
 
+def deprecate_parameter(parameter_name, version=None, help=None,
+                        stacklevel=3):
+    """Trigger a deprecation warning if the wrapped function is called
+       with the provided parameter."""
+
+    message = ("The '{}' parameter has been deprecated since version {} -- {}"
+               .format(parameter_name, version, help))
+
+    def warn():
+        warnings.warn(message, KwantDeprecationWarning,
+                      stacklevel=stacklevel)
+
+    def wrapper(f=None):
+
+        # Instead of being used as a decorator, can be called with
+        # no arguments to just raise the warning.
+        if f is None:
+            warn()
+            return
+
+        sig = inspect.signature(f)
+
+        @functools.wraps(f)
+        def inner(*args, **kwargs):
+            # If the named argument is truthy
+            if sig.bind(*args, **kwargs).arguments.get(parameter_name):
+                warn()
+            return f(*args, **kwargs)
+
+        return inner
+
+    return wrapper
+
+
+# Deprecation for 'args' parameter; defined once to minimize boilerplate,
+# as this parameter is present all over Kwant.
+deprecate_args = deprecate_parameter('args', version=1.4, help=
+    "Instead, provide named parameters as a dictionary via 'params'.")
+
+
 def interleave(seq):
     """Return an iterator that yields pairs of elements from a sequence.
 
diff --git a/kwant/_system.pyx b/kwant/_system.pyx
index d297dbc1961b5b49e8c1af34f283d33f29746aea..45e0da438e6b9a7322307331982ea76d8cf272c7 100644
--- a/kwant/_system.pyx
+++ b/kwant/_system.pyx
@@ -16,6 +16,7 @@ import types
 from .graph.core cimport CGraph, gintArraySlice
 from .graph.defs cimport gint
 from .graph.defs import gint_dtype
+from ._common import deprecate_args
 
 msg = ('Hopping from site {0} to site {1} does not match the '
        'dimensions of onsite Hamiltonians of these sites.')
@@ -253,6 +254,7 @@ def _check_parameters_match(expected_parameters, params):
         raise TypeError(''.join(msg))
 
 
+@deprecate_args
 @cython.binding(True)
 @cython.embedsignature(True)
 def hamiltonian_submatrix(self, args=(), to_sites=None, from_sites=None,
diff --git a/kwant/builder.py b/kwant/builder.py
index c804c3a8f07e9617f0e953f7036c81f467fbb23c..5b89731ae4c14f1c9d585c65b78a40d4f1110310 100644
--- a/kwant/builder.py
+++ b/kwant/builder.py
@@ -22,7 +22,7 @@ from .linalg import lll
 from .operator import Density
 from .physics import DiscreteSymmetry
 from ._common import (ensure_isinstance, get_parameters, reraise_warnings,
-                      interleave)
+                      interleave, deprecate_args)
 
 
 __all__ = ['Builder', 'Site', 'SiteFamily', 'SimpleSiteFamily', 'Symmetry',
@@ -617,6 +617,7 @@ def _ensure_signature(func):
         return func
 
     # function conforming to old API: needs wrapping
+    @deprecate_args
     def wrapper(energy, args=(), *, params=None):
         return func(energy, args)
 
@@ -647,6 +648,7 @@ class SelfEnergyLead(Lead):
         """Trivial finalization: the object is returned itself."""
         return self
 
+    @deprecate_args
     def selfenergy(self, energy, args=(), *, params=None):
         return self.selfenergy_func(energy, args, params=params)
 
@@ -677,9 +679,11 @@ class ModesLead(Lead):
         """Trivial finalization: the object is returned itself."""
         return self
 
+    @deprecate_args
     def modes(self, energy, args=(), *, params=None):
         return self.modes_func(energy, args, params=params)
 
+    @deprecate_args
     def selfenergy(self, energy, args=(), *, params=None):
         stabilized = self.modes(energy, args, params=params)[1]
         return stabilized.selfenergy()
@@ -1894,6 +1898,7 @@ class _FinalizedBuilderMixin:
                 value = herm_conj(value)
         return value
 
+    @deprecate_args
     def discrete_symmetry(self, args=(), *, params=None):
         if self._cons_law is not None:
             eigvals, eigvecs = self._cons_law
diff --git a/kwant/operator.pyx b/kwant/operator.pyx
index 0b1a8e279f5b3a95609e61dd23ddb183474889c5..002eb8b66e188aa95a44b3e5c6248187dfaa6319 100644
--- a/kwant/operator.pyx
+++ b/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, get_parameters
+from ._common import UserCodeError, get_parameters, deprecate_args
 
 
 ################ Generic Utility functions
@@ -499,7 +499,7 @@ 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'.
+            Deprecated in favor of 'params' (and mutually exclusive with it).
         params : dict, optional
             Dictionary of parameter names and their values. Mutually exclusive
             with 'args'.
@@ -514,6 +514,11 @@ cdef class _LocalOperator:
             raise ValueError("Extra arguments are already bound to this "
                              "operator. You should call this operator "
                              "providing neither 'args' nor 'params'.")
+        if args:
+            # deprecate_args does not play nicely with methods of cdef classes,
+            # when used as a decorator, so we manually raise the
+            # deprecation warning here.
+            deprecate_args()
         if args and params:
             raise TypeError("'args' and 'params' are mutually exclusive.")
         if bra is None:
@@ -555,7 +560,8 @@ 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. Mutually exclusive with 'params'.
+            the operator ``onsite`` function.
+            Deprecated in favor of 'params' (and mutually exclusive with it).
         params : dict, optional
             Dictionary of parameter names and their values. Mutually exclusive
             with 'args'.
@@ -568,6 +574,11 @@ cdef class _LocalOperator:
             raise ValueError("Extra arguments are already bound to this "
                              "operator. You should call this operator "
                              "providing neither 'args' nor 'params'.")
+        if args:
+            # deprecate_args does not play nicely with methods of cdef classes,
+            # when used as a decorator, so we manually raise the
+            # deprecation warning here.
+            deprecate_args()
         if args and params:
             raise TypeError("'args' and 'params' are mutually exclusive.")
 
@@ -588,7 +599,15 @@ cdef class _LocalOperator:
 
         Returns a copy of this operator that does not need to be passed extra
         arguments when subsequently called or when using the ``act`` method.
+
+        Providing positional arguments via 'args' is deprecated,
+        instead provide named parameters as a dictionary via 'params'.
         """
+        if args:
+            # deprecate_args does not play nicely with methods of cdef classes,
+            # when used as a decorator, so we manually raise the
+            # deprecation warning here.
+            deprecate_args()
         if args and params:
             raise TypeError("'args' and 'params' are mutually exclusive.")
         # generic creation of new instance
@@ -621,7 +640,8 @@ 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. Mutually exclusive with 'params'.
+            the operator ``onsite`` function.
+            Deprecated in favor of 'params' (and mutually exclusive with it).
         op : operation
             The operation to perform.
             `MAT_ELS`: calculate matrix elements between `bra` and `ket`
@@ -786,7 +806,11 @@ cdef class Density(_LocalOperator):
     @cython.cdivision(True)
     @cython.embedsignature
     def tocoo(self, args=(), *, params=None):
-        """Convert the operator to coordinate format sparse matrix."""
+        """Convert the operator to coordinate format sparse matrix.
+
+        Providing positional arguments via 'args' is deprecated,
+        instead provide named parameters as a dictionary via 'params'.
+        """
         cdef int blk, blk_size, n_blocks, n, k = 0
         cdef int [:, :] offsets, shapes
         cdef int [:] row, col
@@ -794,6 +818,11 @@ cdef class Density(_LocalOperator):
            raise ValueError("Extra arguments are already bound to this "
                             "operator. You should call this operator "
                             "providing neither 'args' nor 'params'.")
+        if args:
+            # deprecate_args does not play nicely with methods of cdef classes,
+            # when used as a decorator, so we manually raise the
+            # deprecation warning here.
+            deprecate_args()
         if args and params:
             raise TypeError("'args' and 'params' are mutually exclusive.")
 
@@ -887,6 +916,9 @@ cdef class Current(_LocalOperator):
 
         Returns a copy of this operator that does not need to be passed extra
         arguments when subsequently called or when using the ``act`` method.
+
+        Providing positional arguments via 'args' is deprecated,
+        instead provide named parameters as a dictionary via 'params'.
         """
         q = super().bind(args, params=params)
         q._bound_hamiltonian = self._eval_hamiltonian(args, params)
@@ -1012,6 +1044,9 @@ cdef class Source(_LocalOperator):
 
         Returns a copy of this operator that does not need to be passed extra
         arguments when subsequently called or when using the ``act`` method.
+
+        Providing positional arguments via 'args' is deprecated,
+        instead provide named parameters as a dictionary via 'params'.
         """
         q = super().bind(args, params=params)
         q._bound_hamiltonian = self._eval_hamiltonian(args, params)
diff --git a/kwant/physics/dispersion.py b/kwant/physics/dispersion.py
index 5c98bf2618536fea4262df0ef8476dcc9181e721..f35a7781695b16902a4d0ec781e2656bebcfe073 100644
--- a/kwant/physics/dispersion.py
+++ b/kwant/physics/dispersion.py
@@ -11,7 +11,7 @@
 import math
 import numpy as np
 from .. import system
-from .._common import ensure_isinstance
+from .._common import ensure_isinstance, deprecate_args
 
 __all__ = ['Bands']
 
@@ -27,7 +27,7 @@ class Bands:
         calculated.
     args : tuple, defaults to empty
         Positional arguments to pass to the ``hamiltonian`` method.
-        Mutually exclusive with 'params'.
+        Deprecated in favor or 'params' (and mutually exclusive with it).
     params : dict, optional
         Dictionary of parameter names and their values. Mutually exclusive
         with 'args'.
@@ -50,6 +50,7 @@ class Bands:
     """
     _crossover_size = 8
 
+    @deprecate_args
     def __init__(self, sys, args=(), *, params=None):
         syst = sys
         ensure_isinstance(syst, system.InfiniteSystem)
diff --git a/kwant/plotter.py b/kwant/plotter.py
index 4652271cf1870568294798445cc9a6495e1958c7..35813dd325eb1febc3e9a573a44d5b68fef13de8 100644
--- a/kwant/plotter.py
+++ b/kwant/plotter.py
@@ -27,6 +27,7 @@ from scipy import spatial, interpolate
 from math import cos, sin, pi, sqrt
 
 from . import system, builder, _common
+from ._common import deprecate_args
 
 
 __all__ = ['plot', 'map', 'bands', 'spectrum', 'current', 'density',
@@ -1348,6 +1349,7 @@ def map(sys, value, colorbar=True, cmap=None, vmin=None, vmax=None, a=None,
     return fig
 
 
+@deprecate_args
 def bands(sys, args=(), momenta=65, file=None, show=True, dpi=None,
           fig_size=None, ax=None, *, params=None):
     """Plot band structure of a translationally invariant 1D system.
@@ -1358,7 +1360,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'.
+        Deprecated in favor of 'params' (and mutually exclusive with it).
     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.
diff --git a/kwant/solvers/common.py b/kwant/solvers/common.py
index 2c3f87a6154fed52fee0153d267f039083a2ef0d..ecf33f80003c6aee65ea5c90c3c75b0f47bc469b 100644
--- a/kwant/solvers/common.py
+++ b/kwant/solvers/common.py
@@ -14,7 +14,7 @@ import abc
 from  numbers import Integral
 import numpy as np
 import scipy.sparse as sp
-from .._common import ensure_isinstance
+from .._common import ensure_isinstance, deprecate_args
 from .. import system
 from functools import reduce
 
@@ -113,7 +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'.
+            Deprecated in favor of 'params' (and mutually exclusive with it).
         check_hermiticity : bool
             Check if Hamiltonian matrices are in fact Hermitian.
             Enables deduction of missing transmission coefficients.
@@ -296,6 +296,7 @@ class SparseSolver(metaclass=abc.ABCMeta):
 
         return LinearSys(lhs, rhs, indices, num_orb), lead_info
 
+    @deprecate_args
     def smatrix(self, sys, energy=0, args=(),
                 out_leads=None, in_leads=None, check_hermiticity=True,
                 *, params=None):
@@ -313,7 +314,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'.
+            Deprecated in favor of 'params' (and mutually exclusive with it).
         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
@@ -390,6 +391,7 @@ class SparseSolver(metaclass=abc.ABCMeta):
 
         return SMatrix(data, lead_info, out_leads, in_leads, check_hermiticity)
 
+    @deprecate_args
     def greens_function(self, sys, energy=0, args=(),
                         out_leads=None, in_leads=None, check_hermiticity=True,
                         *, params=None):
@@ -407,7 +409,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'.
+            Deprecated in favor of 'params' (and mutually exclusive with it).
         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
@@ -488,6 +490,7 @@ class SparseSolver(metaclass=abc.ABCMeta):
         return GreensFunction(data, lead_info, out_leads, in_leads,
                               check_hermiticity)
 
+    @deprecate_args
     def ldos(self, sys, energy=0, args=(), check_hermiticity=True,
              *, params=None):
         """
@@ -504,8 +507,8 @@ 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.  Mutually exclusive
-            with 'params'.
+            evaluate the hamiltonian matrix elements.
+            Deprecated in favor of 'params' (and mutually exclusive with it).
         check_hermiticity : ``bool``
             Check if the Hamiltonian matrices are Hermitian.
         params : dict, optional
@@ -553,6 +556,7 @@ class SparseSolver(metaclass=abc.ABCMeta):
 
         return ldos * (0.5 / np.pi)
 
+    @deprecate_args
     def wave_function(self, sys, energy=0, args=(), check_hermiticity=True,
                       *, params=None):
         r"""
@@ -568,8 +572,8 @@ 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. Mutually exclusive
-            with 'params'.
+            evaluate the hamiltonian matrix elements.
+            Deprecated in favor of 'params' (and mutually exclusive with it).
         check_hermiticity : ``bool``
             Check if the Hamiltonian matrices are Hermitian.
         params : dict, optional
diff --git a/kwant/solvers/tests/test_mumps.py b/kwant/solvers/tests/test_mumps.py
index 9ca695228c0382f6e70b02ffded1956a7c9d615b..b56798e6beb3dea1c13fdbff8753b6f4cb3fe5c4 100644
--- a/kwant/solvers/tests/test_mumps.py
+++ b/kwant/solvers/tests/test_mumps.py
@@ -114,5 +114,8 @@ def test_wavefunc_ldos_consistency():
         _test_sparse.test_wavefunc_ldos_consistency(wave_function, ldos)
 
 
+# We need to keep testing 'args', but we don't want to see
+# all the deprecation warnings in the test logs
+@pytest.mark.filterwarnings("ignore:.*'args' parameter")
 def test_arg_passing():
     _test_sparse.test_arg_passing(wave_function, ldos, smatrix)
diff --git a/kwant/solvers/tests/test_sparse.py b/kwant/solvers/tests/test_sparse.py
index a98f9a44d61859b062b789165422e3e2c8c69649..2ddab9b1cb3c3388b32cb25764d825792670dec2 100644
--- a/kwant/solvers/tests/test_sparse.py
+++ b/kwant/solvers/tests/test_sparse.py
@@ -6,6 +6,8 @@
 # the file AUTHORS.rst at the top-level directory of this distribution and at
 # http://kwant-project.org/authors.
 
+import pytest
+
 from  kwant.solvers.sparse import smatrix, greens_function, ldos, wave_function
 from . import _test_sparse
 
@@ -60,5 +62,9 @@ def test_ldos():
 def test_wavefunc_ldos_consistency():
     _test_sparse.test_wavefunc_ldos_consistency(wave_function, ldos)
 
+
+# We need to keep testing 'args', but we don't want to see
+# all the deprecation warnings in the test logs
+@pytest.mark.filterwarnings("ignore:.*'args' parameter")
 def test_arg_passing():
     _test_sparse.test_arg_passing(wave_function, ldos, smatrix)
diff --git a/kwant/system.py b/kwant/system.py
index 26eea6fa4bc96b38850b9a2a43e8e5ffc2785282..d270691508189866a0d75809917048303443c883 100644
--- a/kwant/system.py
+++ b/kwant/system.py
@@ -14,6 +14,7 @@ import abc
 import warnings
 from copy import copy
 from . import _system
+from ._common  import deprecate_args
 
 
 class System(metaclass=abc.ABCMeta):
@@ -59,12 +60,20 @@ class System(metaclass=abc.ABCMeta):
         if ``i != j``, return the hopping between site ``i`` and ``j``.
 
         Hamiltonians may depend (optionally) on positional and
-        keyword arguments
+        keyword arguments.
+
+        Providing positional arguments via 'args' is deprecated,
+        instead, provide named parameters as a dictionary via 'params'.
         """
         pass
 
+    @deprecate_args
     def discrete_symmetry(self, args, *, params=None):
-        """Return the discrete symmetry of the system."""
+        """Return the discrete symmetry of the system.
+
+        Providing positional arguments via 'args' is deprecated,
+        instead, provide named parameters as a dictionary via 'params'.
+        """
         # Avoid the circular import.
         from .physics import DiscreteSymmetry
         return DiscreteSymmetry()
@@ -131,6 +140,7 @@ class FiniteSystem(System, metaclass=abc.ABCMeta):
 
     """
 
+    @deprecate_args
     def precalculate(self, energy=0, args=(), leads=None,
                      what='modes', *, params=None):
         """
@@ -147,7 +157,7 @@ class FiniteSystem(System, metaclass=abc.ABCMeta):
             evaluated.
         args : sequence
             Additional parameters required for calculating the Hamiltionians.
-            Mutually exclusive with 'params'.
+            Deprecated in favor of 'params' (and mutually exclusive with it).
         leads : sequence of integers or None
             Numbers of the leads to be precalculated. If ``None``, all are
             precalculated.
@@ -194,12 +204,16 @@ class FiniteSystem(System, metaclass=abc.ABCMeta):
         result.leads = new_leads
         return result
 
+    @deprecate_args
     def validate_symmetries(self, args=(), *, params=None):
         """Check that the Hamiltonian satisfies discrete symmetries.
 
         Applies `~kwant.physics.DiscreteSymmetry.validate` to the
         Hamiltonian, see its documentation for details on the return
         format.
+
+        Providing positional arguments via 'args' is deprecated,
+        instead, provide named parameters as a dictionary via 'params'.
         """
         symmetries = self.discrete_symmetry(args=args, params=params)
         ham = self.hamiltonian_submatrix(args, sparse=True, params=params)
@@ -248,19 +262,30 @@ class InfiniteSystem(System, metaclass=abc.ABCMeta):
     exchanged, as well as of site 3 and 4.
 
     """
+    @deprecate_args
     def cell_hamiltonian(self, args=(), sparse=False, *, params=None):
-        """Hamiltonian of a single cell of the infinite system."""
+        """Hamiltonian of a single cell of the infinite system.
+
+        Providing positional arguments via 'args' is deprecated,
+        instead, provide named parameters as a dictionary via 'params'.
+        """
         cell_sites = range(self.cell_size)
         return self.hamiltonian_submatrix(args, cell_sites, cell_sites,
                                           sparse=sparse, params=params)
 
+    @deprecate_args
     def inter_cell_hopping(self, args=(), sparse=False, *, params=None):
-        """Hopping Hamiltonian between two cells of the infinite system."""
+        """Hopping Hamiltonian between two cells of the infinite system.
+
+        Providing positional arguments via 'args' is deprecated,
+        instead, provide named parameters as a dictionary via 'params'.
+        """
         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, params=params)
 
+    @deprecate_args
     def modes(self, energy=0, args=(), *, params=None):
         """Return mode decomposition of the lead
 
@@ -272,6 +297,9 @@ class InfiniteSystem(System, metaclass=abc.ABCMeta):
         freedom on the first ``cell_sites`` sites of the system
         (recall that infinite systems store first the sites in the unit
         cell, then connected sites in the neighboring unit cell).
+
+        Providing positional arguments via 'args' is deprecated,
+        instead, provide named parameters as a dictionary via 'params'.
         """
         from . import physics   # Putting this here avoids a circular import.
         ham = self.cell_hamiltonian(args, params=params)
@@ -301,12 +329,16 @@ class InfiniteSystem(System, metaclass=abc.ABCMeta):
             symmetries.particle_hole = symmetries.chiral = None
         return physics.modes(ham, hop, discrete_symmetry=symmetries)
 
+    @deprecate_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
         ``sum(len(self.hamiltonian(i, i)) for i in range(self.graph.num_nodes -
         self.cell_size))``.
+
+        Providing positional arguments via 'args' is deprecated,
+        instead, provide named parameters as a dictionary via 'params'.
         """
         from . import physics   # Putting this here avoids a circular import.
         ham = self.cell_hamiltonian(args, params=params)
@@ -318,12 +350,16 @@ class InfiniteSystem(System, metaclass=abc.ABCMeta):
         return physics.selfenergy(ham,
                                   self.inter_cell_hopping(args, params=params))
 
+    @deprecate_args
     def validate_symmetries(self, args=(), *, params=None):
         """Check that the Hamiltonian satisfies discrete symmetries.
 
         Returns `~kwant.physics.DiscreteSymmetry.validate` applied
         to the onsite matrix and the hopping. See its documentation for
         details on the return format.
+
+        Providing positional arguments via 'args' is deprecated,
+        instead, provide named parameters as a dictionary via 'params'.
         """
         symmetries = self.discrete_symmetry(args=args, params=params)
         ham = self.cell_hamiltonian(args=args, sparse=True, params=params)
@@ -355,6 +391,7 @@ class PrecalculatedLead:
         # is no parametric dependence anymore
         self.parameters = frozenset()
 
+    @deprecate_args
     def modes(self, energy=0, args=(), *, params=None):
         if self._modes is not None:
             return self._modes
@@ -363,6 +400,7 @@ class PrecalculatedLead:
                              "Consider using precalculate() with "
                              "what='modes' or what='all'")
 
+    @deprecate_args
     def selfenergy(self, energy=0, args=(), *, params=None):
         if self._selfenergy is not None:
             return self._selfenergy
diff --git a/kwant/tests/test_builder.py b/kwant/tests/test_builder.py
index cb2cb9c7f0d5dbe968b8f4e65fc22ffbbc0e5751..c6f0874865968079ae24b1e8c8f373fd674acd71 100644
--- a/kwant/tests/test_builder.py
+++ b/kwant/tests/test_builder.py
@@ -14,6 +14,7 @@ from random import Random
 
 import numpy as np
 import tinyarray as ta
+import pytest
 from pytest import raises, warns
 from numpy.testing import assert_almost_equal
 
@@ -1145,12 +1146,14 @@ def test_discrete_symmetries():
     syst[lat(1)] = np.identity(2)
     syst[lat2(1)] = 1
 
-    sym = syst.finalized().discrete_symmetry(args=[0])
+    params=dict(p=0)
+
+    sym = syst.finalized().discrete_symmetry(params=params)
     for proj, should_be in zip(sym.projectors, np.identity(3)):
         assert np.allclose(proj.toarray(), should_be.reshape((3, 1)))
     assert np.allclose(sym.time_reversal.toarray(), np.identity(3))
     syst.conservation_law = lambda site, p: cons_law[site.family]
-    sym = syst.finalized().discrete_symmetry(args=[0])
+    sym = syst.finalized().discrete_symmetry(params=params)
     for proj, should_be in zip(sym.projectors, np.identity(3)):
         assert np.allclose(proj.toarray(), should_be.reshape((-1, 1)))
 
@@ -1186,8 +1189,10 @@ def test_discrete_symmetries():
     assert np.allclose(proj.toarray(), [[1]])
 
 
+# We need to keep testing 'args', but we don't want to see
+# all the deprecation warnings in the test logs
+@pytest.mark.filterwarnings("ignore:.*'args' parameter")
 def test_argument_passing():
-
     chain = kwant.lattice.chain()
 
     # Test for passing parameters to hamiltonian matrix elements
diff --git a/kwant/tests/test_comprehensive.py b/kwant/tests/test_comprehensive.py
index 088aefbb6cf68856dba6a2e64cc8fdf4616fc4d7..9989bf59fca3d7f9d76a2f0d862e064cdca26316 100644
--- a/kwant/tests/test_comprehensive.py
+++ b/kwant/tests/test_comprehensive.py
@@ -48,8 +48,8 @@ def test_qhe(W=16, L=8):
                                        ((3.2, 3.7), 2, 1e-3),
                                        ((5.2, 5.5), 3, 1e-1)]:
         for r_phi in r_phis:
-            args = (1.0 / r_phi, "")
-            pc = syst.precalculate(1.0, args, what='all')
-            for result in [kwant.smatrix(pc, 1, args),
-                           kwant.solvers.default.greens_function(pc, 1, args)]:
+            params = dict(phi=1.0 / r_phi, salt="")
+            pc = syst.precalculate(1.0, params=params, what='all')
+            for result in [kwant.smatrix(pc, 1, params=params),
+                           kwant.solvers.default.greens_function(pc, 1, params=params)]:
                 assert abs(T_nominal - result.transmission(1, 0)) < max_err
diff --git a/kwant/tests/test_operator.py b/kwant/tests/test_operator.py
index cc76209e3b058243ddcb4f00929942737a91117c..6eb0e631a3d0794846bbd078898a27b4c829cab4 100644
--- a/kwant/tests/test_operator.py
+++ b/kwant/tests/test_operator.py
@@ -156,22 +156,22 @@ def test_operator_construction():
         A(fsyst, sum=True).sum == True
 
 
-def _test(A, bra, ket=None, per_el_val=None, reduced_val=None, args=()):
+def _test(A, bra, ket=None, per_el_val=None, reduced_val=None, params=None):
     if per_el_val is not None:
-        val = A(bra, ket, args=args)
+        val = A(bra, ket, params=params)
         assert np.allclose(val, per_el_val)
         # with bound args
-        val = A.bind(args)(bra, ket)
+        val = A.bind(params=params)(bra, ket)
         assert np.allclose(val, per_el_val)
     # test that inner products give the same thing
     ket = bra if ket is None else ket
-    act_val = np.dot(bra.conj(), A.act(ket, args=args))
-    inner_val = np.sum(A(bra, ket, args=args))
+    act_val = np.dot(bra.conj(), A.act(ket, params=params))
+    inner_val = np.sum(A(bra, ket, params=params))
     # check also when sum is done internally by operator
     try:
         sum_reset = A.sum
         A.sum = True
-        sum_inner_val = A(bra, ket, args=args)
+        sum_inner_val = A(bra, ket, params=params)
         assert inner_val == sum_inner_val
     finally:
         A.sum = sum_reset
@@ -300,22 +300,22 @@ def test_opservables_spin():
     syst.attach_lead(lead)
     syst.attach_lead(lead.reversed())
     fsyst = syst.finalized()
-    args = (0.1,)
-    down, up = kwant.wave_function(fsyst, energy=1., args=args)(0)
+    params = dict(B=0.1)
+    down, up = kwant.wave_function(fsyst, energy=1., params=params)(0)
 
     x_hoppings = kwant.builder.HoppingKind((1,), lat)
     spin_current_z = ops.Current(fsyst, sigmaz, where=x_hoppings(syst))
-    _test(spin_current_z, up, args=args, per_el_val=1)
-    _test(spin_current_z, down, args=args, per_el_val=-1)
+    _test(spin_current_z, up, params=params, per_el_val=1)
+    _test(spin_current_z, down, params=params, per_el_val=-1)
 
     # calculate spin_x torque
     spin_torque_x = ops.Source(fsyst, sigmax, where=[lat(L//2)])
     i = fsyst.id_by_site[lat(L//2)]
     psi = up[2*i:2*(i+1)] + down[2*i:2*(i+1)]
-    H_ii = onsite(None, *args)
+    H_ii = onsite(None, **params)
     K = np.dot(H_ii, sigmax) - np.dot(sigmax, H_ii)
     expect = 1j * ft.reduce(np.dot, (psi.conj(), K, psi))
-    _test(spin_torque_x, up+down, args=args, reduced_val=expect)
+    _test(spin_torque_x, up+down, params=params, reduced_val=expect)
 
 
 def test_opservables_gauged():
@@ -397,17 +397,20 @@ def test_tocoo():
     # No accidental transpose.
     syst = kwant.Builder()
     lat2 = kwant.lattice.chain(norbs=2)
-    syst[lat2(0)] = lambda site, paramerer: np.eye(2)
+    syst[lat2(0)] = lambda site, p: np.eye(2)
     syst = syst.finalized()
     op = ops.Density(syst, [[1, 1], [0, 1]], check_hermiticity=False)
     assert np.all(op.tocoo().toarray() == [[1, 1], [0, 1]])
 
     op = ops.Density(syst, lambda site, p: [[1, 1], [0, 1]],
                      check_hermiticity=False)
-    op = op.bind(args=(1,))
+    op = op.bind(params=dict(p=1))
     raises(ValueError, op.tocoo, [1])
 
 
+# We need to keep testing 'args', but we don't want to see
+# all the deprecation warnings in the test logs
+@pytest.mark.filterwarnings("ignore:.*'args' parameter")
 @pytest.mark.parametrize("A", opservables)
 def test_arg_passing(A):
     lat1 = kwant.lattice.chain(norbs=1)
diff --git a/kwant/tests/test_wraparound.py b/kwant/tests/test_wraparound.py
index 9074e3e0f763d65f64b3b4a539a8ba373ff48761..f9b536ec3ac696e955c179d0a236f5d4d08c68c2 100644
--- a/kwant/tests/test_wraparound.py
+++ b/kwant/tests/test_wraparound.py
@@ -41,12 +41,13 @@ def test_consistence_with_bands(kx=1.9, nkys=31):
         wa_keep_1 = wraparound(syst, keep=1).finalized()
         wa_keep_none = wraparound(syst).finalized()
 
-        bands = kwant.physics.Bands(wa_keep_1, (kx,))
+        bands = kwant.physics.Bands(wa_keep_1, params=dict(k_x=kx))
         energies_a = [bands(ky) for ky in kys]
 
         energies_b = []
         for ky in kys:
-            H = wa_keep_none.hamiltonian_submatrix((kx, ky), sparse=False)
+            params = dict(k_x=kx, k_y=ky)
+            H = wa_keep_none.hamiltonian_submatrix(params=params, sparse=False)
             evs = np.sort(np.linalg.eigvalsh(H).real)
             energies_b.append(evs)
 
@@ -63,10 +64,14 @@ def test_opposite_hoppings():
         syst[lat(-1, 0), lat(-1, -1)] = val
 
         fsyst = wraparound(syst).finalized()
-        np.testing.assert_almost_equal(fsyst.hamiltonian_submatrix([0]), 0)
+        params = dict(k_x=0)
+        np.testing.assert_almost_equal(
+            fsyst.hamiltonian_submatrix(params=params),
+            0)
 
 
 def test_value_types(k=(-1.1, 0.5), E=2, t=1):
+    k = dict(zip(('k_x', 'k_y', 'k_z'), k))
     sym_extents = [1, 2, 3]
     lattices = [kwant.lattice.honeycomb(), kwant.lattice.square()]
     lat_syms = [
@@ -75,7 +80,7 @@ def test_value_types(k=(-1.1, 0.5), E=2, t=1):
     ]
     for lat, sym in lat_syms:
         syst = wraparound(_simple_syst(lat, E, t, sym)).finalized()
-        H = syst.hamiltonian_submatrix(k, sparse=False)
+        H = syst.hamiltonian_submatrix(params=k, sparse=False)
         for E1, t1 in [(float(E), float(t)),
                        (np.array([[E]], float), np.array([[t]], float)),
                        (ta.array([[E]], float), ta.array([[t]], float))]:
@@ -84,7 +89,7 @@ def test_value_types(k=(-1.1, 0.5), E=2, t=1):
             for E2 in [E1, lambda a: E1]:
                 for t2 in [t1, lambda a, b: t1]:
                     syst = wraparound(_simple_syst(lat, E2, t2, sym)).finalized()
-                    H_alt = syst.hamiltonian_submatrix(k, sparse=False)
+                    H_alt = syst.hamiltonian_submatrix(params=k, sparse=False)
                     np.testing.assert_equal(H_alt, H)
             # test when Hamiltonian value functions take extra parameters and
             # have incompatible signatures (must be passed with 'params')
@@ -98,7 +103,7 @@ def test_value_types(k=(-1.1, 0.5), E=2, t=1):
                 lambda a, b, t, E: t,
                 lambda a, b, E, t: t,
             ]
-            params = dict(E=E1, t=t1, k_x=k[0], k_y=k[1])
+            params = dict(E=E1, t=t1, **k)
             for E2, t2 in itertools.product(onsites, hoppings):
                 syst = wraparound(_simple_syst(lat, E2, t2, sym)).finalized()
                 H_alt = syst.hamiltonian_submatrix(params=params, sparse=False)
@@ -275,15 +280,17 @@ def test_fd_mismatch():
 
     def spectrum(syst, keep):
         syst = wraparound(syst, keep=keep).finalized()
+        ks = ('k_x', 'k_y', 'k_z')
         if keep is None:
             def _(*args):
-                return np.linalg.eigvalsh(syst.hamiltonian_submatrix(args=args))
+                params = dict(zip(ks, args))
+                return np.linalg.eigvalsh(
+                    syst.hamiltonian_submatrix(params=params))
         else:
             def _(*args):
-                args = list(args)
-                kext = args.pop(keep)
-                kint = args
-                B = kwant.physics.Bands(syst, args=kint)
+                params = dict(zip(ks, args))
+                kext = params.pop(ks[keep])
+                B = kwant.physics.Bands(syst, params=params)
                 return B(kext)
         return _
 
@@ -348,7 +355,8 @@ def test_fd_mismatch():
     finitewrapped.fill(wrapped, shape, start=np.zeros(3));
 
     sysf = finitewrapped.finalized()
-    spectrum1 = [np.linalg.eigvalsh(sysf.hamiltonian_submatrix(args=(k, 0)))
+    spectrum1 = [np.linalg.eigvalsh(
+                    sysf.hamiltonian_submatrix(params=dict(k_x=k, k_y=0)))
                 for k in np.linspace(-np.pi, np.pi, 5)]
 
     # Second choice: doubled UC with third translation purely in z direction
@@ -367,7 +375,8 @@ def test_fd_mismatch():
     finitewrapped.fill(wrapped, shape, start=np.zeros(3));
 
     sysf = finitewrapped.finalized()
-    spectrum2 = [np.linalg.eigvalsh(sysf.hamiltonian_submatrix(args=(k, 0)))
+    spectrum2 = [np.linalg.eigvalsh(
+                    sysf.hamiltonian_submatrix(params=dict(k_x=k, k_y=0)))
                 for k in np.linspace(-np.pi, np.pi, 5)]
 
     assert np.allclose(spectrum1, spectrum2)