new tutorial dealing with superconductivity

doc/Makefile

@@ -26,7 +26,9 @@ TUTORIAL2C_IMAGES = source/images/tutorial2c_result.png source/images/tutorial2c
 TUTORIAL3A_IMAGES = source/images/tutorial3a_result.png source/images/tutorial3a_result.pdf
 TUTORIAL3B_IMAGES = source/images/tutorial3b_result.png source/images/tutorial3b_result.pdf source/images/tutorial3b_sys.png source/images/tutorial3b_sys.pdf
 TUTORIAL4_IMAGES = source/images/tutorial4_result.png source/images/tutorial4_result.pdf source/images/tutorial4_sys1.png source/images/tutorial4_sys1.pdf source/images/tutorial4_sys2.png source/images/tutorial4_sys2.pdf source/images/tutorial4_bs.png source/images/tutorial4_bs.pdf
-ALL_TUTORIAL_IMAGES = $(TUTORIAL1A_IMAGES)  $(TUTORIAL2A_IMAGES) $(TUTORIAL2B_IMAGES) $(TUTORIAL2C_IMAGES) $(TUTORIAL3A_IMAGES) $(TUTORIAL3B_IMAGES) $(TUTORIAL4_IMAGES)
+TUTORIAL5A_IMAGES = source/images/tutorial5a_result.png source/images/tutorial5a_result.pdf
+TUTORIAL5B_IMAGES = source/images/tutorial5b_result.png source/images/tutorial5b_result.pdf
+ALL_TUTORIAL_IMAGES = $(TUTORIAL1A_IMAGES)  $(TUTORIAL2A_IMAGES) $(TUTORIAL2B_IMAGES) $(TUTORIAL2C_IMAGES) $(TUTORIAL3A_IMAGES) $(TUTORIAL3B_IMAGES) $(TUTORIAL4_IMAGES) $(TUTORIAL5A_IMAGES) $(TUTORIAL5B_IMAGES)
 
 .PHONY: help clean realclean html dirhtml pickle json htmlhelp qthelp latex changes linkcheck doctest
 
@@ -154,3 +156,15 @@ $(TUTORIAL4_IMAGES): source/images/.tutorial4_flag
 source/images/.tutorial4_flag: source/images/tutorial4.py
 	cd source/images/ && python tutorial4.py
 	touch source/images/.tutorial4_flag
+
+$(TUTORIAL5A_IMAGES): source/images/.tutorial5a_flag
+	@:
+source/images/.tutorial5a_flag: source/images/tutorial5a.py
+	cd source/images/ && python tutorial5a.py
+	touch source/images/.tutorial5a_flag
+
+$(TUTORIAL5B_IMAGES): source/images/.tutorial5b_flag
+	@:
+source/images/.tutorial5b_flag: source/images/tutorial5b.py
+	cd source/images/ && python tutorial5b.py
+	touch source/images/.tutorial5b_flag

doc/source/images/tutorial5a.py

+# Physics background
+# ------------------
+#  band structure of a superconducting quantum wire in tight-binding
+#  approximation
+#
+# Kwant features highlighted
+# --------------------------
+#  - Repetition of previously used concepts (band structure calculations,
+#    matrices as values in Builder).
+#  - Main motivation is to contrast to the implementation of superconductivity
+#    in tutorial5b.py
+
+import kwant
+import latex, html
+
+import numpy as np
+from math import pi
+
+# For plotting
+import pylab
+
+tau_x = np.array([[0,1],[1,0]])
+tau_z = np.array([[1,0],[0,-1]])
+
+def make_lead(a=1, t=1.0, mu=0.7, Delta=0.1, W=10):
+    # Start with an empty lead with a single square lattice
+    lat = kwant.lattice.Square(a)
+
+    sym_lead = kwant.TranslationalSymmetry([lat.vec((-1, 0))])
+    lead = kwant.Builder(sym_lead)
+    lead.default_site_group = lat
+
+    # build up one unit cell of the lead, and add the hoppings
+    # to the next unit cell
+    for j in xrange(W):
+        lead[(0, j)] = (4 * t - mu) * tau_z + Delta * tau_x
+
+        if j > 0:
+            lead[(0, j), (0, j-1)] = - t * tau_z
+
+        lead[(1, j), (0, j)] = - t * tau_z
+
+    # return a finalized lead
+    return lead.finalized()
+
+
+def plot_bandstructure(flead, momenta):
+    # Use the method ``energies`` of the finalized lead to compute
+    # the bandstructure
+    energy_list = [flead.energies(k) for k in momenta]
+
+    pylab.plot(momenta, energy_list)
+    pylab.xlabel("momentum [in untis of (lattice constant)^-1]")
+    pylab.ylabel("energy [in units of t]")
+    pylab.ylim([-0.8, 0.8])
+    fig = pylab.gcf()
+    pylab.setp(fig.get_axes()[0].get_xticklabels(),
+               fontsize=latex.mpl_tick_size)
+    pylab.setp(fig.get_axes()[0].get_yticklabels(),
+               fontsize=latex.mpl_tick_size)
+    fig.set_size_inches(latex.mpl_width_in, latex.mpl_width_in*3./4.)
+    fig.subplots_adjust(left=0.15, right=0.95, top=0.95, bottom=0.15)
+    fig.savefig("tutorial5a_result.pdf")
+    fig.savefig("tutorial5a_result.png",
+                dpi=(html.figwidth_px/latex.mpl_width_in))
+
+
+def main():
+    flead = make_lead()
+
+    # list of momenta at which the bands should be computed
+    momenta = np.arange(-1.5, 1.5 + .0001, 0.002 * pi)
+
+    plot_bandstructure(flead, momenta)
+
+
+# Call the main function if the script gets executed (as opposed to imported).
+# See <http://docs.python.org/library/__main__.html>.
+if __name__ == '__main__':
+    main()

