doc/source/tutorial/README

-Note for Kwant developers
--------------------------
-
-Say you have modified SCRIPT.py in this directory.  Make sure that the tutorial
-image generation patches still apply by deleting doc/source/images/SCRIPT.py
-and running ``make doc/source/images/SCRIPT.py`` in doc.  Now examine the newly
-created doc/source/images/SCRIPT.py.  If you do not like the result or the
-patch did not apply, edit doc/source/images/SCRIPT.py until you like it.  You
-can run `make html` during your edits to check things.  Finally, even if you
-did not edit the script, execute generate-diffs.sh in doc/source/images.  If
-the patches applied cleanly the diff files will usually stay unchanged.

doc/source/tutorial/ab_ring.py

-# Tutorial 2.3.3. Nontrivial shapes
-# =================================
-#
-# Physics background
-# ------------------
-#  Flux-dependent transmission through a quantum ring
-#
-# Kwant features highlighted
-# --------------------------
-#  - More complex shapes with lattices
-#  - Allows for discussion of subtleties of `attach_lead` (not in the
-#    example, but in the tutorial main text)
-#  - Modifcations of hoppings/sites after they have been added
-
-from cmath import exp
-from math import pi
-
-import kwant
-
-# For plotting
-from matplotlib import pyplot
-
-
-#HIDDEN_BEGIN_eusz
-def make_system(a=1, t=1.0, W=10, r1=10, r2=20):
-    # 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()
-
-    #### Define the scattering region. ####
-    # Now, we aim for a more complex shape, namely a ring (or annulus)
-    def ring(pos):
-        (x, y) = pos
-        rsq = x ** 2 + y ** 2
-        return (r1 ** 2 < rsq < r2 ** 2)
-#HIDDEN_END_eusz
-
-    # and add the corresponding lattice points using the `shape`-function
-#HIDDEN_BEGIN_lcak
-    syst[lat.shape(ring, (0, r1 + 1))] = 4 * t
-    syst[lat.neighbors()] = -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 the Builder instance repeatedly, we
-    # define the modified hoppings as a function that takes the flux as its
-    # parameter phi.
-#HIDDEN_BEGIN_lvkt
-    def hopping_phase(site1, site2, phi):
-        return -t * exp(1j * phi)
-
-    def crosses_branchcut(hop):
-        ix0, iy0 = hop[0].tag
-
-        # builder.HoppingKind with the argument (1, 0) below
-        # returns hoppings ordered as ((i+1, j), (i, j))
-        return iy0 < 0 and ix0 == 1  # ix1 == 0 then implied
-
-    # Modify only those hopings in x-direction that cross the branch cut
-    def hops_across_cut(syst):
-        for hop in kwant.builder.HoppingKind((1, 0), lat, lat)(syst):
-            if crosses_branchcut(hop):
-                yield hop
-    syst[hops_across_cut] = hopping_phase
-#HIDDEN_END_lvkt
-
-    #### Define the leads. ####
-    # left lead
-#HIDDEN_BEGIN_qwgr
-    sym_lead = kwant.TranslationalSymmetry((-a, 0))
-    lead = kwant.Builder(sym_lead)
-
-    def lead_shape(pos):
-        (x, y) = pos
-        return (-W / 2 < y < W / 2)
-
-    lead[lat.shape(lead_shape, (0, 0))] = 4 * t
-    lead[lat.neighbors()] = -t
-#HIDDEN_END_qwgr
-
-    #### Attach the leads and return the system. ####
-#HIDDEN_BEGIN_skbz
-    syst.attach_lead(lead)
-    syst.attach_lead(lead.reversed())
-#HIDDEN_END_skbz
-
-    return syst
-
-
-def plot_conductance(syst, energy, fluxes):
-    # compute conductance
-
-    normalized_fluxes = [flux / (2 * pi) for flux in fluxes]
-    data = []
-    for flux in fluxes:
-        smatrix = kwant.smatrix(syst, energy, args=[flux])
-        data.append(smatrix.transmission(1, 0))
-
-    pyplot.figure()
-    pyplot.plot(normalized_fluxes, data)
-    pyplot.xlabel("flux [flux quantum]")
-    pyplot.ylabel("conductance [e^2/h]")
-    pyplot.show()
-
-
-def main():
-    syst = make_system()
-
-    # Check that the system looks as intended.
-    kwant.plot(syst)
-
-    # Finalize the system.
-    syst = syst.finalized()
-
-    # We should see a conductance that is periodic with the flux quantum
-    plot_conductance(syst, energy=0.15, fluxes=[0.01 * i * 3 * 2 * pi
-                                                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/tutorial/band_structure.py

-# Tutorial 2.4.1. Band structure calculations
-# ===========================================
-#
-# Physics background
-# ------------------
-#  band structure of a simple quantum wire in tight-binding approximation
-#
-# Kwant features highlighted
-# --------------------------
-#  - Computing the band structure of a finalized lead.
-
-import kwant
-
-# For plotting.
-from matplotlib import pyplot
-
-#HIDDEN_BEGIN_zxip
-def make_lead(a=1, t=1.0, W=10):
-    # Start with an empty lead with a single square lattice
-    lat = kwant.lattice.square(a)
-
-    sym_lead = kwant.TranslationalSymmetry((-a, 0))
-    lead = kwant.Builder(sym_lead)
-
-    # build up one unit cell of the lead, and add the hoppings
-    # to the next unit cell
-    for j in range(W):
-        lead[lat(0, j)] = 4 * t
-
-        if j > 0:
-            lead[lat(0, j), lat(0, j - 1)] = -t
-
-        lead[lat(1, j), lat(0, j)] = -t
-
-    return lead
-#HIDDEN_END_zxip
-
-
-#HIDDEN_BEGIN_pejz
-def main():
-    lead = make_lead().finalized()
-    kwant.plotter.bands(lead, show=False)
-    pyplot.xlabel("momentum [(lattice constant)^-1]")
-    pyplot.ylabel("energy [t]")
-    pyplot.show()
-#HIDDEN_END_pejz
-
-
-# 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/boilerplate.py

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

doc/source/tutorial/closed_system.py

-# Tutorial 2.4.2. Closed systems
-# ==============================
-#
-# Physics background
-# ------------------
-#  Fock-darwin spectrum of a quantum dot (energy spectrum in
-#  as a function of a magnetic field)
-#
-# Kwant features highlighted
-# --------------------------
-#  - Use of `hamiltonian_submatrix` in order to obtain a Hamiltonian
-#    matrix.
-
-from cmath import exp
-import numpy as np
-import kwant
-
-# For eigenvalue computation
-#HIDDEN_BEGIN_tibv
-import scipy.sparse.linalg as sla
-#HIDDEN_END_tibv
-
-# For plotting
-from matplotlib import pyplot
-
-
-def make_system(a=1, t=1.0, r=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).
-
-#HIDDEN_BEGIN_qlyd
-    lat = kwant.lattice.square(a)
-
-    syst = kwant.Builder()
-
-    # Define the quantum dot
-    def circle(pos):
-        (x, y) = pos
-        rsq = x ** 2 + y ** 2
-        return rsq < r ** 2
-
-    def hopx(site1, site2, B=0):
-        # The magnetic field is controlled by the parameter B
-        y = site1.pos[1]
-        return -t * exp(-1j * B * y)
-
-    syst[lat.shape(circle, (0, 0))] = 4 * t
-    # hoppings in x-direction
-    syst[kwant.builder.HoppingKind((1, 0), lat, lat)] = hopx
-    # hoppings in y-directions
-    syst[kwant.builder.HoppingKind((0, 1), lat, lat)] = -t
-
-    # It's a closed system for a change, so no leads
-    return syst
-#HIDDEN_END_qlyd
-
-
-#HIDDEN_BEGIN_yvri
-def plot_spectrum(syst, Bfields):
-
-    # In the following, we compute the spectrum of the quantum dot
-    # using dense matrix methods. This works in this toy example, as
-    # the system is tiny. In a real example, one would want to use
-    # sparse matrix methods
-
-    energies = []
-    for B in Bfields:
-        # Obtain the Hamiltonian as a dense matrix
-        ham_mat = syst.hamiltonian_submatrix(args=[B], sparse=True)
-
-        # we only calculate the 15 lowest eigenvalues
-        ev = sla.eigsh(ham_mat, k=15, which='SM', return_eigenvectors=False)
-
-        energies.append(ev)
-
-    pyplot.figure()
-    pyplot.plot(Bfields, energies)
-    pyplot.xlabel("magnetic field [arbitrary units]")
-    pyplot.ylabel("energy [t]")
-    pyplot.show()
-#HIDDEN_END_yvri
-
-
-#HIDDEN_BEGIN_wave
-def plot_wave_function(syst):
-    # Calculate the wave functions in the system.
-    ham_mat = syst.hamiltonian_submatrix(sparse=True)
-    evecs = sla.eigsh(ham_mat, k=20, which='SM')[1]
-
-    # Plot the probability density of the 10th eigenmode.
-    kwant.plotter.map(syst, np.abs(evecs[:, 9])**2,
-                      colorbar=False, oversampling=1)
-#HIDDEN_END_wave
-
-
-def main():
-    syst = make_system()
-
-    # Check that the system looks as intended.
-    kwant.plot(syst)
-
-    # Finalize the system.
-    syst = syst.finalized()
-
-    # The following try-clause can be removed once SciPy 0.9 becomes uncommon.
-    try:
-        # We should observe energy levels that flow towards Landau
-        # level energies with increasing magnetic field.
-        plot_spectrum(syst, [iB * 0.002 for iB in range(100)])
-
-        # Plot an eigenmode of a circular dot. Here we create a larger system for
-        # better spatial resolution.
-        syst = make_system(r=30).finalized()
-        plot_wave_function(syst)
-    except ValueError as e:
-        if e.message == "Input matrix is not real-valued.":
-            print("The calculation of eigenvalues failed because of a bug in SciPy 0.9.")
-            print("Please upgrade to a newer version of SciPy.")
-        else:
-            raise
-
-
-# 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/discretize.rst

+Discretizing continuous Hamiltonians
+------------------------------------
+
+Introduction
+............
+
+In ":ref:`tutorial_discretization_schrodinger`" we have learnt that Kwant works
+with tight-binding Hamiltonians. Often, however, one will start with a
+continuum model and will subsequently need to discretize it to arrive at a
+tight-binding model.
+Although discretizing a Hamiltonian is usually a simple
+process, it is tedious and repetitive. The situation is further exacerbated
+when one introduces additional on-site degrees of freedom, and tracking all
+the necessary terms becomes a chore.
+The `~kwant.continuum` sub-package aims to be a solution to this problem.
+It is a collection of tools for working with
+continuum models and for discretizing them into tight-binding models.
+
+.. seealso::
+    You can execute the code examples live in your browser by
+    activating thebelab:
+
+    .. thebe-button:: Activate Thebelab
+
+.. seealso::
+    The complete source code of this tutorial can be found in
+    :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:
+
+Discretizing by hand
+....................
+
+As an example, let us consider the following continuum Schrödinger equation
+for a semiconducting heterostructure (using the effective mass approximation):
+
+.. math::
+
+    \left( k_x \frac{\hbar^2}{2 m(x)} k_x \right) \psi(x) = E \, \psi(x).
+
+Replacing the momenta by their corresponding differential operators
+
+.. math::
+    k_\alpha = -i \partial_\alpha,
+
+for :math:`\alpha = x, y` or :math:`z`, and discretizing on a regular lattice of
+points with spacing :math:`a`, we obtain the tight-binding model
+
+.. math::
+
+    H = - \frac{1}{a^2} \sum_i A\left(x+\frac{a}{2}\right)
+            \big(\ket{i}\bra{i+1} + h.c.\big)
+        + \frac{1}{a^2} \sum_i
+            \left( A\left(x+\frac{a}{2}\right) + A\left(x-\frac{a}{2}\right)\right)
+            \ket{i} \bra{i},
+
+with :math:`A(x) = \frac{\hbar^2}{2 m(x)}`.
+
+Using `~kwant.continuum.discretize` to obtain a template
+........................................................
+First we must explicitly import the `kwant.continuum` package:
+
+.. 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
+symmetry that serves as a template.
+(We will see how to use the template to build systems with a particular
+shape later).
+
+.. 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.
+
+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``).
+
+.. admonition:: Technical details
+    :class: dropdown note
+
+    - ``kwant.continuum`` uses ``sympy`` internally to handle symbolic
+      expressions. Strings are converted using `kwant.continuum.sympify`,
+      which essentially applies some Kwant-specific rules (such as treating
+      ``k_x`` and ``x`` as non-commutative) before calling ``sympy.sympify``
+
+    - The builder returned by ``discretize`` will have an N-D
+      translational symmetry, where ``N`` is the number of dimensions that were
+      discretized. This is the case, even if there are expressions in the input
+      (e.g. ``V(x, y)``) which in principle *may not* have this symmetry.  When
+      using the returned builder directly, or when using it as a template to
+      construct systems with different/lower symmetry, it is important to
+      ensure that any functional parameters passed to the system respect the
+      symmetry of the system. Kwant provides no consistency check for this.
+
+    - The discretization process consists of taking input
+      :math:`H(k_x, k_y, k_z)`, multiplying it from the right by
+      :math:`\psi(x, y, z)` and iteratively applying a second-order accurate
+      central derivative approximation for every
+      :math:`k_\alpha=-i\partial_\alpha`:
+
+      .. math::
+         \partial_\alpha \psi(\alpha) =
+            \frac{1}{a} \left( \psi\left(\alpha + \frac{a}{2}\right)
+                              -\psi\left(\alpha - \frac{a}{2}\right)\right).
+
+      This process is done separately for every summand in Hamiltonian.
+      Once all symbols denoting operators are applied internal algorithm is
+      calculating ``gcd`` for hoppings coming from each summand in order to
+      find best possible approximation. Please see source code for details.
+
+    - Instead of using ``discretize`` one can use
+      `~kwant.continuum.discretize_symbolic` to obtain symbolic output.
+      When working interactively in `Jupyter notebooks <https://jupyter.org/>`_
+      it can be useful to use this to see a symbolic representation of
+      the discretized Hamiltonian. This works best when combined with ``sympy``
+      `Pretty Printing <https://docs.sympy.org/latest/tutorials/intro-tutorial/printing.html#setting-up-pretty-printing>`_.
+
+    - The symbolic result of discretization obtained with
+      ``discretize_symbolic`` can be converted into a
+      builder using `~kwant.continuum.build_discretized`.
+      This can be useful if one wants to alter the tight-binding Hamiltonian
+      before building the system.
+
+
+Building a Kwant system from the template
+.........................................
+
+Let us now use the output of ``discretize`` as a template to
+build a system and plot some of its energy eigenstate. For this example the
+Hamiltonian will be
+
+.. math::
+
+    H = k_x^2 + k_y^2 + V(x, y),
+
+where :math:`V(x, y)` is some arbitrary potential.
+
+First, use ``discretize`` to obtain a
+builder that we will use as a 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:
+
+.. 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:
+
+.. jupyter-execute::
+
+    def plot_eigenstate(syst, n=2, Vx=.0003, Vy=.0005):
+
+        def potential(x, y):
+            return Vx * x + Vy * y
+
+        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`
+was used in the initial, symbolic, definition of the Hamiltonian.
+
+In addition, the function passed as ``V`` expects two input parameters ``x``
+and ``y``, the same as in the initial continuum Hamiltonian.
+
+
+.. _discretize-bhz-model:
+
+Models with more structure: Bernevig-Hughes-Zhang
+.................................................
+
+When working with multi-band systems, like the Bernevig-Hughes-Zhang (BHZ)
+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:
+
+.. 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:
+
+.. 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`:
+
+.. jupyter-execute::
+
+    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`:
+
+.. jupyter-execute::
+
+    # get scattering wave functions at E=0
+    wf = kwant.wave_function(syst, energy=0, params=params)
+
+    # 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
+.............................
+
+It is important to remember that the discretization of a continuum
+model is an *approximation* that is only valid in the low-energy
+limit. For example, the quadratic continuum Hamiltonian
+
+.. math::
+
+    H_\textrm{continuous}(k_x) = \frac{\hbar^2}{2m}k_x^2
+
+
+and its discretized approximation
+
+.. math::
+
+    H_\textrm{tight-binding}(k_x) = 2t \big(1 - \cos(k_x a)\big),
+
+
+where :math:`t=\frac{\hbar^2}{2ma^2}`, are only valid in the limit
+:math:`E < t`. The grid spacing :math:`a` must be chosen according
+to how high in energy you need your tight-binding model to be valid.
+
+It is possible to set :math:`a` through the ``grid`` parameter
+to `~kwant.continuum.discretize`, as we will illustrate in the following
+example. Let us start from the continuum Hamiltonian
+
+.. math::
+
+  H(k) = k_x^2 \mathbb{1}_{2\times2} + α k_x \sigma_y.
+
+We start by defining this model as a string and setting the value of the
+:math:`α` parameter:
+
+.. 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)`:
+
+.. 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:
+
+.. 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`:
+
+.. 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
+Brillouin zone also scales with grid spacing: :math:`[-\frac{\pi}{a},
+\frac{\pi}{a}]`.
+
+
+Advanced topics
+...............
+
+The input to `kwant.continuum.discretize` and `kwant.continuum.lambdify` can be
+not only a ``string``, as we saw above, but also a ``sympy`` expression or
+a ``sympy`` matrix.
+This functionality will probably be mostly useful to people who
+are already experienced with ``sympy``.
+
+
+It is possible to use ``identity`` (for identity matrix), ``kron`` (for Kronecker product), as well as Pauli matrices ``sigma_0``,
+``sigma_x``, ``sigma_y``, ``sigma_z`` in the input to
+`~kwant.continuum.lambdify` and `~kwant.continuum.discretize`, in order to simplify
+expressions involving matrices. Matrices can also be provided explicitly using
+square ``[]`` brackets. For example, all following expressions are equivalent:
+
+.. 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),
+    )
+
+    print(e[0] == e[1] == e[2])
+
+We can use the ``locals`` keyword parameter to substitute expressions
+and numerical values:
+
+.. jupyter-execute::
+
+    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.
+
+.. admonition:: Technical details
+    :class: dropdown note
+
+    Because of the way that ``sympy`` handles commutation relations all symbols
+    representing position and momentum operators are set to be non commutative.
+    This means that the order of momentum and position operators in the input
+    expression is preserved.  Note that it is not possible to define individual
+    commutation relations within ``sympy``, even expressions such :math:`x k_y x`
+    will not be simplified, even though mathematically :math:`[x, k_y] = 0`.
+
+
+.. rubric:: References
+
+.. [1] `Science, 314, 1757 (2006)
+    <https://arxiv.org/abs/cond-mat/0611399>`_.
+
+.. [2] `Phys. Rev. B 82, 045122 (2010)
+    <https://arxiv.org/abs/1005.1682>`_.

doc/source/tutorial/index.rst

@@ -3,9 +3,14 @@ Tutorial: learning Kwant through examples
 
 .. toctree::
     introduction
-    tutorial1
-    tutorial2
-    tutorial3
-    tutorial4
-    tutorial5
-    tutorial6
+    first_steps
+    spin_potential_shape
+    spectrum
+    graphene
+    superconductors
+    operators
+    plotting
+    kpm
+    discretize
+    magnetic_field
+    faq

doc/source/tutorial/introduction.rst

@@ -7,8 +7,7 @@ commented extensively. In addition, you will find notes about more subtle,
 technical details at the end of each example. At first reading, these notes may
 be safely skipped.
 
-A scientific article about Kwant is available as well, see `Kwant
-website <http://kwant-project.org/citing.html>`_.
+A scientific article about Kwant is available as well, see :doc:`/pre/citing`.
 
 The article introduces Kwant with a somewhat different focus than the tutorial
 and it is the authors' intention that both texts complement each other.  While
@@ -26,14 +25,14 @@ transport in mesoscopic systems" by Supriyo Datta.
 The Python programming language
 ...............................
 
-Kwant is a library for `Python <http://python.org/>`_.  Care was taken to fit
-well with the spirit of the language and to take advantage of its expressive
-power.  If you do not know Python yet, do not fear: Python is widely regarded
-as one of the most accessible programming languages.  For an introduction we
-recommend the `official Python Tutorial <http://docs.python.org/2/tutorial/>`_.
-The `Beginner's Guide to Python <http://wiki.python.org/moin/BeginnersGuide>`_
-contains a wealth of links to other tutorials, guides and books including some
-for absolute beginners.
+Kwant is a library for `Python <https://www.python.org/>`_.  Care was taken to
+fit well with the spirit of the language and to take advantage of its
+expressive power.  If you do not know Python yet, do not fear: Python is widely
+regarded as one of the most accessible programming languages.  For an
+introduction we recommend the `official Python Tutorial
+<https://docs.python.org/3/tutorial/>`_.  The `Beginner's Guide to Python
+<https://wiki.python.org/moin/BeginnersGuide>`_ contains a wealth of links to
+other tutorials, guides and books including some for absolute beginners.
 
 Kwant
 .....

doc/source/tutorial/plot_zincblende.py

-# Tutorial 2.7.2. 3D example: zincblende structure
-# ================================================
-#
-# Physical background
-# -------------------
-#  - 3D Bravais lattices
-#
-# Kwant features highlighted
-# --------------------------
-#  - demonstrate different ways of plotting in 3D
-
-import kwant
-from matplotlib import pyplot
-
-#HIDDEN_BEGIN_zincblende1
-lat = kwant.lattice.general([(0, 0.5, 0.5), (0.5, 0, 0.5), (0.5, 0.5, 0)],
-                            [(0, 0, 0), (0.25, 0.25, 0.25)])
-a, b = lat.sublattices
-#HIDDEN_END_zincblende1
-
-#HIDDEN_BEGIN_zincblende2
-def make_cuboid(a=15, b=10, c=5):
-    def cuboid_shape(pos):
-        x, y, z = pos
-        return 0 <= x < a and 0 <= y < b and 0 <= z < c
-
-    syst = kwant.Builder()
-    syst[lat.shape(cuboid_shape, (0, 0, 0))] = None
-    syst[lat.neighbors()] = None
-
-    return syst
-#HIDDEN_END_zincblende2
-
-
-def main():
-    # the standard plotting style for 3D is mainly useful for
-    # checking shapes:
-#HIDDEN_BEGIN_plot1
-    syst = make_cuboid()
-
-    kwant.plot(syst)
-#HIDDEN_END_plot1
-
-    # visualize the crystal structure better for a very small system
-#HIDDEN_BEGIN_plot2
-    syst = make_cuboid(a=1.5, b=1.5, c=1.5)
-
-    def family_colors(site):
-        return 'r' if site.family == a else 'g'
-
-    kwant.plot(syst, site_size=0.18, site_lw=0.01, hop_lw=0.05,
-               site_color=family_colors)
-#HIDDEN_END_plot2
-
-
-# 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()