diff --git a/doc/source/tutorial/index.rst b/doc/source/tutorial/index.rst
index 1729a23126f5c0a7cfbd46aaa94e44a64881adf6..275e660e1c10c58d806422e917b77b17896dad5a 100644
--- a/doc/source/tutorial/index.rst
+++ b/doc/source/tutorial/index.rst
@@ -11,3 +11,4 @@ Tutorial: learning Kwant through examples
tutorial6
tutorial7
tutorial8
+ tutorial9
diff --git a/doc/source/tutorial/plot_graphene.py b/doc/source/tutorial/plot_graphene.py
index d9181734ba869a624d03617c565f0852c3c50792..e81e331b2bedf7de7241b857a50611a562d51316 100644
--- a/doc/source/tutorial/plot_graphene.py
+++ b/doc/source/tutorial/plot_graphene.py
@@ -1,4 +1,4 @@
-# Tutorial 2.7.1. 2D example: graphene quantum dot
+# Tutorial 2.8.1. 2D example: graphene quantum dot
# ================================================
#
# Physics background
diff --git a/doc/source/tutorial/plot_zincblende.py b/doc/source/tutorial/plot_zincblende.py
index 467f5767e2d5d2e0e4d08d168c6ba7b2de5f69f0..6b75a54e0df3d4b977ab1ec51a00960ef4f6f262 100644
--- a/doc/source/tutorial/plot_zincblende.py
+++ b/doc/source/tutorial/plot_zincblende.py
@@ -1,4 +1,4 @@
-# Tutorial 2.7.2. 3D example: zincblende structure
+# Tutorial 2.8.2. 3D example: zincblende structure
# ================================================
#
# Physical background
diff --git a/doc/source/tutorial/tutorial6.rst b/doc/source/tutorial/tutorial6.rst
deleted file mode 100644
index dbf267042eb65f06d614fe3a04d18594b5cd3ec2..0000000000000000000000000000000000000000
--- a/doc/source/tutorial/tutorial6.rst
+++ /dev/null
@@ -1,286 +0,0 @@
-Plotting Kwant systems and data in various styles
--------------------------------------------------
-
-The plotting functionality of Kwant has been used extensively (through
-`~kwant.plotter.plot` and `~kwant.plotter.map`) in the previous tutorials. In
-addition to this basic use, `~kwant.plotter.plot` offers many options to change
-the plotting style extensively. It is the goal of this tutorial to show how
-these options can be used to achieve various very different objectives.
-
-2D example: graphene quantum dot
-................................
-
-.. seealso::
- The complete source code of this example can be found in
- :download:`tutorial/plot_graphene.py <../../../tutorial/plot_graphene.py>`
-
-We begin by first considering a circular graphene quantum dot (similar to what
-has been used in parts of the tutorial :ref:`tutorial-graphene`.) In contrast
-to previous examples, we will also use hoppings beyond next-nearest neighbors:
-
-.. literalinclude:: plot_graphene.py
- :start-after: #HIDDEN_BEGIN_makesyst
- :end-before: #HIDDEN_END_makesyst
-
-Note that adding hoppings hoppings to the `n`-th nearest neighbors can be
-simply done by passing `n` as an argument to
-`~kwant.lattice.Polyatomic.neighbors`. Also note that we use the method
-`~kwant.builder.Builder.eradicate_dangling` to get rid of single atoms sticking
-out of the shape. It is necessary to do so *before* adding the
-next-nearest-neighbor hopping [#]_.
-
-Of course, the system can be plotted simply with default settings:
-
-.. literalinclude:: plot_graphene.py
- :start-after: #HIDDEN_BEGIN_plotsyst1
- :end-before: #HIDDEN_END_plotsyst1
-
-However, due to the richer structure of the lattice, this results in a rather
-busy plot:
-
-.. image:: ../images/plot_graphene_syst1.*
-
-A much clearer plot can be obtained by using different colors for both
-sublattices, and by having different line widths for different hoppings. This
-can be achieved by passing a function to the arguments of
-`~kwant.plotter.plot`, instead of a constant. For properties of sites, this
-must be a function taking one site as argument, for hoppings a function taking
-the start end end site of hopping as arguments:
-
-.. literalinclude:: plot_graphene.py
- :start-after: #HIDDEN_BEGIN_plotsyst2
- :end-before: #HIDDEN_END_plotsyst2
-
-Note that since we are using an unfinalized Builder, a ``site`` is really an
-instance of `~kwant.builder.Site`. With these adjustments we arrive at a plot
-that carries the same information, but is much easier to interpret:
-
-.. image:: ../images/plot_graphene_syst2.*
-
-Apart from plotting the *system* itself, `~kwant.plotter.plot` can also be used
-to plot *data* living on the system.
-
-As an example, we now compute the eigenstates of the graphene quantum dot and
-intend to plot the wave function probability in the quantum dot. For aesthetic
-reasons (the wave functions look a bit nicer), we restrict ourselves to
-nearest-neighbor hopping. Computing the wave functions is done in the usual
-way (note that for a large-scale system, one would probably want to use sparse
-linear algebra):
-
-.. literalinclude:: plot_graphene.py
- :start-after: #HIDDEN_BEGIN_plotdata1
- :end-before: #HIDDEN_END_plotdata1
-
-In most cases, to plot the wave function probability, one wouldn't use
-`~kwant.plotter.plot`, but rather `~kwant.plotter.map`. Here, we plot the
-`n`-th wave function using it:
-
-.. literalinclude:: plot_graphene.py
- :start-after: #HIDDEN_BEGIN_plotdata2
- :end-before: #HIDDEN_END_plotdata2
-
-This results in a standard pseudocolor plot, showing in this case (``n=225``) a
-graphene edge state, i.e. a wave function mostly localized at the zigzag edges
-of the quantum dot.
-
-.. image:: ../images/plot_graphene_data1.*
-
-However although in general preferable, `~kwant.plotter.map` has a few
-deficiencies for this small system: For example, there are a few distortions at
-the edge of the dot. (This cannot be avoided in the type of interpolation used
-in `~kwant.plotter.map`). However, we can also use `~kwant.plotter.plot` to
-achieve a similar, but smoother result.
-
-For this note that `~kwant.plotter.plot` can also take an array of floats (or
-function returning floats) as value for the ``site_color`` argument (the same
-holds for the hoppings). Via the colormap specified in ``cmap`` these are mapped
-to color, just as `~kwant.plotter.map` does! In addition, we can also change
-the symbol shape depending on the sublattice. With a triangle pointing up and
-down on the respective sublattice, the symbols used by plot fill the space
-completely:
-
-.. literalinclude:: plot_graphene.py
- :start-after: #HIDDEN_BEGIN_plotdata3
- :end-before: #HIDDEN_END_plotdata3
-
-Note that with ``hop_lw=0`` we deactivate plotting the hoppings (that would not
-serve any purpose here). Moreover, ``site_size=0.5`` guarantees that the two
-different types of triangles touch precisely: By default, `~kwant.plotter.plot`
-takes all sizes in units of the nearest-neighbor spacing. ``site_size=0.5``
-thus means half the distance between neighboring sites (and for the triangles
-this is interpreted as the radius of the inner circle).
-
-Finally, note that since we are dealing with a finalized system now, a site `i`
-is represented by an integer. In order to obtain the original
-`~kwant.builder.Site`, ``syst.sites[i]`` can be used.
-
-With this we arrive at
-
-.. image:: ../images/plot_graphene_data2.*
-
-with the same information as `~kwant.plotter.map`, but with a cleaner look.
-
-The way how data is presented of course influences what features of the data
-are best visible in a given plot. With `~kwant.plotter.plot` one can easily go
-beyond pseudocolor-like plots. For example, we can represent the wave function
-probability using the symbols itself:
-
-.. literalinclude:: plot_graphene.py
- :start-after: #HIDDEN_BEGIN_plotdata4
- :end-before: #HIDDEN_END_plotdata4
-
-Here, we choose the symbol size proportional to the wave function probability,
-while the site color is transparent to also allow for overlapping symbols to be
-visible. The hoppings are also plotted in order to show the underlying lattice.
-
-With this, we arrive at
-
-.. image:: ../images/plot_graphene_data3.*
-
-which shows the edge state nature of the wave function most clearly.
-
-.. rubric:: Footnotes
-
-.. [#] A dangling site is defined as having only one hopping connecting it to
- the rest. With next-nearest-neighbor hopping also all sites that are
- dangling with only nearest-neighbor hopping have more than one hopping.
-
-3D example: zincblende structure
-................................
-
-.. seealso::
- The complete source code of this example can be found in
- :download:`tutorial/plot_zincblende.py <../../../tutorial/plot_zincblende.py>`
-
-Zincblende is a very common crystal structure of semiconductors. It is a
-face-centered cubic crystal with two inequivalent atoms in the unit cell
-(i.e. two different types of atoms, unlike diamond which has the same crystal
-structure, but two equivalent atoms per unit cell).
-
-It is very easily generated in Kwant with `kwant.lattice.general`:
-
-.. literalinclude:: plot_zincblende.py
- :start-after: #HIDDEN_BEGIN_zincblende1
- :end-before: #HIDDEN_END_zincblende1
-
-Note how we keep references to the two different sublattices for later use.
-
-A three-dimensional structure is created as easily as in two dimensions, by
-using the `~kwant.lattice.PolyatomicLattice.shape`-functionality:
-
-.. literalinclude:: plot_zincblende.py
- :start-after: #HIDDEN_BEGIN_zincblende2
- :end-before: #HIDDEN_END_zincblende2
-
-We restrict ourselves here to a simple cuboid, and do not bother to add real
-values for onsite and hopping energies, but only the placeholder ``None`` (in a
-real calculation, several atomic orbitals would have to be considered).
-
-`~kwant.plotter.plot` can plot 3D systems just as easily as its two-dimensional
-counterparts:
-
-.. literalinclude:: plot_zincblende.py
- :start-after: #HIDDEN_BEGIN_plot1
- :end-before: #HIDDEN_END_plot1
-
-resulting in
-
-.. image:: ../images/plot_zincblende_syst1.*
-
-You might notice that the standard options for plotting are quite different in
-3D than in 2D. For example, by default hoppings are not printed, but sites are
-instead represented by little "balls" touching each other (which is achieved by
-a default ``site_size=0.5``). In fact, this style of plotting 3D shows quite
-decently the overall geometry of the system.
-
-When plotting into a window, the 3D plots can also be rotated and scaled
-arbitrarily, allowing for a good inspection of the geometry from all sides.
-
-.. note::
-
- Interactive 3D plots usually do not have the proper aspect ratio, but are a
- bit squashed. This is due to bugs in matplotlib's 3D plotting module that
- does not properly honor the corresponding arguments. By resizing the plot
- window however one can manually adjust the aspect ratio.
-
-Also for 3D it is possible to customize the plot. For example, we
-can explicitly plot the hoppings as lines, and color sites differently
-depending on the sublattice:
-
-.. literalinclude:: plot_zincblende.py
- :start-after: #HIDDEN_BEGIN_plot2
- :end-before: #HIDDEN_END_plot2
-
-which results in a 3D plot that allows to interactively (when plotted
-in a window) explore the crystal structure:
-
-.. image:: ../images/plot_zincblende_syst2.*
-
-Hence, a few lines of code using Kwant allow to explore all the different
-crystal lattices out there!
-
-.. note::
-
- - The 3D plots are in fact only *fake* 3D. For example, sites will always
- be plotted above hoppings (this is due to the limitations of matplotlib's
- 3d module)
- - Plotting hoppings in 3D is inherently much slower than plotting sites.
- Hence, this is not done by default.
-
-
-Interpolated density and current: QPC with disorder
-...................................................
-
-.. seealso::
- The complete source code of this example can be found in
- :download:`tutorial/plot_qpc.py <../../../tutorial/plot_qpc.py>`
-
-In the above examples we saw some useful methods for plotting systems where
-single-site resolution is required. Sometimes, however, having single-site
-precision is a hinderance, rather than a help, and looking at *averaged*
-quantities is more useful. This is particularly important in systems with a
-large number of sites, and systems that are discretizations of continuum
-models.
-
-Here we will show how to plot interpolated quantities using `kwant.plotter.map`
-and `kwant.plotter.current` using the example of a quantum point contact (QPC)
-with a perpendicular magnetic field and disorder:
-
-.. literalinclude:: plot_qpc.py
- :start-after: #HIDDEN_BEGIN_syst
- :end-before: #HIDDEN_END_syst
-
-.. image:: ../images/plot_qpc_syst.*
-
-Now we will compute the density of particles and current due to states
-originating in the left lead with energy 0.15.
-
-.. literalinclude:: plot_qpc.py
- :start-after: #HIDDEN_BEGIN_wf
- :end-before: #HIDDEN_END_wf
-
-We can then plot the density using `~kwant.plotter.map`:
-
-.. literalinclude:: plot_qpc.py
- :start-after: #HIDDEN_BEGIN_density
- :end-before: #HIDDEN_END_density
-
-.. image:: ../images/plot_qpc_density.*
-
-We pass ``method='linear'`` to ``map``, which produces a smoother plot than the
-default style. We see that density is concentrated on the edges of the sample,
-as we expect due to Hall effect induced by the perpendicular magnetic field.
-
-Plotting the current in the system will enable us to make even more sense
-of what is going on:
-
-.. literalinclude:: plot_qpc.py
- :start-after: #HIDDEN_BEGIN_current
- :end-before: #HIDDEN_END_current
-
-.. image:: ../images/plot_qpc_current.*
-
-We can now clearly see that current enters the system from the lower-left edge
-of the system (this matches our intuition, as we calculated the current for
-scattering states originating in the left-hand lead), and is backscattered by
-the restriction of the QPC in the centre.
diff --git a/doc/source/tutorial/tutorial7.rst b/doc/source/tutorial/tutorial7.rst
index 3ac75ba8936c8214ed2488246e045ef09d965263..dbf267042eb65f06d614fe3a04d18594b5cd3ec2 100644
--- a/doc/source/tutorial/tutorial7.rst
+++ b/doc/source/tutorial/tutorial7.rst
@@ -1,289 +1,286 @@
-##############################################################
-Calculating spectral density with the kernel polynomial method
-##############################################################
-
-We have already seen in the ":ref:`closed-systems`" tutorial that we can use
-Kwant simply to build Hamiltonians, which we can then directly diagonalize
-using routines from Scipy.
-
-This already allows us to treat systems with a few thousand sites without too
-many problems. For larger systems one is often not so interested in the exact
-eigenenergies and eigenstates, but more in the *density of states*.
-
-The kernel polynomial method (KPM), is an algorithm to obtain a polynomial
-expansion of the density of states. It can also be used to calculate the
-spectral density of arbitrary operators. Kwant has an implementation of the
-KPM method that is based on the algorithms presented in Ref. [1]_.
-
-Roughly speaking, KPM approximates the density of states (or any other spectral
-density) by expanding the action of the Hamiltonian (and operator of interest)
-on a (small) set of *random vectors* as a sum of Chebyshev polynomials up to
-some order, and then averaging. The accuracy of the method can be tuned by
-modifying the order of the Chebyshev expansion and the number of random
-vectors. See notes on accuracy_ below for details.
+Plotting Kwant systems and data in various styles
+-------------------------------------------------
+
+The plotting functionality of Kwant has been used extensively (through
+`~kwant.plotter.plot` and `~kwant.plotter.map`) in the previous tutorials. In
+addition to this basic use, `~kwant.plotter.plot` offers many options to change
+the plotting style extensively. It is the goal of this tutorial to show how
+these options can be used to achieve various very different objectives.
+
+2D example: graphene quantum dot
+................................
.. seealso::
The complete source code of this example can be found in
- :download:`tutorial/kernel_polynomial_method.py <../../../tutorial/kernel_polynomial_method.py>`
+ :download:`tutorial/plot_graphene.py <../../../tutorial/plot_graphene.py>`
+
+We begin by first considering a circular graphene quantum dot (similar to what
+has been used in parts of the tutorial :ref:`tutorial-graphene`.) In contrast
+to previous examples, we will also use hoppings beyond next-nearest neighbors:
-.. _accuracy:
-.. specialnote:: Performance and accuracy
+.. literalinclude:: plot_graphene.py
+ :start-after: #HIDDEN_BEGIN_makesyst
+ :end-before: #HIDDEN_END_makesyst
- The KPM method is especially well suited for large systems, and in the
- case when one is not interested in individual eigenvalues, but rather
- in obtaining an approximate spectral density.
+Note that adding hoppings hoppings to the `n`-th nearest neighbors can be
+simply done by passing `n` as an argument to
+`~kwant.lattice.Polyatomic.neighbors`. Also note that we use the method
+`~kwant.builder.Builder.eradicate_dangling` to get rid of single atoms sticking
+out of the shape. It is necessary to do so *before* adding the
+next-nearest-neighbor hopping [#]_.
- The accuracy in the energy resolution is dominated by the number of
- moments. The lowest accuracy is at the center of the spectrum, while
- slightly higher accuracy is obtained at the edges of the spectrum.
- If we use the KPM method (with the Jackson kernel, see Ref. [1]_) to
- describe a delta peak at the center of the spectrum, we will obtain a
- function similar to a Gaussian of width :math:`σ=πa/N`, where
- :math:`N` is the number of moments, and :math:`a` is the width of the
- spectrum.
+Of course, the system can be plotted simply with default settings:
- On the other hand, the random vectors will *explore* the range of the
- spectrum, and as the system gets bigger, the number of random vectors
- that are necessary to sample the whole spectrum reduces. Thus, a small
- number of random vectors is in general enough, and increasing this number
- will not result in a visible improvement of the approximation.
+.. literalinclude:: plot_graphene.py
+ :start-after: #HIDDEN_BEGIN_plotsyst1
+ :end-before: #HIDDEN_END_plotsyst1
+However, due to the richer structure of the lattice, this results in a rather
+busy plot:
-Introduction
-************
+.. image:: ../images/plot_graphene_syst1.*
-Our aim is to use the kernel polynomial method to obtain the spectral density
-:math:`Ï_A(E)`, as a function of the energy :math:`E`, of some Hilbert space
-operator :math:`A`. We define
+A much clearer plot can be obtained by using different colors for both
+sublattices, and by having different line widths for different hoppings. This
+can be achieved by passing a function to the arguments of
+`~kwant.plotter.plot`, instead of a constant. For properties of sites, this
+must be a function taking one site as argument, for hoppings a function taking
+the start end end site of hopping as arguments:
-.. math::
+.. literalinclude:: plot_graphene.py
+ :start-after: #HIDDEN_BEGIN_plotsyst2
+ :end-before: #HIDDEN_END_plotsyst2
- Ï_A(E) = Ï(E) A(E),
+Note that since we are using an unfinalized Builder, a ``site`` is really an
+instance of `~kwant.builder.Site`. With these adjustments we arrive at a plot
+that carries the same information, but is much easier to interpret:
-where :math:`A(E)` is the expectation value of :math:`A` for all the
-eigenstates of the Hamiltonian with energy :math:`E`, and the density of
-states is
+.. image:: ../images/plot_graphene_syst2.*
-.. math::
+Apart from plotting the *system* itself, `~kwant.plotter.plot` can also be used
+to plot *data* living on the system.
- Ï(E) = \frac{1}{D} \sum_{k=0}^{D-1} δ(E-E_k),
+As an example, we now compute the eigenstates of the graphene quantum dot and
+intend to plot the wave function probability in the quantum dot. For aesthetic
+reasons (the wave functions look a bit nicer), we restrict ourselves to
+nearest-neighbor hopping. Computing the wave functions is done in the usual
+way (note that for a large-scale system, one would probably want to use sparse
+linear algebra):
-:math:`D` being the Hilbert space dimension, and :math:`E_k` the eigenvalues.
+.. literalinclude:: plot_graphene.py
+ :start-after: #HIDDEN_BEGIN_plotdata1
+ :end-before: #HIDDEN_END_plotdata1
-In the special case when :math:`A` is the identity, then :math:`Ï_A(E)` is
-simply :math:`Ï(E)`, the density of states.
+In most cases, to plot the wave function probability, one wouldn't use
+`~kwant.plotter.plot`, but rather `~kwant.plotter.map`. Here, we plot the
+`n`-th wave function using it:
+.. literalinclude:: plot_graphene.py
+ :start-after: #HIDDEN_BEGIN_plotdata2
+ :end-before: #HIDDEN_END_plotdata2
-Calculating the density of states
-*********************************
+This results in a standard pseudocolor plot, showing in this case (``n=225``) a
+graphene edge state, i.e. a wave function mostly localized at the zigzag edges
+of the quantum dot.
-In the following example, we will use the KPM implementation in Kwant
-to obtain the density of states of a graphene disk.
+.. image:: ../images/plot_graphene_data1.*
-We start by importing kwant and defining our system.
+However although in general preferable, `~kwant.plotter.map` has a few
+deficiencies for this small system: For example, there are a few distortions at
+the edge of the dot. (This cannot be avoided in the type of interpolation used
+in `~kwant.plotter.map`). However, we can also use `~kwant.plotter.plot` to
+achieve a similar, but smoother result.
-.. literalinclude:: kernel_polynomial_method.py
- :start-after: #HIDDEN_BEGIN_sys1
- :end-before: #HIDDEN_END_sys1
+For this note that `~kwant.plotter.plot` can also take an array of floats (or
+function returning floats) as value for the ``site_color`` argument (the same
+holds for the hoppings). Via the colormap specified in ``cmap`` these are mapped
+to color, just as `~kwant.plotter.map` does! In addition, we can also change
+the symbol shape depending on the sublattice. With a triangle pointing up and
+down on the respective sublattice, the symbols used by plot fill the space
+completely:
-After making a system we can then create a `~kwant.kpm.SpectralDensity`
-object that represents the density of states for this system.
+.. literalinclude:: plot_graphene.py
+ :start-after: #HIDDEN_BEGIN_plotdata3
+ :end-before: #HIDDEN_END_plotdata3
-.. literalinclude:: kernel_polynomial_method.py
- :start-after: #HIDDEN_BEGIN_kpm1
- :end-before: #HIDDEN_END_kpm1
+Note that with ``hop_lw=0`` we deactivate plotting the hoppings (that would not
+serve any purpose here). Moreover, ``site_size=0.5`` guarantees that the two
+different types of triangles touch precisely: By default, `~kwant.plotter.plot`
+takes all sizes in units of the nearest-neighbor spacing. ``site_size=0.5``
+thus means half the distance between neighboring sites (and for the triangles
+this is interpreted as the radius of the inner circle).
-The `~kwant.kpm.SpectralDensity` can then be called like a function to obtain a
-sequence of energies in the spectrum of the Hamiltonian, and the corresponding
-density of states at these energies.
+Finally, note that since we are dealing with a finalized system now, a site `i`
+is represented by an integer. In order to obtain the original
+`~kwant.builder.Site`, ``syst.sites[i]`` can be used.
-.. literalinclude:: kernel_polynomial_method.py
- :start-after: #HIDDEN_BEGIN_kpm2
- :end-before: #HIDDEN_END_kpm2
+With this we arrive at
-When called with no arguments, an optimal set of energies is chosen (these are
-not evenly distributed over the spectrum, see Ref. [1]_ for details), however
-it is also possible to provide an explicit sequence of energies at which to
-evaluate the density of states.
+.. image:: ../images/plot_graphene_data2.*
-.. literalinclude:: kernel_polynomial_method.py
- :start-after: #HIDDEN_BEGIN_kpm3
- :end-before: #HIDDEN_END_kpm3
+with the same information as `~kwant.plotter.map`, but with a cleaner look.
-.. image:: ../images/kpm_dos.*
+The way how data is presented of course influences what features of the data
+are best visible in a given plot. With `~kwant.plotter.plot` one can easily go
+beyond pseudocolor-like plots. For example, we can represent the wave function
+probability using the symbols itself:
-In addition to being called like functions, `~kwant.kpm.SpectralDensity`
-objects also have a method `~kwant.kpm.SpectralDensity.average` which can be
-used to integrate the density of states against some distribution function over
-the whole spectrum. If no distribution function is specified, then the uniform
-distribution is used:
+.. literalinclude:: plot_graphene.py
+ :start-after: #HIDDEN_BEGIN_plotdata4
+ :end-before: #HIDDEN_END_plotdata4
-.. literalinclude:: kernel_polynomial_method.py
- :start-after: #HIDDEN_BEGIN_av1
- :end-before: #HIDDEN_END_av1
+Here, we choose the symbol size proportional to the wave function probability,
+while the site color is transparent to also allow for overlapping symbols to be
+visible. The hoppings are also plotted in order to show the underlying lattice.
-.. literalinclude:: ../images/kpm_normalization.txt
+With this, we arrive at
-We see that the integral of the density of states is normalized to 1. If
-we wish to calculate, say, the average number of states populated in
-equilibrium, then we should integrate with respect to a Fermi-Dirac
-distribution and multiply by the total number of available states in
-the system:
+.. image:: ../images/plot_graphene_data3.*
-.. literalinclude:: kernel_polynomial_method.py
- :start-after: #HIDDEN_BEGIN_av2
- :end-before: #HIDDEN_END_av2
+which shows the edge state nature of the wave function most clearly.
-.. literalinclude:: ../images/kpm_total_states.txt
+.. rubric:: Footnotes
-.. specialnote:: Stability and performance: spectral bounds
+.. [#] A dangling site is defined as having only one hopping connecting it to
+ the rest. With next-nearest-neighbor hopping also all sites that are
+ dangling with only nearest-neighbor hopping have more than one hopping.
- The KPM method internally rescales the spectrum of the Hamiltonian to the
- interval ``(-1, 1)`` (see Ref [1]_ for details), which requires calculating
- the boundaries of the spectrum (using ``scipy.sparse.linalg.eigsh``). This
- can be very costly for large systems, so it is possible to pass this
- explicitly as via the ``bounds`` parameter when instantiating the
- `~kwant.kpm.SpectralDensity` (see the class documentation for details).
+3D example: zincblende structure
+................................
- Additionally, `~kwant.kpm.SpectralDensity` accepts a parameter ``epsilon``,
- which ensures that the rescaled Hamiltonian (used internally), always has a
- spectrum strictly contained in the interval ``(-1, 1)``. If bounds are not
- provided then the tolerance on the bounds calculated with
- ``scipy.sparse.linalg.eigsh`` is set to ``epsilon/2``.
+.. seealso::
+ The complete source code of this example can be found in
+ :download:`tutorial/plot_zincblende.py <../../../tutorial/plot_zincblende.py>`
+Zincblende is a very common crystal structure of semiconductors. It is a
+face-centered cubic crystal with two inequivalent atoms in the unit cell
+(i.e. two different types of atoms, unlike diamond which has the same crystal
+structure, but two equivalent atoms per unit cell).
-Increasing the accuracy of the approximation
-********************************************
+It is very easily generated in Kwant with `kwant.lattice.general`:
-`~kwant.kpm.SpectralDensity` has two methods for increasing the accuracy
-of the method, each of which offers different levels of control over what
-exactly is changed.
+.. literalinclude:: plot_zincblende.py
+ :start-after: #HIDDEN_BEGIN_zincblende1
+ :end-before: #HIDDEN_END_zincblende1
-The simplest way to obtain a more accurate solution is to use the
-``add_moments`` method:
+Note how we keep references to the two different sublattices for later use.
-.. literalinclude:: kernel_polynomial_method.py
- :start-after: #HIDDEN_BEGIN_acc1
- :end-before: #HIDDEN_END_acc1
+A three-dimensional structure is created as easily as in two dimensions, by
+using the `~kwant.lattice.PolyatomicLattice.shape`-functionality:
-This will update the number of calculated moments and also the
-number of sampling points such that the maximum distance between successive
-energy points is ``energy_resolution`` (see notes on accuracy_).
+.. literalinclude:: plot_zincblende.py
+ :start-after: #HIDDEN_BEGIN_zincblende2
+ :end-before: #HIDDEN_END_zincblende2
-.. image:: ../images/kpm_dos_acc.*
+We restrict ourselves here to a simple cuboid, and do not bother to add real
+values for onsite and hopping energies, but only the placeholder ``None`` (in a
+real calculation, several atomic orbitals would have to be considered).
-Alternatively, you can directly increase the number of moments
-with ``add_moments``, or the number of random vectors with ``add_vectors``.
+`~kwant.plotter.plot` can plot 3D systems just as easily as its two-dimensional
+counterparts:
-.. literalinclude:: kernel_polynomial_method.py
- :start-after: #HIDDEN_BEGIN_acc2
- :end-before: #HIDDEN_END_acc2
+.. literalinclude:: plot_zincblende.py
+ :start-after: #HIDDEN_BEGIN_plot1
+ :end-before: #HIDDEN_END_plot1
-.. image:: ../images/kpm_dos_r.*
+resulting in
+.. image:: ../images/plot_zincblende_syst1.*
-.. _operator_spectral_density:
+You might notice that the standard options for plotting are quite different in
+3D than in 2D. For example, by default hoppings are not printed, but sites are
+instead represented by little "balls" touching each other (which is achieved by
+a default ``site_size=0.5``). In fact, this style of plotting 3D shows quite
+decently the overall geometry of the system.
-Calculating the spectral density of an operator
-***********************************************
+When plotting into a window, the 3D plots can also be rotated and scaled
+arbitrarily, allowing for a good inspection of the geometry from all sides.
-Above, we saw how to calculate the density of states by creating a
-`~kwant.kpm.SpectralDensity` and passing it a finalized Kwant system.
-When instantiating a `~kwant.kpm.SpectralDensity` we may optionally
-supply an operator in addition to the system. In this case it is
-the spectral density of the given operator that is calculated.
+.. note::
-`~kwant.kpm.SpectralDensity` accepts the operators in a few formats:
+ Interactive 3D plots usually do not have the proper aspect ratio, but are a
+ bit squashed. This is due to bugs in matplotlib's 3D plotting module that
+ does not properly honor the corresponding arguments. By resizing the plot
+ window however one can manually adjust the aspect ratio.
-* *explicit matrices* (numpy array of scipy sparse matrices will work)
-* *operators* from `kwant.operator`
+Also for 3D it is possible to customize the plot. For example, we
+can explicitly plot the hoppings as lines, and color sites differently
+depending on the sublattice:
-If an explicit matrix is provided then it must have the same
-shape as the system Hamiltonian.
+.. literalinclude:: plot_zincblende.py
+ :start-after: #HIDDEN_BEGIN_plot2
+ :end-before: #HIDDEN_END_plot2
-.. literalinclude:: kernel_polynomial_method.py
- :start-after: #HIDDEN_BEGIN_op1
- :end-before: #HIDDEN_END_op1
+which results in a 3D plot that allows to interactively (when plotted
+in a window) explore the crystal structure:
+.. image:: ../images/plot_zincblende_syst2.*
-Or, to do the same calculation using `kwant.operator.Density`:
+Hence, a few lines of code using Kwant allow to explore all the different
+crystal lattices out there!
-.. literalinclude:: kernel_polynomial_method.py
- :start-after: #HIDDEN_BEGIN_op2
- :end-before: #HIDDEN_END_op2
+.. note::
-Using operators from `kwant.operator` allows us to calculate quantities
-such as the *local* density of states by telling the operator not to
-sum over all the sites of the system:
+ - The 3D plots are in fact only *fake* 3D. For example, sites will always
+ be plotted above hoppings (this is due to the limitations of matplotlib's
+ 3d module)
+ - Plotting hoppings in 3D is inherently much slower than plotting sites.
+ Hence, this is not done by default.
-.. literalinclude:: kernel_polynomial_method.py
- :start-after: #HIDDEN_BEGIN_op3
- :end-before: #HIDDEN_END_op3
-`~kwant.kpm.SpectralDensity` will properly handle this vector output,
-which allows us to plot the local density of states at different
-point in the spectrum:
+Interpolated density and current: QPC with disorder
+...................................................
-.. literalinclude:: kernel_polynomial_method.py
- :start-after: #HIDDEN_BEGIN_op4
- :end-before: #HIDDEN_END_op4
+.. seealso::
+ The complete source code of this example can be found in
+ :download:`tutorial/plot_qpc.py <../../../tutorial/plot_qpc.py>`
-.. image:: ../images/kpm_ldos.*
+In the above examples we saw some useful methods for plotting systems where
+single-site resolution is required. Sometimes, however, having single-site
+precision is a hinderance, rather than a help, and looking at *averaged*
+quantities is more useful. This is particularly important in systems with a
+large number of sites, and systems that are discretizations of continuum
+models.
-This nicely illustrates the edge states of the graphene dot at zero
-energy, and the bulk states at higher energy.
+Here we will show how to plot interpolated quantities using `kwant.plotter.map`
+and `kwant.plotter.current` using the example of a quantum point contact (QPC)
+with a perpendicular magnetic field and disorder:
+.. literalinclude:: plot_qpc.py
+ :start-after: #HIDDEN_BEGIN_syst
+ :end-before: #HIDDEN_END_syst
-Advanced topics
-***************
+.. image:: ../images/plot_qpc_syst.*
-Custom distributions for random vectors
-^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
-By default `~kwant.kpm.SpectralDensity` will use random vectors
-whose components are unit complex numbers with phases drawn
-from a uniform distribution. There are several reasons why you may
-wish to make a different choice of distribution for your random vectors,
-for example to enforce certain symmetries or to only use real-valued vectors.
+Now we will compute the density of particles and current due to states
+originating in the left lead with energy 0.15.
-To change how the random vectors are generated, you need only specify a
-function that takes the dimension of the Hilbert space as a single parameter,
-and which returns a vector in that Hilbert space:
+.. literalinclude:: plot_qpc.py
+ :start-after: #HIDDEN_BEGIN_wf
+ :end-before: #HIDDEN_END_wf
-.. literalinclude:: kernel_polynomial_method.py
- :start-after: #HIDDEN_BEGIN_fact1
- :end-before: #HIDDEN_END_fact1
+We can then plot the density using `~kwant.plotter.map`:
-Reproducible calculations
-^^^^^^^^^^^^^^^^^^^^^^^^^
-Because KPM internally uses random vectors, running the same calculation
-twice will not give bit-for-bit the same result. However, similarly to
-the funcions in `~kwant.rmt`, the random number generator can be directly
-manipulated by passing a value to the ``rng`` parameter of
-`~kwant.kpm.SpectralDensity`. ``rng`` can itself be a random number generator,
-or it may simply be a seed to pass to the numpy random number generator
-(that is used internally by default).
+.. literalinclude:: plot_qpc.py
+ :start-after: #HIDDEN_BEGIN_density
+ :end-before: #HIDDEN_END_density
-Defining operators as sesquilinear maps
-^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
-`Above`__, we showed how `~kwant.kpm.SpectralDensity` can calculate the
-spectral density of operators, and how we can define operators by using
-`kwant.operator`. If you need even more flexibility,
-`~kwant.kpm.SpectralDensity` will also accept a *function* as its ``operator``
-parameter. This function must itself take two parameters, ``(bra, ket)`` and
-must return either a scalar or a one-dimensional array. In order to be
-meaningful the function must be a sesquilinear map, i.e. antilinear in its
-first argument, and linear in its second argument. Below, we compare two
-methods for computing the local density of states, one using
-`kwant.operator.Density`, and the other using a custom function.
+.. image:: ../images/plot_qpc_density.*
-.. literalinclude:: kernel_polynomial_method.py
- :start-after: #HIDDEN_BEGIN_blm
- :end-before: #HIDDEN_END_blm
+We pass ``method='linear'`` to ``map``, which produces a smoother plot than the
+default style. We see that density is concentrated on the edges of the sample,
+as we expect due to Hall effect induced by the perpendicular magnetic field.
-__ operator_spectral_density_
+Plotting the current in the system will enable us to make even more sense
+of what is going on:
+.. literalinclude:: plot_qpc.py
+ :start-after: #HIDDEN_BEGIN_current
+ :end-before: #HIDDEN_END_current
-.. rubric:: References
+.. image:: ../images/plot_qpc_current.*
-.. [1] `Rev. Mod. Phys., Vol. 78, No. 1 (2006)
- <https://arxiv.org/abs/cond-mat/0504627>`_.
+We can now clearly see that current enters the system from the lower-left edge
+of the system (this matches our intuition, as we calculated the current for
+scattering states originating in the left-hand lead), and is backscattered by
+the restriction of the QPC in the centre.
diff --git a/doc/source/tutorial/tutorial8.rst b/doc/source/tutorial/tutorial8.rst
index c210d9290fb158bb7f57f1a63bf853505b3debb1..f40d1920231ee7602e9d93c2f472a6e9cb023e22 100644
--- a/doc/source/tutorial/tutorial8.rst
+++ b/doc/source/tutorial/tutorial8.rst
@@ -1,323 +1,289 @@
-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.
+##############################################################
+Calculating spectral density with the kernel polynomial method
+##############################################################
+
+We have already seen in the ":ref:`closed-systems`" tutorial that we can use
+Kwant simply to build Hamiltonians, which we can then directly diagonalize
+using routines from Scipy.
+
+This already allows us to treat systems with a few thousand sites without too
+many problems. For larger systems one is often not so interested in the exact
+eigenenergies and eigenstates, but more in the *density of states*.
+
+The kernel polynomial method (KPM), is an algorithm to obtain a polynomial
+expansion of the density of states. It can also be used to calculate the
+spectral density of arbitrary operators. Kwant has an implementation of the
+KPM method that is based on the algorithms presented in Ref. [1]_.
+
+Roughly speaking, KPM approximates the density of states (or any other spectral
+density) by expanding the action of the Hamiltonian (and operator of interest)
+on a (small) set of *random vectors* as a sum of Chebyshev polynomials up to
+some order, and then averaging. The accuracy of the method can be tuned by
+modifying the order of the Chebyshev expansion and the number of random
+vectors. See notes on accuracy_ below for details.
.. seealso::
- The complete source code of this tutorial can be found in
- :download:`tutorial/continuum_discretizer.py <../../../tutorial/continuum_discretizer.py>`
+ The complete source code of this example can be found in
+ :download:`tutorial/kernel_polynomial_method.py <../../../tutorial/kernel_polynomial_method.py>`
+.. _accuracy:
+.. specialnote:: Performance and accuracy
-.. _tutorial_discretizer_introduction:
+ The KPM method is especially well suited for large systems, and in the
+ case when one is not interested in individual eigenvalues, but rather
+ in obtaining an approximate spectral density.
-Discretizing by hand
-....................
+ The accuracy in the energy resolution is dominated by the number of
+ moments. The lowest accuracy is at the center of the spectrum, while
+ slightly higher accuracy is obtained at the edges of the spectrum.
+ If we use the KPM method (with the Jackson kernel, see Ref. [1]_) to
+ describe a delta peak at the center of the spectrum, we will obtain a
+ function similar to a Gaussian of width :math:`σ=πa/N`, where
+ :math:`N` is the number of moments, and :math:`a` is the width of the
+ spectrum.
-As an example, let us consider the following continuum Schrödinger equation
-for a semiconducting heterostructure (using the effective mass approximation):
+ On the other hand, the random vectors will *explore* the range of the
+ spectrum, and as the system gets bigger, the number of random vectors
+ that are necessary to sample the whole spectrum reduces. Thus, a small
+ number of random vectors is in general enough, and increasing this number
+ will not result in a visible improvement of the approximation.
-.. math::
- \left( k_x \frac{\hbar^2}{2 m(x)} k_x \right) \psi(x) = E \, \psi(x).
+Introduction
+************
-Replacing the momenta by their corresponding differential operators
+Our aim is to use the kernel polynomial method to obtain the spectral density
+:math:`Ï_A(E)`, as a function of the energy :math:`E`, of some Hilbert space
+operator :math:`A`. We define
.. 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::
+ Ï_A(E) = Ï(E) A(E),
- 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
-........................................................
-
-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).
-
-.. literalinclude:: continuum_discretizer.py
- :start-after: #HIDDEN_BEGIN_symbolic_discretization
- :end-before: #HIDDEN_END_symbolic_discretization
-
-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.
-
-Setting the ``verbose`` parameter to ``True`` prints extra information about the
-onsite and hopping functions assigned to the ``Builder`` produced
-by ``discretize``:
-
-.. literalinclude:: ../images/discretizer_intro_verbose.txt
-
-.. specialnote:: Technical details
-
- - ``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 <http://docs.sympy.org/latest/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
+where :math:`A(E)` is the expectation value of :math:`A` for all the
+eigenstates of the Hamiltonian with energy :math:`E`, and the density of
+states is
.. 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:
-
-.. literalinclude:: continuum_discretizer.py
- :start-after: #HIDDEN_BEGIN_template
- :end-before: #HIDDEN_END_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:
+ Ï(E) = \frac{1}{D} \sum_{k=0}^{D-1} δ(E-E_k),
-.. literalinclude:: continuum_discretizer.py
- :start-after: #HIDDEN_BEGIN_fill
- :end-before: #HIDDEN_END_fill
+:math:`D` being the Hilbert space dimension, and :math:`E_k` the eigenvalues.
-After finalizing this system, we can plot one of the system's
-energy eigenstates:
+In the special case when :math:`A` is the identity, then :math:`Ï_A(E)` is
+simply :math:`Ï(E)`, the density of states.
-.. literalinclude:: continuum_discretizer.py
- :start-after: #HIDDEN_BEGIN_plot_eigenstate
- :end-before: #HIDDEN_END_plot_eigenstate
-.. image:: ../images/discretizer_gs.*
+Calculating the density of states
+*********************************
-Note in the above that we provided the function ``V`` to
-``syst.hamiltonian_submatrix`` using ``params=dict(V=potential)``, rather than
-via ``args``.
+In the following example, we will use the KPM implementation in Kwant
+to obtain the density of states of a graphene disk.
-In addition, the function passed as ``V`` expects two input parameters ``x``
-and ``y``, the same as in the initial continuum Hamiltonian.
+We start by importing kwant and defining our system.
+.. literalinclude:: kernel_polynomial_method.py
+ :start-after: #HIDDEN_BEGIN_sys1
+ :end-before: #HIDDEN_END_sys1
-Models with more structure: Bernevig-Hughes-Zhang
-.................................................
+After making a system we can then create a `~kwant.kpm.SpectralDensity`
+object that represents the density of states for this system.
-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:
+.. literalinclude:: kernel_polynomial_method.py
+ :start-after: #HIDDEN_BEGIN_kpm1
+ :end-before: #HIDDEN_END_kpm1
-.. literalinclude:: continuum_discretizer.py
- :start-after: #HIDDEN_BEGIN_define_qsh
- :end-before: #HIDDEN_END_define_qsh
+The `~kwant.kpm.SpectralDensity` can then be called like a function to obtain a
+sequence of energies in the spectrum of the Hamiltonian, and the corresponding
+density of states at these energies.
-We can then make a ribbon out of this template system:
+.. literalinclude:: kernel_polynomial_method.py
+ :start-after: #HIDDEN_BEGIN_kpm2
+ :end-before: #HIDDEN_END_kpm2
-.. literalinclude:: continuum_discretizer.py
- :start-after: #HIDDEN_BEGIN_define_qsh_build
- :end-before: #HIDDEN_END_define_qsh_build
+When called with no arguments, an optimal set of energies is chosen (these are
+not evenly distributed over the spectrum, see Ref. [1]_ for details), however
+it is also possible to provide an explicit sequence of energies at which to
+evaluate the density of states.
-and plot its dispersion using `kwant.plotter.bands`:
+.. literalinclude:: kernel_polynomial_method.py
+ :start-after: #HIDDEN_BEGIN_kpm3
+ :end-before: #HIDDEN_END_kpm3
-.. literalinclude:: continuum_discretizer.py
- :start-after: #HIDDEN_BEGIN_plot_qsh_band
- :end-before: #HIDDEN_END_plot_qsh_band
+.. image:: ../images/kpm_dos.*
-.. image:: ../images/discretizer_qsh_band.*
+In addition to being called like functions, `~kwant.kpm.SpectralDensity`
+objects also have a method `~kwant.kpm.SpectralDensity.average` which can be
+used to integrate the density of states against some distribution function over
+the whole spectrum. If no distribution function is specified, then the uniform
+distribution is used:
-In the above we see the edge states of the quantum spin Hall effect, which
-we can visualize using `kwant.plotter.map`:
+.. literalinclude:: kernel_polynomial_method.py
+ :start-after: #HIDDEN_BEGIN_av1
+ :end-before: #HIDDEN_END_av1
-.. literalinclude:: continuum_discretizer.py
- :start-after: #HIDDEN_BEGIN_plot_qsh_wf
- :end-before: #HIDDEN_END_plot_qsh_wf
+.. literalinclude:: ../images/kpm_normalization.txt
-.. image:: ../images/discretizer_qsh_wf.*
+We see that the integral of the density of states is normalized to 1. If
+we wish to calculate, say, the average number of states populated in
+equilibrium, then we should integrate with respect to a Fermi-Dirac
+distribution and multiply by the total number of available states in
+the system:
+.. literalinclude:: kernel_polynomial_method.py
+ :start-after: #HIDDEN_BEGIN_av2
+ :end-before: #HIDDEN_END_av2
-Limitations of discretization
-.............................
+.. literalinclude:: ../images/kpm_total_states.txt
-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
+.. specialnote:: Stability and performance: spectral bounds
-.. math::
-
- H_\textrm{continuous}(k_x) = \frac{\hbar^2}{2m}k_x^2
-
-
-and its discretized approximation
-
-.. math::
+ The KPM method internally rescales the spectrum of the Hamiltonian to the
+ interval ``(-1, 1)`` (see Ref [1]_ for details), which requires calculating
+ the boundaries of the spectrum (using ``scipy.sparse.linalg.eigsh``). This
+ can be very costly for large systems, so it is possible to pass this
+ explicitly as via the ``bounds`` parameter when instantiating the
+ `~kwant.kpm.SpectralDensity` (see the class documentation for details).
- H_\textrm{tight-binding}(k_x) = 2t \big(1 - \cos(k_x a)\big),
+ Additionally, `~kwant.kpm.SpectralDensity` accepts a parameter ``epsilon``,
+ which ensures that the rescaled Hamiltonian (used internally), always has a
+ spectrum strictly contained in the interval ``(-1, 1)``. If bounds are not
+ provided then the tolerance on the bounds calculated with
+ ``scipy.sparse.linalg.eigsh`` is set to ``epsilon/2``.
-where :math:`t=\frac{\hbar^2}{2ma^2}`, are only valid in the limit
-:math:`E \lt t`. The grid spacing :math:`a` must be chosen according
-to how high in energy you need your tight-binding model to be valid.
+Increasing the accuracy of the approximation
+********************************************
-It is possible to set :math:`a` through the ``grid_spacing`` parameter
-to `~kwant.continuum.discretize`, as we will illustrate in the following
-example. Let us start from the continuum Hamiltonian
+`~kwant.kpm.SpectralDensity` has two methods for increasing the accuracy
+of the method, each of which offers different levels of control over what
+exactly is changed.
-.. math::
+The simplest way to obtain a more accurate solution is to use the
+``add_moments`` method:
- H(k) = k_x^2 \mathbb{1}_{2\times2} + α k_x \sigma_y.
+.. literalinclude:: kernel_polynomial_method.py
+ :start-after: #HIDDEN_BEGIN_acc1
+ :end-before: #HIDDEN_END_acc1
-We start by defining this model as a string and setting the value of the
-:math:`α` parameter:
+This will update the number of calculated moments and also the default
+number of sampling points such that the maximum distance between successive
+energy points is ``energy_resolution`` (see notes on accuracy_).
-.. literalinclude:: continuum_discretizer.py
- :start-after: #HIDDEN_BEGIN_ls_def
- :end-before: #HIDDEN_END_ls_def
+.. image:: ../images/kpm_dos_acc.*
-Now we can use `kwant.continuum.lambdify` to obtain a function that computes
-:math:`H(k)`:
+Alternatively, you can directly increase the number of moments
+with ``add_moments``, or the number of random vectors with ``add_vectors``.
-.. literalinclude:: continuum_discretizer.py
- :start-after: #HIDDEN_BEGIN_ls_hk_cont
- :end-before: #HIDDEN_END_ls_hk_cont
+.. literalinclude:: kernel_polynomial_method.py
+ :start-after: #HIDDEN_BEGIN_acc2
+ :end-before: #HIDDEN_END_acc2
-We can also construct a discretized approximation using
-`kwant.continuum.discretize`, in a similar manner to previous examples:
+.. image:: ../images/kpm_dos_r.*
-.. literalinclude:: continuum_discretizer.py
- :start-after: #HIDDEN_BEGIN_ls_hk_tb
- :end-before: #HIDDEN_END_ls_hk_tb
-Below we can see the continuum and tight-binding dispersions for two
-different values of the discretization grid spacing :math:`a`:
+.. _operator_spectral_density:
-.. image:: ../images/discretizer_lattice_spacing.*
+Calculating the spectral density of an operator
+***********************************************
+Above, we saw how to calculate the density of states by creating a
+`~kwant.kpm.SpectralDensity` and passing it a finalized Kwant system.
+When instantiating a `~kwant.kpm.SpectralDensity` we may optionally
+supply an operator in addition to the system. In this case it is
+the spectral density of the given operator that is calculated.
-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}]`.
+`~kwant.kpm.SpectralDensity` accepts the operators in a few formats:
-.. specialnote:: Note
+* *explicit matrices* (numpy array of scipy sparse matrices will work)
+* *operators* from `kwant.operator`
- The use of `~kwant.wraparound.wraparound` makes it easy to extend this
- example to multidimensional band structures. ``wraparound`` replaces
- the translational symmetry in the ``x`` direction by a parameter ``k_x``.
+If an explicit matrix is provided then it must have the same
+shape as the system Hamiltonian.
-Advanced topics
-...............
+.. literalinclude:: kernel_polynomial_method.py
+ :start-after: #HIDDEN_BEGIN_op1
+ :end-before: #HIDDEN_END_op1
-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``.
+Or, to do the same calculation using `kwant.operator.Density`:
-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:
+.. literalinclude:: kernel_polynomial_method.py
+ :start-after: #HIDDEN_BEGIN_op2
+ :end-before: #HIDDEN_END_op2
-.. literalinclude:: continuum_discretizer.py
- :start-after: #HIDDEN_BEGIN_subs_1
- :end-before: #HIDDEN_END_subs_1
+Using operators from `kwant.operator` allows us to calculate quantities
+such as the *local* density of states by telling the operator not to
+sum over all the sites of the system:
-.. literalinclude:: ../images/discretizer_subs_1.txt
+.. literalinclude:: kernel_polynomial_method.py
+ :start-after: #HIDDEN_BEGIN_op3
+ :end-before: #HIDDEN_END_op3
-We can use the ``subs`` keyword parameter to substitute expressions
-and numerical values:
+`~kwant.kpm.SpectralDensity` will properly handle this vector output,
+which allows us to plot the local density of states at different
+point in the spectrum:
-.. literalinclude:: continuum_discretizer.py
- :start-after: #HIDDEN_BEGIN_subs_2
- :end-before: #HIDDEN_END_subs_2
+.. literalinclude:: kernel_polynomial_method.py
+ :start-after: #HIDDEN_BEGIN_op4
+ :end-before: #HIDDEN_END_op4
-.. literalinclude:: ../images/discretizer_subs_2.txt
+.. image:: ../images/kpm_ldos.*
-Symbolic expressions obtained in this way can be directly passed to all
-``discretizer`` functions.
+This nicely illustrates the edge states of the graphene dot at zero
+energy, and the bulk states at higher energy.
-.. specialnote:: Technical details
- 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`.
+Advanced topics
+***************
+
+Custom distributions for random vectors
+^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
+By default `~kwant.kpm.SpectralDensity` will use random vectors
+whose components are unit complex numbers with phases drawn
+from a uniform distribution. There are several reasons why you may
+wish to make a different choice of distribution for your random vectors,
+for example to enforce certain symmetries or to only use real-valued vectors.
+
+To change how the random vectors are generated, you need only specify a
+function that takes the dimension of the Hilbert space as a single parameter,
+and which returns a vector in that Hilbert space:
+
+.. literalinclude:: kernel_polynomial_method.py
+ :start-after: #HIDDEN_BEGIN_fact1
+ :end-before: #HIDDEN_END_fact1
+
+Reproducible calculations
+^^^^^^^^^^^^^^^^^^^^^^^^^
+Because KPM internally uses random vectors, running the same calculation
+twice will not give bit-for-bit the same result. However, similarly to
+the funcions in `~kwant.rmt`, the random number generator can be directly
+manipulated by passing a value to the ``rng`` parameter of
+`~kwant.kpm.SpectralDensity`. ``rng`` can itself be a random number generator,
+or it may simply be a seed to pass to the numpy random number generator
+(that is used internally by default).
+
+Defining operators as sesquilinear maps
+^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
+`Above`__, we showed how `~kwant.kpm.SpectralDensity` can calculate the
+spectral density of operators, and how we can define operators by using
+`kwant.operator`. If you need even more flexibility,
+`~kwant.kpm.SpectralDensity` will also accept a *function* as its ``operator``
+parameter. This function must itself take two parameters, ``(bra, ket)`` and
+must return either a scalar or a one-dimensional array. In order to be
+meaningful the function must be a sesquilinear map, i.e. antilinear in its
+first argument, and linear in its second argument. Below, we compare two
+methods for computing the local density of states, one using
+`kwant.operator.Density`, and the other using a custom function.
+
+.. literalinclude:: kernel_polynomial_method.py
+ :start-after: #HIDDEN_BEGIN_blm
+ :end-before: #HIDDEN_END_blm
+
+__ operator_spectral_density_
.. 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>`_.
+.. [1] `Rev. Mod. Phys., Vol. 78, No. 1 (2006)
+ <https://arxiv.org/abs/cond-mat/0504627>`_.
diff --git a/doc/source/tutorial/tutorial9.rst b/doc/source/tutorial/tutorial9.rst
new file mode 100644
index 0000000000000000000000000000000000000000..216f24b116f0d4036ac62b0ae77d37eccc302501
--- /dev/null
+++ b/doc/source/tutorial/tutorial9.rst
@@ -0,0 +1,318 @@
+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::
+ The complete source code of this tutorial can be found in
+ :download:`tutorial/continuum_discretizer.py <../../../tutorial/continuum_discretizer.py>`
+
+
+.. _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
+........................................................
+
+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).
+
+.. literalinclude:: continuum_discretizer.py
+ :start-after: #HIDDEN_BEGIN_symbolic_discretization
+ :end-before: #HIDDEN_END_symbolic_discretization
+
+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.
+
+Setting the ``verbose`` parameter to ``True`` prints extra information about the
+onsite and hopping functions assigned to the ``Builder`` produced
+by ``discretize``:
+
+.. literalinclude:: ../images/discretizer_intro_verbose.txt
+
+.. specialnote:: Technical details
+
+ - ``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 <http://docs.sympy.org/latest/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:
+
+.. literalinclude:: continuum_discretizer.py
+ :start-after: #HIDDEN_BEGIN_template
+ :end-before: #HIDDEN_END_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:: continuum_discretizer.py
+ :start-after: #HIDDEN_BEGIN_fill
+ :end-before: #HIDDEN_END_fill
+
+After finalizing this system, we can plot one of the system's
+energy eigenstates:
+
+.. literalinclude:: continuum_discretizer.py
+ :start-after: #HIDDEN_BEGIN_plot_eigenstate
+ :end-before: #HIDDEN_END_plot_eigenstate
+
+.. image:: ../images/discretizer_gs.*
+
+Note in the above that we provided the function ``V`` to
+``syst.hamiltonian_submatrix`` using ``params=dict(V=potential)``, rather than
+via ``args``.
+
+In addition, the function passed as ``V`` expects two input parameters ``x``
+and ``y``, the same as in the initial continuum Hamiltonian.
+
+
+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:
+
+.. literalinclude:: continuum_discretizer.py
+ :start-after: #HIDDEN_BEGIN_define_qsh
+ :end-before: #HIDDEN_END_define_qsh
+
+We can then make a ribbon out of this template system:
+
+.. literalinclude:: continuum_discretizer.py
+ :start-after: #HIDDEN_BEGIN_define_qsh_build
+ :end-before: #HIDDEN_END_define_qsh_build
+
+and plot its dispersion using `kwant.plotter.bands`:
+
+.. literalinclude:: continuum_discretizer.py
+ :start-after: #HIDDEN_BEGIN_plot_qsh_band
+ :end-before: #HIDDEN_END_plot_qsh_band
+
+.. image:: ../images/discretizer_qsh_band.*
+
+In the above we see the edge states of the quantum spin Hall effect, which
+we can visualize using `kwant.plotter.map`:
+
+.. literalinclude:: continuum_discretizer.py
+ :start-after: #HIDDEN_BEGIN_plot_qsh_wf
+ :end-before: #HIDDEN_END_plot_qsh_wf
+
+.. image:: ../images/discretizer_qsh_wf.*
+
+
+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 \lt 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_spacing`` 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:
+
+.. literalinclude:: continuum_discretizer.py
+ :start-after: #HIDDEN_BEGIN_ls_def
+ :end-before: #HIDDEN_END_ls_def
+
+Now we can use `kwant.continuum.lambdify` to obtain a function that computes
+:math:`H(k)`:
+
+.. literalinclude:: continuum_discretizer.py
+ :start-after: #HIDDEN_BEGIN_ls_hk_cont
+ :end-before: #HIDDEN_END_ls_hk_cont
+
+We can also construct a discretized approximation using
+`kwant.continuum.discretize`, in a similar manner to previous examples:
+
+.. literalinclude:: continuum_discretizer.py
+ :start-after: #HIDDEN_BEGIN_ls_hk_tb
+ :end-before: #HIDDEN_END_ls_hk_tb
+
+Below we can see the continuum and tight-binding dispersions for two
+different values of the discretization grid spacing :math:`a`:
+
+.. image:: ../images/discretizer_lattice_spacing.*
+
+
+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:
+
+.. literalinclude:: continuum_discretizer.py
+ :start-after: #HIDDEN_BEGIN_subs_1
+ :end-before: #HIDDEN_END_subs_1
+
+.. literalinclude:: ../images/discretizer_subs_1.txt
+
+We can use the ``substitutions`` keyword parameter to substitute expressions
+and numerical values:
+
+.. literalinclude:: continuum_discretizer.py
+ :start-after: #HIDDEN_BEGIN_subs_2
+ :end-before: #HIDDEN_END_subs_2
+
+.. literalinclude:: ../images/discretizer_subs_2.txt
+
+Symbolic expressions obtained in this way can be directly passed to all
+``discretizer`` functions.
+
+.. specialnote:: Technical details
+
+ 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>`_.