doc/source/code/figure/quantum_wire.py.diff

-@@ -1,121 +1,130 @@
- # Tutorial 2.2.2. Transport through a quantum wire
- # ================================================
- #
- # Physics background
- # ------------------
- #  Conductance of a quantum wire; subbands
- #
- # Kwant features highlighted
- # --------------------------
- #  - Builder for setting up transport systems easily
- #  - Making scattering region and leads
- #  - Using the simple sparse solver for computing Landauer conductance
- 
-+import _defs
- from matplotlib import pyplot
- #HIDDEN_BEGIN_dwhx
- import kwant
- #HIDDEN_END_dwhx
- 
- # First, define the tight-binding system
- 
- #HIDDEN_BEGIN_goiq
- syst = kwant.Builder()
- #HIDDEN_END_goiq
- 
- # Here, we are only working with square lattices
- #HIDDEN_BEGIN_suwo
- a = 1
- lat = kwant.lattice.square(a)
- #HIDDEN_END_suwo
- 
- #HIDDEN_BEGIN_zfvr
- t = 1.0
- W = 10
- L = 30
- 
- # Define the scattering region
- 
- for i in range(L):
-     for j in range(W):
-         # On-site Hamiltonian
-         syst[lat(i, j)] = 4 * t
- 
-         # Hopping in y-direction
-         if j > 0:
-             syst[lat(i, j), lat(i, j - 1)] = -t
- 
-         # Hopping in x-direction
-         if i > 0:
-             syst[lat(i, j), lat(i - 1, j)] = -t
- #HIDDEN_END_zfvr
- 
- # Then, define and attach the leads:
- 
- # First the lead to the left
- # (Note: TranslationalSymmetry takes a real-space vector)
- #HIDDEN_BEGIN_xcmc
- sym_left_lead = kwant.TranslationalSymmetry((-a, 0))
- left_lead = kwant.Builder(sym_left_lead)
- #HIDDEN_END_xcmc
- 
- #HIDDEN_BEGIN_ndez
- for j in range(W):
-     left_lead[lat(0, j)] = 4 * t
-     if j > 0:
-         left_lead[lat(0, j), lat(0, j - 1)] = -t
-     left_lead[lat(1, j), lat(0, j)] = -t
- #HIDDEN_END_ndez
- 
- #HIDDEN_BEGIN_fskr
- syst.attach_lead(left_lead)
- #HIDDEN_END_fskr
- 
- # Then the lead to the right
- #HIDDEN_BEGIN_xhqc
- sym_right_lead = kwant.TranslationalSymmetry((a, 0))
- right_lead = kwant.Builder(sym_right_lead)
- 
- for j in range(W):
-     right_lead[lat(0, j)] = 4 * t
-     if j > 0:
-         right_lead[lat(0, j), lat(0, j - 1)] = -t
-     right_lead[lat(1, j), lat(0, j)] = -t
- 
- syst.attach_lead(right_lead)
- #HIDDEN_END_xhqc
- 
- # Plot it, to make sure it's OK
- #HIDDEN_BEGIN_wsgh
--kwant.plot(syst)
-+size = (_defs.figwidth_in, 0.3 * _defs.figwidth_in)
-+for extension in ('pdf', 'png'):
-+    kwant.plot(syst, file="quantum_wire_syst." + extension,
-+               fig_size=size, dpi=_defs.dpi)
- #HIDDEN_END_wsgh
- 
- # Finalize the system
- #HIDDEN_BEGIN_dngj
- syst = syst.finalized()
- #HIDDEN_END_dngj
- 
- # Now that we have the system, we can compute conductance
- #HIDDEN_BEGIN_buzn
- energies = []
- data = []
- for ie in range(100):
-     energy = ie * 0.01
- 
-     # compute the scattering matrix at a given energy
-     smatrix = kwant.smatrix(syst, energy)
- 
-     # compute the transmission probability from lead 0 to
-     # lead 1
-     energies.append(energy)
-     data.append(smatrix.transmission(1, 0))
- #HIDDEN_END_buzn
- 
- # Use matplotlib to write output
- # We should see conductance steps
- #HIDDEN_BEGIN_lliv
--pyplot.figure()
-+fig = pyplot.figure()
- pyplot.plot(energies, data)
--pyplot.xlabel("energy [t]")
--pyplot.ylabel("conductance [e^2/h]")
--pyplot.show()
-+pyplot.xlabel("energy [t]", fontsize=_defs.mpl_label_size)
-+pyplot.ylabel("conductance [e^2/h]", fontsize=_defs.mpl_label_size)
-+pyplot.setp(fig.get_axes()[0].get_xticklabels(), fontsize=_defs.mpl_tick_size)
-+pyplot.setp(fig.get_axes()[0].get_yticklabels(), 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)
-+for extension in ('pdf', 'png'):
-+    fig.savefig("quantum_wire_result." + extension, dpi=_defs.dpi)
- #HIDDEN_END_lliv

doc/source/code/figure/spin_orbit.py.diff

