add args and kwargs support for value functions

AUTHORS

@@ -16,6 +16,7 @@ alternatives.
 Other people that have contributed to kwant include
 
 * Daniel Jaschke
+* Joseph Weston
 
 -------------------------------------------------------------
 

TODO

@@ -28,8 +28,6 @@ Roughly in order of importance.                                     -*-org-*-
 
 * Allow plotting of infinite systems
 
-* Consider additional optional parameters to value functions
-
 * Use sparse linear algebra to calculate bands
   However, SciPy's sparse eigenvalues don't seem to work well.
 

doc/source/images/ab_ring.py.diff

@@ -8,7 +8,27 @@
  
  
  def make_system(a=1, t=1.0, W=10, r1=10, r2=20):
-@@ -82,6 +83,54 @@
+@@ -38,12 +39,13 @@
+     for hopping in lat.nearest:
+         sys[sys.possible_hoppings(*hopping)] = -t
+ 
+-    # In order to introduce a flux through the ring, we introduce a phase on
+-    # the hoppings on the line cut through one of the arms.  Since we want to
+-    # change the flux without modifying the Builder instance repeatedly, we
+-    # define the modified hoppings as a function that takes the flux as its
+-    # parameter phi.
+-    def fluxphase(site1, site2, phi):
++    # In order to introduce a flux through the ring, we introduce a phase
++    # on the hoppings on the line cut through one of the arms
++
++    # since we want to change the flux without modifying Builder repeatedly,
++    # we define the modified hoppings as a function that takes the flux
++    # through the argument phi.
++    def fluxphase(site1, site2, phi=0):
+         return exp(1j * phi)
+ 
+     def crosses_branchcut(hop):
+@@ -81,6 +83,54 @@
      return sys
  
  
@@ -62,9 +82,9 @@
 +
  def plot_conductance(sys, energy, fluxes):
      # compute conductance
-     # global variable phi controls the flux
-@@ -95,18 +144,29 @@
-         smatrix = kwant.solve(sys, energy)
+ 
+@@ -90,18 +140,29 @@
+         smatrix = kwant.solve(sys, energy, kwargs={'phi': flux})
          data.append(smatrix.transmission(1, 0))
  
 -    pyplot.figure()
@@ -98,7 +118,7 @@
  
      # Finalize the system.
      sys = sys.finalized()
-@@ -115,6 +175,15 @@
+@@ -110,6 +171,15 @@
      plot_conductance(sys, energy=0.15, fluxes=[0.01 * i * 3 * 2 * pi
                                                  for i in xrange(100)])
  

doc/source/images/closed_system.py.diff

@@ -8,7 +8,7 @@
  
  
  def make_system(a=1, t=1.0, r=10):
-@@ -70,32 +71,43 @@
+@@ -66,28 +67,40 @@
  
          energies.append(ev)
  
@@ -19,12 +19,12 @@
 -    pyplot.ylabel("energy [in units of t]")
 -    pyplot.show()
 +    pyplot.xlabel("magnetic field [some arbitrary units]",
-+                 fontsize=_defs.mpl_label_size)
++                  fontsize=_defs.mpl_label_size)
 +    pyplot.ylabel("energy [in units of t]", fontsize=_defs.mpl_label_size)
 +    pyplot.setp(fig.get_axes()[0].get_xticklabels(),
-+               fontsize=_defs.mpl_tick_size)
++                fontsize=_defs.mpl_tick_size)
 +    pyplot.setp(fig.get_axes()[0].get_yticklabels(),
-+               fontsize=_defs.mpl_tick_size)
++                fontsize=_defs.mpl_tick_size)
 +    fig.set_size_inches(_defs.mpl_width_in, _defs.mpl_width_in * 3. / 4.)
 +    fig.subplots_adjust(left=0.15, right=0.95, top=0.95, bottom=0.15)
 +    fig.savefig("closed_system_result.pdf")
@@ -32,13 +32,10 @@
  
  
  def plot_wave_function(sys):
-     # We reset the magnetic field to equal to 0.
-     global B
-     B = 0.
 +    size = (_defs.figwidth_in, _defs.figwidth_in)
- 
++
      # Calculate the wave functions in the system.
-     ham_mat = sys.hamiltonian_submatrix(sparse=True)
+     ham_mat = sys.hamiltonian_submatrix(sparse=True, kwargs={'B': 0})
      evecs = sla.eigsh(ham_mat, k=20, which='SM')[1]
  
      # Plot the probability density of the 10th eigenmode.
@@ -60,7 +57,7 @@
      # Finalize the system.
      sys = sys.finalized()
  