doc/source/images/tutorial5b.py

+# Physics background
+# ------------------
+# - conductance of a NS-junction (Andreev reflection, superconducting gap)
+#
+# Kwant features highlighted
+# --------------------------
+# - Implementing electron and hole ("orbital") degrees of freedom
+#   using different lattices
+
+import kwant
+import latex, html
+
+# For plotting
+import pylab
+
+# For matrix support
+import numpy
+
+def make_system(a=1, W=10, L=10, barrier=1.5, barrierpos=(3, 4),
+                mu=0.4, Delta=0.1, Deltapos=4, t=1.0):
+    # Start with an empty tight-binding system and two square lattices,
+    # corresponding to electron and hole degree of freedom
+    lat_e = kwant.lattice.Square(a)
+    lat_h = kwant.lattice.Square(a)
+
+    sys = kwant.Builder()
+
+    #### Define the scattering region. ####
+    sys[(lat_e(x, y) for x in range(L) for y in range(W))] = 4 * t - mu
+    sys[(lat_h(x, y) for x in range(L) for y in range(W))] = mu - 4 * t
+
+    # the tunnel barrier
+    sys[(lat_e(x, y) for x in range(barrierpos[0], barrierpos[1])
+         for y in range(W))] = 4 * t + barrier - mu
+    sys[(lat_h(x, y) for x in range(barrierpos[0], barrierpos[1])
+         for y in range(W))] = mu - 4 * t - barrier
+
+    # hoppings in x and y-directions, for both electrons and holes
+    sys[sys.possible_hoppings((1, 0), lat_e, lat_e)] = - t
+    sys[sys.possible_hoppings((0, 1), lat_e, lat_e)] = - t
+    sys[sys.possible_hoppings((1, 0), lat_h, lat_h)] = t
+    sys[sys.possible_hoppings((0, 1), lat_h, lat_h)] = t
+
+    # Superconducting order parameter enters as hopping between
+    # electrons and holes
+    sys[((lat_e(x,y), lat_h(x, y)) for x in range(Deltapos, L)
+         for y in range(W))] = Delta
+
+    #### Define the leads. ####
+    # left electron lead
+    sym_lead0 = kwant.TranslationalSymmetry([lat_e.vec((-1, 0))])
+    lead0 = kwant.Builder(sym_lead0)
+
+    lead0[(lat_e(0, j) for j in xrange(W))] = 4 * t - mu
+    # hoppings in x and y-direction
+    lead0[lead0.possible_hoppings((1, 0), lat_e, lat_e)] = - t
+    lead0[lead0.possible_hoppings((0, 1), lat_e, lat_e)] = - t
+
+    # left hole lead
+    sym_lead1 = kwant.TranslationalSymmetry([lat_h.vec((-1, 0))])
+    lead1 = kwant.Builder(sym_lead1)
+
+    lead1[(lat_h(0, j) for j in xrange(W))] = mu - 4 * t
+    # hoppings in x and y-direction
+    lead1[lead1.possible_hoppings((1, 0), lat_h, lat_h)] = t
+    lead1[lead1.possible_hoppings((0, 1), lat_h, lat_h)] = t
+
+    # Then the lead to the right
+    # this one is superconducting and thus is comprised of electrons
+    # AND holes
+    sym_lead2 = kwant.TranslationalSymmetry([lat_e.vec((1, 0))])
+    lead2 = kwant.Builder(sym_lead2)
+
+    lead2[(lat_e(0, j) for j in xrange(W))] = 4 * t - mu
+    lead2[(lat_h(0, j) for j in xrange(W))] = mu - 4 * t
+    # hoppings in x and y-direction
+    lead2[lead2.possible_hoppings((1, 0), lat_e, lat_e)] = - t
+    lead2[lead2.possible_hoppings((0, 1), lat_e, lat_e)] = - t
+    lead2[lead2.possible_hoppings((1, 0), lat_h, lat_h)] = t
+    lead2[lead2.possible_hoppings((0, 1), lat_h, lat_h)] = t
+    lead2[((lat_e(0, j), lat_h(0, j)) for j in xrange(W))] = Delta
+
+    #### Attach the leads and return the finalized system. ####
+    sys.attach_lead(lead0)
+    sys.attach_lead(lead1)
+    sys.attach_lead(lead2)
+
+    return sys.finalized()
+
+def plot_conductance(fsys, energies):
+    # Compute conductance
+    data = []
+    for energy in energies:
+        smatrix = kwant.solve(fsys, energy)
+        # Conductance is N - R_ee + R_he
+        data.append(smatrix.submatrix(0, 0).shape[0] -
+                    smatrix.transmission(0, 0) +
+                    smatrix.transmission(1, 0))
+
+    pylab.plot(energies, data)
+    pylab.xlabel("energy [in units of t]")
+    pylab.ylabel("conductance [in units of e^2/h]")
+    fig = pylab.gcf()
+    pylab.setp(fig.get_axes()[0].get_xticklabels(),
+               fontsize=latex.mpl_tick_size)
+    pylab.setp(fig.get_axes()[0].get_yticklabels(),
+               fontsize=latex.mpl_tick_size)
+    fig.set_size_inches(latex.mpl_width_in, latex.mpl_width_in*3./4.)
+    fig.subplots_adjust(left=0.15, right=0.95, top=0.95, bottom=0.15)
+    fig.savefig("tutorial5b_result.pdf")
+    fig.savefig("tutorial5b_result.png",
+                dpi=(html.figwidth_px/latex.mpl_width_in))
+
+
+
+def main():
+    fsys = make_system()
+
+    plot_conductance(fsys, energies=[0.002 * i for i in xrange(100)])
+
+
+# Call the main function if the script gets executed (as opposed to imported).
+# See <http://docs.python.org/library/__main__.html>.
+if __name__ == '__main__':
+    main()

