From b891b56dc536f3737e6b7ad077eeb0bb8b07f5be Mon Sep 17 00:00:00 2001
From: Joseph Weston <joseph.weston08@gmail.com>
Date: Mon, 18 Feb 2013 18:31:52 +0100
Subject: [PATCH] add args and kwargs support for value functions
---
AUTHORS | 1 +
TODO | 2 -
doc/source/images/ab_ring.py.diff | 30 +++++--
doc/source/images/closed_system.py.diff | 17 ++--
doc/source/images/quantum_well.py.diff | 8 +-
doc/source/tutorial/ab_ring.py | 19 ++--
doc/source/tutorial/closed_system.py | 18 ++--
doc/source/tutorial/quantum_well.py | 23 +----
doc/source/tutorial/tutorial2.rst | 111 ++++--------------------
doc/source/whatsnew/0.3.rst | 52 +++++++++++
kwant/_system.pyx | 41 +++++----
kwant/builder.py | 16 ++--
kwant/solvers/common.py | 46 +++++++---
kwant/system.py | 20 +++--
kwant/tests/test_system.py | 21 +++++
15 files changed, 226 insertions(+), 199 deletions(-)
diff --git a/AUTHORS b/AUTHORS
index 3da5f159..901d5b3d 100644
--- a/AUTHORS
+++ b/AUTHORS
@@ -16,6 +16,7 @@ alternatives.
Other people that have contributed to kwant include
* Daniel Jaschke
+* Joseph Weston
-------------------------------------------------------------
diff --git a/TODO b/TODO
index 9568c819..ad2e1f6b 100644
--- a/TODO
+++ b/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.
diff --git a/doc/source/images/ab_ring.py.diff b/doc/source/images/ab_ring.py.diff
index 80d3ed99..062db080 100644
--- a/doc/source/images/ab_ring.py.diff
+++ b/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)])
diff --git a/doc/source/images/closed_system.py.diff b/doc/source/images/closed_system.py.diff
index 9d89782b..3b402123 100644
--- a/doc/source/images/closed_system.py.diff
+++ b/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)
diff --git a/doc/source/images/quantum_well.py.diff b/doc/source/images/quantum_well.py.diff
index c48d17e7..5efd325b 100644
--- a/doc/source/images/quantum_well.py.diff
+++ b/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()
diff --git a/doc/source/tutorial/ab_ring.py b/doc/source/tutorial/ab_ring.py
index a3706d4d..8f1d991b 100644
--- a/doc/source/tutorial/ab_ring.py
+++ b/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()
diff --git a/doc/source/tutorial/closed_system.py b/doc/source/tutorial/closed_system.py
index dc8263b5..2a3de700 100644
--- a/doc/source/tutorial/closed_system.py
+++ b/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.
diff --git a/doc/source/tutorial/quantum_well.py b/doc/source/tutorial/quantum_well.py
index e9b2465e..0e9eb572 100644
--- a/doc/source/tutorial/quantum_well.py
+++ b/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()
diff --git a/doc/source/tutorial/tutorial2.rst b/doc/source/tutorial/tutorial2.rst
index 1c45bd3a..9e5e9d6b 100644
--- a/doc/source/tutorial/tutorial2.rst
+++ b/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:
diff --git a/doc/source/whatsnew/0.3.rst b/doc/source/whatsnew/0.3.rst
index 411291db..94e31342 100644
--- a/doc/source/whatsnew/0.3.rst
+++ b/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`.
diff --git a/kwant/_system.pyx b/kwant/_system.pyx
index d4821e86..1023dbbe 100644
--- a/kwant/_system.pyx
+++ b/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
diff --git a/kwant/builder.py b/kwant/builder.py
index fac7e003..f31553a7 100644
--- a/kwant/builder.py
+++ b/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
diff --git a/kwant/solvers/common.py b/kwant/solvers/common.py
index a1798636..7452c979 100644
--- a/kwant/solvers/common.py
+++ b/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))
diff --git a/kwant/system.py b/kwant/system.py
index aa56515e..e0a6bf6e 100644
--- a/kwant/system.py
+++ b/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))
diff --git a/kwant/tests/test_system.py b/kwant/tests/test_system.py
index ce721a88..58664bbc 100644
--- a/kwant/tests/test_system.py
+++ b/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)
--
GitLab