-@@ -108,6 +120,7 @@
+@@ -100,6 +113,7 @@
      sys = make_system(r=30).finalized()
      plot_wave_function(sys)
  

doc/source/images/quantum_well.py.diff

@@ -6,10 +6,10 @@
  from matplotlib import pyplot
 +import _defs
  
- # global variable governing the behavior of potential() in
- # make_system()
-@@ -74,19 +75,25 @@
-         smatrix = kwant.solve(sys, energy)
+ def make_system(a=1, t=1.0, W=10, L=30, L_well=10):
+     # Start with an empty tight-binding system and a single square lattice.
+@@ -62,19 +63,25 @@
+         smatrix = kwant.solve(sys, energy, kwargs={'pot': -welldepth})
          data.append(smatrix.transmission(1, 0))
  
 -    pyplot.figure()

doc/source/tutorial/ab_ring.py

@@ -42,14 +42,13 @@ def make_system(a=1, t=1.0, W=10, r1=10, r2=20):
         sys[sys.possible_hoppings(*hopping)] = -t
 #HIDDEN_END_lcak
 
-    # In order to introduce a flux through the ring, we introduce a phase
-    # on the hoppings on the line cut through one of the arms
-
-    # since we want to change the flux without modifying Builder repeatedly,
-    # we define the modified hoppings as a function that takes the flux
-    # through the global variable phi.
+    # In order to introduce a flux through the ring, we introduce a phase on
+    # the hoppings on the line cut through one of the arms.  Since we want to
+    # change the flux without modifying the Builder instance repeatedly, we
+    # define the modified hoppings as a function that takes the flux as its
+    # parameter phi.
 #HIDDEN_BEGIN_lvkt
-    def fluxphase(site1, site2):
+    def fluxphase(site1, site2, phi):
         return exp(1j * phi)
 
     def crosses_branchcut(hop):
@@ -94,15 +93,11 @@ def make_system(a=1, t=1.0, W=10, r1=10, r2=20):
 
 def plot_conductance(sys, energy, fluxes):
     # compute conductance
-    # global variable phi controls the flux
-    global phi
 
     normalized_fluxes = [flux / (2 * pi) for flux in fluxes]
     data = []
     for flux in fluxes:
-        phi = flux
-
-        smatrix = kwant.solve(sys, energy)
+        smatrix = kwant.solve(sys, energy, kwargs={'phi': flux})
         data.append(smatrix.transmission(1, 0))
 
     pyplot.figure()

doc/source/tutorial/closed_system.py

@@ -37,8 +37,8 @@ def make_system(a=1, t=1.0, r=10):
         rsq = x ** 2 + y ** 2
         return rsq < r ** 2
 
-    def hopx(site1, site2):
-        # The magnetic field is controlled by the global variable B
+    def hopx(site1, site2, B):
+        # The magnetic field is controlled by the parameter B
         y = site1.pos[1]
         return -t * exp(-1j * B * y)
 
@@ -55,8 +55,6 @@ def make_system(a=1, t=1.0, r=10):
 
 #HIDDEN_BEGIN_yvri
 def plot_spectrum(sys, Bfields):
-    # global variable B controls the magnetic field
-    global B
 
     # In the following, we compute the spectrum of the quantum dot
     # using dense matrix methods. This works in this toy example, as
@@ -64,11 +62,9 @@ def plot_spectrum(sys, Bfields):
     # sparse matrix methods
 
     energies = []
-    for Bfield in Bfields:
-        B = Bfield
-
+    for B in Bfields:
         # Obtain the Hamiltonian as a dense matrix
-        ham_mat = sys.hamiltonian_submatrix(sparse=True)
+        ham_mat = sys.hamiltonian_submatrix(sparse=True, kwargs={'B': B})
 
         # we only calculate the 15 lowest eigenvalues
         ev = sla.eigsh(ham_mat, k=15, which='SM', return_eigenvectors=False)
@@ -85,12 +81,8 @@ def plot_spectrum(sys, Bfields):
 
 #HIDDEN_BEGIN_wave
 def plot_wave_function(sys):
-    # We reset the magnetic field to equal to 0.
-    global B
-    B = 0.
-
     # Calculate the wave functions in the system.
-    ham_mat = sys.hamiltonian_submatrix(sparse=True)
+    ham_mat = sys.hamiltonian_submatrix(sparse=True, kwargs={'B': 0})
     evecs = sla.eigsh(ham_mat, k=20, which='SM')[1]
 
     # Plot the probability density of the 10th eigenmode.

doc/source/tutorial/quantum_well.py