doc/source/tutorial/index.rst

@@ -13,3 +13,4 @@ these notes maybe safely skipped.
     tutorial2
     tutorial3
     tutorial4
+    tutorial5

doc/source/tutorial/tutorial5.rst

+.. _tutorial-superconductor:
+
+Superconductors: orbital vs lattice degrees of freedom
+------------------------------------------------------
+
+This example deals with superconductivity on the level of the
+Bogoliubov-de Gennes (BdG) equation. In this framework, the Hamiltonian
+is given as
+
+.. math::
+
+    H = \begin{pmatrix} H_0 - \mu& \Delta\\ \Delta^\dagger&\mu-\mathcal{T}H\mathcal{T}^{-1}\end{pmatrix}
+
+where :math:`H_0` is the Hamiltonian of the system without
+superconductivity, :math:`\mu` the chemical potential, :math:`\Delta`
+the superconducting order parameter, and :math:`\mathcal{T}`
+the time-reversal operator. The BdG Hamiltonian introduces
+electron and hole degrees of freedom (an artificial doubling -
+be aware of the fact that electron and hole excitations
+are related!), which we now implement in `kwant`.
+
+For this we restrict ourselves to a simple spin-less system without
+magnetic field, so that :math:`\Delta` is just a number (which we
+choose real), and :math:`\mathcal{T}H\mathcal{T}^{-1}=H_0^*=H_0`.
+
+"Orbital description": Using matrices
+.....................................
+
+We begin by computing the band structure of a superconducting wire.
+The most natural way to implement the BdG Hamiltonian is by using a
+2x2 matrix structure for all Hamiltonian matrix elements:
+
+.. literalinclude:: ../../../examples/tutorial5a.py
+    :lines: 21-45
+
+As you see, the example is syntactically equivalent to our
+:ref:`spin example <tutorial_spinorbit>`, the only difference
+is now that the Pauli matrices act in electron-hole space.
+
+Computing the band structure then yields the result
+
+.. image:: ../images/tutorial5a_result.*
+
+We clearly observe the superconducting gap in the spectrum. That was easy,
+he?
+
+.. seealso::
+    The full source code can be found in
+    :download:`examples/tutorial5a.py <../../../examples/tutorial5a.py>`
+
+
+"Lattice description": Using different lattices
+...............................................
+
+While it seems most natural to implement the BdG Hamiltonian
+using a 2x2 matrix structure for the matrix elements of the Hamiltonian,
+we run into a problem when we want to compute electronic transport in
+a system consisting of a normal and a superconducting lead:
+Since electrons and holes carry charge with opposite sign, we need to
+separate electron and hole degrees of freedom in the scattering matrix.
+In particular, the conductance of a N-S-junction is given as
+
+.. math::
+
+    G = \frac{e^2}{h} (N - R_\text{ee} + R_\text{he})\,,
+
+where :math:`N` is the number of channels in the normal lead, and
+:math:`R_\text{ee}` the total probability of reflection from electrons
+to electrons in the normal lead, and :math:`R_\text{eh}` the total
+probability of reflection from electrons to holes in the normal
+lead. However, the current version of kwant does not allow for an easy
+and elegant partitioning of the scattering matrix in these two degrees
+of freedom (well, there is one since v0.1.3, see the technical notes
+below).
+
+In the following, we will circumvent this problem by introducing
+separate "leads" for electrons and holes, making use of different
+lattices. The system we consider consists of a normal lead on the left,
+a superconductor on the right, and a tunnel barrier inbetween:
+
+.. image:: ../images/tutorial5b_sketch.*
+
+As already mentioned above, we begin by introducing two different
+square lattices representing electron and hole degrees of freedom:
+
+.. literalinclude:: ../../../examples/tutorial5b.py
+    :lines: 18-19,17,23-24
+
+Any diagonal entry (kinetic energy, potentials, ...) in the BdG
+Hamiltonian then corresponds to on-site energies or hoppings within
+the *same* lattice, whereas any off-diagonal entry (essentially, the
+superconducting order parameter :math:`\Delta`) corresponds
+to a hopping between *different* lattices:
+
+.. literalinclude:: ../../../examples/tutorial5b.py
+    :lines: 25-46
+
+Note that the tunnel barrier is added by overwriting previously set
+on-site matrix elements.
+
+Note further, that in the code above, the superconducting order
+parameter is nonzero only in a part of the scattering region.
+Consequently, we have added hoppings between electron and hole
+lattices only in this region, they remain uncoupled in the normal
+part. We use this fact to attach purely electron and hole leads
+(comprised of only electron *or* hole lattices) to the
+system:
+
+.. literalinclude:: ../../../examples/tutorial5b.py
+    :lines: 49-65
+
+This separation into two different leads allows us then later to compute the
+reflection probablities between electrons and holes explicitely.
+
+On the superconducting side, we cannot do this separation, and can
+only define a single lead coupling electrons and holes:
+
+.. literalinclude:: ../../../examples/tutorial5b.py
+    :lines: 70-80
+
+We now have on the left side two leads that are sitting in the same
+spatial position, but in different lattice spaces. This ensures that
+we can still attach all leads as before:
+
+.. literalinclude:: ../../../examples/tutorial5b.py
+    :lines: 83-87
+
+When computing the conductance, we can now extract reflection from
+electrons to electrons as ``smatrix.transmission(0, 0)`` (Don't get
+confused by the fact that it says ``transmission`` -- transmission
+into the same lead is reflection), and reflection from electrons to holes
+as ``smatrix.transmission(1, 0)``, by virtue of our electron and hole leads:
+
+.. literalinclude:: ../../../examples/tutorial5b.py
+    :lines: 89-97
+
+Note that ``smatrix.submatrix(0,0)`` returns the block concerning reflection
+within (electron) lead 0, and from its size we can extract the number of modes
+:math:`N`.
+
+Finally, for the default parameters, we obtain the following result:
+
+.. image:: ../images/tutorial5b_result.*
+
+We a see a conductance that is proportional to the square of the tunneling
+probability within the gap, and proportional to the tunneling probability
+above the gap. At the gap edge, we observe a resonant Andreev reflection.
+
+.. seealso::
+    The full source code can be found in
+    :download:`examples/tutorial5b.py <../../../examples/tutorial5b.py>`
+
+.. specialnote:: Technical details
+
+    - If you are only interested in particle (thermal) currents you do not need
+      to define separate electron and hole leads. In this case, you do not
+      need to distinguish them. Still, separating the leads into electron
+      and hole degrees of freedom makes the lead calculation in the solving
+      phase more efficient.
+
+    - It is in fact possible to separate electron and hole degrees of
+      freedom in the scattering matrix, even if one uses matrices for
+      these degrees of freedom. In the solve step,
+      `~kwant.solvers.sparse.solve` returns an array containing the
+      transverse wave functions of the lead modes, if
+      ``return_modes=True``. By inspecting the wave functions,
+      electron and hole wave functions can be distinguished (they only
+      have entries in either the electron part *or* the hole part. If
+      you encounter modes with entries in both parts, you hit a very
+      unlikely situation in which the standard procedure to compute
+      the modes gave you a superposition of electron and hole
+      modes. That is still OK for computing particle current, but not
+      for electrical current).