-@@ -1,104 +1,110 @@
- # Tutorial 2.3.1. Matrix structure of on-site and hopping elements
- # ================================================================
- #
- # Physics background
- # ------------------
- #  Gaps in quantum wires with spin-orbit coupling and Zeeman splititng,
- #  as theoretically predicted in
- #   http://prl.aps.org/abstract/PRL/v90/i25/e256601
- #  and (supposedly) experimentally oberved in
- #   http://www.nature.com/nphys/journal/v6/n5/abs/nphys1626.html
- #
- # Kwant features highlighted
- # --------------------------
- #  - Numpy matrices as values in Builder
- 
-+import _defs
- import kwant
- 
- # For plotting
- from matplotlib import pyplot
- 
- # For matrix support
- #HIDDEN_BEGIN_xumz
- import tinyarray
- #HIDDEN_END_xumz
- 
- # define Pauli-matrices for convenience
- #HIDDEN_BEGIN_hwbt
- sigma_0 = tinyarray.array([[1, 0], [0, 1]])
- sigma_x = tinyarray.array([[0, 1], [1, 0]])
- sigma_y = tinyarray.array([[0, -1j], [1j, 0]])
- sigma_z = tinyarray.array([[1, 0], [0, -1]])
- #HIDDEN_END_hwbt
- 
- 
- def make_system(t=1.0, alpha=0.5, e_z=0.08, W=10, L=30):
-     # Start with an empty tight-binding system and a single square lattice.
-     # `a` is the lattice constant (by default set to 1 for simplicity).
-     lat = kwant.lattice.square()
- 
-     syst = kwant.Builder()
- 
-     #### Define the scattering region. ####
- #HIDDEN_BEGIN_uxrm
-     syst[(lat(x, y) for x in range(L) for y in range(W))] = \
-         4 * t * sigma_0 + e_z * sigma_z
-     # hoppings in x-direction
-     syst[kwant.builder.HoppingKind((1, 0), lat, lat)] = \
-         -t * sigma_0 + 1j * alpha * sigma_y / 2
-     # hoppings in y-directions
-     syst[kwant.builder.HoppingKind((0, 1), lat, lat)] = \
-         -t * sigma_0 - 1j * alpha * sigma_x / 2
- #HIDDEN_END_uxrm
- 
-     #### Define the left lead. ####
-     lead = kwant.Builder(kwant.TranslationalSymmetry((-1, 0)))
- 
- #HIDDEN_BEGIN_yliu
-     lead[(lat(0, j) for j in range(W))] = 4 * t * sigma_0 + e_z * sigma_z
-     # hoppings in x-direction
-     lead[kwant.builder.HoppingKind((1, 0), lat, lat)] = \
-         -t * sigma_0 + 1j * alpha * sigma_y / 2
-     # hoppings in y-directions
-     lead[kwant.builder.HoppingKind((0, 1), lat, lat)] = \
-         -t * sigma_0 - 1j * alpha * sigma_x / 2
- #HIDDEN_END_yliu
- 
-     #### Attach the leads and return the finalized system. ####
-     syst.attach_lead(lead)
-     syst.attach_lead(lead.reversed())
- 
-     return syst
- 
- 
- def plot_conductance(syst, energies):
-     # Compute conductance
-     data = []
-     for energy in energies:
-         smatrix = kwant.smatrix(syst, energy)
-         data.append(smatrix.transmission(1, 0))
- 
--    pyplot.figure()
-+    fig = pyplot.figure()
-     pyplot.plot(energies, data)
--    pyplot.xlabel("energy [t]")
--    pyplot.ylabel("conductance [e^2/h]")
--    pyplot.show()
-+    pyplot.xlabel("energy [t]", fontsize=_defs.mpl_label_size)
-+    pyplot.ylabel("conductance [e^2/h]",
-+                 fontsize=_defs.mpl_label_size)
-+    pyplot.setp(fig.get_axes()[0].get_xticklabels(),
-+               fontsize=_defs.mpl_tick_size)
-+    pyplot.setp(fig.get_axes()[0].get_yticklabels(),
-+               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)
-+    for extension in ('pdf', 'png'):
-+        fig.savefig("spin_orbit_result." + extension, dpi=_defs.dpi)
- 
- 
- def main():
-     syst = make_system()
- 
--    # Check that the system looks as intended.
--    kwant.plot(syst)
--
-     # Finalize the system.
-     syst = syst.finalized()
- 
-     # We should see non-monotonic conductance steps.
-     plot_conductance(syst, energies=[0.01 * i - 0.3 for i in range(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/code/figure/superconductor.py.diff

-@@ -1,132 +1,141 @@
- # Tutorial 2.6. "Superconductors": orbitals, conservation laws and symmetries
- # ===========================================================================
- #
- # Physics background
- # ------------------
- # - conductance of a NS-junction (Andreev reflection, superconducting gap)
- #
- # Kwant features highlighted
- # --------------------------
- # - Implementing electron and hole ("orbital") degrees of freedom
- #   using conservation laws.
- # - Use of discrete symmetries to relate scattering states.
- 
-+import _defs
- import kwant
- 
- import tinyarray
- import numpy as np
- 
- # For plotting
- from matplotlib import pyplot
-+from contextlib import redirect_stdout
- 
- tau_x = tinyarray.array([[0, 1], [1, 0]])
- tau_y = tinyarray.array([[0, -1j], [1j, 0]])
- tau_z = tinyarray.array([[1, 0], [0, -1]])
- 
- #HIDDEN_BEGIN_nbvn
- 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. On each site, there
-     # are now electron and hole orbitals, so we must specify the
-     # number of orbitals per site. The orbital structure is the same
-     # as in the Hamiltonian.
-     lat = kwant.lattice.square(norbs=2)
-     syst = kwant.Builder()
- 
-     #### Define the scattering region. ####
-     # The superconducting order parameter couples electron and hole orbitals
-     # on each site, and hence enters as an onsite potential.
-     # The pairing is only included beyond the point 'Deltapos' in the scattering region.
-     syst[(lat(x, y) for x in range(Deltapos) for y in range(W))] = (4 * t - mu) * tau_z
-     syst[(lat(x, y) for x in range(Deltapos, L) for y in range(W))] = (4 * t - mu) * tau_z + Delta * tau_x
- 
-     # The tunnel barrier
-     syst[(lat(x, y) for x in range(barrierpos[0], barrierpos[1])
-          for y in range(W))] = (4 * t + barrier - mu) * tau_z
- 
-     # Hoppings
-     syst[lat.neighbors()] = -t * tau_z
- #HIDDEN_END_nbvn
- #HIDDEN_BEGIN_ttth
-     #### Define the leads. ####
-     # Left lead - normal, so the order parameter is zero.
-     sym_left = kwant.TranslationalSymmetry((-a, 0))
-     # Specify the conservation law used to treat electrons and holes separately.
-     # We only do this in the left lead, where the pairing is zero.
-     lead0 = kwant.Builder(sym_left, conservation_law=-tau_z, particle_hole=tau_y)
-     lead0[(lat(0, j) for j in range(W))] = (4 * t - mu) * tau_z
-     lead0[lat.neighbors()] = -t * tau_z
- #HIDDEN_END_ttth
- #HIDDEN_BEGIN_zuuw
-     # Right lead - superconducting, so the order parameter is included.
-     sym_right = kwant.TranslationalSymmetry((a, 0))
-     lead1 = kwant.Builder(sym_right)
-     lead1[(lat(0, j) for j in range(W))] = (4 * t - mu) * tau_z + Delta * tau_x
-     lead1[lat.neighbors()] = -t * tau_z
- 
-     #### Attach the leads and return the system. ####
-     syst.attach_lead(lead0)
-     syst.attach_lead(lead1)
- 
-     return syst
- #HIDDEN_END_zuuw
- 
- #HIDDEN_BEGIN_jbjt
- def plot_conductance(syst, energies):
-     # Compute conductance
-     data = []
-     for energy in energies:
-         smatrix = kwant.smatrix(syst, energy)
-         # Conductance is N - R_ee + R_he
-         data.append(smatrix.submatrix((0, 0), (0, 0)).shape[0] -
-                     smatrix.transmission((0, 0), (0, 0)) +
-                     smatrix.transmission((0, 1), (0, 0)))
- #HIDDEN_END_jbjt
--    pyplot.figure()
-+    fig = pyplot.figure()
-     pyplot.plot(energies, data)
-     pyplot.xlabel("energy [t]")
-     pyplot.ylabel("conductance [e^2/h]")
--    pyplot.show()
-+    pyplot.setp(fig.get_axes()[0].get_xticklabels(),
-+                fontsize=_defs.mpl_tick_size)
-+    pyplot.setp(fig.get_axes()[0].get_yticklabels(),
-+                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)
-+    for extension in ('pdf', 'png'):
-+        fig.savefig("superconductor_transport_result." + extension,
-+                    dpi=_defs.dpi)
- 
- #HIDDEN_BEGIN_pqmp
- def check_PHS(syst):
-     # Scattering matrix
-     s = kwant.smatrix(syst, energy=0)
-     # Electron to electron block
-     s_ee = s.submatrix((0,0), (0,0))
-     # Hole to hole block
-     s_hh = s.submatrix((0,1), (0,1))
-     print('s_ee: \n', np.round(s_ee, 3))
-     print('s_hh: \n', np.round(s_hh[::-1, ::-1], 3))
-     print('s_ee - s_hh^*: \n',
-           np.round(s_ee - s_hh[::-1, ::-1].conj(), 3), '\n')
-     # Electron to hole block
-     s_he = s.submatrix((0,1), (0,0))
-     # Hole to electron block
-     s_eh = s.submatrix((0,0), (0,1))
-     print('s_he: \n', np.round(s_he, 3))
-     print('s_eh: \n', np.round(s_eh[::-1, ::-1], 3))
-     print('s_he + s_eh^*: \n',
-           np.round(s_he + s_eh[::-1, ::-1].conj(), 3))
- #HIDDEN_END_pqmp
- 
- def main():
-     syst = make_system(W=10)
- 
--    # Check that the system looks as intended.
--    kwant.plot(syst)
--
-     # Finalize the system.
-     syst = syst.finalized()
- 
-     # Check particle-hole symmetry of the scattering matrix
--    check_PHS(syst)
-+    with open('check_PHS_out.txt', 'w') as f:
-+        with redirect_stdout(f):
-+            check_PHS(syst)
- 
-     # Compute and plot the conductance
-     plot_conductance(syst, energies=[0.002 * i for i in range(-10, 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/code/include/quantum_wire_revisited.py

-# Tutorial 2.2.3. Building the same system with less code
-# =======================================================
-#
-# Physics background
-# ------------------
-#  Conductance of a quantum wire; subbands
-#
-# Kwant features highlighted
-# --------------------------
-#  - Using iterables and builder.HoppingKind for making systems
-#  - introducing `reversed()` for the leads
-#
-# Note: Does the same as tutorial1a.py, but using other features of Kwant.
-
-#HIDDEN_BEGIN_xkzy
-import kwant
-
-# For plotting
-from matplotlib import pyplot
-
-
-def make_system(a=1, t=1.0, W=10, L=30):
-    # Start with an empty tight-binding system and a single square lattice.
-    # `a` is the lattice constant (by default set to 1 for simplicity.
-    lat = kwant.lattice.square(a)
-
-    syst = kwant.Builder()
-#HIDDEN_END_xkzy
-
-    #### Define the scattering region. ####
-#HIDDEN_BEGIN_vvjt
-    syst[(lat(x, y) for x in range(L) for y in range(W))] = 4 * t
-#HIDDEN_END_vvjt
-#HIDDEN_BEGIN_nooi
-    syst[lat.neighbors()] = -t
-#HIDDEN_END_nooi
-
-    #### Define and attach the leads. ####
-    # Construct the left lead.
-#HIDDEN_BEGIN_iepx
-    lead = kwant.Builder(kwant.TranslationalSymmetry((-a, 0)))
-    lead[(lat(0, j) for j in range(W))] = 4 * t
-    lead[lat.neighbors()] = -t
-#HIDDEN_END_iepx
-
-    # Attach the left lead and its reversed copy.
-#HIDDEN_BEGIN_yxot
-    syst.attach_lead(lead)
-    syst.attach_lead(lead.reversed())
-
-    return syst
-#HIDDEN_END_yxot
-
-
-#HIDDEN_BEGIN_ayuk
-def plot_conductance(syst, energies):
-    # Compute conductance
-    data = []
-    for energy in energies:
-        smatrix = kwant.smatrix(syst, energy)
-        data.append(smatrix.transmission(1, 0))
-
-    pyplot.figure()
-    pyplot.plot(energies, data)
-    pyplot.xlabel("energy [t]")
-    pyplot.ylabel("conductance [e^2/h]")
-    pyplot.show()
-#HIDDEN_END_ayuk
-
-
-#HIDDEN_BEGIN_cjel
-def main():
-    syst = make_system()
-
-    # Check that the system looks as intended.
-    kwant.plot(syst)
-
-    # Finalize the system.
-    syst = syst.finalized()
-
-    # We should see conductance steps.
-    plot_conductance(syst, energies=[0.01 * i for i in range(100)])
-#HIDDEN_END_cjel
-
-
-# Call the main function if the script gets executed (as opposed to imported).
-# See <http://docs.python.org/library/__main__.html>.
-#HIDDEN_BEGIN_ypbj
-if __name__ == '__main__':
-    main()
-#HIDDEN_END_ypbj

doc/source/conf.py

@@ -15,8 +15,15 @@ import sys
 import os
 import string
 from distutils.util import get_platform
-sys.path.insert(0, "../../build/lib.{0}-{1}.{2}".format(
-        get_platform(), *sys.version_info[:2]))
+
+package_path = os.path.abspath(
+    "../../build/lib.{0}-{1}.{2}"
+    .format(get_platform(), *sys.version_info[:2]))
+
+# Insert into sys.path so that we can import kwant here
+sys.path.insert(0, package_path)
+# Insert into PYTHONPATH so that jupyter-sphinx will pick it up
+os.environ['PYTHONPATH'] = ':'.join((package_path, os.environ.get('PYTHONPATH','')))
 
 import kwant
 import kwant.qsymm
@@ -31,7 +38,8 @@ sys.path.insert(0, os.path.abspath('../sphinxext'))
 
 extensions = ['sphinx.ext.autodoc', 'sphinx.ext.autosummary',
               'sphinx.ext.todo', 'sphinx.ext.mathjax', 'numpydoc',
-              'kwantdoc', 'sphinx.ext.linkcode']
+              'kwantdoc', 'sphinx.ext.linkcode', 'jupyter_sphinx.execute',
+              'sphinxcontrib.rsvgconverter']
 
 # Add any paths that contain templates here, relative to this directory.
 templates_path = ['../templates']
@@ -239,7 +247,6 @@ autosummary_generate = True
 autoclass_content = "both"
 autodoc_default_flags = ['show-inheritance']
 
-
 # -- Teach Sphinx to document bound methods like functions ---------------------
 import types
 from sphinx.ext import autodoc

doc/source/pre/whatsnew/1.3.rst

@@ -28,11 +28,43 @@ tight-binding band structures or construct systems with different/lower
 symmetry. For example in just a few lines we can construct a two-band model that exhibits
 a quantum anomalous spin Hall phase:
 
-.. literalinclude:: ../../code/include/plot_qahe.py
-    :start-after: HIDDEN_BEGIN_model
-    :end-before: HIDDEN_END_model
+.. jupyter-kernel::
+    :id: plot_qahe
+
+.. jupyter-execute::
+    :hide-code:
+
+    # Comprehensive example: quantum anomalous Hall effect
+    # ====================================================
+    #
+    # Physics background
+    # ------------------
+    # + Quantum anomalous Hall effect
+    #
+    # Features highlighted
+    # --------------------
+    # + Use of `kwant.continuum` to discretize a continuum Hamiltonian
+    # + Use of `kwant.operator` to compute local current
+    # + Use of `kwant.plotter.current` to plot local current
+
+    import math
+    import matplotlib.pyplot
+    import kwant
+    import kwant.continuum
+
+.. jupyter-execute:: ../../tutorial/boilerplate.py
+    :hide-code:
+
+.. jupyter-execute::
 
-From: :download:`QAHE example script <../../code/download/plot_qahe.py>`
+    def make_model(a):
+        ham = ("alpha * (k_x * sigma_x - k_y * sigma_y)"
+               "+ (m + beta * kk) * sigma_z"
+               "+ (gamma * kk + U) * sigma_0")
+        subs = {"kk": "k_x**2 + k_y**2"}
+        return kwant.continuum.discretize(ham, locals=subs, grid=a)
+
+From: :jupyter-download:script:`plot_qahe`
 
 See the tutorial: :doc:`../../tutorial/discretize`
 
@@ -71,13 +103,47 @@ The example below shows edge states of a quantum anomalous Hall phase
 of the two-band model shown in the `above section
 <#tools-for-continuum-hamiltonians>`_:
 
-.. literalinclude:: ../../code/include/plot_qahe.py
-    :start-after: HIDDEN_BEGIN_current
-    :end-before: HIDDEN_END_current
+.. jupyter-execute::
+    :hide-code:
+
+    def make_system(model, L):
+        def lead_shape(site):
+            x, y = site.pos / L
+            return abs(y) < 0.5
+
+        # QPC shape: a rectangle with 2 gaussians
+        # etched out of the top and bottom edge.
+        def central_shape(site):
+            x, y = site.pos / L
+            return abs(x) < 3/5 and abs(y) < 0.5 - 0.4 * math.exp(-40 * x**2)
+
+        lead = kwant.Builder(kwant.TranslationalSymmetry(
+            model.lattice.vec((-1, 0))))
+        lead.fill(model, lead_shape, (0, 0))
+
+        syst = kwant.Builder()
+        syst.fill(model, central_shape, (0, 0))
+        syst.attach_lead(lead)
+        syst.attach_lead(lead.reversed())
+
+        return syst.finalized()
+
+    # Set up our model and system, and define the model parameters.
+    params = dict(alpha=0.365, beta=0.686, gamma=0.512, m=-0.01, U=0)
+    model = make_model(1)
+    syst = make_system(model, 70)
+
+    # Calculate the scattering states at energy 'm' coming from the left
+    # lead, and the associated particle current.
+    psi = kwant.wave_function(syst, energy=params['m'], params=params)(0)
+
+.. jupyter-execute::
 
-.. image:: ../../code/figure/plot_qahe_current.*
+    J = kwant.operator.Current(syst).bind(params=params)
+    current = sum(J(p) for p in psi)
+    kwant.plotter.current(syst, current);
 
-From: :download:`QAHE example script <../../code/download/plot_qahe.py>`
+From: :jupyter-download:script:`plot_qahe`
 
 Scattering states with discrete symmetries and conservation laws
 ----------------------------------------------------------------

doc/source/pre/whatsnew/1.5.rst

+What's new in Kwant 1.5
+=======================
+
+This article explains the user-visible changes in Kwant 1.5.0.
+
+Improved tutorial building
+--------------------------
+Improving or adding to Kwant's tutorial is now much simpler. Now
+the text and code for each tutorial is kept in the same file, making
+it easy to see where changes need to be made, and images generated by
+the code are inserted directly into the document thanks to the magic of
+`jupyter-sphinx <https://github.com/jupyter-widgets/jupyter-sphinx/>`_.
+It has never been easier to get started contributing to Kwant by
+helping us improve our documentation.

doc/source/pre/whatsnew/index.rst

@@ -2,6 +2,7 @@ What's new in Kwant
 ===================
 
 .. toctree::
+    1.5
     1.4
     1.3
     1.2

doc/source/tutorial/boilerplate.py

+import matplotlib
+import matplotlib.pyplot
+from IPython.display import set_matplotlib_formats
+
+matplotlib.rcParams['figure.figsize'] = matplotlib.pyplot.figaspect(1) * 2
+set_matplotlib_formats('svg')

doc/source/tutorial/discretize.rst

@@ -18,8 +18,32 @@ continuum models and for discretizing them into tight-binding models.
 
 .. seealso::
     The complete source code of this tutorial can be found in
-    :download:`discretize.py </code/download/discretize.py>`
+    :jupyter-download:script:`discretize`
 
+.. jupyter-kernel::
+    :id: discretize
+
+.. jupyter-execute::
+    :hide-code:
+
+    # Tutorial 2.9. Processing continuum Hamiltonians with discretize
+    # ===============================================================
+    #
+    # Physics background
+    # ------------------
+    #  - tight-binding approximation of continuous Hamiltonians
+    #
+    # Kwant features highlighted
+    # --------------------------
+    #  - kwant.continuum.discretize
+
+    import matplotlib as mpl
+    from matplotlib import pyplot
+
+    import kwant
+
+.. jupyter-execute:: boilerplate.py
+    :hide-code:
 
 .. _tutorial_discretizer_introduction:
 
@@ -55,9 +79,16 @@ Using `~kwant.continuum.discretize` to obtain a template
 ........................................................
 First we must explicitly import the `kwant.continuum` package:
 
-.. literalinclude:: /code/include/discretize.py
-    :start-after: #HIDDEN_BEGIN_import
-    :end-before: #HIDDEN_END_import
+.. jupyter-execute::
+
+    import kwant.continuum
+
+.. jupyter-execute::
+    :hide-code:
+
+    import scipy.sparse.linalg
+    import scipy.linalg
+    import numpy as np
 
 The function `kwant.continuum.discretize` takes a symbolic Hamiltonian and
 turns it into a `~kwant.builder.Builder` instance with appropriate spatial
@@ -65,17 +96,16 @@ symmetry that serves as a template.
 (We will see how to use the template to build systems with a particular
 shape later).
 
-.. literalinclude:: /code/include/discretize.py
-    :start-after: #HIDDEN_BEGIN_symbolic_discretization
-    :end-before: #HIDDEN_END_symbolic_discretization
+.. jupyter-execute::
+
+    template = kwant.continuum.discretize('k_x * A(x) * k_x')
+    print(template)
 
 It is worth noting that ``discretize`` treats ``k_x`` and ``x`` as
 non-commuting operators, and so their order is preserved during the
 discretization process.
 
-The builder produced by ``discretize`` may be printed to show the source code of its onsite and hopping functions (this is a special feature of builders returned by ``discretize``):
-
-.. literalinclude:: /code/figure/discretizer_intro_verbose.txt
+Printing the Builder produced by ``discretize`` shows the source code of its onsite and hopping functions (this is a special feature of builders returned by ``discretize``).
 
 .. specialnote:: Technical details
 
@@ -139,26 +169,45 @@ where :math:`V(x, y)` is some arbitrary potential.
 First, use ``discretize`` to obtain a
 builder that we will use as a template:
 
-.. literalinclude:: /code/include/discretize.py
-    :start-after: #HIDDEN_BEGIN_template
-    :end-before: #HIDDEN_END_template
+.. jupyter-execute::
+
+    hamiltonian = "k_x**2 + k_y**2 + V(x, y)"
+    template = kwant.continuum.discretize(hamiltonian)
+    print(template)
 
 We now use this system with the `~kwant.builder.Builder.fill`
 method of `~kwant.builder.Builder` to construct the system we
 want to investigate:
 
-.. literalinclude:: /code/include/discretize.py
-    :start-after: #HIDDEN_BEGIN_fill
-    :end-before: #HIDDEN_END_fill
+.. jupyter-execute::
+
+    def stadium(site):
+        (x, y) = site.pos
+        x = max(abs(x) - 20, 0)
+        return x**2 + y**2 < 30**2
+
+    syst = kwant.Builder()
+    syst.fill(template, stadium, (0, 0));
+    syst = syst.finalized()
 
 After finalizing this system, we can plot one of the system's
 energy eigenstates:
 
-.. literalinclude:: /code/include/discretize.py
-    :start-after: #HIDDEN_BEGIN_plot_eigenstate
-    :end-before: #HIDDEN_END_plot_eigenstate
+.. jupyter-execute::
+
+    def plot_eigenstate(syst, n=2, Vx=.0003, Vy=.0005):
+
+        def potential(x, y):
+            return Vx * x + Vy * y
 
-.. image:: /code/figure/discretizer_gs.*
+        ham = syst.hamiltonian_submatrix(params=dict(V=potential), sparse=True)
+        evecs = scipy.sparse.linalg.eigsh(ham, k=10, which='SM')[1]
+        kwant.plotter.map(syst, abs(evecs[:, n])**2, show=False)
+
+.. jupyter-execute::
+    :hide-code:
+
+    plot_eigenstate(syst)
 
 Note in the above that we pass the spatially varying potential *function*
 to our system via a parameter called ``V``, because the symbol :math:`V`
@@ -176,32 +225,95 @@ model [1]_ [2]_, one can provide matrix input to `~kwant.continuum.discretize`
 using ``identity`` and ``kron``. For example, the definition of the BHZ model can be
 written succinctly as:
 
-.. literalinclude:: /code/include/discretize.py
-    :start-after: #HIDDEN_BEGIN_define_qsh
-    :end-before: #HIDDEN_END_define_qsh
+.. jupyter-execute::
+
+    hamiltonian = """
+       + C * identity(4) + M * kron(sigma_0, sigma_z)
+       - B * (k_x**2 + k_y**2) * kron(sigma_0, sigma_z)
+       - D * (k_x**2 + k_y**2) * kron(sigma_0, sigma_0)
+       + A * k_x * kron(sigma_z, sigma_x)
+       - A * k_y * kron(sigma_0, sigma_y)
+    """
+
+    a = 20
+
+    template = kwant.continuum.discretize(hamiltonian, grid=a)
 
 We can then make a ribbon out of this template system:
 
-.. literalinclude:: /code/include/discretize.py
-    :start-after: #HIDDEN_BEGIN_define_qsh_build
-    :end-before: #HIDDEN_END_define_qsh_build
+.. jupyter-execute::
+
+    L, W = 2000, 1000
+
+    def shape(site):
+        (x, y) = site.pos
+        return (0 <= y < W and 0 <= x < L)
+
+    def lead_shape(site):
+        (x, y) = site.pos
+        return (0 <= y < W)
+
+    syst = kwant.Builder()
+    syst.fill(template, shape, (0, 0))
+
+    lead = kwant.Builder(kwant.TranslationalSymmetry([-a, 0]))
+    lead.fill(template, lead_shape, (0, 0))
+
+    syst.attach_lead(lead)
+    syst.attach_lead(lead.reversed())
+
+    syst = syst.finalized()
 
 and plot its dispersion using `kwant.plotter.bands`:
 
-.. literalinclude:: /code/include/discretize.py
-    :start-after: #HIDDEN_BEGIN_plot_qsh_band
-    :end-before: #HIDDEN_END_plot_qsh_band
+.. jupyter-execute::
 
-.. image:: /code/figure/discretizer_qsh_band.*
+    params = dict(A=3.65, B=-68.6, D=-51.1, M=-0.01, C=0)
+
+    kwant.plotter.bands(syst.leads[0], params=params,
+                        momenta=np.linspace(-0.3, 0.3, 201), show=False)
+
+    pyplot.grid()
+    pyplot.xlim(-.3, 0.3)
+    pyplot.ylim(-0.05, 0.05)
+    pyplot.xlabel('momentum [1/A]')
+    pyplot.ylabel('energy [eV]')
+    pyplot.show()
 
 In the above we see the edge states of the quantum spin Hall effect, which
 we can visualize using `kwant.plotter.map`:
 
-.. literalinclude:: /code/include/discretize.py
-    :start-after: #HIDDEN_BEGIN_plot_qsh_wf
-    :end-before: #HIDDEN_END_plot_qsh_wf
+.. jupyter-execute::
+
+    # get scattering wave functions at E=0
+    wf = kwant.wave_function(syst, energy=0, params=params)
 
-.. image:: /code/figure/discretizer_qsh_wf.*
+    # prepare density operators
+    sigma_z = np.array([[1, 0], [0, -1]])
+    prob_density = kwant.operator.Density(syst, np.eye(4))
+    spin_density = kwant.operator.Density(syst, np.kron(sigma_z, np.eye(2)))
+
+    # calculate expectation values and plot them
+    wf_sqr = sum(prob_density(psi) for psi in wf(0))  # states from left lead
+    rho_sz = sum(spin_density(psi) for psi in wf(0))  # states from left lead
+
+    fig, (ax1, ax2) = pyplot.subplots(1, 2, sharey=True, figsize=(16, 4))
+    kwant.plotter.map(syst, wf_sqr, ax=ax1)
+    kwant.plotter.map(syst, rho_sz, ax=ax2)
+
+    ax = ax1
+    im = [obj for obj in ax.get_children()
+          if isinstance(obj, mpl.image.AxesImage)][0]
+    fig.colorbar(im, ax=ax)
+
+    ax = ax2
+    im = [obj for obj in ax.get_children()
+          if isinstance(obj, mpl.image.AxesImage)][0]
+    fig.colorbar(im, ax=ax)
+
+    ax1.set_title('Probability density')
+    ax2.set_title('Spin density')
+    pyplot.show()
 
 
 Limitations of discretization
@@ -238,29 +350,62 @@ example. Let us start from the continuum Hamiltonian
 We start by defining this model as a string and setting the value of the
 :math:`α` parameter:
 
-.. literalinclude:: /code/include/discretize.py
-    :start-after: #HIDDEN_BEGIN_ls_def
-    :end-before: #HIDDEN_END_ls_def
+.. jupyter-execute::
+
+    hamiltonian = "k_x**2 * identity(2) + alpha * k_x * sigma_y"
+    params = dict(alpha=.5)
 
 Now we can use `kwant.continuum.lambdify` to obtain a function that computes
 :math:`H(k)`:
 
-.. literalinclude:: /code/include/discretize.py
-    :start-after: #HIDDEN_BEGIN_ls_hk_cont
-    :end-before: #HIDDEN_END_ls_hk_cont
+.. jupyter-execute::
+
+    h_k = kwant.continuum.lambdify(hamiltonian, locals=params)
+    k_cont = np.linspace(-4, 4, 201)
+    e_cont = [scipy.linalg.eigvalsh(h_k(k_x=ki)) for ki in k_cont]
 
 We can also construct a discretized approximation using
 `kwant.continuum.discretize`, in a similar manner to previous examples:
 
-.. literalinclude:: /code/include/discretize.py
-    :start-after: #HIDDEN_BEGIN_ls_hk_tb
-    :end-before: #HIDDEN_END_ls_hk_tb
+.. jupyter-execute::
+
+    def plot(ax, a=1):
+        template = kwant.continuum.discretize(hamiltonian, grid=a)
+        syst = kwant.wraparound.wraparound(template).finalized()
+
+        def h_k(k_x):
+            p = dict(k_x=k_x, **params)
+            return syst.hamiltonian_submatrix(params=p)
+
+        k_tb = np.linspace(-np.pi/a, np.pi/a, 201)
+        e_tb = [scipy.linalg.eigvalsh(h_k(k_x=a*ki)) for ki in k_tb]
+
+        ax.plot(k_cont, e_cont, 'r-')
+        ax.plot(k_tb, e_tb, 'k-')
+
+        ax.plot([], [], 'r-', label='continuum')
+        ax.plot([], [], 'k-', label='tight-binding')
+
+        ax.set_xlim(-4, 4)
+        ax.set_ylim(-1, 14)
+        ax.set_title('a={}'.format(a))
+
+        ax.set_xlabel('momentum [a.u.]')
+        ax.set_ylabel('energy [a.u.]')
+        ax.grid()
+        ax.legend()
 
 Below we can see the continuum and tight-binding dispersions for two
 different values of the discretization grid spacing :math:`a`:
 
-.. image:: /code/figure/discretizer_lattice_spacing.*
+.. jupyter-execute::
+    :hide-code:
+
+    _, (ax1, ax2) = pyplot.subplots(1, 2, sharey=True, figsize=(12, 4))
 
+    plot(ax1, a=1)
+    plot(ax2, a=.25)
+    pyplot.show()
 
 We clearly see that the smaller grid spacing is, the better we approximate
 the original continuous dispersion. It is also worth remembering that the
@@ -284,20 +429,26 @@ It is possible to use ``identity`` (for identity matrix), ``kron`` (for Kronecke
 expressions involving matrices. Matrices can also be provided explicitly using
 square ``[]`` brackets. For example, all following expressions are equivalent:
 
-.. literalinclude:: /code/include/discretize.py
-    :start-after: #HIDDEN_BEGIN_subs_1
-    :end-before: #HIDDEN_END_subs_1
+.. jupyter-execute::
+
+    sympify = kwant.continuum.sympify
+    subs = {'sx': [[0, 1], [1, 0]], 'sz': [[1, 0], [0, -1]]}
+
+    e = (
+        sympify('[[k_x**2, alpha * k_x], [k_x * alpha, -k_x**2]]'),
+        sympify('k_x**2 * sigma_z + alpha * k_x * sigma_x'),
+        sympify('k_x**2 * sz + alpha * k_x * sx', locals=subs),
+    )
 
-.. literalinclude:: /code/figure/discretizer_subs_1.txt
+    print(e[0] == e[1] == e[2])
 
 We can use the ``locals`` keyword parameter to substitute expressions
 and numerical values:
 
-.. literalinclude:: /code/include/discretize.py
-    :start-after: #HIDDEN_BEGIN_subs_2
-    :end-before: #HIDDEN_END_subs_2
+.. jupyter-execute::
 
-.. literalinclude:: /code/figure/discretizer_subs_2.txt
+    subs = {'A': 'A(x) + B', 'V': 'V(x) + V_0', 'C': 5}
+    print(sympify('k_x * A * k_x + V + C', locals=subs))
 
 Symbolic expressions obtained in this way can be directly passed to all
 ``discretizer`` functions.

doc/source/tutorial/graphene.rst

@@ -5,7 +5,7 @@ Beyond square lattices: graphene
 
 .. seealso::
     The complete source code of this example can be found in
-    :download:`graphene.py </code/download/graphene.py>`
+    :jupyter-download:script:`graphene`
 
 In the following example, we are going to calculate the
 conductance through a graphene quantum dot with a p-n junction
@@ -18,9 +18,43 @@ We begin by defining the honeycomb lattice of graphene. This is
 in principle already done in `kwant.lattice.honeycomb`, but we do it
 explicitly here to show how to define a new lattice:
 
-.. literalinclude:: /code/include/graphene.py
-    :start-after: #HIDDEN_BEGIN_hnla
-    :end-before: #HIDDEN_END_hnla
+.. jupyter-kernel::
+    :id: graphene
+
+.. jupyter-execute::
+    :hide-code:
+
+    # Tutorial 2.5. Beyond square lattices: graphene
+    # ==============================================
+    #
+    # Physics background
+    # ------------------
+    #  Transport through a graphene quantum dot with a pn-junction
+    #
+    # Kwant features highlighted
+    # --------------------------
+    #  - Application of all the aspects of tutorials 1-3 to a more complicated
+    #    lattice, namely graphene
+
+    from math import pi, sqrt, tanh
+
+    from matplotlib import pyplot
+
+    import kwant
+
+    # For computing eigenvalues
+    import scipy.sparse.linalg as sla
+
+    sin_30, cos_30 = (1 / 2, sqrt(3) / 2)
+
+.. jupyter-execute:: boilerplate.py
+    :hide-code:
+
+.. jupyter-execute::
+
+    graphene = kwant.lattice.general([(1, 0), (sin_30, cos_30)],
+                                     [(0, 0), (0, 1 / sqrt(3))])
+    a, b = graphene.sublattices
 
 The first argument to the `~kwant.lattice.general` function is the list of
 primitive vectors of the lattice; the second one is the coordinates of basis
@@ -31,9 +65,30 @@ itself forms a regular lattice of the same type as well, and those
 In the next step we define the shape of the scattering region (circle again)
 and add all lattice points using the ``shape``-functionality:
 
-.. literalinclude:: /code/include/graphene.py
-    :start-after: #HIDDEN_BEGIN_shzy
-    :end-before: #HIDDEN_END_shzy
+.. jupyter-execute::
+    :hide-code:
+
+    r = 10
+    w = 2.0
+    pot = 0.1
+
+.. jupyter-execute::
+
+    #### Define the scattering region. ####
+    # circular scattering region
+    def circle(pos):
+        x, y = pos
+        return x ** 2 + y ** 2 < r ** 2
+
+    syst = kwant.Builder()
+
+    # w: width and pot: potential maximum of the p-n junction
+    def potential(site):
+        (x, y) = site.pos
+        d = y * cos_30 + x * sin_30
+        return pot * tanh(d / w)
+
+    syst[graphene.shape(circle, (0, 0))] = potential
 
 As you can see, this works exactly the same for any kind of lattice.
 We add the onsite energies using a function describing the p-n junction;
@@ -45,9 +100,11 @@ As a next step we add the hoppings, making use of
 `~kwant.builder.HoppingKind`. For illustration purposes we define
 the hoppings ourselves instead of using ``graphene.neighbors()``:
 
-.. literalinclude:: /code/include/graphene.py
-    :start-after: #HIDDEN_BEGIN_hsmc
-    :end-before: #HIDDEN_END_hsmc
+.. jupyter-execute::
+
+    # specify the hoppings of the graphene lattice in the
+    # format expected by builder.HoppingKind
+    hoppings = (((0, 0), a, b), ((0, 1), a, b), ((-1, 1), a, b))
 
 The nearest-neighbor model for graphene contains only
 hoppings between different basis atoms. For this type of
@@ -63,27 +120,50 @@ respect to the two primitive vectors ``[(1, 0), (sin_30, cos_30)]``.
 
 Adding the hoppings however still works the same way:
 
-.. literalinclude:: /code/include/graphene.py
-    :start-after: #HIDDEN_BEGIN_bfwb
-    :end-before: #HIDDEN_END_bfwb
+.. jupyter-execute::
+
+    syst[[kwant.builder.HoppingKind(*hopping) for hopping in hoppings]] = -1
 
 Modifying the scattering region is also possible as before. Let's
 do something crazy, and remove an atom in sublattice A
 (which removes also the hoppings from/to this site) as well
 as add an additional link:
 
-.. literalinclude:: /code/include/graphene.py
-    :start-after: #HIDDEN_BEGIN_efut
-    :end-before: #HIDDEN_END_efut
+.. jupyter-execute::
+
+    # Modify the scattering region
+    del syst[a(0, 0)]
+    syst[a(-2, 1), b(2, 2)] = -1
 
 Note again that the conversion from a tuple `(i,j)` to site
 is done by the sublattices `a` and `b`.
 
 The leads are defined almost as before:
 
-.. literalinclude:: /code/include/graphene.py
-    :start-after: #HIDDEN_BEGIN_aakh
-    :end-before: #HIDDEN_END_aakh
+.. jupyter-execute::
+
+    #### Define the leads. ####
+    # left lead
+    sym0 = kwant.TranslationalSymmetry(graphene.vec((-1, 0)))
+
+    def lead0_shape(pos):
+        x, y = pos
+        return (-0.4 * r < y < 0.4 * r)
+
+    lead0 = kwant.Builder(sym0)
+    lead0[graphene.shape(lead0_shape, (0, 0))] = -pot
+    lead0[[kwant.builder.HoppingKind(*hopping) for hopping in hoppings]] = -1
+
+    # The second lead, going to the top right
+    sym1 = kwant.TranslationalSymmetry(graphene.vec((0, 1)))
+
+    def lead1_shape(pos):
+        v = pos[1] * sin_30 - pos[0] * cos_30
+        return (-0.4 * r < v < 0.4 * r)
+
+    lead1 = kwant.Builder(sym1)
+    lead1[graphene.shape(lead1_shape, (0, 0))] = pot
+    lead1[[kwant.builder.HoppingKind(*hopping) for hopping in hoppings]] = -1
 
 Note the method `~kwant.lattice.Polyatomic.vec` used in calculating the
 parameter for `~kwant.lattice.TranslationalSymmetry`.  The latter expects a
@@ -98,28 +178,63 @@ in a square lattice -- they follow the non-orthogonal primitive vectors defined
 in the beginning.
 
 Later, we will compute some eigenvalues of the closed scattering region without
-leads. This is why we postpone attaching the leads to the system. Instead,
-we return the scattering region and the leads separately.
+leads. This is why we postpone attaching the leads to the system.
 
-.. literalinclude:: /code/include/graphene.py
-    :start-after: #HIDDEN_BEGIN_kmmw
-    :end-before: #HIDDEN_END_kmmw
 
 The computation of some eigenvalues of the closed system is done
 in the following piece of code:
 
-.. literalinclude:: /code/include/graphene.py
-    :start-after: #HIDDEN_BEGIN_zydk
-    :end-before: #HIDDEN_END_zydk
+.. jupyter-execute::
+
+    def compute_evs(syst):
+        # Compute some eigenvalues of the closed system
+        sparse_mat = syst.hamiltonian_submatrix(sparse=True)
+
+        evs = sla.eigs(sparse_mat, 2)[0]
+        print(evs.real)
 
 The code for computing the band structure and the conductance is identical
 to the previous examples, and needs not be further explained here.
 
-Finally, in the `main` function we make and plot the system:
+Finally, we plot the system:
+
+.. jupyter-execute::
+    :hide-code:
+
+    def plot_conductance(syst, energies):
+        # Compute transmission as a function of energy
+        data = []
+        for energy in energies:
+            smatrix = kwant.smatrix(syst, energy)
+            data.append(smatrix.transmission(0, 1))
+
+        pyplot.figure()
+        pyplot.plot(energies, data)
+        pyplot.xlabel("energy [t]")
+        pyplot.ylabel("conductance [e^2/h]")
+        pyplot.show()
+
+
+    def plot_bandstructure(flead, momenta):
+        bands = kwant.physics.Bands(flead)
+        energies = [bands(k) for k in momenta]
+
+        pyplot.figure()
+        pyplot.plot(momenta, energies)
+        pyplot.xlabel("momentum [(lattice constant)^-1]")
+        pyplot.ylabel("energy [t]")
+        pyplot.show()
+
+
+.. jupyter-execute::
+
+    # To highlight the two sublattices of graphene, we plot one with
+    # a filled, and the other one with an open circle:
+    def family_colors(site):
+        return 0 if site.family == a else 1
 
-.. literalinclude:: /code/include/graphene.py
-    :start-after: #HIDDEN_BEGIN_itkk
-    :end-before: #HIDDEN_END_itkk
+    # Plot the closed system without leads.
+    kwant.plot(syst, site_color=family_colors, site_lw=0.1, colorbar=False);
 
 We customize the plotting: we set the `site_colors` argument of
 `~kwant.plotter.plot` to a function which returns 0 for
@@ -132,27 +247,43 @@ The function `~kwant.plotter.plot` shows these values using a color scale
 (grayscale by default). The symbol `size` is specified in points, and is
 independent on the overall figure size.
 
-Plotting the closed system gives this result:
-
-.. image:: /code/figure/graphene_syst1.*
 
 Computing the eigenvalues of largest magnitude,
 
-.. literalinclude:: /code/include/graphene.py
-    :start-after: #HIDDEN_BEGIN_jmbi
-    :end-before: #HIDDEN_END_jmbi
+.. jupyter-execute::
 
-should yield two eigenvalues equal to ``[ 3.07869311,
+    compute_evs(syst.finalized())
+
+yields two eigenvalues equal to ``[ 3.07869311,
 -3.06233144]``.
 
-The remaining code of `main` attaches the leads to the system and plots it
+The remaining code attaches the leads to the system and plots it
 again:
 
-.. image:: /code/figure/graphene_syst2.*
+.. jupyter-execute::
+
+    # Attach the leads to the system.
+    for lead in [lead0, lead1]:
+        syst.attach_lead(lead)
+
+    # Then, plot the system with leads.
+    kwant.plot(syst, site_color=family_colors, site_lw=0.1,
+               lead_site_lw=0, colorbar=False);
 
-It computes the band structure of one of lead 0:
+We then finalize the system:
 
-.. image:: /code/figure/graphene_bs.*
+.. jupyter-execute::
+
+    syst = syst.finalized()
+
+and compute the band structure of one of lead 0:
+
+.. jupyter-execute::
+
+
+    # Compute the band structure of lead 0.
+    momenta = [-pi + 0.02 * pi * i for i in range(101)]
+    plot_bandstructure(syst.leads[0], momenta)
 
 showing all the features of a zigzag lead, including the flat
 edge state bands (note that the band structure is not symmetric around
@@ -160,7 +291,11 @@ zero energy, due to a potential in the leads).
 
 Finally the transmission through the system is computed,
 
-.. image:: /code/figure/graphene_result.*
+.. jupyter-execute::
+
+    # Plot conductance.
+    energies = [-2 * pot + 4. / 50. * pot * i for i in range(51)]
+    plot_conductance(syst, energies)
 
 showing the typical resonance-like transmission probability through
 an open quantum dot