@@ -11,15 +11,7 @@ import kwant
 # For plotting
 from matplotlib import pyplot
 
-# global variable governing the behavior of potential() in
-# make_system()
-#HIDDEN The following code line is included verbatim in the tutorial text
-#HIDDEN because nested code examples are not supported.  Remember to update
-#HIDDEN the tutorial text when you modify this line.
 #HIDDEN_BEGIN_ehso
-pot = 0
-
-
 def make_system(a=1, t=1.0, W=10, L=30, L_well=10):
     # Start with an empty tight-binding system and a single square lattice.
     # `a` is the lattice constant (by default set to 1 for simplicity.
@@ -29,18 +21,17 @@ def make_system(a=1, t=1.0, W=10, L=30, L_well=10):
 
     #### Define the scattering region. ####
     # Potential profile
-    def potential(site):
+    def potential(site, pot):
         (x, y) = site.pos
         if (L - L_well) / 2 < x < (L + L_well) / 2:
-            # The potential value is provided using a global variable
             return pot
         else:
             return 0
 #HIDDEN_END_ehso
 
 #HIDDEN_BEGIN_coid
-    def onsite(site):
-        return 4 * t + potential(site)
+    def onsite(site, pot=0):
+        return 4 * t + potential(site, pot)
 
     sys[(lat(x, y) for x in range(L) for y in range(W))] = onsite
     for hopping in lat.nearest:
@@ -68,18 +59,12 @@ def make_system(a=1, t=1.0, W=10, L=30, L_well=10):
 
 
 def plot_conductance(sys, energy, welldepths):
-    # We specify that we want to not only read, but also write to a
-    # global variable.
 #HIDDEN_BEGIN_sqvr
-    global pot
 
     # Compute conductance
     data = []
     for welldepth in welldepths:
-        # Set the global variable that defines the potential well depth
-        pot = -welldepth
-
-        smatrix = kwant.solve(sys, energy)
+        smatrix = kwant.solve(sys, energy, kwargs={'pot': -welldepth})
         data.append(smatrix.transmission(1, 0))
 
     pyplot.figure()

doc/source/tutorial/tutorial2.rst

@@ -148,15 +148,12 @@ define the potential profile of a quantum well as:
     :start-after: #HIDDEN_BEGIN_ehso
     :end-before: #HIDDEN_END_ehso
 
-This function takes one argument which is of type
+This function takes two arguments: the first which is of type
 `~kwant.builder.Site`, from which you can get the real-space
-coordinates using ``site.pos``. Note that we use several global
+coordinates using ``site.pos``, and the value of the potential as a
+second argument. Note that we can use global
 variables to define the behavior of `potential()`: `L` and `L_well`
-are variables taken from the namespace of `make_system`, the variable `pot`
-is taken from the global module namespace. By this one can change the
-behavior of `potential()` at another place, for example by setting
-`pot` to a different value. We will use this later to compute
-the transmission as a function of well depth.
+are variables taken from the namespace of `make_system`.
 
 kwant now allows us to pass a function as a value to
 `~kwant.builder.Builder`:
@@ -165,8 +162,11 @@ kwant now allows us to pass a function as a value to
     :start-after: #HIDDEN_BEGIN_coid
     :end-before: #HIDDEN_END_coid
 
-For each lattice point, the corresponding site is then passed to the
-function `onsite()`. Note that we had to define `onsite()`, as it is
+For each lattice point, the corresponding site is then passed as the
+first argument to the function `onsite()`. The values of any additional
+parameters, which can be used to alter the Hamiltonian matrix elements
+at a later stage, are specified later during the call to `solve()`.
+Note that we had to define `onsite()`, as it is
 not possible to mix values and functions as in ``sys[...] = 4 * t +
 potential``.
 
@@ -181,12 +181,11 @@ Finally, we compute the transmission probability:
     :start-after: #HIDDEN_BEGIN_sqvr
     :end-before: #HIDDEN_END_sqvr
 
-Since we change the value of the global variable `pot` to vary the
-well depth, python requires us to write ``global pot`` to `enable
-access to it
-<http://docs.python.org/faq/programming.html#what-are-the-rules-for-local-and-global-variables-in-python>`_.
-Subsequent calls to :func:`kwant.solve <kwant.solvers.sparse.solve>`
-then will use the updated value of pot, and we get the result:
+``kwant.solve`` allows us to specify a dictionary, `kwargs`, that will
+be passed as additional keyword arguments to the function that evaluates the
+Hamiltonian matrix elements.  In this example we are able to solve the system
+for different depths of the potential well by passing the keyword argument
+`pot`. We obtain the result:
 
 .. image:: ../images/quantum_well_result.*
 