doc/source/whatsnew/0.2.rst

@@ -3,6 +3,10 @@ What's New in kwant 0.2
 
 This article explains the user-visible changes in kwant 0.2.
 
+New tutorial dealing with superconductivity
+-------------------------------------------
+:doc:`../tutorial/tutorial5`
+
 `~kwant.plotter.plot` more useful for low level systems
 -------------------------------------------------------
 The behavior of `~kwant.plotter.plot` has been changed when a `low level system

examples/tutorial5a.py

+# Physics background
+# ------------------
+#  band structure of a superconducting quantum wire in tight-binding
+#  approximation
+#
+# Kwant features highlighted
+# --------------------------
+#  - Repetition of previously used concepts (band structure calculations,
+#    matrices as values in Builder).
+#  - Main motivation is to contrast to the implementation of superconductivity
+#    in tutorial5b.py
+
+import kwant
+
+import numpy as np
+from math import pi
+
+# For plotting
+import pylab
+
+tau_x = np.array([[0,1],[1,0]])
+tau_z = np.array([[1,0],[0,-1]])
+
+def make_lead(a=1, t=1.0, mu=0.7, Delta=0.1, W=10):
+    # Start with an empty lead with a single square lattice
+    lat = kwant.lattice.Square(a)
+
+    sym_lead = kwant.TranslationalSymmetry([lat.vec((-1, 0))])
+    lead = kwant.Builder(sym_lead)
+    lead.default_site_group = lat
+
+    # build up one unit cell of the lead, and add the hoppings
+    # to the next unit cell
+    for j in xrange(W):
+        lead[(0, j)] = (4 * t - mu) * tau_z + Delta * tau_x
+
+        if j > 0:
+            lead[(0, j), (0, j-1)] = - t * tau_z
+
+        lead[(1, j), (0, j)] = - t * tau_z
+
+    # return a finalized lead
+    return lead.finalized()
+
+
+def plot_bandstructure(flead, momenta):
+    # Use the method ``energies`` of the finalized lead to compute
+    # the bandstructure
+    energy_list = [flead.energies(k) for k in momenta]
+
+    pylab.plot(momenta, energy_list)
+    pylab.xlabel("momentum [in untis of (lattice constant)^-1]")
+    pylab.ylabel("energy [in units of t]")
+    pylab.ylim([-0.8, 0.8])
+    pylab.show()
+
+
+def main():
+    flead = make_lead()
+
+    # list of momenta at which the bands should be computed
+    momenta = np.arange(-1.5, 1.5 + .0001, 0.002 * pi)
+
+    plot_bandstructure(flead, momenta)
+
+
+# Call the main function if the script gets executed (as opposed to imported).
+# See <http://docs.python.org/library/__main__.html>.
+if __name__ == '__main__':
+    main()

examples/tutorial5b.py

+# Physics background
+# ------------------
+# - conductance of a NS-junction (Andreev reflection, superconducting gap)
+#
+# Kwant features highlighted
+# --------------------------
+# - Implementing electron and hole ("orbital") degrees of freedom
+#   using different lattices
+
+import kwant
+
+# For plotting
+import pylab
+
+# For matrix support
+import numpy
+
+def make_system(a=1, W=10, L=10, barrier=1.5, barrierpos=(3, 4),
+                mu=0.4, Delta=0.1, Deltapos=4, t=1.0):
+    # Start with an empty tight-binding system and two square lattices,
+    # corresponding to electron and hole degree of freedom
+    lat_e = kwant.lattice.Square(a)
+    lat_h = kwant.lattice.Square(a)
+
+    sys = kwant.Builder()
+
+    #### Define the scattering region. ####
+    sys[(lat_e(x, y) for x in range(L) for y in range(W))] = 4 * t - mu
+    sys[(lat_h(x, y) for x in range(L) for y in range(W))] = mu - 4 * t
+
+    # the tunnel barrier
+    sys[(lat_e(x, y) for x in range(barrierpos[0], barrierpos[1])
+         for y in range(W))] = 4 * t + barrier - mu
+    sys[(lat_h(x, y) for x in range(barrierpos[0], barrierpos[1])
+         for y in range(W))] = mu - 4 * t - barrier
+
+    # hoppings in x and y-directions, for both electrons and holes
+    sys[sys.possible_hoppings((1, 0), lat_e, lat_e)] = - t
+    sys[sys.possible_hoppings((0, 1), lat_e, lat_e)] = - t
+    sys[sys.possible_hoppings((1, 0), lat_h, lat_h)] = t
+    sys[sys.possible_hoppings((0, 1), lat_h, lat_h)] = t
+
+    # Superconducting order parameter enters as hopping between
+    # electrons and holes
+    sys[((lat_e(x,y), lat_h(x, y)) for x in range(Deltapos, L)
+         for y in range(W))] = Delta
+
+    #### Define the leads. ####
+    # left electron lead
+    sym_lead0 = kwant.TranslationalSymmetry([lat_e.vec((-1, 0))])
+    lead0 = kwant.Builder(sym_lead0)
+
+    lead0[(lat_e(0, j) for j in xrange(W))] = 4 * t - mu
+    # hoppings in x and y-direction
+    lead0[lead0.possible_hoppings((1, 0), lat_e, lat_e)] = - t
+    lead0[lead0.possible_hoppings((0, 1), lat_e, lat_e)] = - t
+
+    # left hole lead
+    sym_lead1 = kwant.TranslationalSymmetry([lat_h.vec((-1, 0))])
+    lead1 = kwant.Builder(sym_lead1)
+
+    lead1[(lat_h(0, j) for j in xrange(W))] = mu - 4 * t
+    # hoppings in x and y-direction
+    lead1[lead1.possible_hoppings((1, 0), lat_h, lat_h)] = t
+    lead1[lead1.possible_hoppings((0, 1), lat_h, lat_h)] = t
+
+    # Then the lead to the right
+    # this one is superconducting and thus is comprised of electrons
+    # AND holes
+    sym_lead2 = kwant.TranslationalSymmetry([lat_e.vec((1, 0))])
+    lead2 = kwant.Builder(sym_lead2)
+
+    lead2[(lat_e(0, j) for j in xrange(W))] = 4 * t - mu
+    lead2[(lat_h(0, j) for j in xrange(W))] = mu - 4 * t
+    # hoppings in x and y-direction
+    lead2[lead2.possible_hoppings((1, 0), lat_e, lat_e)] = - t
+    lead2[lead2.possible_hoppings((0, 1), lat_e, lat_e)] = - t
+    lead2[lead2.possible_hoppings((1, 0), lat_h, lat_h)] = t
+    lead2[lead2.possible_hoppings((0, 1), lat_h, lat_h)] = t
+    lead2[((lat_e(0, j), lat_h(0, j)) for j in xrange(W))] = Delta
+
+    #### Attach the leads and return the finalized system. ####
+    sys.attach_lead(lead0)
+    sys.attach_lead(lead1)
+    sys.attach_lead(lead2)
+
+    return sys.finalized()
+
+def plot_conductance(fsys, energies):
+    # Compute conductance
+    data = []
+    for energy in energies:
+        smatrix = kwant.solve(fsys, energy)
+        # Conductance is N - R_ee + R_he
+        data.append(smatrix.submatrix(0, 0).shape[0] -
+                    smatrix.transmission(0, 0) +
+                    smatrix.transmission(1, 0))
+
+    pylab.plot(energies, data)
+    pylab.xlabel("energy [in units of t]")
+    pylab.ylabel("conductance [in units of e^2/h]")
+    pylab.show()
+
+
+def main():
+    fsys = make_system()
+
+    # Check that the system looks as intended.
+    kwant.plot(fsys)
+
+    plot_conductance(fsys, energies=[0.002 * i for i in xrange(100)])
+
+
+# Call the main function if the script gets executed (as opposed to imported).
+# See <http://docs.python.org/library/__main__.html>.
+if __name__ == '__main__':
+    main()