@@ -206,88 +205,12 @@ oscillatory transmission behavior through resonances in the quantum well.
 .. specialnote:: Technical details
 
   - Functions can also be used for hoppings. In this case, they take
-    two `~kwant.builder.Site`'s as arguments.
-
-  - In tutorial/quantum_well.py, the line ::
-
-        pot = 0
-
-    is not really necessary. If this line was left out, the
-    global variable `pot` would in fact be created by the
-    first assignment in `plot_conductance()`.
+    two `~kwant.builder.Site`'s as arguments and then an arbitrary number
+    of additional arguments.
 
   - Apart from the real-space position `pos`, `~kwant.builder.Site` has also an
     attribute `tag` containing the lattice indices of the site.
 
-  - Since we use a global variable to change the value of the
-    potential, let us briefly reflect on how python determines
-    which variable to use.
-
-    In our example, the function `potential()` uses the variable
-    `pot` which is not defined in the function itself. In this case,
-    python looks for the variable in the enclosing scopes, i.e.
-    inside the functions/modules/scripts that enclose the
-    corresponding piece of code. For example, in
-
-    >>> def f():
-    ...     def g():
-    ...         print string
-    ...     return g
-    ...
-    >>> g = f()
-    >>> string = "global"
-    >>> g()
-    global
-
-    function `g()` defined inside `f()` uses the global variable
-    `string` (which was actually created only *after* the
-    definition of `g()`!). Note that this only works as long as
-    one only reads `string`; if `g()` was to write to string,
-    it would need to add ``global string`` to `g()`, as we
-    did in `plot_conductance()`.
-
-    Things change if the function `f()` also contains a variable
-    of the same name:
-
-    >>> def f():
-    ...     def g():
-    ...         print string
-    ...     string = "local"
-    ...     return g
-    ...
-    >>> g = f()
-    >>> g()
-    local
-    >>> string = "global"
-    >>> g()
-    local
-
-    In this case, `g()` always uses the local variable inside `f()`
-    (unless we would add ``global string`` in `g()`).
-
-  - `~kwant.builder.Builder` in fact accepts not only functions but any python
-    object which is callable.  We can take advantage of the fact that instances
-    of python classes with a `__call__` method can be called just as if they
-    were functions::
-
-        class Well:
-            def __init__(self, a, b=0):
-                self.a = a
-                self.b = b
-
-            def __call__(self, site):
-                x, y = site.pos
-                return self.a * (x**2 + y**2) + b
-
-        well = Well(3, 4)
-
-        sys[...] = well
-
-        well.a = ...
-
-    This approach allows to avoid the use of global variables.  Parameters can
-    be changed inside the object.
-
 .. _tutorial-abring:
 
 Nontrivial shapes
@@ -348,7 +271,7 @@ the branch cut to `fluxphase`. The rationale
 behind using a function instead of a constant value for the hopping
 is again that we want to vary the flux through the ring, without
 constantly rebuilding the system -- instead the flux is governed
-by the global variable `phi`.
+by the parameter `phi`.
 
 For the leads, we can also use the ``lat.shape()``-functionality:
 

doc/source/whatsnew/0.3.rst

@@ -19,3 +19,55 @@ make two site groups which have the same geometry, but mean different things,
 as for instance in :doc:`../tutorial/tutorial5`, then the `name` argument has
 to be used when creating a lattice, e.g. `a = kwant.lattice.square(name='a');
 b = kwant.lattice.square(name='b')`.
+
+Parameters to Hamiltonian
+-------------------------
+kwant now allows the Hamiltonian matrix elements to be described with functions
+that depend on an arbitrary number of parameters in addition to the sites on
+which they are defined.
+
+Previously, functions defining the Hamiltonian matrix elements had to have the
+following prototypes::
+
+    def onsite(site):
+        ...
+
+    def hopping(site1, site2):
+        ...
+
+If the Hamiltonian elements need to depend on some other external parameters
+(e.g. magnetic field) then those had to be provided by some other means than
+regular function parameters (e.g. global variables).
+
+Now the value functions may accept arbitrary arguments after the `Site`
+arguments.  These extra arguments can be specified when
+`~kwant.solvers.default.solve` is called by setting the keyword arguments:
+
+args
+    A tuple of values to be passed as the positional arguments to the
+    Hamiltonian value functions (not including the `Site` arguments).
+
+kwargs
+    A dictionary of keyword-value pairs to be passed as the keyword
+    arguments to the Hamiltonian value functions.
+
+For example, if the hopping and onsite Hamiltonian value functions have
+the following prototype::
+
+    def onsite(site, t, B, pot):
+        ...
+
+    def hopping(site1, site2, t, B, pot):
+        ...
+
+then the values of `t`, `B` and `pot` for which to solve the system can be
+passed to `solve` like this::
+
+    time = 2
+    magnetic_field = 3
+    potential = 4
+    kwant.solve(sys, energy,
+                args=(time,), kwargs={'B': magnetic_field, 'pot': potential})
+
+Arguments can be passed in an equivalent way in calls to
+`~kwant.solvers.default.wave_func` and `~kwant.solvers.default.ldos`.

kwant/_system.pyx

@@ -20,7 +20,8 @@ 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, CGraph gr, diag, gint [:] from_sites, n_by_to_site,
+def make_sparse(ham, args, kwargs, CGraph gr, diag,
+                gint [:] from_sites, n_by_to_site,
                 gint [:] to_norb, gint [:] to_off,
                 gint [:] from_norb, gint [:] from_off):
     """For internal use by hamiltonian_submatrix."""
@@ -68,7 +69,7 @@ def make_sparse(ham, CGraph gr, diag, gint [:] from_sites, n_by_to_site,
             if ts not in n_by_to_site:
                 continue
             n_ts = n_by_to_site[ts]
-            h = matrix(ham(ts, fs), complex)
+            h = matrix(ham(ts, fs, *args, **kwargs), 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 xrange(h.shape[0]):
@@ -82,7 +83,7 @@ def make_sparse(ham, CGraph gr, diag, gint [:] from_sites, n_by_to_site,
 
 
 @cython.boundscheck(False)
-def make_sparse_full(ham, CGraph gr, diag,
+def make_sparse_full(ham, args, kwargs, CGraph gr, diag,
                      gint [:] to_norb, gint [:] to_off,
                      gint [:] from_norb, gint [:] from_off):
     """For internal use by hamiltonian_submatrix."""
@@ -124,7 +125,7 @@ def make_sparse_full(ham, CGraph gr, diag,
         for ts in nbors.data[:nbors.size]:
             if ts < fs:
                 continue
-            h = matrix(ham(ts, fs), complex)
+            h = matrix(ham(ts, fs, *args, **kwargs), complex)
             if h.shape[0] != to_norb[ts] or h.shape[1] != from_norb[fs]:
                 raise ValueError(msg.format(fs, ts))
             for i in xrange(h.shape[0]):
@@ -139,7 +140,8 @@ def make_sparse_full(ham, CGraph gr, diag,
 
 
 @cython.boundscheck(False)
-def make_dense(ham, CGraph gr, diag, gint [:] from_sites, n_by_to_site,
+def make_dense(ham, args, kwargs, CGraph gr, diag,
+               gint [:] from_sites, n_by_to_site,
                gint [:] to_norb, gint [:] to_off,
                gint [:] from_norb, gint [:] from_off):
     """For internal use by hamiltonian_submatrix."""
@@ -167,7 +169,7 @@ def make_dense(ham, CGraph gr, diag, gint [:] from_sites, n_by_to_site,
             if ts not in n_by_to_site:
                 continue
             n_ts = n_by_to_site[ts]
-            h = matrix(ham(ts, fs), complex)
+            h = matrix(ham(ts, fs, *args, **kwargs), 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],
@@ -176,7 +178,7 @@ def make_dense(ham, CGraph gr, diag, gint [:] from_sites, n_by_to_site,
 
 
 @cython.boundscheck(False)
-def make_dense_full(ham, CGraph gr, diag,
+def make_dense_full(ham, args, kwargs, CGraph gr, diag,
                     gint [:] to_norb, gint [:] to_off,
                     gint [:] from_norb, gint [:] from_off):
     """For internal use by hamiltonian_submatrix."""
@@ -200,7 +202,7 @@ def make_dense_full(ham, CGraph gr, diag,
         for ts in nbors.data[:nbors.size]:
             if ts < fs:
                 continue
-            h = mat = matrix(ham(ts, fs), complex)
+            h = mat = matrix(ham(ts, fs, *args, **kwargs), 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))
@@ -212,7 +214,8 @@ def make_dense_full(ham, CGraph gr, diag,
 
 
 def hamiltonian_submatrix(self, to_sites=None, from_sites=None,
-                          sparse=False, return_norb=False):
+                          sparse=False, return_norb=False,
+                          args=(), kwargs={}):
     """Return a submatrix of the system Hamiltonian.
 
     Parameters
@@ -223,6 +226,10 @@ def hamiltonian_submatrix(self, to_sites=None, from_sites=None,
         Whether to return a sparse or a dense matrix. Defaults to `False`.
     return_norb : bool
         Whether to return arrays of numbers of orbitals.  Defaults to `False`.
+    args : tuple, defaults to empty
+        Positional arguments to pass to the ``hamiltonian`` method.
+    kwargs : dictionary, defaults to empty
+        Keyword arguments to pass to the ``hamiltonian`` method.
 
     Returns
     -------
@@ -253,7 +260,7 @@ def hamiltonian_submatrix(self, to_sites=None, from_sites=None,
         diag = n * [None]
         from_norb = np.empty(n, gint_dtype)
         for site in xrange(n):
-            diag[site] = h = matrix(ham(site, site), complex)
+            diag[site] = h = matrix(ham(site, site, *args, **kwargs), complex)
             from_norb[site] = h.shape[0]
     else:
         diag = len(from_sites) * [None]
@@ -261,7 +268,8 @@ def hamiltonian_submatrix(self, 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), complex)
+            diag[n_site] = h = matrix(ham(site, site, *args, **kwargs),
+                                      complex)
             from_norb[n_site] = h.shape[0]
     from_off = np.empty(from_norb.shape[0] + 1, gint_dtype)
     from_off[0] = 0
@@ -274,14 +282,14 @@ def hamiltonian_submatrix(self, to_sites=None, from_sites=None,
         if to_sites is None:
             to_norb = np.empty(n, gint_dtype)
             for site in xrange(n):
-                h = matrix(ham(site, site), complex)
+                h = matrix(ham(site, site, *args, **kwargs), 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), complex)
+                h = matrix(ham(site, site, *args, **kwargs), complex)
                 to_norb[n_site] = h.shape[0]
         to_off = np.empty(to_norb.shape[0] + 1, gint_dtype)
         to_off[0] = 0
@@ -290,7 +298,8 @@ def hamiltonian_submatrix(self, 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, self.graph, diag, to_norb, to_off, from_norb, from_off)
+        mat = func(ham, args, kwargs, self.graph, diag, to_norb, to_off,
+                   from_norb, from_off)
     else:
         if to_sites is None:
             to_sites = np.arange(n)
@@ -305,6 +314,6 @@ def hamiltonian_submatrix(self, 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, self.graph, diag, from_sites, n_by_to_site,
-                   to_norb, to_off, from_norb, from_off)
+        mat = func(ham, args, kwargs, 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

@@ -1186,11 +1186,12 @@ class FiniteSystem(system.FiniteSystem):
 
     Usable as input for the solvers in `kwant.solvers`.
     """
-    def hamiltonian(self, i, j):
+    def hamiltonian(self, i, j, *args, **kwargs):
         if i == j:
             value = self.onsite_hamiltonians[i]
             if hasattr(value, '__call__'):
-                value = value(self.symmetry.to_fd(self.sites[i]))
+                value = value(self.symmetry.to_fd(self.sites[i]),
+                                                  *args, **kwargs)
             return value
         else:
             edge_id = self.graph.first_edge_id(i, j)
@@ -1203,7 +1204,8 @@ class FiniteSystem(system.FiniteSystem):
             if hasattr(value, '__call__'):
                 site_i = self.sites[i]
                 site_j = self.sites[j]
-                value = value(*self.symmetry.to_fd(site_i, site_j))
+                site_i, site_j = self.symmetry.to_fd(site_i,site_j)
+                value = value(site_i, site_j, *args, **kwargs)
             if conj:
                 value = herm_conj(value)
             return value
@@ -1217,13 +1219,14 @@ class FiniteSystem(system.FiniteSystem):
 
 class InfiniteSystem(system.InfiniteSystem):
     """Finalized infinite system, extracted from a `Builder`."""
-    def hamiltonian(self, i, j):
+    def hamiltonian(self, i, j, *args, **kwargs):
         if i == j:
             if i >= self.slice_size:
                 i -= self.slice_size
             value = self.onsite_hamiltonians[i]
             if hasattr(value, '__call__'):
-                value = value(self.symmetry.to_fd(self.sites[i]))
+                value = value(self.symmetry.to_fd(self.sites[i]),
+                                                  *args, **kwargs)
             return value
         else:
             edge_id = self.graph.first_edge_id(i, j)
@@ -1236,7 +1239,8 @@ class InfiniteSystem(system.InfiniteSystem):
             if hasattr(value, '__call__'):
                 site_i = self.sites[i]
                 site_j = self.sites[j]
-                value = value(*self.symmetry.to_fd(site_i, site_j))
+                site_i, site_j = self.symmetry.to_fd(site_i,site_j)
+                value = value(site_i, site_j, *args, **kwargs)
             if conj:
                 value = herm_conj(value)
             return value

kwant/solvers/common.py

@@ -92,7 +92,8 @@ class SparseSolver(object):
         pass
 
     def _make_linear_sys(self, sys, out_leads, in_leads, energy=0,
-                        force_realspace=False, check_hermiticity=True):
+                        force_realspace=False, check_hermiticity=True,
+                        args=(), kwargs={}):
         """
         Make a sparse linear system of equations defining a scattering
         problem.
@@ -114,6 +115,10 @@ class SparseSolver(object):
             more computationally expensive and less stable.
         check_hermiticity : bool
             Check if Hamiltonian matrices are in fact Hermitian.
+        args : tuple, defaults to empty
+            Positional arguments to pass to the ``hamiltonian`` method.
+        kwargs : dictionary, defaults to empty
+            Keyword arguments to pass to the ``hamiltonian`` method.
 
         Returns
         -------
@@ -144,7 +149,8 @@ class SparseSolver(object):
         if not sys.lead_interfaces:
             raise ValueError('System contains no leads.')
         lhs, norb = sys.hamiltonian_submatrix(sparse=True,
-                                              return_norb=True)[:2]
+                                              return_norb=True,
+                                              args=args, kwargs=kwargs)[:2]
         lhs = getattr(lhs, 'to' + self.lhsformat)()
         lhs = lhs - energy * sp.identity(lhs.shape[0], format=self.lhsformat)
 
@@ -273,7 +279,8 @@ class SparseSolver(object):
         return LinearSys(lhs, rhs, kept_vars), lead_info
 
     def solve(self, sys, energy=0, out_leads=None, in_leads=None,
-              force_realspace=False, check_hermiticity=True):
+              force_realspace=False, check_hermiticity=True,
+              args=(), kwargs={}):
         """
         Compute the scattering matrix or Green's function between leads.
 
@@ -298,6 +305,10 @@ class SparseSolver(object):
             more computationally expensive and less stable.
         check_hermiticity : ``bool``
             Check if the Hamiltonian matrices are Hermitian.
+        args : tuple, defaults to empty
+            Positional arguments to pass to the ``hamiltonian`` method.
+        kwargs : dictionary, defaults to empty
+            Keyword arguments to pass to the ``hamiltonian`` method.
 
         Returns
         -------
@@ -349,7 +360,8 @@ class SparseSolver(object):
 
         linsys, lead_info = self._make_linear_sys(sys, out_leads, in_leads,
                                                   energy, force_realspace,
-                                                  check_hermiticity)
+                                                  check_hermiticity,
+                                                  args, kwargs)
 
         flhs = self._factorized(linsys.lhs)
         data = self._solve_linear_sys(flhs, linsys.rhs, linsys.kept_vars)
@@ -357,7 +369,7 @@ class SparseSolver(object):
         return BlockResult(data, lead_info, out_leads, in_leads)
 
 
-    def ldos(self, fsys, energy=0):
+    def ldos(self, fsys, energy=0, args=(), kwargs={}):
         """
         Calculate the local density of states of a system at a given energy.
 
@@ -368,6 +380,12 @@ class SparseSolver(object):
             scattering region.
         energy : number
             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
+        kwargs : dictionary of keyword, argument pairs or empty dictionary
+            Optional keyword arguments to pass to the function(s) which
+            evaluate the hamiltonian matrix elements
 
         Returns
         -------
@@ -381,7 +399,8 @@ class SparseSolver(object):
                                  "tight binding systems.")
 
         (h, rhs, kept_vars), lead_info = \
-            self._make_linear_sys(fsys, [], xrange(len(fsys.leads)), energy)
+            self._make_linear_sys(fsys, [], xrange(len(fsys.leads)), energy,
+                                  args=args, kwargs=kwargs)
 
         Modes = physics.Modes
         num_extra_vars = sum(li.vecs.shape[1] - li.nmodes
@@ -404,7 +423,7 @@ class SparseSolver(object):
 
         return ldos * (0.5 / np.pi)
 
-    def wave_func(self, sys, energy=0):
+    def wave_func(self, sys, energy=0, args=(), kwargs={}):
         """
         Return a callable object for the computation of the wave function
         inside the scattering region.
@@ -414,6 +433,12 @@ class SparseSolver(object):
         sys : `kwant.system.FiniteSystem`
             The low level system for which the wave functions are to be
             calculated.
+        args : tuple of arguments, or empty tuple
+            Positional arguments to pass to the function(s) which
+            evaluate the hamiltonian matrix elements
+        kwargs : dictionary of keyword, argument pairs or empty dictionary
+            Optional keyword arguments to pass to the function(s) which
+            evaluate the hamiltonian matrix elements
 
         Notes
         -----
@@ -427,18 +452,19 @@ class SparseSolver(object):
         >>> wf = kwant.solvers.default.wave_func(some_sys, some_energy)
         >>> wfs_of_lead_2 = wf(2)
         """
-        return WaveFunc(self, sys, energy)
+        return WaveFunc(self, sys, energy, args, kwargs)
 
 
 class WaveFunc(object):
-    def __init__(self, solver, sys, energy=0):
+    def __init__(self, solver, sys, energy=0, args=(), kwargs={}):
         for lead in sys.leads:
             if not isinstance(lead, system.InfiniteSystem):
                 # TODO: fix this
                 msg = 'All leads must be tight binding systems.'
                 raise ValueError(msg)
         (h, self.rhs, kept_vars), lead_info = \
-            solver._make_linear_sys(sys, [], xrange(len(sys.leads)), energy)
+            solver._make_linear_sys(sys, [], xrange(len(sys.leads)),
+                                    energy, args=args, kwargs=kwargs)
         Modes = physics.Modes
         num_extra_vars = sum(li.vecs.shape[1] - li.nmodes
                              for li in lead_info if isinstance(li, Modes))

kwant/system.py

@@ -45,12 +45,15 @@ class System(object):
         return 1 if np.isscalar(ham) else ham.shape[0]
 
     @abc.abstractmethod
-    def hamiltonian(self, i, j):
+    def hamiltonian(self, i, j, *args, **kwargs):
         """Return the hamiltonian matrix element for sites `i` and `j`.
 
         If ``i == j``, return the on-site Hamiltonian of site `i`.
 
         if ``i != j``, return the hopping between site `i` and `j`.
+
+        Hamiltonians may depend (optionally) on positional and
+        keyword arguments
         """
         pass
 
@@ -129,29 +132,30 @@ class InfiniteSystem(System):
     """
     __metaclass__ = abc.ABCMeta
 
-    def slice_hamiltonian(self, sparse=False):
+    def slice_hamiltonian(self, sparse=False, args=(), kwargs={}):
         """Hamiltonian of a single slice of the infinite system."""
         slice_sites = xrange(self.slice_size)
         return self.hamiltonian_submatrix(slice_sites, slice_sites,
-                                          sparse=sparse)
+                                          sparse=sparse, kwargs=kwargs)
 
-    def inter_slice_hopping(self, sparse=False):
+    def inter_slice_hopping(self, sparse=False, args=(), kwargs={}):
         """Hopping Hamiltonian between two slices of the infinite system."""
         slice_sites = xrange(self.slice_size)
         neighbor_sites = xrange(self.slice_size, self.graph.num_nodes)
         return self.hamiltonian_submatrix(slice_sites, neighbor_sites,
-                                          sparse=sparse)
+                                          sparse=sparse, kwargs=kwargs)
 
-    def self_energy(self, energy):
+    def self_energy(self, energy, args=(), kwargs={}):
         """Return self-energy of a lead.
 
         The returned matrix has the shape (n, n), where n is
         ``sum(self.num_orbitals(i) for i in range(self.slice_size))``.
         """
-        ham = self.slice_hamiltonian()
+        ham = self.slice_hamiltonian(args=args, kwargs=kwargs)
         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.self_energy(ham, self.inter_slice_hopping())
+        return physics.self_energy(ham, self.inter_slice_hopping(args=args,
+                                                                 kwargs=kwargs))

kwant/tests/test_system.py

@@ -65,3 +65,24 @@ def test_hamiltonian_submatrix():
     sys2 = sys.finalized()
     assert_raises(ValueError, sys2.hamiltonian_submatrix)
     assert_raises(ValueError, sys2.hamiltonian_submatrix, None, None, 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
+
+    sys = kwant.Builder()
+    sys[(gr(i) for i in xrange(3))] = onsite
+    sys[((gr(i), gr(i + 1)) for i in xrange(2))] = hopping
+    sys2 = sys.finalized()
+    mat = sys2.hamiltonian_submatrix(args=(2,), kwargs={'p2': 1})
+    mat_should_be = [[5, 1, 0], [1, 4, 1.], [0, 1, 3]]
+
+    # Sorting is required due to unknown compression order of builder.
+    onsite_hamiltonians = mat.flat[::3]
+    perm = np.argsort(onsite_hamiltonians)
+    mat = mat[perm, :]
+    mat = mat[:, perm]
+    np.testing.assert_array_equal(mat, mat_should_be)