diff --git a/.gitignore b/.gitignore
index 72b8ea2fd3736c407ac8b321b17bc6d4937ae0cb..ba94c807d204266b925eb7c25a92bc57bfe42971 100644
--- a/.gitignore
+++ b/.gitignore
@@ -8,13 +8,7 @@
/dist
/doc/build
/doc/source/reference/generated/
-/doc/source/code/include/*.py
-/doc/source/code/figure/*.png
-/doc/source/code/figure/*.pdf
-/doc/source/code/figure/.*_flag
-/doc/source/code/figure/[a-zA-Z]*.py
-/doc/source/code/figure/*.txt
-/doc/source/code/download/
+/doc/source/figure/*.pdf
/build.conf
/kwant.egg-info/
/MANIFEST.in
diff --git a/.gitlab-ci.yml b/.gitlab-ci.yml
index ffffe78d75dda5e855b7acdbf8088e76770ff1c9..a3ef10b22292cf12ec26d224aae75e1409dc8c23 100644
--- a/.gitlab-ci.yml
+++ b/.gitlab-ci.yml
@@ -239,7 +239,10 @@ build documentation:
- build:latest
stage: test
script:
- - make -C doc realclean; make -C doc html SPHINXOPTS='-A website_deploy=True -n -W' SOURCE_LINK_TEMPLATE="$CI_PROJECT_URL"/blob/\$\$r/\$\$f
+ - apt-get update && apt-get install -y librsvg2-bin # for converting svgs to pdfs
+ - pip install git+https://github.com/jupyter-widgets/jupyter-sphinx sphinxcontrib-svg2pdfconverter
+ - python -m ipykernel install --user --name kwant-latest
+ - make -C doc clean; make -C doc html SPHINXOPTS='-A website_deploy=True -n -W -D jupyter_execute_default_kernel=kwant-latest' SOURCE_LINK_TEMPLATE="$CI_PROJECT_URL"/blob/\$\$r/\$\$f
artifacts:
paths:
- doc/build/html/
@@ -251,7 +254,10 @@ build PDF documentation:
- build:latest
stage: test
script:
- - make -C doc latex SPHINXOPTS='-n -W'
+ - apt-get update && apt-get install -y librsvg2-bin # for converting svgs to pdfs
+ - pip install git+https://github.com/jupyter-widgets/jupyter-sphinx sphinxcontrib-svg2pdfconverter
+ - python -m ipykernel install --user --name kwant-latest
+ - make -C doc latex SPHINXOPTS='-n -W -D jupyter_execute_default_kernel=kwant-latest'
- cd doc/build/latex
- make all-pdf
artifacts:
diff --git a/INSTALL.rst b/INSTALL.rst
index 336dd2db0ddbe2fa067559f79884ba910f00fbbc..db5d098791853708de117e989bbb9c80a3fdfa58 100644
--- a/INSTALL.rst
+++ b/INSTALL.rst
@@ -144,9 +144,11 @@ Building the documentation
To build the documentation, the `Sphinx documentation generator
<http://www.sphinx-doc.org/en/stable/>`_ is required with ``numpydoc`` extension
-(version 0.5 or newer). If PDF documentation is to be built, the tools
-from the `libRSVG <https://wiki.gnome.org/action/show/Projects/LibRsvg>`_ (Debian/Ubuntu package
-``librsvg2-bin``) are needed to convert SVG drawings into the PDF format.
+(version 0.5 or newer), as well as ``jupyter-sphinx`` (version 0.2 or newer).
+If PDF documentation is to be built, the tools
+from the `libRSVG <https://wiki.gnome.org/action/show/Projects/LibRsvg>`_
+(Debian/Ubuntu package ``librsvg2-bin``) are needed to convert SVG drawings
+into the PDF format.
As a prerequisite for building the documentation, Kwant must have been built
successfully using ``python3 setup.py build`` as described above (or Kwant must
@@ -160,25 +162,6 @@ Because of some quirks of how Sphinx works, it might be necessary to execute
done, Sphinx may mistakenly use PNG files for PDF output or other problems may
appear.
-When ``make html`` is run, modified tutorial example scripts are executed to
-update any figures that might have changed. The machinery behind this works as
-follows. The canonical source for a tutorial script, say ``graphene.py`` is
-the file ``doc/source/images/graphene.py.diff``. This diff file contains the
-information to recreate two versions of ``graphene.py``: a version that is
-presented in the documentation (``doc/source/tutorial/graphene.py``), and a
-version that is used to generate the figures for the documentation
-(``doc/source/images/graphene.py``). Both versions are related but differ
-e.g. in the details of the plotting. When ``make html`` is run, both versions
-are extracted form the diff file.
-
-The diff file may be modified directly. Another possible way of working is to
-directly modify either the tutorial script or the figure generation script.
-Then ``make html`` will use the command line tool `wiggle
-<https://github.com/neilbrown/wiggle>`_ to propagate the modifications accordingly.
-This will often just work, but may sometimes result in conflicts, in which case
-a message will be printed. The conflicts then have to be resolved much like
-with a version control system.
-
****************************
Hints for specific platforms
****************************
diff --git a/doc/Makefile b/doc/Makefile
index 1bdce83437737749dadea3aa5f4f9823b365e68a..a0b1e46b7fa11f3ec117e2c99d89c44c573b067b 100644
--- a/doc/Makefile
+++ b/doc/Makefile
@@ -1,183 +1,40 @@
-# Makefile for Sphinx documentation
-
-# Copyright 2011-2017 Kwant authors.
+# Minimal makefile for Sphinx documentation
#
-# This file is part of Kwant. It is subject to the license terms in the file
-# LICENSE.rst found in the top-level directory of this distribution and at
-# http://kwant-project.org/license. A list of Kwant authors can be found in
-# the file AUTHORS.rst at the top-level directory of this distribution and at
-# http://kwant-project.org/authors.
# You can set these variables from the command line.
SPHINXOPTS =
-SPHINXBUILD = python3 -c 'import sys, sphinx; sys.exit(sphinx.main(sys.argv))'
-PAPER =
+SPHINXBUILD = sphinx-build
+SOURCEDIR = source
BUILDDIR = build
-# Internal variables.
-PAPEROPT_a4 = -D latex_paper_size=a4
-PAPEROPT_letter = -D latex_paper_size=letter
-ALLSPHINXOPTS = -d $(BUILDDIR)/doctrees $(PAPEROPT_$(PAPER)) $(SPHINXOPTS) source
+# Put it first so that "make" without argument is like "make help".
+help:
+ @$(SPHINXBUILD) -M help "$(SOURCEDIR)" "$(BUILDDIR)" $(SPHINXOPTS) $(O)
+
+.PHONY: help Makefile
# We convert all SVG files to PDF for LaTeX output. For HTML output, we don't
# create PNGs but rather use the SVG files directly.
FIGURESOURCES = $(shell find source -name "*.svg")
GENERATEDPDF = $(patsubst %.svg,%.pdf,$(FIGURESOURCES))
-# Figure generation from patched tutorial scripts
-#
-# As make does not support the generation of multiple targets by a single
-# invocation of a (non-implicit) rule, we use a trick: We pretend to be
-# generating a single (empty) flag file per invocation. The figure files are
-# generated as well, but only as side-effects. Each flag file is used to
-# remember the time at which the corresponding figure-generating script was run.
-# This works perfectly unless the actual output files are deleted without
-# deleting the corresponding flag file.
-FIGSCRIPTS = $(patsubst %.diff,%,$(notdir $(wildcard source/code/figure/*.py.diff)))
-FIGURES = $(patsubst %.py,source/code/figure/.%_flag,$(FIGSCRIPTS))
-SCRIPTS = $(sort $(FIGSCRIPTS) $(notdir $(wildcard source/code/include/*.py)))
-INCLUDES = $(patsubst %,source/code/include/%,$(SCRIPTS))
-DOWNLOADS = $(patsubst %,source/code/download/%,$(SCRIPTS))
-
-.PHONY: help clean realclean html dirhtml pickle json htmlhelp qthelp latex changes linkcheck doctest
-
-help:
- @echo "Please use \`make <target>' where <target> is one of"
- @echo " html to make standalone HTML files"
- @echo " dirhtml to make HTML files named index.html in directories"
- @echo " pickle to make pickle files"
- @echo " json to make JSON files"
- @echo " htmlhelp to make HTML files and a HTML help project"
- @echo " qthelp to make HTML files and a qthelp project"
- @echo " latex to make LaTeX files, you can set PAPER=a4 or PAPER=letter"
- @echo " changes to make an overview of all changed/added/deprecated items"
- @echo " linkcheck to check all external links for integrity"
- @echo " doctest to run all doctests embedded in the documentation (if enabled)"
- @echo
- @echo "Append SPHINXOPTS='-A website_deploy=True' to include web analytics code."
-
-clean:
- -rm -rf $(BUILDDIR)/* $(GENERATEDPDF)
- -rm -rf source/reference/generated
-
-realclean: clean
- -rm -f $(FIGURES)
- -rm -f $(patsubst %,source/code/include/%,$(FIGSCRIPTS))
- -rm -f $(DOWNLOADS)
- -rm -f $(patsubst %,source/code/figure/%,$(FIGSCRIPTS))
- -rm -f $(patsubst %.py,source/code/figure/%_*.png,$(FIGSCRIPTS))
- -rm -f $(patsubst %.py,source/code/figure/%_*.pdf,$(FIGSCRIPTS))
-
-html: $(FIGURES) $(INCLUDES) $(DOWNLOADS)
- $(SPHINXBUILD) -b html $(ALLSPHINXOPTS) $(BUILDDIR)/html
- @echo
- @echo "Build finished. The HTML pages are in $(BUILDDIR)/html."
-
-dirhtml: $(FIGURES) $(INCLUDES) $(DOWNLOADS)
- $(SPHINXBUILD) -b dirhtml $(ALLSPHINXOPTS) $(BUILDDIR)/dirhtml
- @echo
- @echo "Build finished. The HTML pages are in $(BUILDDIR)/dirhtml."
-
-pickle: $(FIGURES) $(INCLUDES) $(DOWNLOADS)
- $(SPHINXBUILD) -b pickle $(ALLSPHINXOPTS) $(BUILDDIR)/pickle
- @echo
- @echo "Build finished; now you can process the pickle files."
-
-json: $(FIGURES) $(INCLUDES) $(DOWNLOADS)
- $(SPHINXBUILD) -b json $(ALLSPHINXOPTS) $(BUILDDIR)/json
- @echo
- @echo "Build finished; now you can process the JSON files."
-
-htmlhelp: $(FIGURES) $(INCLUDES) $(DOWNLOADS)
- $(SPHINXBUILD) -b htmlhelp $(ALLSPHINXOPTS) $(BUILDDIR)/htmlhelp
- @echo
- @echo "Build finished; now you can run HTML Help Workshop with the" \
- ".hhp project file in $(BUILDDIR)/htmlhelp."
-
-qthelp: $(FIGURES) $(INCLUDES) $(DOWNLOADS)
- $(SPHINXBUILD) -b qthelp $(ALLSPHINXOPTS) $(BUILDDIR)/qthelp
- @echo
- @echo "Build finished; now you can run "qcollectiongenerator" with the" \
- ".qhcp project file in $(BUILDDIR)/qthelp, like this:"
- @echo "# qcollectiongenerator $(BUILDDIR)/qthelp/kwant.qhcp"
- @echo "To view the help file:"
- @echo "# assistant -collectionFile $(BUILDDIR)/qthelp/kwant.qhc"
-
-latex: $(GENERATEDPDF) $(FIGURES) $(INCLUDES) $(DOWNLOADS)
- $(SPHINXBUILD) -b latex $(ALLSPHINXOPTS) $(BUILDDIR)/latex
- @echo
- @echo "Build finished; the LaTeX files are in $(BUILDDIR)/latex."
- @echo "Run \`make all-pdf' or \`make all-ps' in that directory to" \
- "run these through (pdf)latex."
-
-changes:
- $(SPHINXBUILD) -b changes $(ALLSPHINXOPTS) $(BUILDDIR)/changes
- @echo
- @echo "The overview file is in $(BUILDDIR)/changes."
-
-linkcheck:
- $(SPHINXBUILD) -b linkcheck $(ALLSPHINXOPTS) $(BUILDDIR)/linkcheck
- @echo
- @echo "Link check complete; look for any errors in the above output " \
- "or in $(BUILDDIR)/linkcheck/output.txt."
-
-doctest:
- $(SPHINXBUILD) -b doctest $(ALLSPHINXOPTS) $(BUILDDIR)/doctest
- @echo "Testing of doctests in the sources finished, look at the " \
- "results in $(BUILDDIR)/doctest/output.txt."
-
%.pdf: %.svg
inkscape --export-pdf=$@ $<
-#### Tutorial and figure script generation machinery ####
-# See source/code/README for an explanation.
-
-# Make tutorial scripts by extracting the (complete!) context of the "patches".
-# We make sure not to use 'wiggle' here.
-.SECONDARY:
-source/code/include/%.py: source/code/figure/%.py.diff
- @sed -n '/^[- ]/ s/^.//p' <$< >$@
- @touch -r $< $@
+# Emtpy target required so that the default target is not triggered
+%.svg:
-source/code/download/%.py: source/code/include/%.py
- @mkdir -p source/code/download
- @grep -v '^#HIDDEN' <$< >$@
-# Make the figure generation scripts by patching tutorial scripts. If the
-# tutorial scripts haven't been modified, don't patch but directly extract the
-# figure generation scripts. This means that 'wiggle' is only needed when the
-# tutorial scripts have been modified.
-.SECONDARY:
-source/code/figure/%.py: source/code/include/%.py
- @if [ $< -nt $@.diff ]; then \
- cp $< $@; \
- rm -f $@.porig; \
- if ! wiggle --replace $@ $@.diff; then \
- command -v wiggle >/dev/null 2>&1 && \
- echo "Resolve conflicts by editing the files named below"; \
- touch -d@0 $@; \
- exit 1; \
- fi \
- else \
- sed -n '/^[+ ]/ s/^.//p' <$@.diff >$@; \
- touch -r $@.diff $@; \
- fi
-
-# Make the figure generation scripts also depend on the diffs.
-define makedep
-source/code/figure/$(1): source/code/figure/$(1).diff
-endef
-$(foreach name,$(FIGSCRIPTS),$(eval $(call makedep,$(name))))
-
-# Run an figure generation script. When successful, and if the script is newer
-# than the corresponding diff, recreate the latter. Note that the
-# corresponding tutorial script cannot be newer, since if it is, the figure
-# generation script is generated from it by patching.
-source/code/figure/.%_flag: source/code/figure/%.py
- cd $(dir $<) && python3 $(notdir $<)
- @if [ ! -f $<.diff -o $< -nt $<.diff ]; then \
- wiggle --diff --lines source/code/include/$(notdir $<) $< >$<.diff; \
- touch -r $< $<.diff; \
- fi
- @rm -f $<.porig
- @touch $@
+clean:
+ rm -f $(GENERATEDPDF)
+ @$(SPHINXBUILD) -M $@ "$(SOURCEDIR)" "$(BUILDDIR)" $(SPHINXOPTS) $(O)
+
+latex: Makefile $(GENERATEDPDF)
+ cd .. ; python setup.py build ; cd -
+ @$(SPHINXBUILD) -M $@ "$(SOURCEDIR)" "$(BUILDDIR)" $(SPHINXOPTS) $(O)
+
+# Catch-all target: route all unknown targets to Sphinx using the new
+# "make mode" option. $(O) is meant as a shortcut for $(SPHINXOPTS).
+%: Makefile
+ cd .. ; python setup.py build ; cd -
+ @$(SPHINXBUILD) -M $@ "$(SOURCEDIR)" "$(BUILDDIR)" $(SPHINXOPTS) $(O)
diff --git a/doc/source/code/README b/doc/source/code/README
deleted file mode 100644
index 896e658aca4cc6ec6d2d757dab099beca7bea2d1..0000000000000000000000000000000000000000
--- a/doc/source/code/README
+++ /dev/null
@@ -1,42 +0,0 @@
-This directory contains the code examples from the documentation.
-Most scripts are present in three related but different versions that
-correspond to three different usages.
-
-* Subdirectory 'figure': scripts used for figure generation. Figures
- are not displayed but saved to disk.
-
-* Subdirectory 'include': scripts that display figures on screen.
- They contain commented marks for including snippets in the
- documentation.
-
-* Subdirectory 'download': complete scripts to be offered for download
- by readers. Like 'include' but with the include marks removed.
-
-Most scripts are extracted from corresponding '*.py.diff' files inside
-'figure/'. These are patches from the 'include' version to the
-'figure' version. The patches include complete context and as such
-can be used to recreate both files. It's these patches that are kept
-under version control.
-
-running 'make html' or 'make latex' inside '/doc' will automatically
-update all these scripts according to the following scheme:
-
- ---->------------->------
- / \
- / download/x.py \
-figure/x.py.diff ^ \
- ^ \ | \
- | -> include/x.py ---(patch)---> figure/x.py
- | | |
- | | |
- \ v /
- ----<----------(diff)--------------<--------
-
-Thus, it is possible to update figure/x.py.diff, include/x.py or
-figure/x.py and any changes will be propagated automatically when
-'make' is run. (Only download/x.py is a dead end.) The user will be
-informed about any conflicts. The makefile will only update files
-that are older than their sources and is careful to propagate time
-stamps in order to avoid infinite loops.
-
-Editing only figure/x.py.diff is a sure way to avoid any conflicts.
diff --git a/doc/source/code/figure/_defs.py b/doc/source/code/figure/_defs.py
deleted file mode 100644
index 985e2e9a3cc8f0057e53ea9726ff252975ad7343..0000000000000000000000000000000000000000
--- a/doc/source/code/figure/_defs.py
+++ /dev/null
@@ -1,34 +0,0 @@
-################################################################
-# Make matplotlib work without X11
-################################################################
-import matplotlib
-matplotlib.use('Agg')
-
-################################################################
-# Prepend Kwant's build directory to sys.path
-################################################################
-import sys
-from distutils.util import get_platform
-sys.path.insert(0, "../../../../build/lib.{0}-{1}.{2}".format(
- get_platform(), *sys.version_info[:2]))
-
-################################################################
-# Define constants for plotting
-################################################################
-pt_to_in = 1. / 72.
-
-# Default width of figures in pts
-figwidth_pt = 600
-figwidth_in = figwidth_pt * pt_to_in
-
-# Width for smaller figures
-figwidth_small_pt = 400
-figwidth_small_in = figwidth_small_pt * pt_to_in
-
-# Sizes for matplotlib figures
-mpl_width_in = figwidth_pt * pt_to_in
-mpl_label_size = None # font sizes in points
-mpl_tick_size = None
-
-# dpi for conversion from inches
-dpi = 90
diff --git a/doc/source/code/figure/ab_ring.py.diff b/doc/source/code/figure/ab_ring.py.diff
deleted file mode 100644
index c4fb178b462033ee9f0cdc906c4cfb5f349e403b..0000000000000000000000000000000000000000
--- a/doc/source/code/figure/ab_ring.py.diff
+++ /dev/null
@@ -1,202 +0,0 @@
-@@ -1,127 +1,196 @@
- # 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
-
-+import _defs
- 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 make_system_note1(a=1, t=1.0, W=10, r1=10, r2=20):
-+ lat = kwant.lattice.square(a)
-+ syst = kwant.Builder()
-+ def ring(pos):
-+ (x, y) = pos
-+ rsq = x**2 + y**2
-+ return ( r1**2 < rsq < r2**2)
-+ syst[lat.shape(ring, (0, 11))] = 4 * t
-+ syst[lat.neighbors()] = -t
-+ sym_lead0 = kwant.TranslationalSymmetry((-a, 0))
-+ lead0 = kwant.Builder(sym_lead0)
-+ def lead_shape(pos):
-+ (x, y) = pos
-+ return (-1 < x < 1) and ( 0.5 * W < y < 1.5 * W )
-+ lead0[lat.shape(lead_shape, (0, W))] = 4 * t
-+ lead0[lat.neighbors()] = -t
-+ lead1 = lead0.reversed()
-+ syst.attach_lead(lead0)
-+ syst.attach_lead(lead1)
-+ return syst
-+
-+
-+def make_system_note2(a=1, t=1.0, W=10, r1=10, r2=20):
-+ lat = kwant.lattice.square(a)
-+ syst = kwant.Builder()
-+ def ring(pos):
-+ (x, y) = pos
-+ rsq = x**2 + y**2
-+ return ( r1**2 < rsq < r2**2)
-+ syst[lat.shape(ring, (0, 11))] = 4 * t
-+ syst[lat.neighbors()] = -t
-+ sym_lead0 = kwant.TranslationalSymmetry((-a, 0))
-+ lead0 = kwant.Builder(sym_lead0)
-+ def lead_shape(pos):
-+ (x, y) = pos
-+ return (-1 < x < 1) and ( -W/2 < y < W/2 )
-+ lead0[lat.shape(lead_shape, (0, 0))] = 4 * t
-+ lead0[lat.neighbors()] = -t
-+ lead1 = lead0.reversed()
-+ syst.attach_lead(lead0)
-+ syst.attach_lead(lead1, lat(0, 0))
-+ 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, params=dict(phi=flux))
- data.append(smatrix.transmission(1, 0))
-
-- pyplot.figure()
-+ fig = pyplot.figure()
- pyplot.plot(normalized_fluxes, data)
-- pyplot.xlabel("flux [flux quantum]")
-- pyplot.ylabel("conductance [e^2/h]")
-- pyplot.show()
-+ pyplot.xlabel("flux [flux quantum]",
-+ fontsize=_defs.mpl_label_size)
-+ pyplot.ylabel("conductance [e^2/h]",
-+ fontsize=_defs.mpl_label_size)
-+ pyplot.setp(fig.get_axes()[0].get_xticklabels(),
-+ fontsize=_defs.mpl_tick_size)
-+ pyplot.setp(fig.get_axes()[0].get_yticklabels(),
-+ fontsize=_defs.mpl_tick_size)
-+ fig.set_size_inches(_defs.mpl_width_in, _defs.mpl_width_in * 3. / 4.)
-+ fig.subplots_adjust(left=0.15, right=0.95, top=0.95, bottom=0.15)
-+ fig.savefig("ab_ring_result.pdf")
-+ fig.savefig("ab_ring_result.png", dpi=_defs.dpi)
-
-
- def main():
- syst = make_system()
-
- # Check that the system looks as intended.
-- kwant.plot(syst)
-+ size = (_defs.figwidth_in, _defs.figwidth_in)
-+ for extension in ('pdf', 'png'):
-+ kwant.plot(syst, file="ab_ring_syst." + extension,
-+ fig_size=size, dpi=_defs.dpi)
-+
-
- # 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)])
-
-
-+ # Finally, some plots needed for the notes
-+ syst = make_system_note1()
-+ for extension in ('pdf', 'png'):
-+ kwant.plot(syst, file="ab_ring_note1." + extension,
-+ fig_size=size, dpi=_defs.dpi)
-+ syst = make_system_note2()
-+ for extension in ('pdf', 'png'):
-+ kwant.plot(syst, file="ab_ring_note2." + extension,
-+ fig_size=size, dpi=_defs.dpi)
-+
-+
- # 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()
diff --git a/doc/source/code/figure/band_structure.py.diff b/doc/source/code/figure/band_structure.py.diff
deleted file mode 100644
index 8a2652072aea317f4a6e5090497c37b43afdef15..0000000000000000000000000000000000000000
--- a/doc/source/code/figure/band_structure.py.diff
+++ /dev/null
@@ -1,67 +0,0 @@
-@@ -1,52 +1,62 @@
- # 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 _defs
- 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()
-+ fig = kwant.plotter.bands(lead, show=False)
-+ pyplot.xlabel("momentum [(lattice constant)^-1]",
-+ fontsize=_defs.mpl_label_size)
-+ pyplot.ylabel("energy [t]", fontsize=_defs.mpl_label_size)
-+ pyplot.setp(fig.get_axes()[0].get_xticklabels(),
-+ fontsize=_defs.mpl_tick_size)
-+ pyplot.setp(fig.get_axes()[0].get_yticklabels(),
-+ fontsize=_defs.mpl_tick_size)
-+ fig.set_size_inches(_defs.mpl_width_in, _defs.mpl_width_in * 3. / 4.)
-+ fig.subplots_adjust(left=0.15, right=0.95, top=0.95, bottom=0.15)
-+ for extension in ('pdf', 'png'):
-+ fig.savefig("band_structure_result." + extension, dpi=_defs.dpi)
-+
- #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()
diff --git a/doc/source/code/figure/closed_system.py.diff b/doc/source/code/figure/closed_system.py.diff
deleted file mode 100644
index 748e9eca48693b42d4f54faf64e988e034096101..0000000000000000000000000000000000000000
--- a/doc/source/code/figure/closed_system.py.diff
+++ /dev/null
@@ -1,168 +0,0 @@
-@@ -1,144 +1,161 @@
- # 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.
-
-+import _defs
- 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, norbs=1)
-
- 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):
- # 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):
-
- energies = []
- for B in Bfields:
- # Obtain the Hamiltonian as a sparse matrix
- ham_mat = syst.hamiltonian_submatrix(params=dict(B=B), sparse=True)
-
- # we only calculate the 15 lowest eigenvalues
- ev = sla.eigsh(ham_mat.tocsc(), k=15, sigma=0,
- return_eigenvectors=False)
-
- energies.append(ev)
-
-- pyplot.figure()
-+ fig = pyplot.figure()
- pyplot.plot(Bfields, energies)
-- pyplot.xlabel("magnetic field [arbitrary units]")
-- pyplot.ylabel("energy [t]")
-- pyplot.show()
-+ pyplot.xlabel("magnetic field [arbitrary units]",
-+ fontsize=_defs.mpl_label_size)
-+ pyplot.ylabel("energy [t]", fontsize=_defs.mpl_label_size)
-+ pyplot.setp(fig.get_axes()[0].get_xticklabels(),
-+ fontsize=_defs.mpl_tick_size)
-+ pyplot.setp(fig.get_axes()[0].get_yticklabels(),
-+ fontsize=_defs.mpl_tick_size)
-+ fig.set_size_inches(_defs.mpl_width_in, _defs.mpl_width_in * 3. / 4.)
-+ fig.subplots_adjust(left=0.15, right=0.95, top=0.95, bottom=0.15)
-+ for extension in ('pdf', 'png'):
-+ fig.savefig("closed_system_result." + extension, dpi=_defs.dpi)
- #HIDDEN_END_yvri
-
- def sorted_eigs(ev):
- evals, evecs = ev
- evals, evecs = map(np.array, zip(*sorted(zip(evals, evecs.transpose()))))
- return evals, evecs.transpose()
-
- #HIDDEN_BEGIN_wave
- def plot_wave_function(syst, B=0.001):
-+ size = (_defs.figwidth_in, _defs.figwidth_in)
-+
- # Calculate the wave functions in the system.
- ham_mat = syst.hamiltonian_submatrix(sparse=True, params=dict(B=B))
- evals, evecs = sorted_eigs(sla.eigsh(ham_mat.tocsc(), k=20, sigma=0))
-
- # Plot the probability density of the 10th eigenmode.
-- kwant.plotter.map(syst, np.abs(evecs[:, 9])**2,
-- colorbar=False, oversampling=1)
-+ for extension in ('pdf', 'png'):
-+ kwant.plotter.map(
-+ syst, np.abs(evecs[:, 9])**2, colorbar=False, oversampling=1,
-+ file="closed_system_eigenvector." + extension,
-+ fig_size=size, dpi=_defs.dpi)
- #HIDDEN_END_wave
-
-
- #HIDDEN_BEGIN_current
- def plot_current(syst, B=0.001):
-+ size = (_defs.figwidth_in, _defs.figwidth_in)
-+
- # Calculate the wave functions in the system.
- ham_mat = syst.hamiltonian_submatrix(sparse=True, params=dict(B=B))
- evals, evecs = sorted_eigs(sla.eigsh(ham_mat.tocsc(), k=20, sigma=0))
-
- # Calculate and plot the local current of the 10th eigenmode.
- J = kwant.operator.Current(syst)
- current = J(evecs[:, 9], params=dict(B=B))
-- kwant.plotter.current(syst, current, colorbar=False)
-+ for extension in ('pdf', 'png'):
-+ kwant.plotter.current(
-+ syst, current, colorbar=False,
-+ file="closed_system_current." + extension,
-+ fig_size=size, dpi=_defs.dpi)
- #HIDDEN_END_current
-
-
- 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)
- plot_current(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()
diff --git a/doc/source/code/figure/discretize.py.diff b/doc/source/code/figure/discretize.py.diff
deleted file mode 100644
index aebaf3e03c10d291d2b25789c1145b90d6210a15..0000000000000000000000000000000000000000
--- a/doc/source/code/figure/discretize.py.diff
+++ /dev/null
@@ -1,248 +0,0 @@
-@@ -1,225 +1,239 @@
- # Tutorial 2.9. Processing continuum Hamiltonians with discretize
- # ===============================================================
- #
- # Physics background
- # ------------------
- # - tight-binding approximation of continuous Hamiltonians
- #
- # Kwant features highlighted
- # --------------------------
- # - kwant.continuum.discretize
-
-+import _defs
-
- import kwant
- #HIDDEN_BEGIN_import
- import kwant.continuum
- #HIDDEN_END_import
- import scipy.sparse.linalg
- import scipy.linalg
- import numpy as np
-
- # For plotting
- import matplotlib as mpl
- from matplotlib import pyplot as plt
-
-
-+def save_figure(file_name):
-+ if not file_name:
-+ return
-+ for extension in ('pdf', 'png'):
-+ plt.savefig('.'.join((file_name,extension)),
-+ dpi=_defs.dpi, bbox_inches='tight')
-+
-+
- def stadium_system(L=20, W=20):
- #HIDDEN_BEGIN_template
- hamiltonian = "k_x**2 + k_y**2 + V(x, y)"
- template = kwant.continuum.discretize(hamiltonian)
-- print(template)
-+ with open('discretizer_verbose.txt', 'w') as f:
-+ print(template, file=f)
- #HIDDEN_END_template
-
- #HIDDEN_BEGIN_fill
- 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()
- #HIDDEN_END_fill
- return syst
-
-
- #HIDDEN_BEGIN_plot_eigenstate
- 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)
- #HIDDEN_END_plot_eigenstate
-- plt.show()
-+ save_figure('discretizer_gs')
-
-
- def qsh_system(a=20, L=2000, W=1000):
- #HIDDEN_BEGIN_define_qsh
- 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)
- """
-
- template = kwant.continuum.discretize(hamiltonian, grid=a)
- #HIDDEN_END_define_qsh
-
- #HIDDEN_BEGIN_define_qsh_build
- 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()
- #HIDDEN_END_define_qsh_build
- return syst
-
-
- def analyze_qsh():
- params = dict(A=3.65, B=-68.6, D=-51.1, M=-0.01, C=0)
- syst = qsh_system()
-
- #HIDDEN_BEGIN_plot_qsh_band
- kwant.plotter.bands(syst.leads[0], params=params,
- momenta=np.linspace(-0.3, 0.3, 201), show=False)
- #HIDDEN_END_plot_qsh_band
-
- plt.grid()
- plt.xlim(-.3, 0.3)
- plt.ylim(-0.05, 0.05)
- plt.xlabel('momentum [1/A]')
- plt.ylabel('energy [eV]')
-- plt.show()
-+ save_figure('discretizer_qsh_band')
- #HIDDEN_BEGIN_plot_qsh_wf
-+
- # 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) = plt.subplots(1, 2, sharey=True, figsize=(16, 4))
- kwant.plotter.map(syst, wf_sqr, ax=ax1)
- kwant.plotter.map(syst, rho_sz, ax=ax2)
- #HIDDEN_END_plot_qsh_wf
- 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')
-- plt.show()
-+ save_figure('discretizer_qsh_wf')
-
-
- def lattice_spacing():
- #HIDDEN_BEGIN_ls_def
- hamiltonian = "k_x**2 * identity(2) + alpha * k_x * sigma_y"
- params = dict(alpha=.5)
- #HIDDEN_END_ls_def
-
- def plot(ax, a=1):
- #HIDDEN_BEGIN_ls_hk_cont
- 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]
- #HIDDEN_END_ls_hk_cont
-
- #HIDDEN_BEGIN_ls_hk_tb
- 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]
- #HIDDEN_END_ls_hk_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()
-
- _, (ax1, ax2) = plt.subplots(1, 2, sharey=True, figsize=(12, 4))
-
- plot(ax1, a=1)
- plot(ax2, a=.25)
-- plt.show()
-+ save_figure('discretizer_lattice_spacing')
-
-
- def substitutions():
- #HIDDEN_BEGIN_subs_1
- 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])
-+ with open('discretizer_subs_1.txt', 'w') as f:
-+ print(e[0] == e[1] == e[2], file=f)
- #HIDDEN_END_subs_1
-
- #HIDDEN_BEGIN_subs_2
- subs = {'A': 'A(x) + B', 'V': 'V(x) + V_0', 'C': 5}
-- print(sympify('k_x * A * k_x + V + C', locals=subs))
-+ with open('discretizer_subs_2.txt', 'w') as f:
-+ print(sympify('k_x * A * k_x + V + C', locals=subs), file=f)
- #HIDDEN_END_subs_2
-
-
- def main():
- #HIDDEN_BEGIN_symbolic_discretization
- template = kwant.continuum.discretize('k_x * A(x) * k_x')
-- print(template)
-+ with open('discretizer_intro_verbose.txt', 'w') as f:
-+ print(template, file=f)
- #HIDDEN_END_symbolic_discretization
-
- syst = stadium_system()
- plot_eigenstate(syst)
-
- analyze_qsh()
- lattice_spacing()
- substitutions()
-
- # 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()
diff --git a/doc/source/code/figure/faq.py.diff b/doc/source/code/figure/faq.py.diff
deleted file mode 100644
index 7c8889dceec701572fde707e7c5aed3d5cf9f427..0000000000000000000000000000000000000000
--- a/doc/source/code/figure/faq.py.diff
+++ /dev/null
@@ -1,483 +0,0 @@
-@@ -1,450 +1,479 @@
- # Frequently asked questions
- #
- # This script is a disorganized collection of code snippets. As a whole, it is
- # not meant as an example of good programming practice.
-
-+import _defs
- import kwant
- import numpy as np
- import tinyarray
- import matplotlib
- from matplotlib import pyplot as plt
- matplotlib.rcParams['figure.figsize'] = (3.5, 3.5)
-
-
-+def save_figure(file_name):
-+ for extension in ('pdf', 'png'):
-+ plt.savefig('.'.join((file_name, extension)),
-+ dpi=_defs.dpi, bbox_inches='tight')
-+
-+
- #### What is a Site?
-
- #HIDDEN_BEGIN_site
- a = 1
- lat = kwant.lattice.square(a)
- syst = kwant.Builder()
-
- syst[lat(1, 0)] = 4
- syst[lat(1, 1)] = 4
-
- kwant.plot(syst)
-+save_figure("faq_site")
- #HIDDEN_END_site
-
-
- #### What is a hopping?
-
- a = 1
- lat = kwant.lattice.square(a)
- syst = kwant.Builder()
-
- syst[lat(1, 0)] = 4
- syst[lat(1, 1)] = 4
- #HIDDEN_BEGIN_hopping
- syst[(lat(1, 0), lat(1, 1))] = 1j
- #HIDDEN_END_hopping
-
- kwant.plot(syst)
-+save_figure("faq_hopping")
-
-
- #### What is a lattice?
-
- #HIDDEN_BEGIN_lattice_monatomic
- # Two monatomic lattices
- primitive_vectors = [(1, 0), (0, 1)]
- lat_a = kwant.lattice.Monatomic(primitive_vectors, offset=(0, 0))
- lat_b = kwant.lattice.Monatomic(primitive_vectors, offset=(0.5, 0.5))
- # lat1 is equivalent to kwant.lattice.square()
-
- syst = kwant.Builder()
-
- syst[lat_a(0, 0)] = 4
- syst[lat_b(0, 0)] = 4
-
- kwant.plot(syst)
-+save_figure("faq_lattice")
- #HIDDEN_END_lattice_monatomic
-
- #HIDDEN_BEGIN_lattice_polyatomic
- # One polyatomic lattice containing two sublattices
- lat = kwant.lattice.Polyatomic([(1, 0), (0, 1)], [(0, 0), (0.5, 0.5)])
- sub_a, sub_b = lat.sublattices
- #HIDDEN_END_lattice_polyatomic
-
- #### How to make a hole in a system?
-
- #HIDDEN_BEGIN_hole
- # Define the lattice and the (empty) system
- a = 2
- lat = kwant.lattice.cubic(a)
- syst = kwant.Builder()
-
- L = 10
- W = 10
- H = 2
-
- # Add sites to the system in a cuboid
-
- syst[(lat(i, j, k) for i in range(L) for j in range(W) for k in range(H))] = 4
- kwant.plot(syst)
-+save_figure("faq_hole1")
-
- # Delete sites to create a hole
-
- def in_hole(site):
- x, y, z = site.pos / a - (L/2, W/2, H/2) # position relative to centre
- return abs(x) < L / 4 and abs(y) < W / 4
-
- for site in filter(in_hole, list(syst.sites())):
- del syst[site]
-
- kwant.plot(syst)
-+save_figure("faq_hole2")
- #HIDDEN_END_hole
-
-
- #### How can we get access to the sites of our system?
-
- builder = kwant.Builder()
- lat = kwant.lattice.square()
- builder[(lat(i, j) for i in range(3) for j in range(3))] = 4
- #HIDDEN_BEGIN_sites1
- # Before finalizing the system
-
- sites = list(builder.sites()) # sites() doe *not* return a list
- #HIDDEN_END_sites1
- #HIDDEN_BEGIN_sites2
- # After finalizing the system
- syst = builder.finalized()
- sites = syst.sites # syst.sites is an actual list
- #HIDDEN_END_sites2
- #HIDDEN_BEGIN_sites3
- i = syst.id_by_site[lat(0, 2)] # we want the id of the site lat(0, 2)
- #HIDDEN_END_sites3
-
-
- #### How to plot a polyatomic lattice with different colors?
-
- #HIDDEN_BEGIN_colors1
- lat = kwant.lattice.kagome()
- syst = kwant.Builder()
-
- a, b, c = lat.sublattices # The kagome lattice has 3 sublattices
- #HIDDEN_END_colors1
-
- #HIDDEN_BEGIN_colors2
- # Plot sites from different families in different colors
- def family_color(site):
- if site.family == a:
- return 'red'
- if site.family == b:
- return 'green'
- else:
- return 'blue'
-
- def plot_system(syst):
- kwant.plot(syst, site_lw=0.1, site_color=family_color)
-
- ## Add sites and hoppings.
- for i in range(4):
- for j in range (4):
- syst[a(i, j)] = 4
- syst[b(i, j)] = 4
- syst[c(i, j)] = 4
-
- syst[lat.neighbors()] = -1
-
- ## Plot the system.
- plot_system(syst)
-+save_figure("faq_colors")
- #HIDDEN_END_colors2
-
--
- #### How to create all hoppings in a given direction using Hoppingkind?
-
- # Monatomic lattice
-
- #HIDDEN_BEGIN_direction1
- # Create hopping between neighbors with HoppingKind
- a = 1
- syst = kwant.Builder()
- lat = kwant.lattice.square(a)
- syst[ (lat(i, j) for i in range(5) for j in range(5)) ] = 4
-
- syst[kwant.builder.HoppingKind((1, 0), lat)] = -1
- kwant.plot(syst)
-+save_figure("faq_direction1")
- #HIDDEN_END_direction1
-
- # Polyatomic lattice
-
- lat = kwant.lattice.kagome()
- syst = kwant.Builder()
-
- a, b, c = lat.sublattices # The kagome lattice has 3 sublattices
-
- def family_color(site):
- if site.family == a:
- return 'red'
- if site.family == b:
- return 'green'
- else:
- return 'blue'
-
- def plot_system(syst):
- kwant.plot(syst, site_size=0.15, site_lw=0.05, site_color=family_color)
-
-
- for i in range(4):
- for j in range (4):
- syst[a(i, j)] = 4
- syst[b(i, j)] = 4
- syst[c(i, j)] = 4
-
-
- #HIDDEN_BEGIN_direction2
- # equivalent to syst[kwant.builder.HoppingKind((0, 1), b)] = -1
- syst[kwant.builder.HoppingKind((0, 1), b, b)] = -1
- #HIDDEN_END_direction2
- plot_system(syst)
-+save_figure("faq_direction2")
- # Delete the hoppings previously created
- del syst[kwant.builder.HoppingKind((0, 1), b, b)]
- #HIDDEN_BEGIN_direction3
- syst[kwant.builder.HoppingKind((0, 0), a, b)] = -1
- syst[kwant.builder.HoppingKind((0, 0), a, c)] = -1
- syst[kwant.builder.HoppingKind((0, 0), c, b)] = -1
- #HIDDEN_END_direction3
- plot_system(syst)
-+save_figure("faq_direction3")
-
-
- #### How to create the hoppings between adjacent sites?
-
- # Monatomic lattice
-
- #HIDDEN_BEGIN_adjacent1
- # Create hoppings with lat.neighbors()
- syst = kwant.Builder()
- lat = kwant.lattice.square()
- syst[(lat(i, j) for i in range(3) for j in range(3))] = 4
-
- syst[lat.neighbors()] = -1 # Equivalent to lat.neighbors(1)
- kwant.plot(syst)
-+save_figure("faq_adjacent1")
-
- del syst[lat.neighbors()] # Delete all nearest-neighbor hoppings
- syst[lat.neighbors(2)] = -1
-
- kwant.plot(syst)
-+save_figure("faq_adjacent2")
- #HIDDEN_END_adjacent1
-
- # Polyatomic lattice
-
- #HIDDEN_BEGIN_FAQ6
-
- # Hoppings using .neighbors()
- #HIDDEN_BEGIN_adjacent2
- # Create the system
- lat = kwant.lattice.kagome()
- syst = kwant.Builder()
- a, b, c = lat.sublattices # The kagome lattice has 3 sublattices
-
- for i in range(4):
- for j in range (4):
- syst[a(i, j)] = 4 # red
- syst[b(i, j)] = 4 # green
- syst[c(i, j)] = 4 # blue
-
- syst[lat.neighbors()] = -1
- #HIDDEN_END_adjacent2
- plot_system(syst)
-+save_figure("faq_adjacent3")
- del syst[lat.neighbors()] # Delete the hoppings previously created
- #HIDDEN_BEGIN_adjacent3
- syst[a.neighbors()] = -1
- #HIDDEN_END_adjacent3
- plot_system(syst)
-+save_figure("faq_adjacent4")
- del syst[a.neighbors()] # Delete the hoppings previously created
-
-
- syst[lat.neighbors(2)] = -1
- plot_system(syst)
- del syst[lat.neighbors(2)]
-
-
- #### How to create a lead with a lattice different from the scattering region?
-
- # Plot sites from different families in different colors
-
- def plot_system(syst):
-
- def family_color(site):
- if site.family == subA:
- return 'blue'
- if site.family == subB:
- return 'yellow'
- else:
- return 'green'
-
- kwant.plot(syst, site_lw=0.1, site_color=family_color)
-
-
- #HIDDEN_BEGIN_different_lattice1
- # Define the scattering Region
- L = 5
- W = 5
-
- lat = kwant.lattice.honeycomb()
- subA, subB = lat.sublattices
-
- syst = kwant.Builder()
- syst[(subA(i, j) for i in range(L) for j in range(W))] = 4
- syst[(subB(i, j) for i in range(L) for j in range(W))] = 4
- syst[lat.neighbors()] = -1
- #HIDDEN_END_different_lattice1
- plot_system(syst)
-+save_figure("faq_different_lattice1")
-
- #HIDDEN_BEGIN_different_lattice2
- # Create a lead
- lat_lead = kwant.lattice.square()
- sym_lead1 = kwant.TranslationalSymmetry((0, 1))
-
- lead1 = kwant.Builder(sym_lead1)
- lead1[(lat_lead(i, 0) for i in range(2, 7))] = 4
- lead1[lat_lead.neighbors()] = -1
- #HIDDEN_END_different_lattice2
- plot_system(lead1)
-+save_figure("faq_different_lattice2")
-
- #HIDDEN_BEGIN_different_lattice3
- syst[(lat_lead(i, 5) for i in range(2, 7))] = 4
- syst[lat_lead.neighbors()] = -1
-
- # Manually attach sites from graphene to square lattice
- syst[((lat_lead(i+2, 5), subB(i, 4)) for i in range(5))] = -1
- #HIDDEN_END_different_lattice3
- plot_system(syst)
-+save_figure("faq_different_lattice3")
-
- #HIDDEN_BEGIN_different_lattice4
- syst.attach_lead(lead1)
- #HIDDEN_END_different_lattice4
- plot_system(syst)
-+save_figure("faq_different_lattice4")
-
-
- #### How to cut a finite system out of a system with translationnal symmetries?
-
- #HIDDEN_BEGIN_fill1
- # Create 3d model.
- cubic = kwant.lattice.cubic()
- sym_3d = kwant.TranslationalSymmetry([1, 0, 0], [0, 1, 0], [0, 0, 1])
- model = kwant.Builder(sym_3d)
- model[cubic(0, 0, 0)] = 4
- model[cubic.neighbors()] = -1
- #HIDDEN_END_fill1
-
- #HIDDEN_BEGIN_fill2
- # Build scattering region (white).
- def cuboid_shape(site):
- x, y, z = abs(site.pos)
- return x < 4 and y < 10 and z < 3
-
- cuboid = kwant.Builder()
- cuboid.fill(model, cuboid_shape, (0, 0, 0));
- #HIDDEN_END_fill2
- kwant.plot(cuboid);
-+save_figure("faq_fill2")
-
- #HIDDEN_BEGIN_fill3
- # Build electrode (black).
- def electrode_shape(site):
- x, y, z = site.pos - (0, 5, 2)
- return y**2 + z**2 < 2.3**2
-
- electrode = kwant.Builder(kwant.TranslationalSymmetry([1, 0, 0]))
- electrode.fill(model, electrode_shape, (0, 5, 2)) # lead
-
- # Scattering region
- cuboid.fill(electrode, lambda s: abs(s.pos[0]) < 7, (0, 5, 4))
-
- cuboid.attach_lead(electrode)
- #HIDDEN_END_fill3
- kwant.plot(cuboid);
-+save_figure("faq_fill3")
-
-
- #### How does Kwant order the propagating modes of a lead?
-
- #HIDDEN_BEGIN_pm
- lat = kwant.lattice.square()
-
- lead = kwant.Builder(kwant.TranslationalSymmetry((-1, 0)))
- lead[(lat(0, i) for i in range(3))] = 4
- lead[lat.neighbors()] = -1
-
- flead = lead.finalized()
-
- E = 2.5
- prop_modes, _ = flead.modes(energy=E)
- #HIDDEN_END_pm
-
- def plot_and_label_modes(lead, E):
- # Plot the different modes
- pmodes, _ = lead.modes(energy=E)
- kwant.plotter.bands(lead, show=False)
- for i, k in enumerate(pmodes.momenta):
- plt.plot(k, E, 'ko')
- plt.annotate(str(i), xy=(k, E), xytext=(-5, 8),
- textcoords='offset points',
- bbox=dict(boxstyle='round,pad=0.1',fc='white', alpha=0.7))
- plt.plot([-3, 3], [E, E], 'r--')
- plt.ylim(E-1, E+1)
- plt.xlim(-2, 2)
- plt.xlabel("momentum")
- plt.ylabel("energy")
- plt.show()
-
- plot_and_label_modes(flead, E)
-+save_figure('faq_pm1')
-
- # More involved example
-
- s0 = np.eye(2)
- sz = np.array([[1, 0], [0, -1]])
-
- lead2 = kwant.Builder(kwant.TranslationalSymmetry((-1, 0)))
-
- lead2[(lat(0, i) for i in range(2))] = np.diag([1.8, -1])
- lead2[lat.neighbors()] = -1 * sz
-
- flead2 = lead2.finalized()
-
- plot_and_label_modes(flead2, 1)
-+save_figure('faq_pm2')
-
-
- #### How does Kwant order components of an individual wavefunction?
- def circle(R):
- return lambda r: np.linalg.norm(r) < R
-
-
- def make_system(lat):
- norbs = lat.norbs
- syst = kwant.Builder()
- syst[lat.shape(circle(3), (0, 0))] = 4 * np.eye(norbs)
- syst[lat.neighbors()] = -1 * np.eye(norbs)
-
- lead = kwant.Builder(kwant.TranslationalSymmetry((-1, 0)))
- lead[(lat(0, i) for i in range(-1, 2))] = 4 * np.eye(norbs)
- lead[lat.neighbors()] = -1 * np.eye(norbs)
-
- syst.attach_lead(lead)
- syst.attach_lead(lead.reversed())
-
- return syst.finalized()
-
-
- #HIDDEN_BEGIN_ord1
- lat = kwant.lattice.square(norbs=1)
- syst = make_system(lat)
- scattering_states = kwant.wave_function(syst, energy=1)
- wf = scattering_states(0)[0] # scattering state from lead 0 incoming in mode 0
-
- idx = syst.id_by_site[lat(0, 0)] # look up index of site
-
--print('wavefunction on lat(0, 0): ', wf[idx])
-+with open('faq_ord1.txt', 'w') as f:
-+ print('wavefunction on lat(0, 0): ', wf[idx], file=f)
- #HIDDEN_END_ord1
-
- #HIDDEN_BEGIN_ord2
- lat = kwant.lattice.square(norbs=2)
- syst = make_system(lat)
- scattering_states = kwant.wave_function(syst, energy=1)
- wf = scattering_states(0)[0] # scattering state from lead 0 incoming in mode 0
-
- idx = syst.id_by_site[lat(0, 0)] # look up index of site
-
- # Group consecutive degrees of freedom from 'wf' together; these correspond
- # to degrees of freedom on the same site.
- wf = wf.reshape(-1, 2)
-
--print('wavefunction on lat(0, 0): ', wf[idx])
-+with open('faq_ord2.txt', 'w') as f:
-+ print('wavefunction on lat(0, 0): ', wf[idx], file=f)
- #HIDDEN_END_ord2
diff --git a/doc/source/code/figure/graphene.py.diff b/doc/source/code/figure/graphene.py.diff
deleted file mode 100644
index 1352ed270df555a9cd1aa5c6ef51632590cd7337..0000000000000000000000000000000000000000
--- a/doc/source/code/figure/graphene.py.diff
+++ /dev/null
@@ -1,217 +0,0 @@
-@@ -1,179 +1,203 @@
- # Tutorial 2.5. Beyond square lattices: graphene
- # ==============================================
- #
- # Physics background
- # ------------------
- # Transport through a graphene quantum dot with a pn-junction
- #
- # Kwant features highlighted
- # --------------------------
- # - Application of all the aspects of tutorials 1-3 to a more complicated
- # lattice, namely graphene
-
-+import _defs
- from math import pi, sqrt, tanh
-
- import kwant
-
- # For computing eigenvalues
- import scipy.sparse.linalg as sla
-
- # For plotting
- from matplotlib import pyplot
-
-
- # Define the graphene lattice
- sin_30, cos_30 = (1 / 2, sqrt(3) / 2)
- #HIDDEN_BEGIN_hnla
- graphene = kwant.lattice.general([(1, 0), (sin_30, cos_30)],
- [(0, 0), (0, 1 / sqrt(3))])
- a, b = graphene.sublattices
- #HIDDEN_END_hnla
-
-
- #HIDDEN_BEGIN_shzy
- def make_system(r=10, w=2.0, pot=0.1):
-
- #### Define the scattering region. ####
- # circular scattering region
- def circle(pos):
- x, y = pos
- return x ** 2 + y ** 2 < r ** 2
-
- syst = kwant.Builder()
-
- # w: width and pot: potential maximum of the p-n junction
- def potential(site):
- (x, y) = site.pos
- d = y * cos_30 + x * sin_30
- return pot * tanh(d / w)
-
- syst[graphene.shape(circle, (0, 0))] = potential
- #HIDDEN_END_shzy
-
- # specify the hoppings of the graphene lattice in the
- # format expected by builder.HoppingKind
- #HIDDEN_BEGIN_hsmc
- hoppings = (((0, 0), a, b), ((0, 1), a, b), ((-1, 1), a, b))
- #HIDDEN_END_hsmc
- #HIDDEN_BEGIN_bfwb
- syst[[kwant.builder.HoppingKind(*hopping) for hopping in hoppings]] = -1
- #HIDDEN_END_bfwb
-
- # Modify the scattering region
- #HIDDEN_BEGIN_efut
- del syst[a(0, 0)]
- syst[a(-2, 1), b(2, 2)] = -1
- #HIDDEN_END_efut
-
- #### Define the leads. ####
- #HIDDEN_BEGIN_aakh
- # left lead
- sym0 = kwant.TranslationalSymmetry(graphene.vec((-1, 0)))
-
- def lead0_shape(pos):
- x, y = pos
- return (-0.4 * r < y < 0.4 * r)
-
- lead0 = kwant.Builder(sym0)
- lead0[graphene.shape(lead0_shape, (0, 0))] = -pot
- lead0[[kwant.builder.HoppingKind(*hopping) for hopping in hoppings]] = -1
-
- # The second lead, going to the top right
- sym1 = kwant.TranslationalSymmetry(graphene.vec((0, 1)))
-
- def lead1_shape(pos):
- v = pos[1] * sin_30 - pos[0] * cos_30
- return (-0.4 * r < v < 0.4 * r)
-
- lead1 = kwant.Builder(sym1)
- lead1[graphene.shape(lead1_shape, (0, 0))] = pot
- lead1[[kwant.builder.HoppingKind(*hopping) for hopping in hoppings]] = -1
- #HIDDEN_END_aakh
-
- #HIDDEN_BEGIN_kmmw
- return syst, [lead0, lead1]
- #HIDDEN_END_kmmw
-
-
- #HIDDEN_BEGIN_zydk
- def compute_evs(syst):
- # Compute some eigenvalues of the closed system
- sparse_mat = syst.hamiltonian_submatrix(sparse=True)
-
- evs = sla.eigs(sparse_mat, 2)[0]
- print(evs.real)
- #HIDDEN_END_zydk
-
-
- def plot_conductance(syst, energies):
- # Compute transmission as a function of energy
- data = []
- for energy in energies:
- smatrix = kwant.smatrix(syst, energy)
- data.append(smatrix.transmission(0, 1))
-
-- pyplot.figure()
-+ fig = pyplot.figure()
- pyplot.plot(energies, data)
-- pyplot.xlabel("energy [t]")
-- pyplot.ylabel("conductance [e^2/h]")
-- pyplot.show()
-+ pyplot.xlabel("energy [t]",
-+ fontsize=_defs.mpl_label_size)
-+ pyplot.ylabel("conductance [e^2/h]",
-+ fontsize=_defs.mpl_label_size)
-+ pyplot.setp(fig.get_axes()[0].get_xticklabels(),
-+ fontsize=_defs.mpl_tick_size)
-+ pyplot.setp(fig.get_axes()[0].get_yticklabels(),
-+ fontsize=_defs.mpl_tick_size)
-+ fig.set_size_inches(_defs.mpl_width_in, _defs.mpl_width_in * 3. / 4.)
-+ fig.subplots_adjust(left=0.15, right=0.95, top=0.95, bottom=0.15)
-+ for extension in ('pdf', 'png'):
-+ fig.savefig("graphene_result." + extension, dpi=_defs.dpi)
-
-
- def plot_bandstructure(flead, momenta):
- bands = kwant.physics.Bands(flead)
- energies = [bands(k) for k in momenta]
-
-- pyplot.figure()
-+ fig = pyplot.figure()
- pyplot.plot(momenta, energies)
-- pyplot.xlabel("momentum [(lattice constant)^-1]")
-- pyplot.ylabel("energy [t]")
-- pyplot.show()
-+ pyplot.xlabel("momentum [(lattice constant)^-1]",
-+ fontsize=_defs.mpl_label_size)
-+ pyplot.ylabel("energy [t]",
-+ fontsize=_defs.mpl_label_size)
-+ pyplot.setp(fig.get_axes()[0].get_xticklabels(),
-+ fontsize=_defs.mpl_tick_size)
-+ pyplot.setp(fig.get_axes()[0].get_yticklabels(),
-+ fontsize=_defs.mpl_tick_size)
-+ fig.set_size_inches(_defs.mpl_width_in, _defs.mpl_width_in * 3. / 4.)
-+ fig.subplots_adjust(left=0.15, right=0.95, top=0.95, bottom=0.15)
-+ for extension in ('pdf', 'png'):
-+ fig.savefig("graphene_bs." + extension, dpi=_defs.dpi)
-
-
- #HIDDEN The part of the following code block which begins with family_colors
- #HIDDEN is included verbatim in the tutorial text because nested code examples
- #HIDDEN are not supported. Remember to update the tutorial text when you
- #HIDDEN modify this block.
- #HIDDEN_BEGIN_itkk
- def main():
- pot = 0.1
- syst, leads = make_system(pot=pot)
-
- # To highlight the two sublattices of graphene, we plot one with
- # a filled, and the other one with an open circle:
- def family_colors(site):
- return 0 if site.family == a else 1
-
-- # Plot the closed system without leads.
-- kwant.plot(syst, site_color=family_colors, site_lw=0.1, colorbar=False)
-+ size = (_defs.figwidth_in, _defs.figwidth_in)
-+ for extension in ('pdf', 'png'):
-+ kwant.plot(syst, site_color=family_colors, site_lw=0.1, colorbar=False,
-+ file="graphene_syst1." + extension,
-+ fig_size=size, dpi=_defs.dpi)
- #HIDDEN_END_itkk
-
- # Compute some eigenvalues.
- #HIDDEN_BEGIN_jmbi
- compute_evs(syst.finalized())
- #HIDDEN_END_jmbi
-
- # Attach the leads to the system.
- for lead in leads:
- syst.attach_lead(lead)
-
-- # Then, plot the system with leads.
-- kwant.plot(syst, site_color=family_colors, site_lw=0.1,
-- lead_site_lw=0, colorbar=False)
-+ size = (_defs.figwidth_in, 0.9 * _defs.figwidth_in)
-+ for extension in ('pdf', 'png'):
-+ kwant.plot(syst, site_color=family_colors, colorbar=False, site_lw=0.1,
-+ file="graphene_syst2." + extension,
-+ fig_size=size, dpi=_defs.dpi, lead_site_lw=0)
-
- # Finalize the system.
- syst = syst.finalized()
-
- # Compute the band structure of lead 0.
- momenta = [-pi + 0.02 * pi * i for i in range(101)]
- plot_bandstructure(syst.leads[0], momenta)
-
- # Plot conductance.
- energies = [-2 * pot + 4. / 50. * pot * i for i in range(51)]
- plot_conductance(syst, energies)
-
-
- # 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()
diff --git a/doc/source/code/figure/kernel_polynomial_method.py.diff b/doc/source/code/figure/kernel_polynomial_method.py.diff
deleted file mode 100644
index 3ac532512ea1d3add0d5604c897aa0e9ceb49399..0000000000000000000000000000000000000000
--- a/doc/source/code/figure/kernel_polynomial_method.py.diff
+++ /dev/null
@@ -1,374 +0,0 @@
-@@ -1,219 +1,242 @@
- # Tutorial 2.8. Calculating spectral density with the Kernel Polynomial Method
- # ============================================================================
- #
- # Physics background
- # ------------------
- # - Chebyshev polynomials, random trace approximation, spectral densities.
- #
- # Kwant features highlighted
- # --------------------------
- # - kpm module,kwant operators.
-
- import scipy
-
-+import _defs
-+from contextlib import redirect_stdout
-+
- # For plotting
- from matplotlib import pyplot as plt
-
- #HIDDEN_BEGIN_sys1
- # necessary imports
- import kwant
- import numpy as np
-
-
- # define the system
- def make_syst(r=30, t=-1, a=1):
- syst = kwant.Builder()
- lat = kwant.lattice.honeycomb(a, norbs=1)
-
- def circle(pos):
- x, y = pos
- return x ** 2 + y ** 2 < r ** 2
-
- syst[lat.shape(circle, (0, 0))] = 0.
- syst[lat.neighbors()] = t
- syst.eradicate_dangling()
-
- return syst
- #HIDDEN_END_sys1
-
- #HIDDEN_BEGIN_sys2
- # define a Haldane system
- def make_syst_topo(r=30, a=1, t=1, t2=0.5):
- syst = kwant.Builder()
- lat = kwant.lattice.honeycomb(a, norbs=1, name=['a', 'b'])
-
- def circle(pos):
- x, y = pos
- return x ** 2 + y ** 2 < r ** 2
-
- syst[lat.shape(circle, (0, 0))] = 0.
- syst[lat.neighbors()] = t
- # add second neighbours hoppings
- syst[lat.a.neighbors()] = 1j * t2
- syst[lat.b.neighbors()] = -1j * t2
- syst.eradicate_dangling()
-
- return lat, syst.finalized()
- #HIDDEN_END_sys2
-
-
- #HIDDEN_BEGIN_sys3
- # define the system
- def make_syst_staggered(r=30, t=-1, a=1, m=0.1):
- syst = kwant.Builder()
- lat = kwant.lattice.honeycomb(a, norbs=1)
-
- def circle(pos):
- x, y = pos
- return x ** 2 + y ** 2 < r ** 2
-
- syst[lat.a.shape(circle, (0, 0))] = m
- syst[lat.b.shape(circle, (0, 0))] = -m
- syst[lat.neighbors()] = t
- syst.eradicate_dangling()
-
- return syst
- #HIDDEN_END_sys3
-
- # Plot several density of states curves on the same axes.
--def plot_dos(labels_to_data):
-+def plot_dos(labels_to_data, file_name=None, ylabel="DoS [a.u.]"):
- plt.figure(figsize=(5,4))
- for label, (x, y) in labels_to_data:
- plt.plot(x, y.real, label=label, linewidth=2)
- plt.legend(loc=2, framealpha=0.5)
- plt.xlabel("energy [t]")
- plt.ylabel(ylabel)
-- plt.show()
-+ save_figure(file_name)
- plt.clf()
-
-
- # Plot fill density of states plus curves on the same axes.
--def plot_dos_and_curves(dos labels_to_data):
-+def plot_dos_and_curves(dos, labels_to_data, file_name=None, ylabel="DoS [a.u.]"):
- plt.figure(figsize=(5,4))
- plt.fill_between(dos[0], dos[1], label="DoS [a.u.]",
- alpha=0.5, color='gray')
- for label, (x, y) in labels_to_data:
- plt.plot(x, y, label=label, linewidth=2)
- plt.legend(loc=2, framealpha=0.5)
- plt.xlabel("energy [t]")
- plt.ylabel(ylabel)
-- plt.show()
-+ save_figure(file_name)
- plt.clf()
-
- def site_size_conversion(densities):
- return 3 * np.abs(densities) / max(densities)
-
-
- # Plot several local density of states maps in different subplots
- def plot_ldos(syst, densities, file_name=None):
- fig, axes = plt.subplots(1, len(densities), figsize=(7*len(densities), 7))
- for ax, (title, rho) in zip(axes, densities):
- kwant.plotter.density(syst, rho.real, ax=ax)
- ax.set_title(title)
- ax.set(adjustable='box', aspect='equal')
-- plt.show()
-+ save_figure(file_name)
- plt.clf()
-
-+def save_figure(file_name):
-+ if not file_name:
-+ return
-+ for extension in ('pdf', 'png'):
-+ plt.savefig('.'.join((file_name,extension)),
-+ dpi=_defs.dpi, bbox_inches='tight')
-+
-+
- def simple_dos_example():
- #HIDDEN_BEGIN_kpm1
- fsyst = make_syst().finalized()
-
- spectrum = kwant.kpm.SpectralDensity(fsyst)
- #HIDDEN_END_kpm1
-
- #HIDDEN_BEGIN_kpm2
- energies, densities = spectrum()
- #HIDDEN_END_kpm2
-
- #HIDDEN_BEGIN_kpm3
- energy_subset = np.linspace(0, 2)
- density_subset = spectrum(energy_subset)
- #HIDDEN_END_kpm3
-
- plot_dos([
- ('densities', (energies, densities)),
- ('density subset', (energy_subset, density_subset)),
-- ])
-+ ],
-+ file_name='kpm_dos'
-+ )
-
-
- def dos_integrating_example(fsyst):
- spectrum = kwant.kpm.SpectralDensity(fsyst)
-
- #HIDDEN_BEGIN_int1
-- print('identity resolution:', spectrum.integrate())
-+ with open('kpm_normalization.txt', 'w') as f:
-+ with redirect_stdout(f):
-+ print('identity resolution:', spectrum.integrate())
- #HIDDEN_END_int1
-
- #HIDDEN_BEGIN_int2
- # Fermi energy 0.1 and temperature 0.2
- fermi = lambda E: 1 / (np.exp((E - 0.1) / 0.2) + 1)
-
-- print('number of filled states:', spectrum.integrate(fermi))
-+ with open('kpm_total_states.txt', 'w') as f:
-+ with redirect_stdout(f):
-+ print('number of filled states:', spectrum.integrate(fermi))
- #HIDDEN_END_int2
-
-
- def increasing_accuracy_example(fsyst):
- spectrum = kwant.kpm.SpectralDensity(fsyst)
- original_dos = spectrum() # get unaltered DoS
-
- #HIDDEN_BEGIN_acc1
- spectrum.add_moments(energy_resolution=0.03)
- #HIDDEN_END_acc1
-
- increased_resolution_dos = spectrum()
-
- plot_dos([
- ('density', original_dos),
- ('higher energy resolution', increased_resolution_dos),
-- ])
-+ ],
-+ file_name='kpm_dos_acc'
-+ )
-
- #HIDDEN_BEGIN_acc2
- spectrum.add_moments(100)
- spectrum.add_vectors(5)
- #HIDDEN_END_acc2
-
- increased_moments_dos = spectrum()
-
- plot_dos([
- ('density', original_dos),
- ('higher number of moments', increased_moments_dos),
-- ])
-+ ],
-+ file_name='kpm_dos_r'
-+ )
-
-
- def operator_example(fsyst):
- #HIDDEN_BEGIN_op1
- # identity matrix
- matrix_op = scipy.sparse.eye(len(fsyst.sites))
- matrix_spectrum = kwant.kpm.SpectralDensity(fsyst, operator=matrix_op)
- #HIDDEN_END_op1
-
- #HIDDEN_BEGIN_op2
- # 'sum=True' means we sum over all the sites
- kwant_op = kwant.operator.Density(fsyst, sum=True)
- operator_spectrum = kwant.kpm.SpectralDensity(fsyst, operator=kwant_op)
- #HIDDEN_END_op2
-
- plot_dos([
- ('identity matrix', matrix_spectrum()),
- ('kwant.operator.Density', operator_spectrum()),
- ])
-
-
- def ldos_example(fsyst):
- #HIDDEN_BEGIN_op3
- # 'sum=False' is the default, but we include it explicitly here for clarity.
- kwant_op = kwant.operator.Density(fsyst, sum=False)
- local_dos = kwant.kpm.SpectralDensity(fsyst, operator=kwant_op)
- #HIDDEN_END_op3
-
- #HIDDEN_BEGIN_op4
- zero_energy_ldos = local_dos(energy=0)
- finite_energy_ldos = local_dos(energy=1)
- plot_ldos(fsyst, [
- ('energy = 0', zero_energy_ldos),
- ('energy = 1', finite_energy_ldos)
-- ])
-+ ],
-+ file_name='kpm_ldos'
-+ )
- #HIDDEN_END_op4
-
-
- def ldos_sites_example():
- fsyst = make_syst_staggered().finalized()
- #HIDDEN_BEGIN_op5
- # find 'A' and 'B' sites in the unit cell at the center of the disk
- center_tag = np.array([0, 0])
- where = lambda s: s.tag == center_tag
- # make local vectors
- vector_factory = kwant.kpm.LocalVectors(fsyst, where)
- #HIDDEN_END_op5
-
- #HIDDEN_BEGIN_op6
- # 'num_vectors' can be unspecified when using 'LocalVectors'
- local_dos = kwant.kpm.SpectralDensity(fsyst, num_vectors=None,
- vector_factory=vector_factory,
- mean=False)
- energies, densities = local_dos()
- plot_dos([
- ('A sublattice', (energies, densities[:, 0])),
- ('B sublattice', (energies, densities[:, 1])),
-- ])
-+ ],
-+ file_name='kpm_ldos_sites'
-+ )
- #HIDDEN_END_op6
-
-
- def vector_factory_example(fsyst):
- spectrum = kwant.kpm.SpectralDensity(fsyst)
- #HIDDEN_BEGIN_fact1
- # construct a generator of vectors with n random elements -1 or +1.
- n = fsyst.hamiltonian_submatrix(sparse=True).shape[0]
- def binary_vectors():
- while True:
- yield np.rint(np.random.random_sample(n)) * 2 - 1
-
- custom_factory = kwant.kpm.SpectralDensity(fsyst,
- vector_factory=binary_vectors())
- #HIDDEN_END_fact1
- plot_dos([
- ('default vector factory', spectrum()),
- ('binary vector factory', custom_factory()),
- ])
-
-
- def bilinear_map_operator_example(fsyst):
- #HIDDEN_BEGIN_blm
- rho = kwant.operator.Density(fsyst, sum=True)
-
- # sesquilinear map that does the same thing as `rho`
- def rho_alt(bra, ket):
- return np.vdot(bra, ket)
-
- rho_spectrum = kwant.kpm.SpectralDensity(fsyst, operator=rho)
- rho_alt_spectrum = kwant.kpm.SpectralDensity(fsyst, operator=rho_alt)
- #HIDDEN_END_blm
-
- plot_dos([
- ('kwant.operator.Density', rho_spectrum()),
- ('bilinear operator', rho_alt_spectrum()),
- ])
-
- def conductivity_example():
- #HIDDEN_BEGIN_cond
- # construct the Haldane model
- lat, fsyst = make_syst_topo()
- # find 'A' and 'B' sites in the unit cell at the center of the disk
- where = lambda s: np.linalg.norm(s.pos) < 1
-
- # component 'xx'
- s_factory = kwant.kpm.LocalVectors(fsyst, where)
- cond_xx = kwant.kpm.conductivity(fsyst, alpha='x', beta='x', mean=True,
- num_vectors=None, vector_factory=s_factory)
- # component 'xy'
- s_factory = kwant.kpm.LocalVectors(fsyst, where)
- cond_xy = kwant.kpm.conductivity(fsyst, alpha='x', beta='y', mean=True,
- num_vectors=None, vector_factory=s_factory)
-
- energies = cond_xx.energies
- cond_array_xx = np.array([cond_xx(e, temperature=0.01) for e in energies])
- cond_array_xy = np.array([cond_xy(e, temperature=0.01) for e in energies])
-
- # area of the unit cell per site
- area_per_site = np.abs(np.cross(*lat.prim_vecs)) / len(lat.sublattices)
- cond_array_xx /= area_per_site
- cond_array_xy /= area_per_site
- #HIDDEN_END_cond
- # ldos
- s_factory = kwant.kpm.LocalVectors(fsyst, where)
- spectrum = kwant.kpm.SpectralDensity(fsyst, num_vectors=None,
- vector_factory=s_factory)
-
- plot_dos_and_curves(
- (spectrum.energies, spectrum.densities * 8),
- [
- (r'Longitudinal conductivity $\sigma_{xx} / 4$',
- (energies, cond_array_xx / 4)),
- (r'Hall conductivity $\sigma_{xy}$',
- (energies, cond_array_xy))],
- ylabel=r'$\sigma [e^2/h]$',
- file_name='kpm_cond'
- )
-
-
- def main():
- simple_dos_example()
-
- fsyst = make_syst().finalized()
-
- dos_integrating_example(fsyst)
- increasing_accuracy_example(fsyst)
- operator_example(fsyst)
- ldos_example(fsyst)
- ldos_sites_example()
- vector_factory_example(fsyst)
- bilinear_map_operator_example(fsyst)
- conductivity_example()
-
-
- # 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()
diff --git a/doc/source/code/figure/magnetic_texture.py.diff b/doc/source/code/figure/magnetic_texture.py.diff
deleted file mode 100644
index a29a458e43ea5314b09e93426bf907872228f8f9..0000000000000000000000000000000000000000
--- a/doc/source/code/figure/magnetic_texture.py.diff
+++ /dev/null
@@ -1,286 +0,0 @@
-@@ -1,245 +1,268 @@
- # Tutorial 2.7. Spin textures
- # ===========================
- #
- # Physics background
- # ------------------
- # - Spin textures
- # - Skyrmions
- #
- # Kwant features highlighted
- # --------------------------
- # - operators
- # - plotting vector fields
-
-+import _defs
-+from contextlib import redirect_stdout
-+
- from math import sin, cos, tanh, pi
- import itertools
- import numpy as np
- import tinyarray as ta
- import matplotlib.pyplot as plt
-
- import kwant
-
- sigma_0 = ta.array([[1, 0], [0, 1]])
- sigma_x = ta.array([[0, 1], [1, 0]])
- sigma_y = ta.array([[0, -1j], [1j, 0]])
- sigma_z = ta.array([[1, 0], [0, -1]])
-
- # vector of Pauli matrices σ_αiβ where greek
- # letters denote spinor indices
- sigma = np.rollaxis(np.array([sigma_x, sigma_y, sigma_z]), 1)
-
- #HIDDEN_BEGIN_model
- def field_direction(pos, r0, delta):
- x, y = pos
- r = np.linalg.norm(pos)
- r_tilde = (r - r0) / delta
- theta = (tanh(r_tilde) - 1) * (pi / 2)
-
- if r == 0:
- m_i = [0, 0, -1]
- else:
- m_i = [
- (x / r) * sin(theta),
- (y / r) * sin(theta),
- cos(theta),
- ]
-
- return np.array(m_i)
-
-
- def scattering_onsite(site, r0, delta, J):
- m_i = field_direction(site.pos, r0, delta)
- return J * np.dot(m_i, sigma)
-
-
- def lead_onsite(site, J):
- return J * sigma_z
- #HIDDEN_END_model
-
-
- #HIDDEN_BEGIN_syst
- lat = kwant.lattice.square(norbs=2)
-
- def make_system(L=80):
-
- syst = kwant.Builder()
-
- def square(pos):
- return all(-L/2 < p < L/2 for p in pos)
-
- syst[lat.shape(square, (0, 0))] = scattering_onsite
- syst[lat.neighbors()] = -sigma_0
-
- lead = kwant.Builder(kwant.TranslationalSymmetry((-1, 0)),
- conservation_law=-sigma_z)
-
- lead[lat.shape(square, (0, 0))] = lead_onsite
- lead[lat.neighbors()] = -sigma_0
-
- syst.attach_lead(lead)
- syst.attach_lead(lead.reversed())
-
- return syst
- #HIDDEN_END_syst
-
-
--def plot_vector_field(syst, params):
-+def save_plot(fname):
-+ for extension in ('.pdf', '.png'):
-+ plt.savefig(fname + extension,
-+ dpi=_defs.dpi, bbox_inches='tight')
-+
-+
-+def plot_vector_field(syst, params, fname=None):
- xmin, ymin = min(s.tag for s in syst.sites)
- xmax, ymax = max(s.tag for s in syst.sites)
- x, y = np.meshgrid(np.arange(xmin, xmax+1), np.arange(ymin, ymax+1))
-
- m_i = [field_direction(p, **params) for p in zip(x.flat, y.flat)]
- m_i = np.reshape(m_i, x.shape + (3,))
- m_i = np.rollaxis(m_i, 2, 0)
-
-- fig, ax = plt.subplots(1, 1)
-+ fig, ax = plt.subplots(1, 1, figsize=(9, 7))
- im = ax.quiver(x, y, *m_i, pivot='mid', scale=75)
- fig.colorbar(im)
-- plt.show()
-+ if fname:
-+ save_plot(fname)
-+ else:
-+ plt.show()
-
-
--def plot_densities(syst, densities):
-- fig, axes = plt.subplots(1, len(densities))
-+def plot_densities(syst, densities, fname=None):
-+ fig, axes = plt.subplots(1, len(densities), figsize=(7*len(densities), 7))
- for ax, (title, rho) in zip(axes, densities):
- kwant.plotter.map(syst, rho, ax=ax, a=4)
- ax.set_title(title)
-- plt.show()
-+ if fname:
-+ save_plot(fname)
-+ else:
-+ plt.show()
-
-
--def plot_currents(syst, currents):
-- fig, axes = plt.subplots(1, len(currents))
-+
-+def plot_currents(syst, currents, fname=None):
-+ fig, axes = plt.subplots(1, len(currents), figsize=(7*len(currents), 7))
- if not hasattr(axes, '__len__'):
- axes = (axes,)
- for ax, (title, current) in zip(axes, currents):
- kwant.plotter.current(syst, current, ax=ax, colorbar=False)
- ax.set_title(title)
-- plt.show()
-+ if fname:
-+ save_plot(fname)
-+ else:
-+ plt.show()
-
-
- def main():
- syst = make_system().finalized()
-
- #HIDDEN_BEGIN_wavefunction
- params = dict(r0=20, delta=10, J=1)
- wf = kwant.wave_function(syst, energy=-1, params=params)
- psi = wf(0)[0]
- #HIDDEN_END_wavefunction
-
-- plot_vector_field(syst, dict(r0=20, delta=10))
-+ plot_vector_field(syst, dict(r0=20, delta=10), fname='mag_field_direction')
-
- #HIDDEN_BEGIN_ldos
- # even (odd) indices correspond to spin up (down)
- up, down = psi[::2], psi[1::2]
- density = np.abs(up)**2 + np.abs(down)**2
- #HIDDEN_END_ldos
-
- #HIDDEN_BEGIN_lsdz
- # spin down components have a minus sign
- spin_z = np.abs(up)**2 - np.abs(down)**2
- #HIDDEN_END_lsdz
-
- #HIDDEN_BEGIN_lsdy
- # spin down components have a minus sign
- spin_y = 1j * (down.conjugate() * up - up.conjugate() * down)
- #HIDDEN_END_lsdy
-
- #HIDDEN_BEGIN_lden
- rho = kwant.operator.Density(syst)
- rho_sz = kwant.operator.Density(syst, sigma_z)
- rho_sy = kwant.operator.Density(syst, sigma_y)
-
- # calculate the expectation values of the operators with 'psi'
- density = rho(psi)
- spin_z = rho_sz(psi)
- spin_y = rho_sy(psi)
- #HIDDEN_END_lden
-
- plot_densities(syst, [
- ('$σ_0$', density),
- ('$σ_z$', spin_z),
- ('$σ_y$', spin_y),
-- ])
-+ ], fname='spin_densities')
-
- #HIDDEN_BEGIN_current
- J_0 = kwant.operator.Current(syst)
- J_z = kwant.operator.Current(syst, sigma_z)
- J_y = kwant.operator.Current(syst, sigma_y)
-
- # calculate the expectation values of the operators with 'psi'
- current = J_0(psi)
- spin_z_current = J_z(psi)
- spin_y_current = J_y(psi)
- #HIDDEN_END_current
-
- plot_currents(syst, [
- ('$J_{σ_0}$', current),
- ('$J_{σ_z}$', spin_z_current),
- ('$J_{σ_y}$', spin_y_current),
-- ])
-+ ], fname='spin_currents')
-
- #HIDDEN_BEGIN_following
- def following_m_i(site, r0, delta):
- m_i = field_direction(site.pos, r0, delta)
- return np.dot(m_i, sigma)
-
- J_m = kwant.operator.Current(syst, following_m_i)
-
- # evaluate the operator
- m_current = J_m(psi, params=dict(r0=25, delta=10))
- #HIDDEN_END_following
-
- plot_currents(syst, [
- (r'$J_{\mathbf{m}_i}$', m_current),
- ('$J_{σ_z}$', spin_z_current),
-- ])
-+ ], fname='spin_current_comparison')
-
-
- #HIDDEN_BEGIN_density_cut
- def circle(site):
- return np.linalg.norm(site.pos) < 20
-
- rho_circle = kwant.operator.Density(syst, where=circle, sum=True)
-
- all_states = np.vstack((wf(0), wf(1)))
- dos_in_circle = sum(rho_circle(p) for p in all_states) / (2 * pi)
-- print('density of states in circle:', dos_in_circle)
-+ with open('circle_dos.txt', 'w') as f:
-+ with redirect_stdout(f):
-+ print('density of states in circle:', dos_in_circle)
- #HIDDEN_END_density_cut
-
- #HIDDEN_BEGIN_current_cut
- def left_cut(site_to, site_from):
- return site_from.pos[0] <= -39 and site_to.pos[0] > -39
-
- def right_cut(site_to, site_from):
- return site_from.pos[0] < 39 and site_to.pos[0] >= 39
-
- J_left = kwant.operator.Current(syst, where=left_cut, sum=True)
- J_right = kwant.operator.Current(syst, where=right_cut, sum=True)
-
- Jz_left = kwant.operator.Current(syst, sigma_z, where=left_cut, sum=True)
- Jz_right = kwant.operator.Current(syst, sigma_z, where=right_cut, sum=True)
-
-- print('J_left:', J_left(psi), ' J_right:', J_right(psi))
-- print('Jz_left:', Jz_left(psi), ' Jz_right:', Jz_right(psi))
-+ with open('current_cut.txt', 'w') as f:
-+ with redirect_stdout(f):
-+ print('J_left:', J_left(psi), ' J_right:', J_right(psi))
-+ print('Jz_left:', Jz_left(psi), ' Jz_right:', Jz_right(psi))
- #HIDDEN_END_current_cut
-
- #HIDDEN_BEGIN_bind
- J_m = kwant.operator.Current(syst, following_m_i)
- J_z = kwant.operator.Current(syst, sigma_z)
-
- J_m_bound = J_m.bind(params=dict(r0=25, delta=10, J=1))
- J_z_bound = J_z.bind(params=dict(r0=25, delta=10, J=1))
-
- # Sum current local from all scattering states on the left at energy=-1
- wf_left = wf(0)
- J_m_left = sum(J_m_bound(p) for p in wf_left)
- J_z_left = sum(J_z_bound(p) for p in wf_left)
- #HIDDEN_END_bind
-
- plot_currents(syst, [
- (r'$J_{\mathbf{m}_i}$ (from left)', J_m_left),
- (r'$J_{σ_z}$ (from left)', J_z_left),
-- ])
-+ ], fname='bound_current')
-
-
- if __name__ == '__main__':
- main()
diff --git a/doc/source/code/figure/plot_graphene.py.diff b/doc/source/code/figure/plot_graphene.py.diff
deleted file mode 100644
index 82c7555424f6a0206c68cbce0235da3e6fa7ef73..0000000000000000000000000000000000000000
--- a/doc/source/code/figure/plot_graphene.py.diff
+++ /dev/null
@@ -1,141 +0,0 @@
-@@ -1,112 +1,130 @@
- # Tutorial 2.8.1. 2D example: graphene quantum dot
- # ================================================
- #
- # Physics background
- # ------------------
- # - graphene edge states
- #
- # Kwant features highlighted
- # --------------------------
- # - demonstrate different ways of plotting
-
-+import _defs
- import kwant
- from matplotlib import pyplot
-
- #HIDDEN_BEGIN_makesyst
- lat = kwant.lattice.honeycomb()
- a, b = lat.sublattices
-
- def make_system(r=8, t=-1, tp=-0.1):
-
- def circle(pos):
- x, y = pos
- return x**2 + y**2 < r**2
-
- syst = kwant.Builder()
-- syst[lat.shape(circle, (0, 0))] = 0
-+ syst[lat.shape(circle, (0,0))] = 0
- syst[lat.neighbors()] = t
- syst.eradicate_dangling()
- if tp:
- syst[lat.neighbors(2)] = tp
-
- return syst
- #HIDDEN_END_makesyst
-
-
- #HIDDEN_BEGIN_plotsyst1
- def plot_system(syst):
-- kwant.plot(syst)
- #HIDDEN_END_plotsyst1
-- # the standard plot is ok, but not very intelligible. One can do
-- # better by playing wioth colors and linewidths
-+ # standard plot - not very intelligible for this particular situation
-+ size = (_defs.figwidth_in, _defs.figwidth_in)
-+ for extension in ('pdf', 'png'):
-+ kwant.plot(syst, file="plot_graphene_syst1." + extension,
-+ fig_size=size, dpi=_defs.dpi)
-
- # use color and linewidths to get a better plot
- #HIDDEN_BEGIN_plotsyst2
- def family_color(site):
- return 'black' if site.family == a else 'white'
-
- def hopping_lw(site1, site2):
- return 0.04 if site1.family == site2.family else 0.1
-
-- kwant.plot(syst, site_lw=0.1, site_color=family_color, hop_lw=hopping_lw)
-+ size = (_defs.figwidth_in, _defs.figwidth_in)
-+ for extension in ('pdf', 'png'):
-+ kwant.plot(syst, site_lw=0.1, site_color=family_color,
-+ hop_lw=hopping_lw, file="plot_graphene_syst2." + extension,
-+ fig_size=size, dpi=_defs.dpi)
- #HIDDEN_END_plotsyst2
-
-
- #HIDDEN_BEGIN_plotdata1
- def plot_data(syst, n):
- import scipy.linalg as la
-
- syst = syst.finalized()
- ham = syst.hamiltonian_submatrix()
- evecs = la.eigh(ham)[1]
-
- wf = abs(evecs[:, n])**2
- #HIDDEN_END_plotdata1
-
- # the usual - works great in general, looks just a bit crufty for
- # small systems
-
- #HIDDEN_BEGIN_plotdata2
-- kwant.plotter.map(syst, wf, oversampling=10, cmap='gist_heat_r')
-+ size = (_defs.figwidth_in, _defs.figwidth_in)
-+ for extension in ('pdf', 'png'):
-+ kwant.plotter.map(syst, wf, oversampling=10, cmap='gist_heat_r',
-+ file="plot_graphene_data1." + extension,
-+ fig_size=size, dpi=_defs.dpi)
- #HIDDEN_END_plotdata2
-
- # use two different sort of triangles to cleanly fill the space
- #HIDDEN_BEGIN_plotdata3
- def family_shape(i):
- site = syst.sites[i]
- return ('p', 3, 180) if site.family == a else ('p', 3, 0)
-
- def family_color(i):
- return 'black' if syst.sites[i].family == a else 'white'
-
-- kwant.plot(syst, site_color=wf, site_symbol=family_shape,
-- site_size=0.5, hop_lw=0, cmap='gist_heat_r')
-+ size = (_defs.figwidth_in, _defs.figwidth_in)
-+ for extension in ('pdf', 'png'):
-+ kwant.plot(syst, site_color=wf, site_symbol=family_shape,
-+ site_size=0.5, hop_lw=0, cmap='gist_heat_r',
-+ file="plot_graphene_data2." + extension,
-+ fig_size=size, dpi=_defs.dpi)
- #HIDDEN_END_plotdata3
-
- # plot by changing the symbols itself
- #HIDDEN_BEGIN_plotdata4
- def site_size(i):
- return 3 * wf[i] / wf.max()
-
-- kwant.plot(syst, site_size=site_size, site_color=(0, 0, 1, 0.3),
-- hop_lw=0.1)
-+ size = (_defs.figwidth_in, _defs.figwidth_in)
-+ for extension in ('pdf', 'png'):
-+ kwant.plot(syst, site_size=site_size, site_color=(0,0,1,0.3),
-+ hop_lw=0.1, file="plot_graphene_data3." + extension,
-+ fig_size=size, dpi=_defs.dpi)
- #HIDDEN_END_plotdata4
-
-
- def main():
- # plot the graphene system in different styles
- syst = make_system()
-
- plot_system(syst)
-
- # compute a wavefunction (number 225) and plot it in different
- # styles
- syst = make_system(tp=0)
-
- plot_data(syst, 225)
-
-
- # 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()
diff --git a/doc/source/code/figure/plot_qahe.py.diff b/doc/source/code/figure/plot_qahe.py.diff
deleted file mode 100644
index 1788f677df857d1b75b8706aa0907fad4b9a467d..0000000000000000000000000000000000000000
--- a/doc/source/code/figure/plot_qahe.py.diff
+++ /dev/null
@@ -1,78 +0,0 @@
-@@ -1,72 +1,76 @@
- # Comprehensive example: quantum anomalous Hall effect
- # ====================================================
- #
- # Physics background
- # ------------------
- # + Quantum anomalous Hall effect
- #
- # Features highlighted
- # --------------------
- # + Use of `kwant.continuum` to discretize a continuum Hamiltonian
- # + Use of `kwant.operator` to compute local current
- # + Use of `kwant.plotter.current` to plot local current
-
-+import _defs
- import math
- import matplotlib.pyplot
- import kwant
- import kwant.continuum
-
-
- # 2 band model exhibiting quantum anomalous Hall effect
- #HIDDEN_BEGIN_model
- def make_model(a):
- ham = ("alpha * (k_x * sigma_x - k_y * sigma_y)"
- "+ (m + beta * kk) * sigma_z"
- "+ (gamma * kk + U) * sigma_0")
- subs = {"kk": "k_x**2 + k_y**2"}
- return kwant.continuum.discretize(ham, locals=subs, grid=a)
- #HIDDEN_END_model
-
-
- def make_system(model, L):
- def lead_shape(site):
- x, y = site.pos / L
- return abs(y) < 0.5
-
- # QPC shape: a rectangle with 2 gaussians
- # etched out of the top and bottom edge.
- def central_shape(site):
- x, y = site.pos / L
- return abs(x) < 3/5 and abs(y) < 0.5 - 0.4 * math.exp(-40 * x**2)
-
- lead = kwant.Builder(kwant.TranslationalSymmetry(
- model.lattice.vec((-1, 0))))
- lead.fill(model, lead_shape, (0, 0))
-
- syst = kwant.Builder()
- syst.fill(model, central_shape, (0, 0))
- syst.attach_lead(lead)
- syst.attach_lead(lead.reversed())
-
- return syst.finalized()
-
-
- def main():
- # Set up our model and system, and define the model parameters.
- params = dict(alpha=0.365, beta=0.686, gamma=0.512, m=-0.01, U=0)
- model = make_model(1)
- syst = make_system(model, 70)
- kwant.plot(syst)
-
- # Calculate the scattering states at energy 'm' coming from the left
- # lead, and the associated particle current.
- psi = kwant.wave_function(syst, energy=params['m'], params=params)(0)
- #HIDDEN_BEGIN_current
- J = kwant.operator.Current(syst).bind(params=params)
- current = sum(J(p) for p in psi)
-- kwant.plotter.current(syst, current)
-+ for extension in ('pdf', 'png'):
-+ kwant.plotter.current(syst, current,
-+ file="plot_qahe_current." + extension,
-+ dpi=_defs.dpi)
- #HIDDEN_END_current
-
-
- if __name__ == '__main__':
- main()
diff --git a/doc/source/code/figure/plot_zincblende.py.diff b/doc/source/code/figure/plot_zincblende.py.diff
deleted file mode 100644
index ad5007871331dd28bf8653ab5be63e102d58205b..0000000000000000000000000000000000000000
--- a/doc/source/code/figure/plot_zincblende.py.diff
+++ /dev/null
@@ -1,71 +0,0 @@
-@@ -1,59 +1,67 @@
- # Tutorial 2.8.2. 3D example: zincblende structure
- # ================================================
- #
- # Physical background
- # -------------------
- # - 3D Bravais lattices
- #
- # Kwant features highlighted
- # --------------------------
- # - demonstrate different ways of plotting in 3D
-
-+import _defs
- 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)
-+ size = (_defs.figwidth_in, _defs.figwidth_in)
-+ for extension in ('pdf', 'png'):
-+ kwant.plot(syst, file="plot_zincblende_syst1." + extension,
-+ fig_size=size, dpi=_defs.dpi)
- #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)
-+ size = (_defs.figwidth_in, _defs.figwidth_in)
-+ for extension in ('pdf', 'png'):
-+ kwant.plot(syst, site_size=0.18, site_lw=0.01, hop_lw=0.05,
-+ site_color=family_colors,
-+ file="plot_zincblende_syst2." + extension,
-+ fig_size=size, dpi=_defs.dpi)
- #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()
diff --git a/doc/source/code/figure/quantum_well.py.diff b/doc/source/code/figure/quantum_well.py.diff
deleted file mode 100644
index 8b1be4bb17cd8513816464e1782381842961ab88..0000000000000000000000000000000000000000
--- a/doc/source/code/figure/quantum_well.py.diff
+++ /dev/null
@@ -1,103 +0,0 @@
-@@ -1,88 +1,95 @@
- # Tutorial 2.3.2. Spatially dependent values through functions
- # ============================================================
- #
- # Physics background
- # ------------------
- # transmission through a quantum well
- #
- # Kwant features highlighted
- # --------------------------
- # - Functions as values in Builder
-
-+import _defs
- import kwant
-
- # For plotting
- from matplotlib import pyplot
-
-
- #HIDDEN_BEGIN_ehso
- def make_system(a=1, t=1.0, W=10, L=30, L_well=10):
- # Start with an empty tight-binding system and a single square lattice.
- # `a` is the lattice constant (by default set to 1 for simplicity.
- lat = kwant.lattice.square(a)
-
- syst = kwant.Builder()
-
- #### Define the scattering region. ####
- # Potential profile
- def potential(site, pot):
- (x, y) = site.pos
- if (L - L_well) / 2 < x < (L + L_well) / 2:
- return pot
- else:
- return 0
- #HIDDEN_END_ehso
-
- #HIDDEN_BEGIN_coid
- def onsite(site, pot):
- return 4 * t + potential(site, pot)
-
- syst[(lat(x, y) for x in range(L) for y in range(W))] = onsite
- syst[lat.neighbors()] = -t
- #HIDDEN_END_coid
-
- #### Define and attach the leads. ####
- lead = kwant.Builder(kwant.TranslationalSymmetry((-a, 0)))
- lead[(lat(0, j) for j in range(W))] = 4 * t
- lead[lat.neighbors()] = -t
- syst.attach_lead(lead)
- syst.attach_lead(lead.reversed())
-
- return syst
-
-
- def plot_conductance(syst, energy, welldepths):
- #HIDDEN_BEGIN_sqvr
-
- # Compute conductance
- data = []
- for welldepth in welldepths:
- smatrix = kwant.smatrix(syst, energy, params=dict(pot=-welldepth))
- data.append(smatrix.transmission(1, 0))
-
-- pyplot.figure()
-+ fig = pyplot.figure()
- pyplot.plot(welldepths, data)
-- pyplot.xlabel("well depth [t]")
-- pyplot.ylabel("conductance [e^2/h]")
-- pyplot.show()
-+ pyplot.xlabel("well depth [t]",
-+ fontsize=_defs.mpl_label_size)
-+ pyplot.ylabel("conductance [e^2/h]",
-+ fontsize=_defs.mpl_label_size)
-+ pyplot.setp(fig.get_axes()[0].get_xticklabels(),
-+ fontsize=_defs.mpl_tick_size)
-+ pyplot.setp(fig.get_axes()[0].get_yticklabels(),
-+ fontsize=_defs.mpl_tick_size)
-+ fig.set_size_inches(_defs.mpl_width_in, _defs.mpl_width_in * 3. / 4.)
-+ fig.subplots_adjust(left=0.15, right=0.95, top=0.95, bottom=0.15)
-+ for extension in ('pdf', 'png'):
-+ fig.savefig("quantum_well_result." + extension, dpi=_defs.dpi)
- #HIDDEN_END_sqvr
-
-
- def main():
- syst = make_system()
-
-- # Check that the system looks as intended.
-- kwant.plot(syst)
--
- # Finalize the system.
- syst = syst.finalized()
-
- # We should see conductance steps.
- plot_conductance(syst, energy=0.2,
- welldepths=[0.01 * i 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()
diff --git a/doc/source/code/figure/quantum_wire.py.diff b/doc/source/code/figure/quantum_wire.py.diff
deleted file mode 100644
index 396ab255359afa8c53f7790a1c53b7b9fe2de178..0000000000000000000000000000000000000000
--- a/doc/source/code/figure/quantum_wire.py.diff
+++ /dev/null
@@ -1,136 +0,0 @@
-@@ -1,121 +1,130 @@
- # Tutorial 2.2.2. Transport through a quantum wire
- # ================================================
- #
- # Physics background
- # ------------------
- # Conductance of a quantum wire; subbands
- #
- # Kwant features highlighted
- # --------------------------
- # - Builder for setting up transport systems easily
- # - Making scattering region and leads
- # - Using the simple sparse solver for computing Landauer conductance
-
-+import _defs
- from matplotlib import pyplot
- #HIDDEN_BEGIN_dwhx
- import kwant
- #HIDDEN_END_dwhx
-
- # First, define the tight-binding system
-
- #HIDDEN_BEGIN_goiq
- syst = kwant.Builder()
- #HIDDEN_END_goiq
-
- # Here, we are only working with square lattices
- #HIDDEN_BEGIN_suwo
- a = 1
- lat = kwant.lattice.square(a)
- #HIDDEN_END_suwo
-
- #HIDDEN_BEGIN_zfvr
- t = 1.0
- W = 10
- L = 30
-
- # Define the scattering region
-
- for i in range(L):
- for j in range(W):
- # On-site Hamiltonian
- syst[lat(i, j)] = 4 * t
-
- # Hopping in y-direction
- if j > 0:
- syst[lat(i, j), lat(i, j - 1)] = -t
-
- # Hopping in x-direction
- if i > 0:
- syst[lat(i, j), lat(i - 1, j)] = -t
- #HIDDEN_END_zfvr
-
- # Then, define and attach the leads:
-
- # First the lead to the left
- # (Note: TranslationalSymmetry takes a real-space vector)
- #HIDDEN_BEGIN_xcmc
- sym_left_lead = kwant.TranslationalSymmetry((-a, 0))
- left_lead = kwant.Builder(sym_left_lead)
- #HIDDEN_END_xcmc
-
- #HIDDEN_BEGIN_ndez
- for j in range(W):
- left_lead[lat(0, j)] = 4 * t
- if j > 0:
- left_lead[lat(0, j), lat(0, j - 1)] = -t
- left_lead[lat(1, j), lat(0, j)] = -t
- #HIDDEN_END_ndez
-
- #HIDDEN_BEGIN_fskr
- syst.attach_lead(left_lead)
- #HIDDEN_END_fskr
-
- # Then the lead to the right
- #HIDDEN_BEGIN_xhqc
- sym_right_lead = kwant.TranslationalSymmetry((a, 0))
- right_lead = kwant.Builder(sym_right_lead)
-
- for j in range(W):
- right_lead[lat(0, j)] = 4 * t
- if j > 0:
- right_lead[lat(0, j), lat(0, j - 1)] = -t
- right_lead[lat(1, j), lat(0, j)] = -t
-
- syst.attach_lead(right_lead)
- #HIDDEN_END_xhqc
-
- # Plot it, to make sure it's OK
- #HIDDEN_BEGIN_wsgh
--kwant.plot(syst)
-+size = (_defs.figwidth_in, 0.3 * _defs.figwidth_in)
-+for extension in ('pdf', 'png'):
-+ kwant.plot(syst, file="quantum_wire_syst." + extension,
-+ fig_size=size, dpi=_defs.dpi)
- #HIDDEN_END_wsgh
-
- # Finalize the system
- #HIDDEN_BEGIN_dngj
- syst = syst.finalized()
- #HIDDEN_END_dngj
-
- # Now that we have the system, we can compute conductance
- #HIDDEN_BEGIN_buzn
- energies = []
- data = []
- for ie in range(100):
- energy = ie * 0.01
-
- # compute the scattering matrix at a given energy
- smatrix = kwant.smatrix(syst, energy)
-
- # compute the transmission probability from lead 0 to
- # lead 1
- energies.append(energy)
- data.append(smatrix.transmission(1, 0))
- #HIDDEN_END_buzn
-
- # Use matplotlib to write output
- # We should see conductance steps
- #HIDDEN_BEGIN_lliv
--pyplot.figure()
-+fig = pyplot.figure()
- pyplot.plot(energies, data)
--pyplot.xlabel("energy [t]")
--pyplot.ylabel("conductance [e^2/h]")
--pyplot.show()
-+pyplot.xlabel("energy [t]", fontsize=_defs.mpl_label_size)
-+pyplot.ylabel("conductance [e^2/h]", fontsize=_defs.mpl_label_size)
-+pyplot.setp(fig.get_axes()[0].get_xticklabels(), fontsize=_defs.mpl_tick_size)
-+pyplot.setp(fig.get_axes()[0].get_yticklabels(), fontsize=_defs.mpl_tick_size)
-+fig.set_size_inches(_defs.mpl_width_in, _defs.mpl_width_in * 3. / 4.)
-+fig.subplots_adjust(left=0.15, right=0.95, top=0.95, bottom=0.15)
-+for extension in ('pdf', 'png'):
-+ fig.savefig("quantum_wire_result." + extension, dpi=_defs.dpi)
- #HIDDEN_END_lliv
diff --git a/doc/source/code/figure/spin_orbit.py.diff b/doc/source/code/figure/spin_orbit.py.diff
deleted file mode 100644
index 794c1b5b755e22ea1a175c40a9078bd3e5b3e12b..0000000000000000000000000000000000000000
--- a/doc/source/code/figure/spin_orbit.py.diff
+++ /dev/null
@@ -1,118 +0,0 @@
-@@ -1,104 +1,110 @@
- # Tutorial 2.3.1. Matrix structure of on-site and hopping elements
- # ================================================================
- #
- # Physics background
- # ------------------
- # Gaps in quantum wires with spin-orbit coupling and Zeeman splititng,
- # as theoretically predicted in
- # http://prl.aps.org/abstract/PRL/v90/i25/e256601
- # and (supposedly) experimentally oberved in
- # http://www.nature.com/nphys/journal/v6/n5/abs/nphys1626.html
- #
- # Kwant features highlighted
- # --------------------------
- # - Numpy matrices as values in Builder
-
-+import _defs
- import kwant
-
- # For plotting
- from matplotlib import pyplot
-
- # For matrix support
- #HIDDEN_BEGIN_xumz
- import tinyarray
- #HIDDEN_END_xumz
-
- # define Pauli-matrices for convenience
- #HIDDEN_BEGIN_hwbt
- sigma_0 = tinyarray.array([[1, 0], [0, 1]])
- sigma_x = tinyarray.array([[0, 1], [1, 0]])
- sigma_y = tinyarray.array([[0, -1j], [1j, 0]])
- sigma_z = tinyarray.array([[1, 0], [0, -1]])
- #HIDDEN_END_hwbt
-
-
- def make_system(t=1.0, alpha=0.5, e_z=0.08, W=10, L=30):
- # Start with an empty tight-binding system and a single square lattice.
- # `a` is the lattice constant (by default set to 1 for simplicity).
- lat = kwant.lattice.square()
-
- syst = kwant.Builder()
-
- #### Define the scattering region. ####
- #HIDDEN_BEGIN_uxrm
- syst[(lat(x, y) for x in range(L) for y in range(W))] = \
- 4 * t * sigma_0 + e_z * sigma_z
- # hoppings in x-direction
- syst[kwant.builder.HoppingKind((1, 0), lat, lat)] = \
- -t * sigma_0 + 1j * alpha * sigma_y / 2
- # hoppings in y-directions
- syst[kwant.builder.HoppingKind((0, 1), lat, lat)] = \
- -t * sigma_0 - 1j * alpha * sigma_x / 2
- #HIDDEN_END_uxrm
-
- #### Define the left lead. ####
- lead = kwant.Builder(kwant.TranslationalSymmetry((-1, 0)))
-
- #HIDDEN_BEGIN_yliu
- lead[(lat(0, j) for j in range(W))] = 4 * t * sigma_0 + e_z * sigma_z
- # hoppings in x-direction
- lead[kwant.builder.HoppingKind((1, 0), lat, lat)] = \
- -t * sigma_0 + 1j * alpha * sigma_y / 2
- # hoppings in y-directions
- lead[kwant.builder.HoppingKind((0, 1), lat, lat)] = \
- -t * sigma_0 - 1j * alpha * sigma_x / 2
- #HIDDEN_END_yliu
-
- #### Attach the leads and return the finalized system. ####
- syst.attach_lead(lead)
- syst.attach_lead(lead.reversed())
-
- return syst
-
-
- def plot_conductance(syst, energies):
- # Compute conductance
- data = []
- for energy in energies:
- smatrix = kwant.smatrix(syst, energy)
- data.append(smatrix.transmission(1, 0))
-
-- pyplot.figure()
-+ fig = pyplot.figure()
- pyplot.plot(energies, data)
-- pyplot.xlabel("energy [t]")
-- pyplot.ylabel("conductance [e^2/h]")
-- pyplot.show()
-+ pyplot.xlabel("energy [t]", fontsize=_defs.mpl_label_size)
-+ pyplot.ylabel("conductance [e^2/h]",
-+ fontsize=_defs.mpl_label_size)
-+ pyplot.setp(fig.get_axes()[0].get_xticklabels(),
-+ fontsize=_defs.mpl_tick_size)
-+ pyplot.setp(fig.get_axes()[0].get_yticklabels(),
-+ fontsize=_defs.mpl_tick_size)
-+ fig.set_size_inches(_defs.mpl_width_in, _defs.mpl_width_in * 3. / 4.)
-+ fig.subplots_adjust(left=0.15, right=0.95, top=0.95, bottom=0.15)
-+ for extension in ('pdf', 'png'):
-+ fig.savefig("spin_orbit_result." + extension, dpi=_defs.dpi)
-
-
- def main():
- syst = make_system()
-
-- # Check that the system looks as intended.
-- kwant.plot(syst)
--
- # Finalize the system.
- syst = syst.finalized()
-
- # We should see non-monotonic conductance steps.
- plot_conductance(syst, energies=[0.01 * i - 0.3 for i in range(100)])
-
-
- # Call the main function if the script gets executed (as opposed to imported).
- # See <http://docs.python.org/library/__main__.html>.
- if __name__ == '__main__':
- main()
diff --git a/doc/source/code/figure/superconductor.py.diff b/doc/source/code/figure/superconductor.py.diff
deleted file mode 100644
index 75d9b061f233f5fb1a45c874977e3156cf7f1296..0000000000000000000000000000000000000000
--- a/doc/source/code/figure/superconductor.py.diff
+++ /dev/null
@@ -1,148 +0,0 @@
-@@ -1,132 +1,141 @@
- # Tutorial 2.6. "Superconductors": orbitals, conservation laws and symmetries
- # ===========================================================================
- #
- # Physics background
- # ------------------
- # - conductance of a NS-junction (Andreev reflection, superconducting gap)
- #
- # Kwant features highlighted
- # --------------------------
- # - Implementing electron and hole ("orbital") degrees of freedom
- # using conservation laws.
- # - Use of discrete symmetries to relate scattering states.
-
-+import _defs
- import kwant
-
- import tinyarray
- import numpy as np
-
- # For plotting
- from matplotlib import pyplot
-+from contextlib import redirect_stdout
-
- tau_x = tinyarray.array([[0, 1], [1, 0]])
- tau_y = tinyarray.array([[0, -1j], [1j, 0]])
- tau_z = tinyarray.array([[1, 0], [0, -1]])
-
- #HIDDEN_BEGIN_nbvn
- def make_system(a=1, W=10, L=10, barrier=1.5, barrierpos=(3, 4),
- mu=0.4, Delta=0.1, Deltapos=4, t=1.0):
- # Start with an empty tight-binding system. On each site, there
- # are now electron and hole orbitals, so we must specify the
- # number of orbitals per site. The orbital structure is the same
- # as in the Hamiltonian.
- lat = kwant.lattice.square(norbs=2)
- syst = kwant.Builder()
-
- #### Define the scattering region. ####
- # The superconducting order parameter couples electron and hole orbitals
- # on each site, and hence enters as an onsite potential.
- # The pairing is only included beyond the point 'Deltapos' in the scattering region.
- syst[(lat(x, y) for x in range(Deltapos) for y in range(W))] = (4 * t - mu) * tau_z
- syst[(lat(x, y) for x in range(Deltapos, L) for y in range(W))] = (4 * t - mu) * tau_z + Delta * tau_x
-
- # The tunnel barrier
- syst[(lat(x, y) for x in range(barrierpos[0], barrierpos[1])
- for y in range(W))] = (4 * t + barrier - mu) * tau_z
-
- # Hoppings
- syst[lat.neighbors()] = -t * tau_z
- #HIDDEN_END_nbvn
- #HIDDEN_BEGIN_ttth
- #### Define the leads. ####
- # Left lead - normal, so the order parameter is zero.
- sym_left = kwant.TranslationalSymmetry((-a, 0))
- # Specify the conservation law used to treat electrons and holes separately.
- # We only do this in the left lead, where the pairing is zero.
- lead0 = kwant.Builder(sym_left, conservation_law=-tau_z, particle_hole=tau_y)
- lead0[(lat(0, j) for j in range(W))] = (4 * t - mu) * tau_z
- lead0[lat.neighbors()] = -t * tau_z
- #HIDDEN_END_ttth
- #HIDDEN_BEGIN_zuuw
- # Right lead - superconducting, so the order parameter is included.
- sym_right = kwant.TranslationalSymmetry((a, 0))
- lead1 = kwant.Builder(sym_right)
- lead1[(lat(0, j) for j in range(W))] = (4 * t - mu) * tau_z + Delta * tau_x
- lead1[lat.neighbors()] = -t * tau_z
-
- #### Attach the leads and return the system. ####
- syst.attach_lead(lead0)
- syst.attach_lead(lead1)
-
- return syst
- #HIDDEN_END_zuuw
-
- #HIDDEN_BEGIN_jbjt
- def plot_conductance(syst, energies):
- # Compute conductance
- data = []
- for energy in energies:
- smatrix = kwant.smatrix(syst, energy)
- # Conductance is N - R_ee + R_he
- data.append(smatrix.submatrix((0, 0), (0, 0)).shape[0] -
- smatrix.transmission((0, 0), (0, 0)) +
- smatrix.transmission((0, 1), (0, 0)))
- #HIDDEN_END_jbjt
-- pyplot.figure()
-+ fig = pyplot.figure()
- pyplot.plot(energies, data)
- pyplot.xlabel("energy [t]")
- pyplot.ylabel("conductance [e^2/h]")
-- pyplot.show()
-+ pyplot.setp(fig.get_axes()[0].get_xticklabels(),
-+ fontsize=_defs.mpl_tick_size)
-+ pyplot.setp(fig.get_axes()[0].get_yticklabels(),
-+ fontsize=_defs.mpl_tick_size)
-+ fig.set_size_inches(_defs.mpl_width_in, _defs.mpl_width_in * 3. / 4.)
-+ fig.subplots_adjust(left=0.15, right=0.95, top=0.95, bottom=0.15)
-+ for extension in ('pdf', 'png'):
-+ fig.savefig("superconductor_transport_result." + extension,
-+ dpi=_defs.dpi)
-
- #HIDDEN_BEGIN_pqmp
- def check_PHS(syst):
- # Scattering matrix
- s = kwant.smatrix(syst, energy=0)
- # Electron to electron block
- s_ee = s.submatrix((0,0), (0,0))
- # Hole to hole block
- s_hh = s.submatrix((0,1), (0,1))
- print('s_ee: \n', np.round(s_ee, 3))
- print('s_hh: \n', np.round(s_hh[::-1, ::-1], 3))
- print('s_ee - s_hh^*: \n',
- np.round(s_ee - s_hh[::-1, ::-1].conj(), 3), '\n')
- # Electron to hole block
- s_he = s.submatrix((0,1), (0,0))
- # Hole to electron block
- s_eh = s.submatrix((0,0), (0,1))
- print('s_he: \n', np.round(s_he, 3))
- print('s_eh: \n', np.round(s_eh[::-1, ::-1], 3))
- print('s_he + s_eh^*: \n',
- np.round(s_he + s_eh[::-1, ::-1].conj(), 3))
- #HIDDEN_END_pqmp
-
- def main():
- syst = make_system(W=10)
-
-- # Check that the system looks as intended.
-- kwant.plot(syst)
--
- # Finalize the system.
- syst = syst.finalized()
-
- # Check particle-hole symmetry of the scattering matrix
-- check_PHS(syst)
-+ with open('check_PHS_out.txt', 'w') as f:
-+ with redirect_stdout(f):
-+ check_PHS(syst)
-
- # Compute and plot the conductance
- plot_conductance(syst, energies=[0.002 * i for i in range(-10, 100)])
-
-
- # Call the main function if the script gets executed (as opposed to imported).
- # See <http://docs.python.org/library/__main__.html>.
- if __name__ == '__main__':
- main()
diff --git a/doc/source/code/include/quantum_wire_revisited.py b/doc/source/code/include/quantum_wire_revisited.py
deleted file mode 100644
index 60add47870d7c6690e30cd4a9e46adc69d8d4d6f..0000000000000000000000000000000000000000
--- a/doc/source/code/include/quantum_wire_revisited.py
+++ /dev/null
@@ -1,91 +0,0 @@
-# Tutorial 2.2.3. Building the same system with less code
-# =======================================================
-#
-# Physics background
-# ------------------
-# Conductance of a quantum wire; subbands
-#
-# Kwant features highlighted
-# --------------------------
-# - Using iterables and builder.HoppingKind for making systems
-# - introducing `reversed()` for the leads
-#
-# Note: Does the same as tutorial1a.py, but using other features of Kwant.
-
-#HIDDEN_BEGIN_xkzy
-import kwant
-
-# For plotting
-from matplotlib import pyplot
-
-
-def make_system(a=1, t=1.0, W=10, L=30):
- # Start with an empty tight-binding system and a single square lattice.
- # `a` is the lattice constant (by default set to 1 for simplicity.
- lat = kwant.lattice.square(a)
-
- syst = kwant.Builder()
-#HIDDEN_END_xkzy
-
- #### Define the scattering region. ####
-#HIDDEN_BEGIN_vvjt
- syst[(lat(x, y) for x in range(L) for y in range(W))] = 4 * t
-#HIDDEN_END_vvjt
-#HIDDEN_BEGIN_nooi
- syst[lat.neighbors()] = -t
-#HIDDEN_END_nooi
-
- #### Define and attach the leads. ####
- # Construct the left lead.
-#HIDDEN_BEGIN_iepx
- lead = kwant.Builder(kwant.TranslationalSymmetry((-a, 0)))
- lead[(lat(0, j) for j in range(W))] = 4 * t
- lead[lat.neighbors()] = -t
-#HIDDEN_END_iepx
-
- # Attach the left lead and its reversed copy.
-#HIDDEN_BEGIN_yxot
- syst.attach_lead(lead)
- syst.attach_lead(lead.reversed())
-
- return syst
-#HIDDEN_END_yxot
-
-
-#HIDDEN_BEGIN_ayuk
-def plot_conductance(syst, energies):
- # Compute conductance
- data = []
- for energy in energies:
- smatrix = kwant.smatrix(syst, energy)
- data.append(smatrix.transmission(1, 0))
-
- pyplot.figure()
- pyplot.plot(energies, data)
- pyplot.xlabel("energy [t]")
- pyplot.ylabel("conductance [e^2/h]")
- pyplot.show()
-#HIDDEN_END_ayuk
-
-
-#HIDDEN_BEGIN_cjel
-def main():
- syst = make_system()
-
- # Check that the system looks as intended.
- kwant.plot(syst)
-
- # Finalize the system.
- syst = syst.finalized()
-
- # We should see conductance steps.
- plot_conductance(syst, energies=[0.01 * i for i in range(100)])
-#HIDDEN_END_cjel
-
-
-# Call the main function if the script gets executed (as opposed to imported).
-# See <http://docs.python.org/library/__main__.html>.
-#HIDDEN_BEGIN_ypbj
-if __name__ == '__main__':
- main()
-#HIDDEN_END_ypbj
diff --git a/doc/source/conf.py b/doc/source/conf.py
index d3d7fceb5022796e9fe464cde6def683778bed46..080a71ff0039aaf8528ddf692c45e1569050697a 100644
--- a/doc/source/conf.py
+++ b/doc/source/conf.py
@@ -15,8 +15,15 @@ import sys
import os
import string
from distutils.util import get_platform
-sys.path.insert(0, "../../build/lib.{0}-{1}.{2}".format(
- get_platform(), *sys.version_info[:2]))
+
+package_path = os.path.abspath(
+ "../../build/lib.{0}-{1}.{2}"
+ .format(get_platform(), *sys.version_info[:2]))
+
+# Insert into sys.path so that we can import kwant here
+sys.path.insert(0, package_path)
+# Insert into PYTHONPATH so that jupyter-sphinx will pick it up
+os.environ['PYTHONPATH'] = ':'.join((package_path, os.environ.get('PYTHONPATH','')))
import kwant
import kwant.qsymm
@@ -31,7 +38,8 @@ sys.path.insert(0, os.path.abspath('../sphinxext'))
extensions = ['sphinx.ext.autodoc', 'sphinx.ext.autosummary',
'sphinx.ext.todo', 'sphinx.ext.mathjax', 'numpydoc',
- 'kwantdoc', 'sphinx.ext.linkcode']
+ 'kwantdoc', 'sphinx.ext.linkcode', 'jupyter_sphinx.execute',
+ 'sphinxcontrib.rsvgconverter']
# Add any paths that contain templates here, relative to this directory.
templates_path = ['../templates']
@@ -239,7 +247,6 @@ autosummary_generate = True
autoclass_content = "both"
autodoc_default_flags = ['show-inheritance']
-
# -- Teach Sphinx to document bound methods like functions ---------------------
import types
from sphinx.ext import autodoc
diff --git a/doc/source/code/figure/ab_ring_sketch.svg b/doc/source/figure/ab_ring_sketch.svg
similarity index 100%
rename from doc/source/code/figure/ab_ring_sketch.svg
rename to doc/source/figure/ab_ring_sketch.svg
diff --git a/doc/source/code/figure/ab_ring_sketch2.svg b/doc/source/figure/ab_ring_sketch2.svg
similarity index 100%
rename from doc/source/code/figure/ab_ring_sketch2.svg
rename to doc/source/figure/ab_ring_sketch2.svg
diff --git a/doc/source/code/figure/superconductor_transport_sketch.svg b/doc/source/figure/superconductor_transport_sketch.svg
similarity index 100%
rename from doc/source/code/figure/superconductor_transport_sketch.svg
rename to doc/source/figure/superconductor_transport_sketch.svg
diff --git a/doc/source/pre/whatsnew/1.3.rst b/doc/source/pre/whatsnew/1.3.rst
index f7dbd8fc0ef4373499d976fc3b908f04548e629e..8336e99d9c66dc1fb7628414fe28cf648c18319c 100644
--- a/doc/source/pre/whatsnew/1.3.rst
+++ b/doc/source/pre/whatsnew/1.3.rst
@@ -28,11 +28,43 @@ tight-binding band structures or construct systems with different/lower
symmetry. For example in just a few lines we can construct a two-band model that exhibits
a quantum anomalous spin Hall phase:
-.. literalinclude:: ../../code/include/plot_qahe.py
- :start-after: HIDDEN_BEGIN_model
- :end-before: HIDDEN_END_model
+.. jupyter-kernel::
+ :id: plot_qahe
+
+.. jupyter-execute::
+ :hide-code:
+
+ # Comprehensive example: quantum anomalous Hall effect
+ # ====================================================
+ #
+ # Physics background
+ # ------------------
+ # + Quantum anomalous Hall effect
+ #
+ # Features highlighted
+ # --------------------
+ # + Use of `kwant.continuum` to discretize a continuum Hamiltonian
+ # + Use of `kwant.operator` to compute local current
+ # + Use of `kwant.plotter.current` to plot local current
+
+ import math
+ import matplotlib.pyplot
+ import kwant
+ import kwant.continuum
+
+.. jupyter-execute:: ../../tutorial/boilerplate.py
+ :hide-code:
+
+.. jupyter-execute::
-From: :download:`QAHE example script <../../code/download/plot_qahe.py>`
+ def make_model(a):
+ ham = ("alpha * (k_x * sigma_x - k_y * sigma_y)"
+ "+ (m + beta * kk) * sigma_z"
+ "+ (gamma * kk + U) * sigma_0")
+ subs = {"kk": "k_x**2 + k_y**2"}
+ return kwant.continuum.discretize(ham, locals=subs, grid=a)
+
+From: :jupyter-download:script:`plot_qahe`
See the tutorial: :doc:`../../tutorial/discretize`
@@ -71,13 +103,47 @@ The example below shows edge states of a quantum anomalous Hall phase
of the two-band model shown in the `above section
<#tools-for-continuum-hamiltonians>`_:
-.. literalinclude:: ../../code/include/plot_qahe.py
- :start-after: HIDDEN_BEGIN_current
- :end-before: HIDDEN_END_current
+.. jupyter-execute::
+ :hide-code:
+
+ def make_system(model, L):
+ def lead_shape(site):
+ x, y = site.pos / L
+ return abs(y) < 0.5
+
+ # QPC shape: a rectangle with 2 gaussians
+ # etched out of the top and bottom edge.
+ def central_shape(site):
+ x, y = site.pos / L
+ return abs(x) < 3/5 and abs(y) < 0.5 - 0.4 * math.exp(-40 * x**2)
+
+ lead = kwant.Builder(kwant.TranslationalSymmetry(
+ model.lattice.vec((-1, 0))))
+ lead.fill(model, lead_shape, (0, 0))
+
+ syst = kwant.Builder()
+ syst.fill(model, central_shape, (0, 0))
+ syst.attach_lead(lead)
+ syst.attach_lead(lead.reversed())
+
+ return syst.finalized()
+
+ # Set up our model and system, and define the model parameters.
+ params = dict(alpha=0.365, beta=0.686, gamma=0.512, m=-0.01, U=0)
+ model = make_model(1)
+ syst = make_system(model, 70)
+
+ # Calculate the scattering states at energy 'm' coming from the left
+ # lead, and the associated particle current.
+ psi = kwant.wave_function(syst, energy=params['m'], params=params)(0)
+
+.. jupyter-execute::
-.. image:: ../../code/figure/plot_qahe_current.*
+ J = kwant.operator.Current(syst).bind(params=params)
+ current = sum(J(p) for p in psi)
+ kwant.plotter.current(syst, current);
-From: :download:`QAHE example script <../../code/download/plot_qahe.py>`
+From: :jupyter-download:script:`plot_qahe`
Scattering states with discrete symmetries and conservation laws
----------------------------------------------------------------
diff --git a/doc/source/pre/whatsnew/1.5.rst b/doc/source/pre/whatsnew/1.5.rst
index b985489ec6a45dd48c04bcab293d8454f76abab9..9a9c085c8d924c0cda45d84dc75a5b1a0b893016 100644
--- a/doc/source/pre/whatsnew/1.5.rst
+++ b/doc/source/pre/whatsnew/1.5.rst
@@ -2,3 +2,13 @@ What's new in Kwant 1.5
=======================
This article explains the user-visible changes in Kwant 1.5.0.
+
+Improved tutorial building
+--------------------------
+Improving or adding to Kwant's tutorial is now much simpler. Now
+the text and code for each tutorial is kept in the same file, making
+it easy to see where changes need to be made, and images generated by
+the code are inserted directly into the document thanks to the magic of
+`jupyter-sphinx <https://github.com/jupyter-widgets/jupyter-sphinx/>`_.
+It has never been easier to get started contributing to Kwant by
+helping us improve our documentation.
diff --git a/doc/source/tutorial/boilerplate.py b/doc/source/tutorial/boilerplate.py
new file mode 100644
index 0000000000000000000000000000000000000000..db2fba37a6f8546b66d4ca8f06a8eda668438ca7
--- /dev/null
+++ b/doc/source/tutorial/boilerplate.py
@@ -0,0 +1,6 @@
+import matplotlib
+import matplotlib.pyplot
+from IPython.display import set_matplotlib_formats
+
+matplotlib.rcParams['figure.figsize'] = matplotlib.pyplot.figaspect(1) * 2
+set_matplotlib_formats('svg')
diff --git a/doc/source/tutorial/discretize.rst b/doc/source/tutorial/discretize.rst
index b737b00ca074cbfd88a7e863be67ebc03063d902..31c4b3bb3ca217888c3c0e1b53e444a270906705 100644
--- a/doc/source/tutorial/discretize.rst
+++ b/doc/source/tutorial/discretize.rst
@@ -18,8 +18,32 @@ continuum models and for discretizing them into tight-binding models.
.. seealso::
The complete source code of this tutorial can be found in
- :download:`discretize.py </code/download/discretize.py>`
+ :jupyter-download:script:`discretize`
+.. jupyter-kernel::
+ :id: discretize
+
+.. jupyter-execute::
+ :hide-code:
+
+ # Tutorial 2.9. Processing continuum Hamiltonians with discretize
+ # ===============================================================
+ #
+ # Physics background
+ # ------------------
+ # - tight-binding approximation of continuous Hamiltonians
+ #
+ # Kwant features highlighted
+ # --------------------------
+ # - kwant.continuum.discretize
+
+ import matplotlib as mpl
+ from matplotlib import pyplot
+
+ import kwant
+
+.. jupyter-execute:: boilerplate.py
+ :hide-code:
.. _tutorial_discretizer_introduction:
@@ -55,9 +79,16 @@ Using `~kwant.continuum.discretize` to obtain a template
........................................................
First we must explicitly import the `kwant.continuum` package:
-.. literalinclude:: /code/include/discretize.py
- :start-after: #HIDDEN_BEGIN_import
- :end-before: #HIDDEN_END_import
+.. jupyter-execute::
+
+ import kwant.continuum
+
+.. jupyter-execute::
+ :hide-code:
+
+ import scipy.sparse.linalg
+ import scipy.linalg
+ import numpy as np
The function `kwant.continuum.discretize` takes a symbolic Hamiltonian and
turns it into a `~kwant.builder.Builder` instance with appropriate spatial
@@ -65,17 +96,16 @@ symmetry that serves as a template.
(We will see how to use the template to build systems with a particular
shape later).
-.. literalinclude:: /code/include/discretize.py
- :start-after: #HIDDEN_BEGIN_symbolic_discretization
- :end-before: #HIDDEN_END_symbolic_discretization
+.. jupyter-execute::
+
+ template = kwant.continuum.discretize('k_x * A(x) * k_x')
+ print(template)
It is worth noting that ``discretize`` treats ``k_x`` and ``x`` as
non-commuting operators, and so their order is preserved during the
discretization process.
-The builder produced by ``discretize`` may be printed to show the source code of its onsite and hopping functions (this is a special feature of builders returned by ``discretize``):
-
-.. literalinclude:: /code/figure/discretizer_intro_verbose.txt
+Printing the Builder produced by ``discretize`` shows the source code of its onsite and hopping functions (this is a special feature of builders returned by ``discretize``).
.. specialnote:: Technical details
@@ -139,26 +169,45 @@ where :math:`V(x, y)` is some arbitrary potential.
First, use ``discretize`` to obtain a
builder that we will use as a template:
-.. literalinclude:: /code/include/discretize.py
- :start-after: #HIDDEN_BEGIN_template
- :end-before: #HIDDEN_END_template
+.. jupyter-execute::
+
+ hamiltonian = "k_x**2 + k_y**2 + V(x, y)"
+ template = kwant.continuum.discretize(hamiltonian)
+ print(template)
We now use this system with the `~kwant.builder.Builder.fill`
method of `~kwant.builder.Builder` to construct the system we
want to investigate:
-.. literalinclude:: /code/include/discretize.py
- :start-after: #HIDDEN_BEGIN_fill
- :end-before: #HIDDEN_END_fill
+.. jupyter-execute::
+
+ def stadium(site):
+ (x, y) = site.pos
+ x = max(abs(x) - 20, 0)
+ return x**2 + y**2 < 30**2
+
+ syst = kwant.Builder()
+ syst.fill(template, stadium, (0, 0));
+ syst = syst.finalized()
After finalizing this system, we can plot one of the system's
energy eigenstates:
-.. literalinclude:: /code/include/discretize.py
- :start-after: #HIDDEN_BEGIN_plot_eigenstate
- :end-before: #HIDDEN_END_plot_eigenstate
+.. jupyter-execute::
+
+ def plot_eigenstate(syst, n=2, Vx=.0003, Vy=.0005):
+
+ def potential(x, y):
+ return Vx * x + Vy * y
-.. image:: /code/figure/discretizer_gs.*
+ ham = syst.hamiltonian_submatrix(params=dict(V=potential), sparse=True)
+ evecs = scipy.sparse.linalg.eigsh(ham, k=10, which='SM')[1]
+ kwant.plotter.map(syst, abs(evecs[:, n])**2, show=False)
+
+.. jupyter-execute::
+ :hide-code:
+
+ plot_eigenstate(syst)
Note in the above that we pass the spatially varying potential *function*
to our system via a parameter called ``V``, because the symbol :math:`V`
@@ -176,32 +225,95 @@ model [1]_ [2]_, one can provide matrix input to `~kwant.continuum.discretize`
using ``identity`` and ``kron``. For example, the definition of the BHZ model can be
written succinctly as:
-.. literalinclude:: /code/include/discretize.py
- :start-after: #HIDDEN_BEGIN_define_qsh
- :end-before: #HIDDEN_END_define_qsh
+.. jupyter-execute::
+
+ hamiltonian = """
+ + C * identity(4) + M * kron(sigma_0, sigma_z)
+ - B * (k_x**2 + k_y**2) * kron(sigma_0, sigma_z)
+ - D * (k_x**2 + k_y**2) * kron(sigma_0, sigma_0)
+ + A * k_x * kron(sigma_z, sigma_x)
+ - A * k_y * kron(sigma_0, sigma_y)
+ """
+
+ a = 20
+
+ template = kwant.continuum.discretize(hamiltonian, grid=a)
We can then make a ribbon out of this template system:
-.. literalinclude:: /code/include/discretize.py
- :start-after: #HIDDEN_BEGIN_define_qsh_build
- :end-before: #HIDDEN_END_define_qsh_build
+.. jupyter-execute::
+
+ L, W = 2000, 1000
+
+ def shape(site):
+ (x, y) = site.pos
+ return (0 <= y < W and 0 <= x < L)
+
+ def lead_shape(site):
+ (x, y) = site.pos
+ return (0 <= y < W)
+
+ syst = kwant.Builder()
+ syst.fill(template, shape, (0, 0))
+
+ lead = kwant.Builder(kwant.TranslationalSymmetry([-a, 0]))
+ lead.fill(template, lead_shape, (0, 0))
+
+ syst.attach_lead(lead)
+ syst.attach_lead(lead.reversed())
+
+ syst = syst.finalized()
and plot its dispersion using `kwant.plotter.bands`:
-.. literalinclude:: /code/include/discretize.py
- :start-after: #HIDDEN_BEGIN_plot_qsh_band
- :end-before: #HIDDEN_END_plot_qsh_band
+.. jupyter-execute::
-.. image:: /code/figure/discretizer_qsh_band.*
+ params = dict(A=3.65, B=-68.6, D=-51.1, M=-0.01, C=0)
+
+ kwant.plotter.bands(syst.leads[0], params=params,
+ momenta=np.linspace(-0.3, 0.3, 201), show=False)
+
+ pyplot.grid()
+ pyplot.xlim(-.3, 0.3)
+ pyplot.ylim(-0.05, 0.05)
+ pyplot.xlabel('momentum [1/A]')
+ pyplot.ylabel('energy [eV]')
+ pyplot.show()
In the above we see the edge states of the quantum spin Hall effect, which
we can visualize using `kwant.plotter.map`:
-.. literalinclude:: /code/include/discretize.py
- :start-after: #HIDDEN_BEGIN_plot_qsh_wf
- :end-before: #HIDDEN_END_plot_qsh_wf
+.. jupyter-execute::
+
+ # get scattering wave functions at E=0
+ wf = kwant.wave_function(syst, energy=0, params=params)
-.. image:: /code/figure/discretizer_qsh_wf.*
+ # prepare density operators
+ sigma_z = np.array([[1, 0], [0, -1]])
+ prob_density = kwant.operator.Density(syst, np.eye(4))
+ spin_density = kwant.operator.Density(syst, np.kron(sigma_z, np.eye(2)))
+
+ # calculate expectation values and plot them
+ wf_sqr = sum(prob_density(psi) for psi in wf(0)) # states from left lead
+ rho_sz = sum(spin_density(psi) for psi in wf(0)) # states from left lead
+
+ fig, (ax1, ax2) = pyplot.subplots(1, 2, sharey=True, figsize=(16, 4))
+ kwant.plotter.map(syst, wf_sqr, ax=ax1)
+ kwant.plotter.map(syst, rho_sz, ax=ax2)
+
+ ax = ax1
+ im = [obj for obj in ax.get_children()
+ if isinstance(obj, mpl.image.AxesImage)][0]
+ fig.colorbar(im, ax=ax)
+
+ ax = ax2
+ im = [obj for obj in ax.get_children()
+ if isinstance(obj, mpl.image.AxesImage)][0]
+ fig.colorbar(im, ax=ax)
+
+ ax1.set_title('Probability density')
+ ax2.set_title('Spin density')
+ pyplot.show()
Limitations of discretization
@@ -238,29 +350,62 @@ example. Let us start from the continuum Hamiltonian
We start by defining this model as a string and setting the value of the
:math:`α` parameter:
-.. literalinclude:: /code/include/discretize.py
- :start-after: #HIDDEN_BEGIN_ls_def
- :end-before: #HIDDEN_END_ls_def
+.. jupyter-execute::
+
+ hamiltonian = "k_x**2 * identity(2) + alpha * k_x * sigma_y"
+ params = dict(alpha=.5)
Now we can use `kwant.continuum.lambdify` to obtain a function that computes
:math:`H(k)`:
-.. literalinclude:: /code/include/discretize.py
- :start-after: #HIDDEN_BEGIN_ls_hk_cont
- :end-before: #HIDDEN_END_ls_hk_cont
+.. jupyter-execute::
+
+ h_k = kwant.continuum.lambdify(hamiltonian, locals=params)
+ k_cont = np.linspace(-4, 4, 201)
+ e_cont = [scipy.linalg.eigvalsh(h_k(k_x=ki)) for ki in k_cont]
We can also construct a discretized approximation using
`kwant.continuum.discretize`, in a similar manner to previous examples:
-.. literalinclude:: /code/include/discretize.py
- :start-after: #HIDDEN_BEGIN_ls_hk_tb
- :end-before: #HIDDEN_END_ls_hk_tb
+.. jupyter-execute::
+
+ def plot(ax, a=1):
+ template = kwant.continuum.discretize(hamiltonian, grid=a)
+ syst = kwant.wraparound.wraparound(template).finalized()
+
+ def h_k(k_x):
+ p = dict(k_x=k_x, **params)
+ return syst.hamiltonian_submatrix(params=p)
+
+ k_tb = np.linspace(-np.pi/a, np.pi/a, 201)
+ e_tb = [scipy.linalg.eigvalsh(h_k(k_x=a*ki)) for ki in k_tb]
+
+ ax.plot(k_cont, e_cont, 'r-')
+ ax.plot(k_tb, e_tb, 'k-')
+
+ ax.plot([], [], 'r-', label='continuum')
+ ax.plot([], [], 'k-', label='tight-binding')
+
+ ax.set_xlim(-4, 4)
+ ax.set_ylim(-1, 14)
+ ax.set_title('a={}'.format(a))
+
+ ax.set_xlabel('momentum [a.u.]')
+ ax.set_ylabel('energy [a.u.]')
+ ax.grid()
+ ax.legend()
Below we can see the continuum and tight-binding dispersions for two
different values of the discretization grid spacing :math:`a`:
-.. image:: /code/figure/discretizer_lattice_spacing.*
+.. jupyter-execute::
+ :hide-code:
+
+ _, (ax1, ax2) = pyplot.subplots(1, 2, sharey=True, figsize=(12, 4))
+ plot(ax1, a=1)
+ plot(ax2, a=.25)
+ pyplot.show()
We clearly see that the smaller grid spacing is, the better we approximate
the original continuous dispersion. It is also worth remembering that the
@@ -284,20 +429,26 @@ It is possible to use ``identity`` (for identity matrix), ``kron`` (for Kronecke
expressions involving matrices. Matrices can also be provided explicitly using
square ``[]`` brackets. For example, all following expressions are equivalent:
-.. literalinclude:: /code/include/discretize.py
- :start-after: #HIDDEN_BEGIN_subs_1
- :end-before: #HIDDEN_END_subs_1
+.. jupyter-execute::
+
+ sympify = kwant.continuum.sympify
+ subs = {'sx': [[0, 1], [1, 0]], 'sz': [[1, 0], [0, -1]]}
+
+ e = (
+ sympify('[[k_x**2, alpha * k_x], [k_x * alpha, -k_x**2]]'),
+ sympify('k_x**2 * sigma_z + alpha * k_x * sigma_x'),
+ sympify('k_x**2 * sz + alpha * k_x * sx', locals=subs),
+ )
-.. literalinclude:: /code/figure/discretizer_subs_1.txt
+ print(e[0] == e[1] == e[2])
We can use the ``locals`` keyword parameter to substitute expressions
and numerical values:
-.. literalinclude:: /code/include/discretize.py
- :start-after: #HIDDEN_BEGIN_subs_2
- :end-before: #HIDDEN_END_subs_2
+.. jupyter-execute::
-.. literalinclude:: /code/figure/discretizer_subs_2.txt
+ subs = {'A': 'A(x) + B', 'V': 'V(x) + V_0', 'C': 5}
+ print(sympify('k_x * A * k_x + V + C', locals=subs))
Symbolic expressions obtained in this way can be directly passed to all
``discretizer`` functions.
diff --git a/doc/source/tutorial/faq.rst b/doc/source/tutorial/faq.rst
index c2ba0700dbc8d3ef4aaf3c2665fbf4ed77a944eb..ec52448d7afaf0d3ee4a2440088829a8dfc1b921 100644
--- a/doc/source/tutorial/faq.rst
+++ b/doc/source/tutorial/faq.rst
@@ -5,6 +5,17 @@ questions that are discussed there. The `Kwant paper
<https://downloads.kwant-project.org/doc/kwant-paper.pdf>`_ also digs deeper
into Kwant's structure.
+.. jupyter-execute::
+ :hide-code:
+
+ import kwant
+ import numpy as np
+ import tinyarray
+ import matplotlib
+ from matplotlib import pyplot as plt
+
+.. jupyter-execute:: boilerplate.py
+ :hide-code:
What is a system, and what is a builder?
========================================
@@ -48,11 +59,16 @@ combination of family and tag uniquely defines a site.
For example let us create an empty tight binding system and add two sites:
-.. literalinclude:: /code/include/faq.py
- :start-after: #HIDDEN_BEGIN_site
- :end-before: #HIDDEN_END_site
+.. jupyter-execute::
+
+ a = 1
+ lat = kwant.lattice.square(a)
+ syst = kwant.Builder()
-.. image:: /code/figure/faq_site.*
+ syst[lat(1, 0)] = 4
+ syst[lat(1, 1)] = 4
+
+ kwant.plot(syst);
In the above snippet we added 2 sites: ``lat(1, 0)`` and ``lat(0, 1)``. Both
of these sites belong to the same family, ``lat``, but have different tags:
@@ -79,11 +95,11 @@ tuples, for example lists, are not treated as hoppings.
Starting from the example code from `What is a site?`_, we can add a hopping
to our system in the following way:
-.. literalinclude:: /code/include/faq.py
- :start-after: #HIDDEN_BEGIN_hopping
- :end-before: #HIDDEN_END_hopping
+.. jupyter-execute::
+
+ syst[(lat(1, 0), lat(1, 1))] = 1j
-.. image:: /code/figure/faq_hopping.*
+ kwant.plot(syst);
Visually, a hopping is represented as a line that joins two sites.
@@ -139,18 +155,29 @@ Let's create two monatomic lattices (``lat_a`` and ``lat_b``). ``(1, 0)``
and ``(0, 1)`` will be the primitive vectors and ``(0, 0)`` and ``(0.5, 0.5)``
the origins of the two lattices:
-.. literalinclude:: /code/include/faq.py
- :start-after: #HIDDEN_BEGIN_lattice_monatomic
- :end-before: #HIDDEN_END_lattice_monatomic
+.. jupyter-execute::
+
+ # Two monatomic lattices
+ primitive_vectors = [(1, 0), (0, 1)]
+ lat_a = kwant.lattice.Monatomic(primitive_vectors, offset=(0, 0))
+ lat_b = kwant.lattice.Monatomic(primitive_vectors, offset=(0.5, 0.5))
+ # lat1 is equivalent to kwant.lattice.square()
-.. image:: /code/figure/faq_lattice.*
+ syst = kwant.Builder()
+
+ syst[lat_a(0, 0)] = 4
+ syst[lat_b(0, 0)] = 4
+
+ kwant.plot(syst);
We can also create a ``Polyatomic`` lattice with the same primitive vectors and
two sites in the basis:
-.. literalinclude:: /code/include/faq.py
- :start-after: #HIDDEN_BEGIN_lattice_polyatomic
- :end-before: #HIDDEN_END_lattice_polyatomic
+.. jupyter-execute::
+
+ # One polyatomic lattice containing two sublattices
+ lat = kwant.lattice.Polyatomic([(1, 0), (0, 1)], [(0, 0), (0.5, 0.5)])
+ sub_a, sub_b = lat.sublattices
The two sublattices ``sub_a`` and ``sub_b`` are nothing else than ``Monatomic``
instances, and are equivalent to ``lat_a`` and ``lat_b`` that we created
@@ -169,17 +196,40 @@ When plotting, how to color the different sublattices differently?
==================================================================
In the following example we shall use a kagome lattice, which has three sublattices.
-.. literalinclude:: /code/include/faq.py
- :start-after: #HIDDEN_BEGIN_colors1
- :end-before: #HIDDEN_END_colors1
+.. jupyter-execute::
+
+ lat = kwant.lattice.kagome()
+ syst = kwant.Builder()
+
+ a, b, c = lat.sublattices # The kagome lattice has 3 sublattices
As we can see below, we create a new plotting function that assigns a color for each family, and a different size for the hoppings depending on the family of the two sites. Finally we add sites and hoppings to our system and plot it with the new function.
-.. literalinclude:: /code/include/faq.py
- :start-after: #HIDDEN_BEGIN_colors2
- :end-before: #HIDDEN_END_colors2
+.. jupyter-execute::
-.. image:: /code/figure/faq_colors.*
+ # Plot sites from different families in different colors
+ def family_color(site):
+ if site.family == a:
+ return 'red'
+ if site.family == b:
+ return 'green'
+ else:
+ return 'blue'
+
+ def plot_system(syst):
+ kwant.plot(syst, site_lw=0.1, site_color=family_color)
+
+ ## Add sites and hoppings.
+ for i in range(4):
+ for j in range (4):
+ syst[a(i, j)] = 4
+ syst[b(i, j)] = 4
+ syst[c(i, j)] = 4
+
+ syst[lat.neighbors()] = -1
+
+ ## Plot the system.
+ plot_system(syst)
How to create many similar hoppings in one go?
@@ -197,31 +247,64 @@ integers).
The following example shows how this can be used:
-.. literalinclude:: /code/include/faq.py
- :start-after: #HIDDEN_BEGIN_direction1
- :end-before: #HIDDEN_END_direction1
+.. jupyter-execute::
+
+ # Create hopping between neighbors with HoppingKind
+ a = 1
+ syst = kwant.Builder()
+ lat = kwant.lattice.square(a)
+ syst[ (lat(i, j) for i in range(5) for j in range(5)) ] = 4
-.. image:: /code/figure/faq_direction1.*
+ syst[kwant.builder.HoppingKind((1, 0), lat)] = -1
+ kwant.plot(syst);
Note that ``HoppingKind`` only works with site families so you cannot use
them directly with ``Polyatomic`` lattices; you have to explicitly specify
the sublattices when creating a ``HoppingKind``:
-.. literalinclude:: /code/include/faq.py
- :start-after: #HIDDEN_BEGIN_direction2
- :end-before: #HIDDEN_END_direction2
+.. jupyter-execute::
+ :hide-code:
+
+ lat = kwant.lattice.kagome()
+ syst = kwant.Builder()
+
+ a, b, c = lat.sublattices # The kagome lattice has 3 sublattices
+
+ for i in range(4):
+ for j in range (4):
+ syst[a(i, j)] = 4
+ syst[b(i, j)] = 4
+ syst[c(i, j)] = 4
+
+.. jupyter-execute::
+
+ # equivalent to syst[kwant.builder.HoppingKind((0, 1), b)] = -1
+ syst[kwant.builder.HoppingKind((0, 1), b, b)] = -1
Here, we want the hoppings between the sites from sublattice b with a direction of (0,1) in the lattice coordinates.
-.. image:: /code/figure/faq_direction2.*
+.. jupyter-execute::
+ :hide-code:
-.. literalinclude:: /code/include/faq.py
- :start-after: #HIDDEN_BEGIN_direction3
- :end-before: #HIDDEN_END_direction3
+ plot_system(syst)
+
+.. jupyter-execute::
+ :hide-code:
+
+ # Delete the hoppings previously created
+ del syst[kwant.builder.HoppingKind((0, 1), b, b)]
+
+.. jupyter-execute::
+
+ syst[kwant.builder.HoppingKind((0, 0), a, b)] = -1
+ syst[kwant.builder.HoppingKind((0, 0), a, c)] = -1
+ syst[kwant.builder.HoppingKind((0, 0), c, b)] = -1
Here, we create hoppings between the sites of the same lattice coordinates but from different families.
-.. image:: /code/figure/faq_direction3.*
+.. jupyter-execute::
+
+ plot_system(syst)
How to set the hoppings between adjacent sites?
@@ -230,32 +313,53 @@ How to set the hoppings between adjacent sites?
that returns a list of ``HoppingKind`` instances that connect sites with their
(n-nearest) neighors:
-.. literalinclude:: /code/include/faq.py
- :start-after: #HIDDEN_BEGIN_adjacent1
- :end-before: #HIDDEN_END_adjacent1
+.. jupyter-execute::
+
+ # Create hoppings with lat.neighbors()
+ syst = kwant.Builder()
+ lat = kwant.lattice.square()
+ syst[(lat(i, j) for i in range(3) for j in range(3))] = 4
-.. image:: /code/figure/faq_adjacent1.*
-.. image:: /code/figure/faq_adjacent2.*
+ syst[lat.neighbors()] = -1 # Equivalent to lat.neighbors(1)
+ kwant.plot(syst);
-As we can see in the figure above, ``lat.neighbors()`` (on the left) returns the hoppings between the first nearest neighbors and ``lat.neighbors(2)`` (on the right) returns the hoppings between the second nearest neighbors.
+.. jupyter-execute::
+
+ del syst[lat.neighbors()] # Delete all nearest-neighbor hoppings
+ syst[lat.neighbors(2)] = -1
+
+ kwant.plot(syst);
+
+As we can see in the figures above, ``lat.neighbors()`` returns the hoppings between the first nearest neighbors and ``lat.neighbors(2)`` returns the hoppings between the second nearest neighbors.
When using a ``Polyatomic`` lattice ``neighbors()`` knows about the different
sublattices:
-.. literalinclude:: /code/include/faq.py
- :start-after: #HIDDEN_BEGIN_adjacent2
- :end-before: #HIDDEN_END_adjacent2
+.. jupyter-execute::
+
+ # Create the system
+ lat = kwant.lattice.kagome()
+ syst = kwant.Builder()
+ a, b, c = lat.sublattices # The kagome lattice has 3 sublattices
+
+ for i in range(4):
+ for j in range (4):
+ syst[a(i, j)] = 4 # red
+ syst[b(i, j)] = 4 # green
+ syst[c(i, j)] = 4 # blue
-.. image:: /code/figure/faq_adjacent3.*
+ syst[lat.neighbors()] = -1
+
+ plot_system(syst)
However, if we use the ``neighbors()`` method of a single sublattice, we will
only get the neighbors *on that sublattice*:
-.. literalinclude:: /code/include/faq.py
- :start-after: #HIDDEN_BEGIN_adjacent3
- :end-before: #HIDDEN_END_adjacent3
+.. jupyter-execute::
-.. image:: /code/figure/faq_adjacent4.*
+ del syst[lat.neighbors()] # Delete the hoppings previously created
+ syst[a.neighbors()] = -1
+ plot_system(syst)
Note in the above that there are *only* hoppings between the red sites. This
is an artifact of the visualisation: the blue and green sites just happen to lie
@@ -268,12 +372,35 @@ To make a hole in the system, use ``del syst[site]``, just like with any other
mapping. In the following example we remove all sites inside some "hole"
region:
-.. literalinclude:: /code/include/faq.py
- :start-after: #HIDDEN_BEGIN_hole
- :end-before: #HIDDEN_END_hole
+.. jupyter-execute::
+
+ # Define the lattice and the (empty) system
+ a = 2
+ lat = kwant.lattice.cubic(a)
+ syst = kwant.Builder()
+
+ L = 10
+ W = 10
+ H = 2
+
+ # Add sites to the system in a cuboid
+
+ syst[(lat(i, j, k) for i in range(L) for j in range(W) for k in range(H))] = 4
+ kwant.plot(syst);
+
+.. jupyter-execute::
+
+ # Delete sites to create a hole
+
+ def in_hole(site):
+ x, y, z = site.pos / a - (L/2, W/2, H/2) # position relative to centre
+ return abs(x) < L / 4 and abs(y) < W / 4
+
+ for site in filter(in_hole, list(syst.sites())):
+ del syst[site]
+
+ kwant.plot(syst);
-.. image:: /code/figure/faq_hole1.*
-.. image:: /code/figure/faq_hole2.*
``del syst[site]`` also works after hoppings have been added to the system.
If a site is deleted, then all the hoppings to/from that site are also deleted.
@@ -287,9 +414,18 @@ what is a builder?`_ if you do not know the difference).
We can access the sites of a ``Builder`` by using its `~kwant.builder.Builder.sites` method:
-.. literalinclude:: /code/include/faq.py
- :start-after: #HIDDEN_BEGIN_sites1
- :end-before: #HIDDEN_END_sites1
+.. jupyter-execute::
+ :hide-code:
+
+ builder = kwant.Builder()
+ lat = kwant.lattice.square()
+ builder[(lat(i, j) for i in range(3) for j in range(3))] = 4
+
+.. jupyter-execute::
+
+ # Before finalizing the system
+
+ sites = list(builder.sites()) # sites() doe *not* return a list
The ``sites()`` method returns an *iterator* over the system sites, and in the
above example we create a list from the contents of this iterator, which
@@ -300,9 +436,11 @@ well be returned in a different order.
After finalization, when we are dealing with a ``System``, the sites themselves
are stored in a list, which can be accessed via the ``sites`` attribute:
-.. literalinclude:: /code/include/faq.py
- :start-after: #HIDDEN_BEGIN_sites2
- :end-before: #HIDDEN_END_sites2
+.. jupyter-execute::
+
+ # After finalizing the system
+ syst = builder.finalized()
+ sites = syst.sites # syst.sites is an actual list
The order of sites in a ``System`` is fixed, and also defines the ordering of
the system Hamiltonian, system wavefunctions etc. (see `How does Kwant order components of an individual wavefunction?`_ for details).
@@ -310,9 +448,9 @@ the system Hamiltonian, system wavefunctions etc. (see `How does Kwant order com
``System`` also contains the inverse mapping, ``id_by_site`` which gives us
the index of a given site within the system:
-.. literalinclude:: /code/include/faq.py
- :start-after: #HIDDEN_BEGIN_sites3
- :end-before: #HIDDEN_END_sites3
+.. jupyter-execute::
+
+ i = syst.id_by_site[lat(0, 2)] # we want the id of the site lat(0, 2)
How to use different lattices for the scattering region and a lead?
@@ -322,19 +460,35 @@ which we want to connect to leads that contain sites from a square lattice.
First we construct the central system:
-.. literalinclude:: /code/include/faq.py
- :start-after: #HIDDEN_BEGIN_different_lattice1
- :end-before: #HIDDEN_END_different_lattice1
+.. jupyter-execute::
-.. image:: /code/figure/faq_different_lattice1.*
+ # Define the scattering Region
+ L = 5
+ W = 5
+
+ lat = kwant.lattice.honeycomb()
+ subA, subB = lat.sublattices
+
+ syst = kwant.Builder()
+ syst[(subA(i, j) for i in range(L) for j in range(W))] = 4
+ syst[(subB(i, j) for i in range(L) for j in range(W))] = 4
+ syst[lat.neighbors()] = -1
+
+ kwant.plot(syst);
and the lead:
-.. literalinclude:: /code/include/faq.py
- :start-after: #HIDDEN_BEGIN_different_lattice2
- :end-before: #HIDDEN_END_different_lattice2
+.. jupyter-execute::
+
+ # Create a lead
+ lat_lead = kwant.lattice.square()
+ sym_lead1 = kwant.TranslationalSymmetry((0, 1))
+
+ lead1 = kwant.Builder(sym_lead1)
+ lead1[(lat_lead(i, 0) for i in range(2, 7))] = 4
+ lead1[lat_lead.neighbors()] = -1
-.. image:: /code/figure/faq_different_lattice2.*
+ kwant.plot(syst);
We cannot simply use `~kwant.builder.Builder.attach_lead` to attach this lead to the
system with the honeycomb lattice because Kwant does not know how sites from
@@ -343,19 +497,22 @@ these two lattices should be connected.
We must first add a layer of sites from the square lattice to the system and manually
add the hoppings from these sites to the sites from the honeycomb lattice:
-.. literalinclude:: /code/include/faq.py
- :start-after: #HIDDEN_BEGIN_different_lattice3
- :end-before: #HIDDEN_END_different_lattice3
+.. jupyter-execute::
-.. image:: /code/figure/faq_different_lattice3.*
+ syst[(lat_lead(i, 5) for i in range(2, 7))] = 4
+ syst[lat_lead.neighbors()] = -1
+
+ # Manually attach sites from graphene to square lattice
+ syst[((lat_lead(i+2, 5), subB(i, 4)) for i in range(5))] = -1
+
+ kwant.plot(syst);
``attach_lead()`` will now be able to attach the lead:
-.. literalinclude:: /code/include/faq.py
- :start-after: #HIDDEN_BEGIN_different_lattice4
- :end-before: #HIDDEN_END_different_lattice4
+.. jupyter-execute::
-.. image:: /code/figure/faq_different_lattice4.*
+ syst.attach_lead(lead1)
+ kwant.plot(syst);
How to cut a finite system out of a system with translational symmetries?
@@ -369,30 +526,50 @@ will be repeated in the desired shape to create the final system.
For example, say we want to create a simple model on a cubic lattice:
-.. literalinclude:: /code/include/faq.py
- :start-after: #HIDDEN_BEGIN_fill1
- :end-before: #HIDDEN_END_fill1
+.. jupyter-execute::
+
+ # Create 3d model.
+ cubic = kwant.lattice.cubic()
+ sym_3d = kwant.TranslationalSymmetry([1, 0, 0], [0, 1, 0], [0, 0, 1])
+ model = kwant.Builder(sym_3d)
+ model[cubic(0, 0, 0)] = 4
+ model[cubic.neighbors()] = -1
We have now created our "template" ``Builder`` which has 3 translational
symmetries. Next we will fill a system with no translational symmetries with
sites and hoppings from the template inside a cuboid:
-.. literalinclude:: /code/include/faq.py
- :start-after: #HIDDEN_BEGIN_fill2
- :end-before: #HIDDEN_END_fill2
+.. jupyter-execute::
-.. image:: /code/figure/faq_fill2.*
+ # Build scattering region (white).
+ def cuboid_shape(site):
+ x, y, z = abs(site.pos)
+ return x < 4 and y < 10 and z < 3
+
+ cuboid = kwant.Builder()
+ cuboid.fill(model, cuboid_shape, (0, 0, 0));
+ kwant.plot(cuboid);
We can then use the original template to create a lead, which has 1 translational
symmetry. We can then use this lead as a template to fill another section of
the system with a cylinder of sites and hoppings:
-.. literalinclude:: /code/include/faq.py
- :start-after: #HIDDEN_BEGIN_fill3
- :end-before: #HIDDEN_END_fill3
+.. jupyter-execute::
+
+ # Build electrode (black).
+ def electrode_shape(site):
+ x, y, z = site.pos - (0, 5, 2)
+ return y**2 + z**2 < 2.3**2
-.. image:: /code/figure/faq_fill3.*
+ electrode = kwant.Builder(kwant.TranslationalSymmetry([1, 0, 0]))
+ electrode.fill(model, electrode_shape, (0, 5, 2)) # lead
+ # Scattering region
+ cuboid.fill(electrode, lambda s: abs(s.pos[0]) < 7, (0, 5, 4))
+
+ cuboid.attach_lead(electrode)
+
+ kwant.plot(cuboid);
How does Kwant order the propagating modes of a lead?
=====================================================
@@ -402,16 +579,44 @@ achieved with the `~kwant.system.InfiniteSystem.modes` method, which returns a
pair of objects, the first of which contains the propagating modes of the
system in a `~kwant.physics.PropagatingModes` object:
-.. literalinclude:: /code/include/faq.py
- :start-after: #HIDDEN_BEGIN_pm
- :end-before: #HIDDEN_END_pm
+.. jupyter-execute::
+
+ lat = kwant.lattice.square()
+
+ lead = kwant.Builder(kwant.TranslationalSymmetry((-1, 0)))
+ lead[(lat(0, i) for i in range(3))] = 4
+ lead[lat.neighbors()] = -1
+
+ flead = lead.finalized()
+
+ E = 2.5
+ prop_modes, _ = flead.modes(energy=E)
``PropagatingModes`` contains the wavefunctions, velocities and momenta of the
modes at the requested energy (2.5 in this example). In order to understand
the order in which these quantities are returned it is often useful to look at
the a section of the band structure for the system in question:
-.. image:: /code/figure/faq_pm1.*
+.. jupyter-execute::
+ :hide-code:
+
+ def plot_and_label_modes(lead, E):
+ # Plot the different modes
+ pmodes, _ = lead.modes(energy=E)
+ kwant.plotter.bands(lead, show=False)
+ for i, k in enumerate(pmodes.momenta):
+ plt.plot(k, E, 'ko')
+ plt.annotate(str(i), xy=(k, E), xytext=(-5, 8),
+ textcoords='offset points',
+ bbox=dict(boxstyle='round,pad=0.1',fc='white', alpha=0.7))
+ plt.plot([-3, 3], [E, E], 'r--')
+ plt.ylim(E-1, E+1)
+ plt.xlim(-2, 2)
+ plt.xlabel("momentum")
+ plt.ylabel("energy")
+ plt.show()
+
+ plot_and_label_modes(flead, E)
On the above band structure we have labelled the 4 modes in the order
that they appear in the output of ``modes()`` at energy 2.5. Note that
@@ -425,7 +630,20 @@ the modes are sorted in the following way:
For more complicated systems and band structures this can lead to some
unintuitive orderings:
-.. image:: /code/figure/faq_pm2.*
+.. jupyter-execute::
+ :hide-code:
+
+ s0 = np.eye(2)
+ sz = np.array([[1, 0], [0, -1]])
+
+ lead2 = kwant.Builder(kwant.TranslationalSymmetry((-1, 0)))
+
+ lead2[(lat(0, i) for i in range(2))] = np.diag([1.8, -1])
+ lead2[lat.neighbors()] = -1 * sz
+
+ flead2 = lead2.finalized()
+
+ plot_and_label_modes(flead2, 1)
How does Kwant order scattering states?
@@ -452,10 +670,38 @@ When all the site families present in a system have only 1 degree of freedom
per site (i.e. the all the onsites are scalars) then the index into a
wavefunction defined over the system is exactly the site index:
-.. literalinclude:: /code/include/faq.py
- :start-after: #HIDDEN_BEGIN_ord1
- :end-before: #HIDDEN_END_ord1
-.. literalinclude:: /code/figure/faq_ord1.txt
+.. jupyter-execute::
+ :hide-code:
+
+ def circle(R):
+ return lambda r: np.linalg.norm(r) < R
+
+
+ def make_system(lat):
+ norbs = lat.norbs
+ syst = kwant.Builder()
+ syst[lat.shape(circle(3), (0, 0))] = 4 * np.eye(norbs)
+ syst[lat.neighbors()] = -1 * np.eye(norbs)
+
+ lead = kwant.Builder(kwant.TranslationalSymmetry((-1, 0)))
+ lead[(lat(0, i) for i in range(-1, 2))] = 4 * np.eye(norbs)
+ lead[lat.neighbors()] = -1 * np.eye(norbs)
+
+ syst.attach_lead(lead)
+ syst.attach_lead(lead.reversed())
+
+ return syst.finalized()
+
+.. jupyter-execute::
+
+ lat = kwant.lattice.square(norbs=1)
+ syst = make_system(lat)
+ scattering_states = kwant.wave_function(syst, energy=1)
+ wf = scattering_states(0)[0] # scattering state from lead 0 incoming in mode 0
+
+ idx = syst.id_by_site[lat(0, 0)] # look up index of site
+
+ print('wavefunction on lat(0, 0): ', wf[idx])
We see that the wavefunction on a single site is a single complex number, as
expected.
@@ -467,10 +713,20 @@ to one another. In the case where all site families in the system have the
wavefunction into a matrix, where the row number indexes the site, and the
column number indexes the degree of freedom on that site:
-.. literalinclude:: /code/include/faq.py
- :start-after: #HIDDEN_BEGIN_ord2
- :end-before: #HIDDEN_END_ord2
-.. literalinclude:: /code/figure/faq_ord2.txt
+.. jupyter-execute::
+
+ lat = kwant.lattice.square(norbs=2)
+ syst = make_system(lat)
+ scattering_states = kwant.wave_function(syst, energy=1)
+ wf = scattering_states(0)[0] # scattering state from lead 0 incoming in mode 0
+
+ idx = syst.id_by_site[lat(0, 0)] # look up index of site
+
+ # Group consecutive degrees of freedom from 'wf' together; these correspond
+ # to degrees of freedom on the same site.
+ wf = wf.reshape(-1, 2)
+
+ print('wavefunction on lat(0, 0): ', wf[idx])
We see that the wavefunction on a single site is a *vector* of 2 complex numbers,
as we expect.
diff --git a/doc/source/tutorial/first_steps.rst b/doc/source/tutorial/first_steps.rst
index 1ddce048c4da222856413471aeb452390ed3e7c7..f509aaa0c262be17a1eaa53b8551812e06e88f22 100644
--- a/doc/source/tutorial/first_steps.rst
+++ b/doc/source/tutorial/first_steps.rst
@@ -13,7 +13,7 @@ Schrödinger equation
.. math::
H = \frac{-\hbar^2}{2m}(\partial_x^2 + \partial_y^2) + V(y)
-with a hard-wall confinement :math:`V(y)` in y-direction.
+with a hard-wall confinement :math:`V(y)` in the y-direction.
To be able to implement the quantum wire with Kwant, the continuous Hamiltonian
:math:`H` has to be discretized thus turning it into a tight-binding
@@ -62,23 +62,50 @@ Transport through a quantum wire
.. seealso::
The complete source code of this example can be found in
- :download:`quantum_wire.py </code/download/quantum_wire.py>`
+ :jupyter-download:script:`quantum_wire`
In order to use Kwant, we need to import it:
-.. literalinclude:: /code/include/quantum_wire.py
- :start-after: #HIDDEN_BEGIN_dwhx
- :end-before: #HIDDEN_END_dwhx
+.. jupyter-kernel::
+ :id: quantum_wire
+
+.. jupyter-execute::
+ :hide-code:
+
+ # Tutorial 2.2.2. Transport through a quantum wire
+ # ================================================
+ #
+ # Physics background
+ # ------------------
+ # Conductance of a quantum wire; subbands
+ #
+ # Kwant features highlighted
+ # --------------------------
+ # - Builder for setting up transport systems easily
+ # - Making scattering region and leads
+ # - Using the simple sparse solver for computing Landauer conductance
+
+ # For plotting
+ from matplotlib import pyplot
+
+.. jupyter-execute:: boilerplate.py
+ :hide-code:
+
+.. jupyter-execute::
+
+ import kwant
Enabling Kwant is as easy as this [#]_ !
-The first step is now the definition of the system with scattering region and
+The first step is now to define the system with scattering region and
leads. For this we make use of the `~kwant.builder.Builder` type that allows to
define a system in a convenient way. We need to create an instance of it:
-.. literalinclude:: /code/include/quantum_wire.py
- :start-after: #HIDDEN_BEGIN_goiq
- :end-before: #HIDDEN_END_goiq
+.. jupyter-execute::
+
+ # First define the tight-binding system
+
+ syst = kwant.Builder()
Observe that we just accessed `~kwant.builder.Builder` by the name
``kwant.Builder``. We could have just as well written
@@ -92,9 +119,10 @@ Apart from `~kwant.builder.Builder` we also need to specify
what kind of sites we want to add to the system. Here we work with
a square lattice. For simplicity, we set the lattice constant to unity:
-.. literalinclude:: /code/include/quantum_wire.py
- :start-after: #HIDDEN_BEGIN_suwo
- :end-before: #HIDDEN_END_suwo
+.. jupyter-execute::
+
+ a = 1
+ lat = kwant.lattice.square(a)
Since we work with a square lattice, we label the points with two
integer coordinates `(i, j)`. `~kwant.builder.Builder` then
@@ -111,9 +139,25 @@ needed in Builder (more about that in the technical details below).
We now build a rectangular scattering region that is `W`
lattice points wide and `L` lattice points long:
-.. literalinclude:: /code/include/quantum_wire.py
- :start-after: #HIDDEN_BEGIN_zfvr
- :end-before: #HIDDEN_END_zfvr
+.. jupyter-execute::
+
+ t = 1.0
+ W, L = 10, 30
+
+ # Define the scattering region
+
+ for i in range(L):
+ for j in range(W):
+ # On-site Hamiltonian
+ syst[lat(i, j)] = 4 * t
+
+ # Hopping in y-direction
+ if j > 0:
+ syst[lat(i, j), lat(i, j - 1)] = -t
+
+ # Hopping in x-direction
+ if i > 0:
+ syst[lat(i, j), lat(i - 1, j)] = -t
Observe how the above code corresponds directly to the terms of the discretized
Hamiltonian:
@@ -140,9 +184,14 @@ Next, we define the leads. Leads are also constructed using
`~kwant.builder.Builder`, but in this case, the
system must have a translational symmetry:
-.. literalinclude:: /code/include/quantum_wire.py
- :start-after: #HIDDEN_BEGIN_xcmc
- :end-before: #HIDDEN_END_xcmc
+.. jupyter-execute::
+
+ # Then, define and attach the leads:
+
+ # First the lead to the left
+ # (Note: TranslationalSymmetry takes a real-space vector)
+ sym_left_lead = kwant.TranslationalSymmetry((-a, 0))
+ left_lead = kwant.Builder(sym_left_lead)
Here, the `~kwant.builder.Builder` takes a
`~kwant.lattice.TranslationalSymmetry` as the optional parameter. Note that the
@@ -155,17 +204,22 @@ as the hoppings inside one unit cell and to the next unit cell of the lead.
For a square lattice, and a lead in y-direction the unit cell is
simply a vertical line of points:
-.. literalinclude:: /code/include/quantum_wire.py
- :start-after: #HIDDEN_BEGIN_ndez
- :end-before: #HIDDEN_END_ndez
+.. jupyter-execute::
+
+ for j in range(W):
+ left_lead[lat(0, j)] = 4 * t
+ if j > 0:
+ left_lead[lat(0, j), lat(0, j - 1)] = -t
+ left_lead[lat(1, j), lat(0, j)] = -t
Note that here it doesn't matter if you add the hoppings to the next or the
previous unit cell -- the translational symmetry takes care of that. The
isolated, infinite is attached at the correct position using
-.. literalinclude:: /code/include/quantum_wire.py
- :start-after: #HIDDEN_BEGIN_fskr
- :end-before: #HIDDEN_END_fskr
+.. jupyter-execute::
+ :hide-output:
+
+ syst.attach_lead(left_lead)
This call returns the lead number which will be used to refer to the lead when
computing transmissions (further down in this tutorial). More details about
@@ -175,9 +229,21 @@ We also want to add a lead on the right side. The only difference to
the left lead is that the vector of the translational
symmetry must point to the right, the remaining code is the same:
-.. literalinclude:: /code/include/quantum_wire.py
- :start-after: #HIDDEN_BEGIN_xhqc
- :end-before: #HIDDEN_END_xhqc
+.. jupyter-execute::
+ :hide-output:
+
+ # Then the lead to the right
+
+ sym_right_lead = kwant.TranslationalSymmetry((a, 0))
+ right_lead = kwant.Builder(sym_right_lead)
+
+ for j in range(W):
+ right_lead[lat(0, j)] = 4 * t
+ if j > 0:
+ right_lead[lat(0, j), lat(0, j - 1)] = -t
+ right_lead[lat(1, j), lat(0, j)] = -t
+
+ syst.attach_lead(right_lead)
Note that here we added points with x-coordinate 0, just as for the left lead.
You might object that the right lead should be placed `L`
@@ -187,13 +253,9 @@ you do not need to worry about that.
Now we have finished building our system! We plot it, to make sure we didn't
make any mistakes:
-.. literalinclude:: /code/include/quantum_wire.py
- :start-after: #HIDDEN_BEGIN_wsgh
- :end-before: #HIDDEN_END_wsgh
+.. jupyter-execute::
-This should bring up this picture:
-
-.. image:: /code/figure/quantum_wire_syst.*
+ kwant.plot(syst);
The system is represented in the usual way for tight-binding systems:
dots represent the lattice points `(i, j)`, and for every
@@ -203,16 +265,29 @@ fading color.
In order to use our system for a transport calculation, we need to finalize it
-.. literalinclude:: /code/include/quantum_wire.py
- :start-after: #HIDDEN_BEGIN_dngj
- :end-before: #HIDDEN_END_dngj
+.. jupyter-execute::
+
+ # Finalize the system
+ syst = syst.finalized()
Having successfully created a system, we now can immediately start to compute
its conductance as a function of energy:
-.. literalinclude:: /code/include/quantum_wire.py
- :start-after: #HIDDEN_BEGIN_buzn
- :end-before: #HIDDEN_END_buzn
+.. jupyter-execute::
+
+ # Now that we have the system, we can compute conductance
+ energies = []
+ data = []
+ for ie in range(100):
+ energy = ie * 0.01
+
+ # compute the scattering matrix at a given energy
+ smatrix = kwant.smatrix(syst, energy)
+
+ # compute the transmission probability from lead 0 to
+ # lead 1
+ energies.append(energy)
+ data.append(smatrix.transmission(1, 0))
We use ``kwant.smatrix`` which is a short name for
`kwant.solvers.default.smatrix` of the default solver module
@@ -227,13 +302,15 @@ Finally we can use ``matplotlib`` to make a plot of the computed data
(although writing to file and using an external viewer such as
gnuplot or xmgrace is just as viable)
-.. literalinclude:: /code/include/quantum_wire.py
- :start-after: #HIDDEN_BEGIN_lliv
- :end-before: #HIDDEN_END_lliv
+.. jupyter-execute::
-This should yield the result
-
-.. image:: /code/figure/quantum_wire_result.*
+ # Use matplotlib to write output
+ # We should see conductance steps
+ pyplot.figure()
+ pyplot.plot(energies, data)
+ pyplot.xlabel("energy [t]")
+ pyplot.ylabel("conductance [e^2/h]")
+ pyplot.show()
We see a conductance quantized in units of :math:`e^2/h`,
increasing in steps as the energy is increased. The
@@ -333,7 +410,12 @@ Building the same system with less code
.. seealso::
The complete source code of this example can be found in
- :download:`quantum_wire_revisited.py </code/download/quantum_wire_revisited.py>`
+ :jupyter-download:script:`quantum_wire_revisited`
+
+
+.. jupyter-kernel::
+ :id: quantum_wire_revisited
+
Kwant allows for more than one way to build a system. The reason is that
`~kwant.builder.Builder` is essentially just a container that can be filled in
@@ -347,19 +429,56 @@ the code into separate entities. In this example we therefore also aim at more
structure.
We begin the program collecting all imports in the beginning of the
-file and put the build-up of the system into a separate function
-`make_system`:
+file and defining the a square lattice and empty scattering region.
+
+.. jupyter-execute::
+ :hide-code:
+
+ # Tutorial 2.2.3. Building the same system with less code
+ # =======================================================
+ #
+ # Physics background
+ # ------------------
+ # Conductance of a quantum wire; subbands
+ #
+ # Kwant features highlighted
+ # --------------------------
+ # - Using iterables and builder.HoppingKind for making systems
+ # - introducing `reversed()` for the leads
+ #
+ # Note: Does the same as tutorial1a.py, but using other features of Kwant.
+
+.. jupyter-execute::
+
+ import kwant
+
+ # For plotting
+ from matplotlib import pyplot
+
+.. jupyter-execute:: boilerplate.py
+ :hide-code:
+
+.. jupyter-execute::
+
+ a = 1
+ t = 1.0
+ W, L = 10, 30
+
+ # Start with an empty tight-binding system and a single square lattice.
+ # `a` is the lattice constant (by default set to 1 for simplicity).
+ lat = kwant.lattice.square(a)
+
+ syst = kwant.Builder()
-.. literalinclude:: /code/include/quantum_wire_revisited.py
- :start-after: #HIDDEN_BEGIN_xkzy
- :end-before: #HIDDEN_END_xkzy
Previously, the scattering region was build using two ``for``-loops.
Instead, we now write:
-.. literalinclude:: /code/include/quantum_wire_revisited.py
- :start-after: #HIDDEN_BEGIN_vvjt
- :end-before: #HIDDEN_END_vvjt
+
+.. jupyter-execute::
+
+ syst[(lat(x, y) for x in range(L) for y in range(W))] = 4 * t
+
Here, all lattice points are added at once in the first line. The
construct ``((i, j) for i in range(L) for j in range(W))`` is a
@@ -375,9 +494,9 @@ hoppings. In this case, an iterable like for the lattice
points becomes a bit cumbersome, and we use instead another
feature of Kwant:
-.. literalinclude:: /code/include/quantum_wire_revisited.py
- :start-after: #HIDDEN_BEGIN_nooi
- :end-before: #HIDDEN_END_nooi
+.. jupyter-execute::
+
+ syst[lat.neighbors()] = -t
In regular lattices, hoppings form large groups such that hoppings within a
group can be transformed into one another by lattice translations. In order to
@@ -397,49 +516,48 @@ values for all the nth-nearest neighbors at once, one can similarly use
The left lead is constructed in an analogous way:
-.. literalinclude:: /code/include/quantum_wire_revisited.py
- :start-after: #HIDDEN_BEGIN_iepx
- :end-before: #HIDDEN_END_iepx
+.. jupyter-execute::
+
+ lead = kwant.Builder(kwant.TranslationalSymmetry((-a, 0)))
+ lead[(lat(0, j) for j in range(W))] = 4 * t
+ lead[lat.neighbors()] = -t
The previous example duplicated almost identical code for the left and
the right lead. The only difference was the direction of the translational
symmetry vector. Here, we only construct the left lead, and use the method
`~kwant.builder.Builder.reversed` of `~kwant.builder.Builder` to obtain a copy
of a lead pointing in the opposite direction. Both leads are attached as
-before and the finished system returned:
+before:
-.. literalinclude:: /code/include/quantum_wire_revisited.py
- :start-after: #HIDDEN_BEGIN_yxot
- :end-before: #HIDDEN_END_yxot
+.. jupyter-execute::
+ :hide-output:
-The remainder of the script has been organized into two functions. One for the
-plotting of the conductance.
+ syst.attach_lead(lead)
+ syst.attach_lead(lead.reversed())
-.. literalinclude:: /code/include/quantum_wire_revisited.py
- :start-after: #HIDDEN_BEGIN_ayuk
- :end-before: #HIDDEN_END_ayuk
+The remainder of the script proceeds identically. We first finalize the system:
-And one ``main`` function.
+.. jupyter-execute::
-.. literalinclude:: /code/include/quantum_wire_revisited.py
- :start-after: #HIDDEN_BEGIN_cjel
- :end-before: #HIDDEN_END_cjel
+ syst = syst.finalized()
-Finally, we use the following standard Python construct [#]_ to execute
-``main`` if the program is used as a script (i.e. executed as
-``python quantum_wire_revisited.py``):
+and then calculate the transmission and plot:
-.. literalinclude:: /code/include/quantum_wire_revisited.py
- :start-after: #HIDDEN_BEGIN_ypbj
- :end-before: #HIDDEN_END_ypbj
+.. jupyter-execute::
-If the example, however, is imported inside Python using ``import
-quantum_wire_revisted as qw``, ``main`` is not executed automatically.
-Instead, you can execute it manually using ``qw.main()``. On the other
-hand, you also have access to the other functions, ``make_system`` and
-``plot_conductance``, and can thus play with the parameters.
+ energies = []
+ data = []
+ for ie in range(100):
+ energy = ie * 0.01
+ smatrix = kwant.smatrix(syst, energy)
+ energies.append(energy)
+ data.append(smatrix.transmission(1, 0))
-The result of the example should be identical to the previous one.
+ pyplot.figure()
+ pyplot.plot(energies, data)
+ pyplot.xlabel("energy [t]")
+ pyplot.ylabel("conductance [e^2/h]")
+ pyplot.show()
.. specialnote:: Technical details
@@ -453,11 +571,9 @@ The result of the example should be identical to the previous one.
For technical reasons it is not possible to add several points
using a tuple of sites. Hence it is worth noting
- a subtle detail in
+ a subtle detail in:
- .. literalinclude:: /code/include/quantum_wire_revisited.py
- :start-after: #HIDDEN_BEGIN_vvjt
- :end-before: #HIDDEN_END_vvjt
+ >>> syst[(lat(x, y) for x in range(L) for y in range(W))] = 4 * t
Note that ``(lat(x, y) for x in range(L) for y in range(W))`` is not
a tuple, but a generator.
@@ -512,6 +628,132 @@ The result of the example should be identical to the previous one.
it would be impossible to distinguish whether one would like to add two
separate sites, or one hopping.
+
+Tips for organizing your simulation scripts
+...........................................
+
+.. seealso::
+ The complete source code of this example can be found in
+ :jupyter-download:script:`quantum_wire_organized`
+
+
+.. jupyter-kernel::
+ :id: quantum_wire_organized
+
+.. jupyter-execute::
+ :hide-code:
+
+ # Tutorial 2.2.4. Organizing a simulation script
+ # ==============================================
+ #
+ # Physics background
+ # ------------------
+ # Conductance of a quantum wire; subbands
+ #
+ # Note: Does the same as quantum_write_revisited.py, but features
+ # better code organization
+
+
+The above two examples illustrate some of the core features of Kwant, however
+the code was presented in a style which is good for exposition, but which is
+bad for making your code understandable and reusable. In this example we will
+lay out some best practices for writing your own simulation scripts.
+
+In the above examples we constructed a single Kwant system, using global variables
+for parameters such as the lattice constant and the length and width of the system.
+Instead, it is preferable to create a *function* that you can call, and which will
+return a Kwant ``Builder``:
+
+.. jupyter-execute::
+
+ from matplotlib import pyplot
+ import kwant
+
+.. jupyter-execute:: boilerplate.py
+ :hide-code:
+
+.. jupyter-execute::
+
+ def make_system(L, W, a=1, t=1.0):
+ lat = kwant.lattice.square(a)
+
+ syst = kwant.Builder()
+ syst[(lat(i, j) for i in range(L) for j in range(W))] = 4 * t
+ syst[lat.neighbors()] = -t
+
+ lead = kwant.Builder(kwant.TranslationalSymmetry((-a, 0)))
+ lead[(lat(0, j) for j in range(W))] = 4 * t
+ lead[lat.neighbors()] = -t
+
+ syst.attach_lead(lead)
+ syst.attach_lead(lead.reversed())
+
+ return syst
+
+
+By encapsulating system creation within ``make_system`` we *document* our code
+by telling readers that *this* is how we create a system, and that creating a system
+depends on *these* parameters (the length and width of the system, in this case, as well
+as the lattice constant and the value for the hopping parameter). By defining a function
+we also ensure that we can consistently create different systems (e.g. of different sizes)
+of the same type (rectangular slab).
+
+We similarly encapsulate the part of the script that does computation and plotting into
+a function ``plot_conductance``:
+
+.. jupyter-execute::
+
+ def plot_conductance(syst, energies):
+ # Compute conductance
+ data = []
+ for energy in energies:
+ smatrix = kwant.smatrix(syst, energy)
+ data.append(smatrix.transmission(1, 0))
+
+ pyplot.figure()
+ pyplot.plot(energies, data)
+ pyplot.xlabel("energy [t]")
+ pyplot.ylabel("conductance [e^2/h]")
+ pyplot.show()
+
+And the ``main`` function that glues together the components that we previously defined:
+
+.. jupyter-execute::
+
+ def main():
+ syst = make_system(W=10, L=30)
+
+ # Check that the system looks as intended.
+ kwant.plot(syst)
+
+ # Finalize the system.
+ fsyst = syst.finalized()
+
+ # We should see conductance steps.
+ plot_conductance(fsyst, energies=[0.01 * i for i in range(100)])
+
+
+Finally, we use the following standard Python construct [#]_ to execute
+``main`` if the program is used as a script (i.e. executed as
+``python quantum_wire_organized.py``):
+
+
+.. jupyter-execute::
+
+ # 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()
+
+If the example, however, is imported inside Python using ``import
+quantum_wire_organized as qw``, ``main`` is not executed automatically.
+Instead, you can execute it manually using ``qw.main()``. On the other
+hand, you also have access to the other functions, ``make_system`` and
+``plot_conductance``, and can thus play with the parameters.
+
+The result of this example should be identical to the previous one.
+
+
.. rubric:: Footnotes
.. [#] https://docs.python.org/3/library/__main__.html
diff --git a/doc/source/tutorial/graphene.rst b/doc/source/tutorial/graphene.rst
index 63bcf5f29859ea09f11b84410c64871bf49fe9b5..6eaa93d6ed142f430bd078399b2850cdea61cf32 100644
--- a/doc/source/tutorial/graphene.rst
+++ b/doc/source/tutorial/graphene.rst
@@ -5,7 +5,7 @@ Beyond square lattices: graphene
.. seealso::
The complete source code of this example can be found in
- :download:`graphene.py </code/download/graphene.py>`
+ :jupyter-download:script:`graphene`
In the following example, we are going to calculate the
conductance through a graphene quantum dot with a p-n junction
@@ -18,9 +18,43 @@ We begin by defining the honeycomb lattice of graphene. This is
in principle already done in `kwant.lattice.honeycomb`, but we do it
explicitly here to show how to define a new lattice:
-.. literalinclude:: /code/include/graphene.py
- :start-after: #HIDDEN_BEGIN_hnla
- :end-before: #HIDDEN_END_hnla
+.. jupyter-kernel::
+ :id: graphene
+
+.. jupyter-execute::
+ :hide-code:
+
+ # Tutorial 2.5. Beyond square lattices: graphene
+ # ==============================================
+ #
+ # Physics background
+ # ------------------
+ # Transport through a graphene quantum dot with a pn-junction
+ #
+ # Kwant features highlighted
+ # --------------------------
+ # - Application of all the aspects of tutorials 1-3 to a more complicated
+ # lattice, namely graphene
+
+ from math import pi, sqrt, tanh
+
+ from matplotlib import pyplot
+
+ import kwant
+
+ # For computing eigenvalues
+ import scipy.sparse.linalg as sla
+
+ sin_30, cos_30 = (1 / 2, sqrt(3) / 2)
+
+.. jupyter-execute:: boilerplate.py
+ :hide-code:
+
+.. jupyter-execute::
+
+ graphene = kwant.lattice.general([(1, 0), (sin_30, cos_30)],
+ [(0, 0), (0, 1 / sqrt(3))])
+ a, b = graphene.sublattices
The first argument to the `~kwant.lattice.general` function is the list of
primitive vectors of the lattice; the second one is the coordinates of basis
@@ -31,9 +65,30 @@ itself forms a regular lattice of the same type as well, and those
In the next step we define the shape of the scattering region (circle again)
and add all lattice points using the ``shape``-functionality:
-.. literalinclude:: /code/include/graphene.py
- :start-after: #HIDDEN_BEGIN_shzy
- :end-before: #HIDDEN_END_shzy
+.. jupyter-execute::
+ :hide-code:
+
+ r = 10
+ w = 2.0
+ pot = 0.1
+
+.. jupyter-execute::
+
+ #### Define the scattering region. ####
+ # circular scattering region
+ def circle(pos):
+ x, y = pos
+ return x ** 2 + y ** 2 < r ** 2
+
+ syst = kwant.Builder()
+
+ # w: width and pot: potential maximum of the p-n junction
+ def potential(site):
+ (x, y) = site.pos
+ d = y * cos_30 + x * sin_30
+ return pot * tanh(d / w)
+
+ syst[graphene.shape(circle, (0, 0))] = potential
As you can see, this works exactly the same for any kind of lattice.
We add the onsite energies using a function describing the p-n junction;
@@ -45,9 +100,11 @@ As a next step we add the hoppings, making use of
`~kwant.builder.HoppingKind`. For illustration purposes we define
the hoppings ourselves instead of using ``graphene.neighbors()``:
-.. literalinclude:: /code/include/graphene.py
- :start-after: #HIDDEN_BEGIN_hsmc
- :end-before: #HIDDEN_END_hsmc
+.. jupyter-execute::
+
+ # specify the hoppings of the graphene lattice in the
+ # format expected by builder.HoppingKind
+ hoppings = (((0, 0), a, b), ((0, 1), a, b), ((-1, 1), a, b))
The nearest-neighbor model for graphene contains only
hoppings between different basis atoms. For this type of
@@ -63,27 +120,50 @@ respect to the two primitive vectors ``[(1, 0), (sin_30, cos_30)]``.
Adding the hoppings however still works the same way:
-.. literalinclude:: /code/include/graphene.py
- :start-after: #HIDDEN_BEGIN_bfwb
- :end-before: #HIDDEN_END_bfwb
+.. jupyter-execute::
+
+ syst[[kwant.builder.HoppingKind(*hopping) for hopping in hoppings]] = -1
Modifying the scattering region is also possible as before. Let's
do something crazy, and remove an atom in sublattice A
(which removes also the hoppings from/to this site) as well
as add an additional link:
-.. literalinclude:: /code/include/graphene.py
- :start-after: #HIDDEN_BEGIN_efut
- :end-before: #HIDDEN_END_efut
+.. jupyter-execute::
+
+ # Modify the scattering region
+ del syst[a(0, 0)]
+ syst[a(-2, 1), b(2, 2)] = -1
Note again that the conversion from a tuple `(i,j)` to site
is done by the sublattices `a` and `b`.
The leads are defined almost as before:
-.. literalinclude:: /code/include/graphene.py
- :start-after: #HIDDEN_BEGIN_aakh
- :end-before: #HIDDEN_END_aakh
+.. jupyter-execute::
+
+ #### Define the leads. ####
+ # left lead
+ sym0 = kwant.TranslationalSymmetry(graphene.vec((-1, 0)))
+
+ def lead0_shape(pos):
+ x, y = pos
+ return (-0.4 * r < y < 0.4 * r)
+
+ lead0 = kwant.Builder(sym0)
+ lead0[graphene.shape(lead0_shape, (0, 0))] = -pot
+ lead0[[kwant.builder.HoppingKind(*hopping) for hopping in hoppings]] = -1
+
+ # The second lead, going to the top right
+ sym1 = kwant.TranslationalSymmetry(graphene.vec((0, 1)))
+
+ def lead1_shape(pos):
+ v = pos[1] * sin_30 - pos[0] * cos_30
+ return (-0.4 * r < v < 0.4 * r)
+
+ lead1 = kwant.Builder(sym1)
+ lead1[graphene.shape(lead1_shape, (0, 0))] = pot
+ lead1[[kwant.builder.HoppingKind(*hopping) for hopping in hoppings]] = -1
Note the method `~kwant.lattice.Polyatomic.vec` used in calculating the
parameter for `~kwant.lattice.TranslationalSymmetry`. The latter expects a
@@ -98,28 +178,63 @@ in a square lattice -- they follow the non-orthogonal primitive vectors defined
in the beginning.
Later, we will compute some eigenvalues of the closed scattering region without
-leads. This is why we postpone attaching the leads to the system. Instead,
-we return the scattering region and the leads separately.
+leads. This is why we postpone attaching the leads to the system.
-.. literalinclude:: /code/include/graphene.py
- :start-after: #HIDDEN_BEGIN_kmmw
- :end-before: #HIDDEN_END_kmmw
The computation of some eigenvalues of the closed system is done
in the following piece of code:
-.. literalinclude:: /code/include/graphene.py
- :start-after: #HIDDEN_BEGIN_zydk
- :end-before: #HIDDEN_END_zydk
+.. jupyter-execute::
+
+ def compute_evs(syst):
+ # Compute some eigenvalues of the closed system
+ sparse_mat = syst.hamiltonian_submatrix(sparse=True)
+
+ evs = sla.eigs(sparse_mat, 2)[0]
+ print(evs.real)
The code for computing the band structure and the conductance is identical
to the previous examples, and needs not be further explained here.
-Finally, in the `main` function we make and plot the system:
+Finally, we plot the system:
+
+.. jupyter-execute::
+ :hide-code:
+
+ def plot_conductance(syst, energies):
+ # Compute transmission as a function of energy
+ data = []
+ for energy in energies:
+ smatrix = kwant.smatrix(syst, energy)
+ data.append(smatrix.transmission(0, 1))
+
+ pyplot.figure()
+ pyplot.plot(energies, data)
+ pyplot.xlabel("energy [t]")
+ pyplot.ylabel("conductance [e^2/h]")
+ pyplot.show()
+
+
+ def plot_bandstructure(flead, momenta):
+ bands = kwant.physics.Bands(flead)
+ energies = [bands(k) for k in momenta]
+
+ pyplot.figure()
+ pyplot.plot(momenta, energies)
+ pyplot.xlabel("momentum [(lattice constant)^-1]")
+ pyplot.ylabel("energy [t]")
+ pyplot.show()
+
+
+.. jupyter-execute::
+
+ # To highlight the two sublattices of graphene, we plot one with
+ # a filled, and the other one with an open circle:
+ def family_colors(site):
+ return 0 if site.family == a else 1
-.. literalinclude:: /code/include/graphene.py
- :start-after: #HIDDEN_BEGIN_itkk
- :end-before: #HIDDEN_END_itkk
+ # Plot the closed system without leads.
+ kwant.plot(syst, site_color=family_colors, site_lw=0.1, colorbar=False);
We customize the plotting: we set the `site_colors` argument of
`~kwant.plotter.plot` to a function which returns 0 for
@@ -132,27 +247,43 @@ The function `~kwant.plotter.plot` shows these values using a color scale
(grayscale by default). The symbol `size` is specified in points, and is
independent on the overall figure size.
-Plotting the closed system gives this result:
-
-.. image:: /code/figure/graphene_syst1.*
Computing the eigenvalues of largest magnitude,
-.. literalinclude:: /code/include/graphene.py
- :start-after: #HIDDEN_BEGIN_jmbi
- :end-before: #HIDDEN_END_jmbi
+.. jupyter-execute::
-should yield two eigenvalues equal to ``[ 3.07869311,
+ compute_evs(syst.finalized())
+
+yields two eigenvalues equal to ``[ 3.07869311,
-3.06233144]``.
-The remaining code of `main` attaches the leads to the system and plots it
+The remaining code attaches the leads to the system and plots it
again:
-.. image:: /code/figure/graphene_syst2.*
+.. jupyter-execute::
+
+ # Attach the leads to the system.
+ for lead in [lead0, lead1]:
+ syst.attach_lead(lead)
+
+ # Then, plot the system with leads.
+ kwant.plot(syst, site_color=family_colors, site_lw=0.1,
+ lead_site_lw=0, colorbar=False);
-It computes the band structure of one of lead 0:
+We then finalize the system:
-.. image:: /code/figure/graphene_bs.*
+.. jupyter-execute::
+
+ syst = syst.finalized()
+
+and compute the band structure of one of lead 0:
+
+.. jupyter-execute::
+
+
+ # Compute the band structure of lead 0.
+ momenta = [-pi + 0.02 * pi * i for i in range(101)]
+ plot_bandstructure(syst.leads[0], momenta)
showing all the features of a zigzag lead, including the flat
edge state bands (note that the band structure is not symmetric around
@@ -160,7 +291,11 @@ zero energy, due to a potential in the leads).
Finally the transmission through the system is computed,
-.. image:: /code/figure/graphene_result.*
+.. jupyter-execute::
+
+ # Plot conductance.
+ energies = [-2 * pot + 4. / 50. * pot * i for i in range(51)]
+ plot_conductance(syst, energies)
showing the typical resonance-like transmission probability through
an open quantum dot
diff --git a/doc/source/tutorial/kpm.rst b/doc/source/tutorial/kpm.rst
index 38c84891cbd84ea044e8600ee56c3e3d0a1cd638..7e1232581335aa1a77969899f59344e34c93eba7 100644
--- a/doc/source/tutorial/kpm.rst
+++ b/doc/source/tutorial/kpm.rst
@@ -17,8 +17,30 @@ KPM method `kwant.kpm`, that is based on the algorithms presented in Ref. [1]_.
.. seealso::
The complete source code of this example can be found in
- :download:`kernel_polynomial_method.py </code/download/kernel_polynomial_method.py>`
+ :jupyter-download:script:`kernel_polynomial_method`
+.. jupyter-kernel::
+ :id: kernel_polynomial_method
+
+.. jupyter-execute::
+ :hide-code:
+
+ # Tutorial 2.8. Calculating spectral density with the Kernel Polynomial Method
+ # ============================================================================
+ #
+ # Physics background
+ # ------------------
+ # - Chebyshev polynomials, random trace approximation, spectral densities.
+ #
+ # Kwant features highlighted
+ # --------------------------
+ # - kpm module,kwant operators.
+
+ import scipy
+ from matplotlib import pyplot
+
+.. jupyter-execute:: boilerplate.py
+ :hide-code:
Introduction
************
@@ -103,35 +125,104 @@ to obtain the (global) density of states of a graphene disk.
We start by importing kwant and defining our system.
-.. literalinclude:: /code/include/kernel_polynomial_method.py
- :start-after: #HIDDEN_BEGIN_sys1
- :end-before: #HIDDEN_END_sys1
+.. jupyter-execute::
+
+ # necessary imports
+ import kwant
+ import numpy as np
+
+
+ # define the system
+ def make_syst(r=30, t=-1, a=1):
+ syst = kwant.Builder()
+ lat = kwant.lattice.honeycomb(a, norbs=1)
+
+ def circle(pos):
+ x, y = pos
+ return x ** 2 + y ** 2 < r ** 2
+
+ syst[lat.shape(circle, (0, 0))] = 0.
+ syst[lat.neighbors()] = t
+ syst.eradicate_dangling()
+
+ return syst
+
+.. jupyter-execute::
+ :hide-code:
+
+ ## common plotting routines ##
+
+ # Plot several density of states curves on the same axes.
+ def plot_dos(labels_to_data):
+ pyplot.figure()
+ for label, (x, y) in labels_to_data:
+ pyplot.plot(x, y.real, label=label, linewidth=2)
+ pyplot.legend(loc=2, framealpha=0.5)
+ pyplot.xlabel("energy [t]")
+ pyplot.ylabel("DoS [a.u.]")
+ pyplot.show()
+
+
+ # Plot fill density of states plus curves on the same axes.
+ def plot_dos_and_curves(dos, labels_to_data):
+ pyplot.figure()
+ pyplot.fill_between(dos[0], dos[1], label="DoS [a.u.]",
+ alpha=0.5, color='gray')
+ for label, (x, y) in labels_to_data:
+ pyplot.plot(x, y, label=label, linewidth=2)
+ pyplot.legend(loc=2, framealpha=0.5)
+ pyplot.xlabel("energy [t]")
+ pyplot.ylabel("$σ [e^2/h]$")
+ pyplot.show()
+
+
+ def site_size_conversion(densities):
+ return 3 * np.abs(densities) / max(densities)
+
+
+ # Plot several local density of states maps in different subplots
+ def plot_ldos(syst, densities):
+ fig, axes = pyplot.subplots(1, len(densities), figsize=(7*len(densities), 7))
+ for ax, (title, rho) in zip(axes, densities):
+ kwant.plotter.density(syst, rho.real, ax=ax)
+ ax.set_title(title)
+ ax.set(adjustable='box', aspect='equal')
+ pyplot.show()
After making a system we can then create a `~kwant.kpm.SpectralDensity`
object that represents the density of states for this system.
-.. literalinclude:: /code/include/kernel_polynomial_method.py
- :start-after: #HIDDEN_BEGIN_kpm1
- :end-before: #HIDDEN_END_kpm1
+.. jupyter-execute::
+
+ fsyst = make_syst().finalized()
+
+ spectrum = kwant.kpm.SpectralDensity(fsyst)
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.
-.. literalinclude:: /code/include/kernel_polynomial_method.py
- :start-after: #HIDDEN_BEGIN_kpm2
- :end-before: #HIDDEN_END_kpm2
+.. jupyter-execute::
+
+ energies, densities = spectrum()
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.
-.. literalinclude:: /code/include/kernel_polynomial_method.py
- :start-after: #HIDDEN_BEGIN_kpm3
- :end-before: #HIDDEN_END_kpm3
+.. jupyter-execute::
+
+ energy_subset = np.linspace(0, 2)
+ density_subset = spectrum(energy_subset)
+
+.. jupyter-execute::
+ :hide-code:
-.. image:: /code/figure/kpm_dos.*
+ plot_dos([
+ ('densities', (energies, densities)),
+ ('density subset', (energy_subset, density_subset)),
+ ])
In addition to being called like functions, `~kwant.kpm.SpectralDensity`
objects also have a method `~kwant.kpm.SpectralDensity.integrate` which can be
@@ -139,22 +230,21 @@ 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:: /code/include/kernel_polynomial_method.py
- :start-after: #HIDDEN_BEGIN_int1
- :end-before: #HIDDEN_END_int1
+.. jupyter-execute::
-.. literalinclude:: /code/figure/kpm_normalization.txt
+ print('identity resolution:', spectrum.integrate())
We see that the integral of the density of states is normalized to the total
number of available states in the system. If we wish to calculate, say, the
number of states populated in equilibrium, then we should integrate with
respect to a Fermi-Dirac distribution:
-.. literalinclude:: /code/include/kernel_polynomial_method.py
- :start-after: #HIDDEN_BEGIN_int2
- :end-before: #HIDDEN_END_int2
+.. jupyter-execute::
-.. literalinclude:: /code/figure/kpm_total_states.txt
+ # Fermi energy 0.1 and temperature 0.2
+ fermi = lambda E: 1 / (np.exp((E - 0.1) / 0.2) + 1)
+
+ print('number of filled states:', spectrum.integrate(fermi))
.. specialnote:: Stability and performance: spectral bounds
@@ -194,23 +284,51 @@ In the following example, we compute the local density of states at the center
of the graphene disk, and we add a staggering potential between the two
sublattices.
-.. literalinclude:: /code/include/kernel_polynomial_method.py
- :start-after: #HIDDEN_BEGIN_sys3
- :end-before: #HIDDEN_END_sys3
+.. jupyter-execute::
+
+ def make_syst_staggered(r=30, t=-1, a=1, m=0.1):
+ syst = kwant.Builder()
+ lat = kwant.lattice.honeycomb(a, norbs=1)
+
+ def circle(pos):
+ x, y = pos
+ return x ** 2 + y ** 2 < r ** 2
+
+ syst[lat.a.shape(circle, (0, 0))] = m
+ syst[lat.b.shape(circle, (0, 0))] = -m
+ syst[lat.neighbors()] = t
+ syst.eradicate_dangling()
+
+ return syst
Next, we choose one site of each sublattice "A" and "B",
-.. literalinclude:: /code/include/kernel_polynomial_method.py
- :start-after: #HIDDEN_BEGIN_op5
- :end-before: #HIDDEN_END_op5
+.. jupyter-execute::
+
+ fsyst_staggered = make_syst_staggered().finalized()
+ # find 'A' and 'B' sites in the unit cell at the center of the disk
+ center_tag = np.array([0, 0])
+ where = lambda s: s.tag == center_tag
+ # make local vectors
+ vector_factory = kwant.kpm.LocalVectors(fsyst_staggered, where)
and plot their respective local density of states.
-.. literalinclude:: /code/include/kernel_polynomial_method.py
- :start-after: #HIDDEN_BEGIN_op6
- :end-before: #HIDDEN_END_op6
+.. jupyter-execute::
+
+ # 'num_vectors' can be unspecified when using 'LocalVectors'
+ local_dos = kwant.kpm.SpectralDensity(fsyst_staggered, num_vectors=None,
+ vector_factory=vector_factory,
+ mean=False)
+ energies, densities = local_dos()
-.. image:: /code/figure/kpm_ldos_sites.*
+.. jupyter-execute::
+ :hide-code:
+
+ plot_dos([
+ ('A sublattice', (energies, densities[:, 0])),
+ ('B sublattice', (energies, densities[:, 1])),
+ ])
Note that there is no noise comming from the random vectors.
@@ -225,9 +343,15 @@ exactly is changed.
The simplest way to obtain a more accurate solution is to use the
``add_moments`` method:
-.. literalinclude:: /code/include/kernel_polynomial_method.py
- :start-after: #HIDDEN_BEGIN_acc1
- :end-before: #HIDDEN_END_acc1
+.. jupyter-execute::
+ :hide-code:
+
+ spectrum = kwant.kpm.SpectralDensity(fsyst)
+ original_dos = spectrum()
+
+.. jupyter-execute::
+
+ spectrum.add_moments(energy_resolution=0.03)
This will update the number of calculated moments and also the default
number of sampling points such that the maximum distance between successive
@@ -237,12 +361,19 @@ Alternatively, you can directly increase the number of moments
with ``add_moments``, or the number of random vectors with ``add_vectors``.
The added vectors will be generated with the ``vector_factory``.
-.. literalinclude:: /code/include/kernel_polynomial_method.py
- :start-after: #HIDDEN_BEGIN_acc2
- :end-before: #HIDDEN_END_acc2
+.. jupyter-execute::
+
+ spectrum.add_moments(100)
+ spectrum.add_vectors(5)
-.. image:: /code/figure/kpm_dos_r.*
+.. jupyter-execute::
+ :hide-code:
+ increased_moments_dos = spectrum()
+ plot_dos([
+ ('density', original_dos),
+ ('higher number of moments', increased_moments_dos),
+ ])
.. _operator_spectral_density:
@@ -264,17 +395,19 @@ the spectral density of the given operator that is calculated.
If an explicit matrix is provided then it must have the same
shape as the system Hamiltonian.
-.. literalinclude:: /code/include/kernel_polynomial_method.py
- :start-after: #HIDDEN_BEGIN_op1
- :end-before: #HIDDEN_END_op1
+.. jupyter-execute::
+ # identity matrix
+ matrix_op = scipy.sparse.eye(len(fsyst.sites))
+ matrix_spectrum = kwant.kpm.SpectralDensity(fsyst, operator=matrix_op)
Or, to do the same calculation using `kwant.operator.Density`:
-.. literalinclude:: /code/include/kernel_polynomial_method.py
- :start-after: #HIDDEN_BEGIN_op2
- :end-before: #HIDDEN_END_op2
+.. jupyter-execute::
+ # 'sum=True' means we sum over all the sites
+ kwant_op = kwant.operator.Density(fsyst, sum=True)
+ operator_spectrum = kwant.kpm.SpectralDensity(fsyst, operator=kwant_op)
Spectral density with random vectors
^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
@@ -283,9 +416,11 @@ 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:: /code/include/kernel_polynomial_method.py
- :start-after: #HIDDEN_BEGIN_op3
- :end-before: #HIDDEN_END_op3
+.. jupyter-execute::
+
+ # 'sum=False' is the default, but we include it explicitly here for clarity.
+ kwant_op = kwant.operator.Density(fsyst, sum=False)
+ local_dos = kwant.kpm.SpectralDensity(fsyst, operator=kwant_op)
`~kwant.kpm.SpectralDensity` will properly handle this vector output,
and will average the local density obtained with random vectors.
@@ -294,12 +429,18 @@ The accuracy of this approximation depends on the number of random vectors used.
This allows us to plot an approximate local density of states at different
points in the spectrum:
-.. literalinclude:: /code/include/kernel_polynomial_method.py
- :start-after: #HIDDEN_BEGIN_op4
- :end-before: #HIDDEN_END_op4
+.. jupyter-execute::
-.. image:: /code/figure/kpm_ldos.*
+ zero_energy_ldos = local_dos(energy=0)
+ finite_energy_ldos = local_dos(energy=1)
+.. jupyter-execute::
+ :hide-code:
+
+ plot_ldos(fsyst, [
+ ('energy = 0', zero_energy_ldos),
+ ('energy = 1', finite_energy_ldos)
+ ])
Calculating Kubo conductivity
*****************************
@@ -325,17 +466,53 @@ with nearest neigbours hoppings. To turn it into a topological
insulator we add imaginary second nearest neigbours hoppings, where
the sign changes for each sublattice.
-.. literalinclude:: /code/include/kernel_polynomial_method.py
- :start-after: #HIDDEN_BEGIN_sys2
- :end-before: #HIDDEN_END_sys2
+.. jupyter-execute::
+
+ def make_syst_topo(r=30, a=1, t=1, t2=0.5):
+ syst = kwant.Builder()
+ lat = kwant.lattice.honeycomb(a, norbs=1, name=['a', 'b'])
+
+ def circle(pos):
+ x, y = pos
+ return x ** 2 + y ** 2 < r ** 2
+
+ syst[lat.shape(circle, (0, 0))] = 0.
+ syst[lat.neighbors()] = t
+ # add second neighbours hoppings
+ syst[lat.a.neighbors()] = 1j * t2
+ syst[lat.b.neighbors()] = -1j * t2
+ syst.eradicate_dangling()
+
+ return lat, syst.finalized()
To calculate the bulk conductivity, we will select sites in the unit cell
in the middle of the sample, and create a vector factory that outputs local
vectors
-.. literalinclude:: /code/include/kernel_polynomial_method.py
- :start-after: #HIDDEN_BEGIN_cond
- :end-before: #HIDDEN_END_cond
+.. jupyter-execute::
+
+ # construct the Haldane model
+ lat, fsyst_topo = make_syst_topo()
+ # find 'A' and 'B' sites in the unit cell at the center of the disk
+ where = lambda s: np.linalg.norm(s.pos) < 1
+
+ # component 'xx'
+ s_factory = kwant.kpm.LocalVectors(fsyst_topo, where)
+ cond_xx = kwant.kpm.conductivity(fsyst_topo, alpha='x', beta='x', mean=True,
+ num_vectors=None, vector_factory=s_factory)
+ # component 'xy'
+ s_factory = kwant.kpm.LocalVectors(fsyst_topo, where)
+ cond_xy = kwant.kpm.conductivity(fsyst_topo, alpha='x', beta='y', mean=True,
+ num_vectors=None, vector_factory=s_factory)
+
+ energies = cond_xx.energies
+ cond_array_xx = np.array([cond_xx(e, temperature=0.01) for e in energies])
+ cond_array_xy = np.array([cond_xy(e, temperature=0.01) for e in energies])
+
+ # area of the unit cell per site
+ area_per_site = np.abs(np.cross(*lat.prim_vecs)) / len(lat.sublattices)
+ cond_array_xx /= area_per_site
+ cond_array_xy /= area_per_site
Note that the Kubo conductivity must be normalized with the area covered
by the vectors. In this case, each local vector represents a site, and
@@ -345,8 +522,22 @@ value of the conductivity over large parts of the system. In this
case, the area that normalizes the result, is the area covered by
the random vectors.
-.. image:: /code/figure/kpm_cond.*
+.. jupyter-execute::
+ :hide-code:
+
+ s_factory = kwant.kpm.LocalVectors(fsyst_topo, where)
+ spectrum = kwant.kpm.SpectralDensity(fsyst_topo, num_vectors=None,
+ vector_factory=s_factory)
+
+ plot_dos_and_curves(
+ (spectrum.energies, spectrum.densities * 8),
+ [
+ (r'Longitudinal conductivity $\sigma_{xx} / 4$',
+ (energies, cond_array_xx / 4)),
+ (r'Hall conductivity $\sigma_{xy}$',
+ (energies, cond_array_xy))],
+ )
.. _advanced_topics:
@@ -370,9 +561,16 @@ 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:: /code/include/kernel_polynomial_method.py
- :start-after: #HIDDEN_BEGIN_fact1
- :end-before: #HIDDEN_END_fact1
+.. jupyter-execute::
+
+ # construct a generator of vectors with n random elements -1 or +1.
+ n = fsyst.hamiltonian_submatrix(sparse=True).shape[0]
+ def binary_vectors():
+ while True:
+ yield np.rint(np.random.random_sample(n)) * 2 - 1
+
+ custom_factory = kwant.kpm.SpectralDensity(fsyst,
+ vector_factory=binary_vectors())
Aditionally, a `~kwant.kpm.LocalVectors` generator is also available, that
returns local vectors that correspond to the sites passed. Note that
@@ -407,9 +605,16 @@ 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:: /code/include/kernel_polynomial_method.py
- :start-after: #HIDDEN_BEGIN_blm
- :end-before: #HIDDEN_END_blm
+.. jupyter-execute::
+
+ rho = kwant.operator.Density(fsyst, sum=True)
+
+ # sesquilinear map that does the same thing as `rho`
+ def rho_alt(bra, ket):
+ return np.vdot(bra, ket)
+
+ rho_spectrum = kwant.kpm.SpectralDensity(fsyst, operator=rho)
+ rho_alt_spectrum = kwant.kpm.SpectralDensity(fsyst, operator=rho_alt)
__ operator_spectral_density_
diff --git a/doc/source/tutorial/operators.rst b/doc/source/tutorial/operators.rst
index ae397e3f55ce58656b7910724c5d78dadd8542d9..1646ae86e2cd5a47747843f6dd952144b0a92175 100644
--- a/doc/source/tutorial/operators.rst
+++ b/doc/source/tutorial/operators.rst
@@ -13,8 +13,46 @@ texture.
.. seealso::
The complete source code of this example can be found in
- :download:`magnetic_texture.py </code/download/magnetic_texture.py>`
-
+ :jupyter-download:script:`magnetic_texture`
+
+.. jupyter-kernel::
+ :id: magnetic_texture
+
+.. jupyter-execute::
+ :hide-code:
+
+ # Tutorial 2.7. Spin textures
+ # ===========================
+ #
+ # Physics background
+ # ------------------
+ # - Spin textures
+ # - Skyrmions
+ #
+ # Kwant features highlighted
+ # --------------------------
+ # - operators
+ # - plotting vector fields
+
+ from math import sin, cos, tanh, pi
+ import itertools
+ import numpy as np
+ import tinyarray as ta
+ import matplotlib.pyplot as plt
+
+ import kwant
+
+ sigma_0 = ta.array([[1, 0], [0, 1]])
+ sigma_x = ta.array([[0, 1], [1, 0]])
+ sigma_y = ta.array([[0, -1j], [1j, 0]])
+ sigma_z = ta.array([[1, 0], [0, -1]])
+
+ # vector of Pauli matrices σ_αiβ where greek
+ # letters denote spinor indices
+ sigma = np.rollaxis(np.array([sigma_x, sigma_y, sigma_z]), 1)
+
+.. jupyter-execute:: boilerplate.py
+ :hide-code:
Introduction
------------
@@ -47,28 +85,119 @@ of site :math:`i`, and :math:`r_i = \sqrt{x_i^2 + y_i^2}`.
To define this model in Kwant we start as usual by defining functions that
depend on the model parameters:
-.. literalinclude:: /code/include/magnetic_texture.py
- :start-after: #HIDDEN_BEGIN_model
- :end-before: #HIDDEN_END_model
+.. jupyter-execute::
+
+ def field_direction(pos, r0, delta):
+ x, y = pos
+ r = np.linalg.norm(pos)
+ r_tilde = (r - r0) / delta
+ theta = (tanh(r_tilde) - 1) * (pi / 2)
+
+ if r == 0:
+ m_i = [0, 0, -1]
+ else:
+ m_i = [
+ (x / r) * sin(theta),
+ (y / r) * sin(theta),
+ cos(theta),
+ ]
+
+ return np.array(m_i)
+
+
+ def scattering_onsite(site, r0, delta, J):
+ m_i = field_direction(site.pos, r0, delta)
+ return J * np.dot(m_i, sigma)
+
+
+ def lead_onsite(site, J):
+ return J * sigma_z
and define our system as a square shape on a square lattice with two orbitals
per site, with leads attached on the left and right:
-.. literalinclude:: /code/include/magnetic_texture.py
- :start-after: #HIDDEN_BEGIN_syst
- :end-before: #HIDDEN_END_syst
+.. jupyter-execute::
+
+ lat = kwant.lattice.square(norbs=2)
+
+ def make_system(L=80):
+
+ syst = kwant.Builder()
+
+ def square(pos):
+ return all(-L/2 < p < L/2 for p in pos)
+
+ syst[lat.shape(square, (0, 0))] = scattering_onsite
+ syst[lat.neighbors()] = -sigma_0
+
+ lead = kwant.Builder(kwant.TranslationalSymmetry((-1, 0)),
+ conservation_law=-sigma_z)
+
+ lead[lat.shape(square, (0, 0))] = lead_onsite
+ lead[lat.neighbors()] = -sigma_0
+
+ syst.attach_lead(lead)
+ syst.attach_lead(lead.reversed())
+
+ return syst
Below is a plot of a projection of :math:`\mathbf{m}_i` onto the x-y plane
inside the scattering region. The z component is shown by the color scale:
-.. image:: /code/figure/mag_field_direction.*
+.. jupyter-execute::
+ :hide-code:
+
+ def plot_vector_field(syst, params):
+ xmin, ymin = min(s.tag for s in syst.sites)
+ xmax, ymax = max(s.tag for s in syst.sites)
+ x, y = np.meshgrid(np.arange(xmin, xmax+1), np.arange(ymin, ymax+1))
+
+ m_i = [field_direction(p, **params) for p in zip(x.flat, y.flat)]
+ m_i = np.reshape(m_i, x.shape + (3,))
+ m_i = np.rollaxis(m_i, 2, 0)
+
+ fig, ax = plt.subplots(1, 1)
+ im = ax.quiver(x, y, *m_i, pivot='mid', scale=75)
+ fig.colorbar(im)
+ plt.show()
+
+
+ def plot_densities(syst, densities):
+ fig, axes = plt.subplots(1, len(densities), figsize=(13, 10))
+ for ax, (title, rho) in zip(axes, densities):
+ kwant.plotter.map(syst, rho, ax=ax, a=4)
+ ax.set_title(title)
+ plt.show()
+
+
+ def plot_currents(syst, currents):
+ fig, axes = plt.subplots(1, len(currents), figsize=(13, 10))
+ if not hasattr(axes, '__len__'):
+ axes = (axes,)
+ for ax, (title, current) in zip(axes, currents):
+ kwant.plotter.current(syst, current, ax=ax, colorbar=False,
+ fig_size=(13, 10))
+ ax.set_title(title)
+ plt.show()
+
+.. jupyter-execute::
+ :hide-code:
+
+ syst = make_system().finalized()
+
+.. jupyter-execute::
+ :hide-code:
+
+ plot_vector_field(syst, dict(r0=20, delta=10))
We will now be interested in analyzing the form of the scattering states
that originate from the left lead:
-.. literalinclude:: /code/include/magnetic_texture.py
- :start-after: #HIDDEN_BEGIN_wavefunction
- :end-before: #HIDDEN_END_wavefunction
+.. jupyter-execute::
+
+ params = dict(r0=20, delta=10, J=1)
+ wf = kwant.wave_function(syst, energy=-1, params=params)
+ psi = wf(0)[0]
Local densities
---------------
@@ -82,34 +211,45 @@ When there are multiple degrees of freedom per site, however, one has to be
more careful. In the present case with two (spin) degrees of freedom per site
one could calculate the per-site density like:
-.. literalinclude:: /code/include/magnetic_texture.py
- :start-after: #HIDDEN_BEGIN_ldos
- :end-before: #HIDDEN_END_ldos
+.. jupyter-execute::
+
+ # even (odd) indices correspond to spin up (down)
+ up, down = psi[::2], psi[1::2]
+ density = np.abs(up)**2 + np.abs(down)**2
With more than one degree of freedom per site we have more freedom as to what
local quantities we can meaningfully compute. For example, we may wish to
calculate the local z-projected spin density. We could calculate
this in the following way:
-.. literalinclude:: /code/include/magnetic_texture.py
- :start-after: #HIDDEN_BEGIN_lsdz
- :end-before: #HIDDEN_END_lsdz
+.. jupyter-execute::
+
+ # spin down components have a minus sign
+ spin_z = np.abs(up)**2 - np.abs(down)**2
If we wanted instead to calculate the local y-projected spin density, we would
need to use an even more complicated expression:
-.. literalinclude:: /code/include/magnetic_texture.py
- :start-after: #HIDDEN_BEGIN_lsdy
- :end-before: #HIDDEN_END_lsdy
+.. jupyter-execute::
+
+ # spin down components have a minus sign
+ spin_y = 1j * (down.conjugate() * up - up.conjugate() * down)
The `kwant.operator` module aims to alleviate somewhat this tedious
book-keeping by providing a simple interface for defining operators that act on
wavefunctions. To calculate the above quantities we would use the
`~kwant.operator.Density` operator like so:
-.. literalinclude:: /code/include/magnetic_texture.py
- :start-after: #HIDDEN_BEGIN_lden
- :end-before: #HIDDEN_END_lden
+.. jupyter-execute::
+
+ rho = kwant.operator.Density(syst)
+ rho_sz = kwant.operator.Density(syst, sigma_z)
+ rho_sy = kwant.operator.Density(syst, sigma_y)
+
+ # calculate the expectation values of the operators with 'psi'
+ density = rho(psi)
+ spin_z = rho_sz(psi)
+ spin_y = rho_sy(psi)
`~kwant.operator.Density` takes a `~kwant.system.System` as its first parameter
as well as (optionally) a square matrix that defines the quantity that you wish
@@ -126,8 +266,14 @@ Below we can see colorplots of the above-calculated quantities. The array that
is returned by evaluating a `~kwant.operator.Density` can be used directly with
`kwant.plotter.density`:
-.. image:: /code/figure/spin_densities.*
+.. jupyter-execute::
+ :hide-code:
+ plot_densities(syst, [
+ ('$σ_0$', density),
+ ('$σ_z$', spin_z),
+ ('$σ_y$', spin_y),
+ ])
.. specialnote:: Technical Details
@@ -180,14 +326,27 @@ where :math:`\mathbf{H}_{ab}` is the hopping matrix from site :math:`b` to site
:math:`a`. For example, to calculate the local current and
spin current:
-.. literalinclude:: /code/include/magnetic_texture.py
- :start-after: #HIDDEN_BEGIN_current
- :end-before: #HIDDEN_END_current
+.. jupyter-execute::
+
+ J_0 = kwant.operator.Current(syst)
+ J_z = kwant.operator.Current(syst, sigma_z)
+ J_y = kwant.operator.Current(syst, sigma_y)
+
+ # calculate the expectation values of the operators with 'psi'
+ current = J_0(psi)
+ spin_z_current = J_z(psi)
+ spin_y_current = J_y(psi)
Evaluating a `~kwant.operator.Current` operator on a wavefunction returns a
1D array of values that can be directly used with `kwant.plotter.current`:
-.. image:: /code/figure/spin_currents.*
+.. jupyter-execute::
+
+ plot_currents(syst, [
+ ('$J_{σ_0}$', current),
+ ('$J_{σ_z}$', spin_z_current),
+ ('$J_{σ_y}$', spin_y_current),
+ ])
.. note::
@@ -262,9 +421,16 @@ scattering region.
Doing this is as simple as passing a *function* when instantiating
the `~kwant.operator.Current`, instead of a constant matrix:
-.. literalinclude:: /code/include/magnetic_texture.py
- :start-after: #HIDDEN_BEGIN_following
- :end-before: #HIDDEN_END_following
+.. jupyter-execute::
+
+ def following_m_i(site, r0, delta):
+ m_i = field_direction(site.pos, r0, delta)
+ return np.dot(m_i, sigma)
+
+ J_m = kwant.operator.Current(syst, following_m_i)
+
+ # evaluate the operator
+ m_current = J_m(psi, params=dict(r0=25, delta=10))
The function must take a `~kwant.builder.Site` as its first parameter,
and may optionally take other parameters (i.e. it must have the same
@@ -274,7 +440,7 @@ matrix that defines the operator we wish to calculate.
.. note::
In the above example we had to pass the extra parameters needed by the
- ``following_operator`` function via the ``param`` keyword argument. In
+ ``following_operator`` function via the ``params`` keyword argument. In
general you must pass all the parameters needed by the Hamiltonian via
``params`` (as you would when calling `~kwant.solvers.default.smatrix` or
`~kwant.solvers.default.wave_function`). In the previous examples,
@@ -286,7 +452,13 @@ Using this we can see that the spin current is essentially oriented along
the direction of :math:`m_i` in the present regime where the onsite term
in the Hamiltonian is dominant:
-.. image:: /code/figure/spin_current_comparison.*
+.. jupyter-execute::
+ :hide-code:
+
+ plot_currents(syst, [
+ (r'$J_{\mathbf{m}_i}$', m_current),
+ ('$J_{σ_z}$', spin_z_current),
+ ])
.. note:: Although this example used exclusively `~kwant.operator.Current`,
you can do the same with `~kwant.operator.Density`.
@@ -308,11 +480,16 @@ Density of states in a circle
To calculate the density of states inside a circle of radius
20 we can simply do:
-.. literalinclude:: /code/include/magnetic_texture.py
- :start-after: #HIDDEN_BEGIN_density_cut
- :end-before: #HIDDEN_END_density_cut
+.. jupyter-execute::
+
+ def circle(site):
+ return np.linalg.norm(site.pos) < 20
-.. literalinclude:: /code/figure/circle_dos.txt
+ rho_circle = kwant.operator.Density(syst, where=circle, sum=True)
+
+ all_states = np.vstack((wf(0), wf(1)))
+ dos_in_circle = sum(rho_circle(p) for p in all_states) / (2 * pi)
+ print('density of states in circle:', dos_in_circle)
note that we also provide ``sum=True``, which means that evaluating the
operator on a wavefunction will produce a single scalar. This is semantically
@@ -325,11 +502,22 @@ Current flowing through a cut
Below we calculate the probability current and z-projected spin current near
the interfaces with the left and right leads.
-.. literalinclude:: /code/include/magnetic_texture.py
- :start-after: #HIDDEN_BEGIN_current_cut
- :end-before: #HIDDEN_END_current_cut
+.. jupyter-execute::
+
+ def left_cut(site_to, site_from):
+ return site_from.pos[0] <= -39 and site_to.pos[0] > -39
+
+ def right_cut(site_to, site_from):
+ return site_from.pos[0] < 39 and site_to.pos[0] >= 39
-.. literalinclude:: /code/figure/current_cut.txt
+ J_left = kwant.operator.Current(syst, where=left_cut, sum=True)
+ J_right = kwant.operator.Current(syst, where=right_cut, sum=True)
+
+ Jz_left = kwant.operator.Current(syst, sigma_z, where=left_cut, sum=True)
+ Jz_right = kwant.operator.Current(syst, sigma_z, where=right_cut, sum=True)
+
+ print('J_left:', J_left(psi), ' J_right:', J_right(psi))
+ print('Jz_left:', Jz_left(psi), ' Jz_right:', Jz_right(psi))
We see that the probability current is conserved across the scattering region,
but the z-projected spin current is not due to the fact that the Hamiltonian
@@ -360,8 +548,23 @@ This can be achieved with `~kwant.operator.Current.bind`:
*different* set of parameters. This will almost certainly give
incorrect results.
-.. literalinclude:: /code/include/magnetic_texture.py
- :start-after: #HIDDEN_BEGIN_bind
- :end-before: #HIDDEN_END_bind
+.. jupyter-execute::
+
+ J_m = kwant.operator.Current(syst, following_m_i)
+ J_z = kwant.operator.Current(syst, sigma_z)
+
+ J_m_bound = J_m.bind(params=dict(r0=25, delta=10, J=1))
+ J_z_bound = J_z.bind(params=dict(r0=25, delta=10, J=1))
+
+ # Sum current local from all scattering states on the left at energy=-1
+ wf_left = wf(0)
+ J_m_left = sum(J_m_bound(p) for p in wf_left)
+ J_z_left = sum(J_z_bound(p) for p in wf_left)
+
+.. jupyter-execute::
+ :hide-code:
-.. image:: /code/figure/bound_current.*
+ plot_currents(syst, [
+ (r'$J_{\mathbf{m}_i}$ (from left)', J_m_left),
+ (r'$J_{σ_z}$ (from left)', J_z_left),
+ ])
diff --git a/doc/source/tutorial/plotting.rst b/doc/source/tutorial/plotting.rst
index f8bb6586a1c2b7d6cd30cf084d983b5eef8b3ce5..ea04094a7381839540b5c5dfe82ea35258ff4767 100644
--- a/doc/source/tutorial/plotting.rst
+++ b/doc/source/tutorial/plotting.rst
@@ -12,15 +12,55 @@ these options can be used to achieve various very different objectives.
.. seealso::
The complete source code of this example can be found in
- :download:`plot_graphene.py </code/download/plot_graphene.py>`
+ :jupyter-download:script:`plot_graphene`
+
+.. jupyter-kernel::
+ :id: plot_graphene
+
+.. jupyter-execute::
+ :hide-code:
+
+ # Tutorial 2.8.1. 2D example: graphene quantum dot
+ # ================================================
+ #
+ # Physics background
+ # ------------------
+ # - graphene edge states
+ #
+ # Kwant features highlighted
+ # --------------------------
+ # - demonstrate different ways of plotting
+
+ from matplotlib import pyplot
+
+ import kwant
+
+.. jupyter-execute:: boilerplate.py
+ :hide-code:
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:: /code/include/plot_graphene.py
- :start-after: #HIDDEN_BEGIN_makesyst
- :end-before: #HIDDEN_END_makesyst
+.. jupyter-execute::
+
+ lat = kwant.lattice.honeycomb()
+ a, b = lat.sublattices
+
+ def make_system(r=8, t=-1, tp=-0.1):
+
+ def circle(pos):
+ x, y = pos
+ return x**2 + y**2 < r**2
+
+ syst = kwant.Builder()
+ syst[lat.shape(circle, (0, 0))] = 0
+ syst[lat.neighbors()] = t
+ syst.eradicate_dangling()
+ if tp:
+ syst[lat.neighbors(2)] = tp
+
+ return syst
Note that adding hoppings hoppings to the `n`-th nearest neighbors can be
simply done by passing `n` as an argument to
@@ -29,16 +69,14 @@ simply done by passing `n` as an argument to
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:: /code/include/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
+Of course, the system can be plotted simply with default settings,
+however, due to the richer structure of the lattice, this results in a rather
busy plot:
-.. image:: /code/figure/plot_graphene_syst1.*
+.. jupyter-execute::
+
+ syst = make_system()
+ kwant.plot(syst);
A much clearer plot can be obtained by using different colors for both
sublattices, and by having different line widths for different hoppings. This
@@ -47,15 +85,19 @@ can be achieved by passing a function to the arguments of
must be a function taking one site as argument, for hoppings a function taking
the start end end site of hopping as arguments:
-.. literalinclude:: /code/include/plot_graphene.py
- :start-after: #HIDDEN_BEGIN_plotsyst2
- :end-before: #HIDDEN_END_plotsyst2
+.. jupyter-execute::
+
+ def family_color(site):
+ return 'black' if site.family == a else 'white'
+
+ def hopping_lw(site1, site2):
+ return 0.04 if site1.family == site2.family else 0.1
+
+ kwant.plot(syst, site_lw=0.1, site_color=family_color, hop_lw=hopping_lw);
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:: /code/figure/plot_graphene_syst2.*
+that carries the same information, but is much easier to interpret.
Apart from plotting the *system* itself, `~kwant.plotter.plot` can also be used
to plot *data* living on the system.
@@ -67,23 +109,27 @@ 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:: /code/include/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:
+.. jupyter-execute::
+
+ import scipy.linalg as la
+
+ syst = make_system(tp=0).finalized()
+ ham = syst.hamiltonian_submatrix()
+ evecs = la.eigh(ham)[1]
-.. literalinclude:: /code/include/plot_graphene.py
- :start-after: #HIDDEN_BEGIN_plotdata2
- :end-before: #HIDDEN_END_plotdata2
+ wf = abs(evecs[:, 225])**2
+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.
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:: /code/figure/plot_graphene_data1.*
+.. jupyter-execute::
+
+ kwant.plotter.map(syst, wf, oversampling=10, cmap='gist_heat_r');
However although in general preferable, `~kwant.plotter.map` has a few
deficiencies for this small system: For example, there are a few distortions at
@@ -99,9 +145,17 @@ 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:: /code/include/plot_graphene.py
- :start-after: #HIDDEN_BEGIN_plotdata3
- :end-before: #HIDDEN_END_plotdata3
+.. jupyter-execute::
+
+ def family_shape(i):
+ site = syst.sites[i]
+ return ('p', 3, 180) if site.family == a else ('p', 3, 0)
+
+ def family_color(i):
+ return 'black' if syst.sites[i].family == a else 'white'
+
+ kwant.plot(syst, site_color=wf, site_symbol=family_shape,
+ site_size=0.5, hop_lw=0, cmap='gist_heat_r');
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
@@ -114,31 +168,23 @@ 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:: /code/figure/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:: /code/include/plot_graphene.py
- :start-after: #HIDDEN_BEGIN_plotdata4
- :end-before: #HIDDEN_END_plotdata4
+.. jupyter-execute::
+
+ def site_size(i):
+ return 3 * wf[i] / wf.max()
+
+ kwant.plot(syst, site_size=site_size, site_color=(0, 0, 1, 0.3),
+ hop_lw=0.1);
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:: /code/figure/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
@@ -150,7 +196,31 @@ which shows the edge state nature of the wave function most clearly.
.. seealso::
The complete source code of this example can be found in
- :download:`plot_zincblende.py </code/download/plot_zincblende.py>`
+ :jupyter-download:script:`plot_zincblende`
+
+.. jupyter-kernel::
+ :id: plot_zincblende
+
+.. jupyter-execute::
+ :hide-code:
+
+ # Tutorial 2.8.2. 3D example: zincblende structure
+ # ================================================
+ #
+ # Physical background
+ # -------------------
+ # - 3D Bravais lattices
+ #
+ # Kwant features highlighted
+ # --------------------------
+ # - demonstrate different ways of plotting in 3D
+
+ from matplotlib import pyplot
+
+ import kwant
+
+.. jupyter-execute:: boilerplate.py
+ :hide-code:
Zincblende is a very common crystal structure of semiconductors. It is a
face-centered cubic crystal with two inequivalent atoms in the unit cell
@@ -159,18 +229,29 @@ structure, but two equivalent atoms per unit cell).
It is very easily generated in Kwant with `kwant.lattice.general`:
-.. literalinclude:: /code/include/plot_zincblende.py
- :start-after: #HIDDEN_BEGIN_zincblende1
- :end-before: #HIDDEN_END_zincblende1
+.. jupyter-execute::
+
+ 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
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:: /code/include/plot_zincblende.py
- :start-after: #HIDDEN_BEGIN_zincblende2
- :end-before: #HIDDEN_END_zincblende2
+.. jupyter-execute::
+
+ 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
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
@@ -179,13 +260,11 @@ 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:: /code/include/plot_zincblende.py
- :start-after: #HIDDEN_BEGIN_plot1
- :end-before: #HIDDEN_END_plot1
+.. jupyter-execute::
-resulting in
+ syst = make_cuboid()
-.. image:: /code/figure/plot_zincblende_syst1.*
+ kwant.plot(syst);
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
@@ -207,14 +286,18 @@ 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:: /code/include/plot_zincblende.py
- :start-after: #HIDDEN_BEGIN_plot2
- :end-before: #HIDDEN_END_plot2
+.. jupyter-execute::
-which results in a 3D plot that allows to interactively (when plotted
-in a window) explore the crystal structure:
+ syst = make_cuboid(a=1.5, b=1.5, c=1.5)
-.. image:: /code/figure/plot_zincblende_syst2.*
+ 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);
+
+which results in a 3D plot that allows to interactively (when plotted
+in a window) explore the crystal structure.
Hence, a few lines of code using Kwant allow to explore all the different
crystal lattices out there!
diff --git a/doc/source/tutorial/spectrum.rst b/doc/source/tutorial/spectrum.rst
index 734f5a1d0e161d2434789bc9ce3f3a880e122419..b2cbc51d3bef2b378891baca292857360ca6d95d 100644
--- a/doc/source/tutorial/spectrum.rst
+++ b/doc/source/tutorial/spectrum.rst
@@ -6,7 +6,32 @@ Band structure calculations
.. seealso::
The complete source code of this example can be found in
- :download:`band_structure.py </code/download/band_structure.py>`
+ :jupyter-download:script:`band_structure`
+
+.. jupyter-kernel::
+ :id: band_structure
+
+.. jupyter-execute::
+ :hide-code:
+
+ # 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
+
+.. jupyter-execute:: boilerplate.py
+ :hide-code:
When doing transport simulations, one also often needs to know the band
structure of the leads, i.e. the energies of the propagating plane waves in the
@@ -19,9 +44,26 @@ tight-binding wire.
Computing band structures in Kwant is easy. Just define a lead in the
usual way:
-.. literalinclude:: /code/include/band_structure.py
- :start-after: #HIDDEN_BEGIN_zxip
- :end-before: #HIDDEN_END_zxip
+.. jupyter-execute::
+
+ 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
"Usual way" means defining a translational symmetry vector, as well
as one unit cell of the lead, and the hoppings to neighboring
@@ -42,13 +84,24 @@ do something more profound with the dispersion relation these energies may be
calculated directly using `~kwant.physics.Bands`. For now we just plot the
bandstructure:
-.. literalinclude:: /code/include/band_structure.py
- :start-after: #HIDDEN_BEGIN_pejz
- :end-before: #HIDDEN_END_pejz
+.. jupyter-execute::
+
+ def main():
+ lead = make_lead().finalized()
+ kwant.plotter.bands(lead, show=False)
+ pyplot.xlabel("momentum [(lattice constant)^-1]")
+ pyplot.ylabel("energy [t]")
+ pyplot.show()
This gives the result:
-.. image:: /code/figure/band_structure_result.*
+.. jupyter-execute::
+ :hide-code:
+
+ # 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()
where we observe the cosine-like dispersion of the square lattice. Close
to ``k=0`` this agrees well with the quadratic dispersion this tight-binding
@@ -61,7 +114,34 @@ Closed systems
.. seealso::
The complete source code of this example can be found in
- :download:`closed_system.py </code/download/closed_system.py>`
+ :jupyter-download:script:`closed_system`
+
+.. jupyter-kernel::
+ :id: closed_system
+
+.. jupyter-execute::
+ :hide-code:
+
+ # 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
+ from matplotlib import pyplot
+ import kwant
+
+.. jupyter-execute:: boilerplate.py
+ :hide-code:
Although Kwant is (currently) mainly aimed towards transport problems, it
can also easily be used to compute properties of closed systems -- after
@@ -74,16 +154,49 @@ To compute the eigenenergies and eigenstates, we will make use of the sparse
linear algebra functionality of `SciPy <https://www.scipy.org>`_, which
interfaces the ARPACK package:
-.. literalinclude:: /code/include/closed_system.py
- :start-after: #HIDDEN_BEGIN_tibv
- :end-before: #HIDDEN_END_tibv
+
+.. jupyter-execute::
+
+ # For eigenvalue computation
+ import scipy.sparse.linalg as sla
We set up the system using the `shape`-function as in
:ref:`tutorial-abring`, but do not add any leads:
-.. literalinclude:: /code/include/closed_system.py
- :start-after: #HIDDEN_BEGIN_qlyd
- :end-before: #HIDDEN_END_qlyd
+.. jupyter-execute::
+ :hide-code:
+
+ a = 1
+ t = 1.0
+ r = 10
+
+.. jupyter-execute::
+
+ def make_system(a=1, t=1.0, r=10):
+
+ lat = kwant.lattice.square(a, norbs=1)
+
+ 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):
+ # 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
We add the magnetic field using a function and a global variable as we
did in the two previous tutorial. (Here, the gauge is chosen such that
@@ -93,14 +206,40 @@ The spectrum can be obtained by diagonalizing the Hamiltonian of the
system, which in turn can be obtained from the finalized
system using `~kwant.system.System.hamiltonian_submatrix`:
-.. literalinclude:: /code/include/closed_system.py
- :start-after: #HIDDEN_BEGIN_yvri
- :end-before: #HIDDEN_END_yvri
+.. jupyter-execute::
+
+ def plot_spectrum(syst, Bfields):
+
+ energies = []
+ for B in Bfields:
+ # Obtain the Hamiltonian as a sparse matrix
+ ham_mat = syst.hamiltonian_submatrix(params=dict(B=B), sparse=True)
+
+ # we only calculate the 15 lowest eigenvalues
+ ev = sla.eigsh(ham_mat.tocsc(), k=15, sigma=0,
+ return_eigenvectors=False)
+
+ energies.append(ev)
+
+ pyplot.figure()
+ pyplot.plot(Bfields, energies)
+ pyplot.xlabel("magnetic field [arbitrary units]")
+ pyplot.ylabel("energy [t]")
+ pyplot.show()
Note that we use sparse linear algebra to efficiently calculate only a
few lowest eigenvalues. Finally, we obtain the result:
-.. image:: /code/figure/closed_system_result.*
+.. jupyter-execute::
+ :hide-code:
+
+ syst = make_system()
+
+ syst = syst.finalized()
+
+ # 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)])
At zero magnetic field several energy levels are degenerate (since our
quantum dot is rather symmetric). These degeneracies are split
@@ -110,11 +249,34 @@ Landau level energies at higher magnetic fields [#]_.
The eigenvectors are obtained very similarly, and can be plotted directly
using `~kwant.plotter.map`:
-.. literalinclude:: /code/include/closed_system.py
- :start-after: #HIDDEN_BEGIN_wave
- :end-before: #HIDDEN_END_wave
+.. jupyter-execute::
+ :hide-code:
+
+ def sorted_eigs(ev):
+ evals, evecs = ev
+ evals, evecs = map(np.array, zip(*sorted(zip(evals, evecs.transpose()))))
+ return evals, evecs.transpose()
+
+.. jupyter-execute::
+
+ def plot_wave_function(syst, B=0.001):
+ # Calculate the wave functions in the system.
+ ham_mat = syst.hamiltonian_submatrix(sparse=True, params=dict(B=B))
+ evals, evecs = sorted_eigs(sla.eigsh(ham_mat.tocsc(), k=20, sigma=0))
+
+ # Plot the probability density of the 10th eigenmode.
+ kwant.plotter.map(syst, np.abs(evecs[:, 9])**2,
+ colorbar=False, oversampling=1)
-.. image:: /code/figure/closed_system_eigenvector.*
+.. jupyter-execute::
+ :hide-code:
+
+ syst = make_system(r=30)
+
+ # 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);
The last two arguments to `~kwant.plotter.map` are optional. The first prevents
a colorbar from appearing. The second, ``oversampling=1``, makes the image look
@@ -127,11 +289,22 @@ that they can have *non-zero* local current. We can calculate the local
current due to a state by using `kwant.operator.Current` and plotting
it using `kwant.plotter.current`:
-.. literalinclude:: /code/include/closed_system.py
- :start-after: #HIDDEN_BEGIN_current
- :end-before: #HIDDEN_END_current
+.. jupyter-execute::
+
+ def plot_current(syst, B=0.001):
+ # Calculate the wave functions in the system.
+ ham_mat = syst.hamiltonian_submatrix(sparse=True, params=dict(B=B))
+ evals, evecs = sorted_eigs(sla.eigsh(ham_mat.tocsc(), k=20, sigma=0))
+
+ # Calculate and plot the local current of the 10th eigenmode.
+ J = kwant.operator.Current(syst)
+ current = J(evecs[:, 9], params=dict(B=B))
+ kwant.plotter.current(syst, current, colorbar=False)
+
+.. jupyter-execute::
+ :hide-code:
-.. image:: /code/figure/closed_system_current.*
+ plot_current(syst);
.. specialnote:: Technical details
diff --git a/doc/source/tutorial/spin_potential_shape.rst b/doc/source/tutorial/spin_potential_shape.rst
index 02ac1fbad566ad8daa570f26677e86c3d83fd01d..770ba4a55f06f0c4a8bf0b89669cb66bb2fad407 100644
--- a/doc/source/tutorial/spin_potential_shape.rst
+++ b/doc/source/tutorial/spin_potential_shape.rst
@@ -11,7 +11,10 @@ Matrix structure of on-site and hopping elements
.. seealso::
The complete source code of this example can be found in
- :download:`spin_orbit.py </code/download/spin_orbit.py>`
+ :jupyter-download:script:`spin_orbit`
+
+.. jupyter-kernel::
+ :id: spin_orbit
We begin by extending the simple 2DEG-Hamiltonian by a Rashba spin-orbit
coupling and a Zeeman splitting due to an external magnetic field:
@@ -42,24 +45,78 @@ use matrices in our program, we import the Tinyarray package. (`NumPy
<http://www.numpy.org/>`_ would work as well, but Tinyarray is much faster
for small arrays.)
-.. literalinclude:: /code/include/spin_orbit.py
- :start-after: #HIDDEN_BEGIN_xumz
- :end-before: #HIDDEN_END_xumz
+.. jupyter-execute::
+ :hide-code:
+
+ # Tutorial 2.3.1. Matrix structure of on-site and hopping elements
+ # ================================================================
+ #
+ # Physics background
+ # ------------------
+ # Gaps in quantum wires with spin-orbit coupling and Zeeman splititng,
+ # as theoretically predicted in
+ # http://prl.aps.org/abstract/PRL/v90/i25/e256601
+ # and (supposedly) experimentally oberved in
+ # http://www.nature.com/nphys/journal/v6/n5/abs/nphys1626.html
+ #
+ # Kwant features highlighted
+ # --------------------------
+ # - Numpy matrices as values in Builder
+
+ import kwant
+
+ # For plotting
+ from matplotlib import pyplot
+
+.. jupyter-execute:: boilerplate.py
+ :hide-code:
+
+.. jupyter-execute::
+
+ # For matrix support
+ import tinyarray
For convenience, we define the Pauli-matrices first (with :math:`\sigma_0` the
unit matrix):
-.. literalinclude:: /code/include/spin_orbit.py
- :start-after: #HIDDEN_BEGIN_hwbt
- :end-before: #HIDDEN_END_hwbt
+.. jupyter-execute::
+
+ # define Pauli-matrices for convenience
+ sigma_0 = tinyarray.array([[1, 0], [0, 1]])
+ sigma_x = tinyarray.array([[0, 1], [1, 0]])
+ sigma_y = tinyarray.array([[0, -1j], [1j, 0]])
+ sigma_z = tinyarray.array([[1, 0], [0, -1]])
+
+and we also define some other parameters useful for constructing our system:
+
+.. jupyter-execute::
+
+ t = 1.0
+ alpha = 0.5
+ e_z = 0.08
+ W, L = 10, 30
Previously, we used numbers as the values of our matrix elements.
However, `~kwant.builder.Builder` also accepts matrices as values, and
we can simply write:
-.. literalinclude:: /code/include/spin_orbit.py
- :start-after: #HIDDEN_BEGIN_uxrm
- :end-before: #HIDDEN_END_uxrm
+.. jupyter-execute::
+ :hide-code:
+
+ lat = kwant.lattice.square()
+ syst = kwant.Builder()
+
+.. jupyter-execute::
+
+ #### Define the scattering region. ####
+ syst[(lat(x, y) for x in range(L) for y in range(W))] = \
+ 4 * t * sigma_0 + e_z * sigma_z
+ # hoppings in x-direction
+ syst[kwant.builder.HoppingKind((1, 0), lat, lat)] = \
+ -t * sigma_0 + 1j * alpha * sigma_y / 2
+ # hoppings in y-directions
+ syst[kwant.builder.HoppingKind((0, 1), lat, lat)] = \
+ -t * sigma_0 - 1j * alpha * sigma_x / 2
Note that the Zeeman energy adds to the onsite term, whereas the Rashba
spin-orbit term adds to the hoppings (due to the derivative operator).
@@ -85,14 +142,49 @@ when specifying `(1, 0)` it is not necessary to specify `(-1, 0)`),
The leads also allow for a matrix structure,
-.. literalinclude:: /code/include/spin_orbit.py
- :start-after: #HIDDEN_BEGIN_yliu
- :end-before: #HIDDEN_END_yliu
+
+.. jupyter-execute::
+ :hide-code:
+
+ #### Define the left lead. ####
+ lead = kwant.Builder(kwant.TranslationalSymmetry((-1, 0)))
+
+.. jupyter-execute::
+
+ lead[(lat(0, j) for j in range(W))] = 4 * t * sigma_0 + e_z * sigma_z
+ # hoppings in x-direction
+ lead[kwant.builder.HoppingKind((1, 0), lat, lat)] = \
+ -t * sigma_0 + 1j * alpha * sigma_y / 2
+ # hoppings in y-directions
+ lead[kwant.builder.HoppingKind((0, 1), lat, lat)] = \
+ -t * sigma_0 - 1j * alpha * sigma_x / 2
+
+.. jupyter-execute::
+ :hide-code:
+
+ #### Attach the leads and finalize the system. ####
+ syst.attach_lead(lead)
+ syst.attach_lead(lead.reversed())
+ syst = syst.finalized()
The remainder of the code is unchanged, and as a result we should obtain
the following, clearly non-monotonic conductance steps:
-.. image:: /code/figure/spin_orbit_result.*
+.. jupyter-execute::
+ :hide-code:
+
+ # Compute conductance
+ energies=[0.01 * i - 0.3 for i in range(100)]
+ data = []
+ for energy in energies:
+ smatrix = kwant.smatrix(syst, energy)
+ data.append(smatrix.transmission(1, 0))
+
+ pyplot.figure()
+ pyplot.plot(energies, data)
+ pyplot.xlabel("energy [t]")
+ pyplot.ylabel("conductance [e^2/h]")
+ pyplot.show()
.. specialnote:: Technical details
@@ -129,7 +221,32 @@ Spatially dependent values through functions
.. seealso::
The complete source code of this example can be found in
- :download:`quantum_well.py </code/download/quantum_well.py>`
+ :jupyter-download:script:`quantum_well`
+
+.. jupyter-kernel::
+ :id: quantum_well
+
+.. jupyter-execute::
+ :hide-code:
+
+ # Tutorial 2.3.2. Spatially dependent values through functions
+ # ============================================================
+ #
+ # Physics background
+ # ------------------
+ # transmission through a quantum well
+ #
+ # Kwant features highlighted
+ # --------------------------
+ # - Functions as values in Builder
+
+ import kwant
+
+ # For plotting
+ from matplotlib import pyplot
+
+.. jupyter-execute:: boilerplate.py
+ :hide-code:
Up to now, all examples had position-independent matrix-elements
(and thus translational invariance along the wire, which
@@ -148,22 +265,50 @@ changing the potential then implies the need to build up the system again.
Instead, we use a python *function* to define the onsite energies. We
define the potential profile of a quantum well as:
-.. literalinclude:: /code/include/quantum_well.py
- :start-after: #HIDDEN_BEGIN_ehso
- :end-before: #HIDDEN_END_ehso
+.. jupyter-execute::
+
+ W, L, L_well = 10, 30, 10
+
+ def potential(site, pot):
+ (x, y) = site.pos
+ if (L - L_well) / 2 < x < (L + L_well) / 2:
+ return pot
+ else:
+ return 0
This function takes two arguments: the first of type `~kwant.builder.Site`,
from which you can get the real-space coordinates using ``site.pos``, and the
value of the potential as the second. Note that in `potential` we can access
-variables of the surrounding function: `L` and `L_well` are taken from the
-namespace of `make_system`.
+variables `L` and `L_well` that are defined globally.
Kwant now allows us to pass a function as a value to
`~kwant.builder.Builder`:
-.. literalinclude:: /code/include/quantum_well.py
- :start-after: #HIDDEN_BEGIN_coid
- :end-before: #HIDDEN_END_coid
+.. jupyter-execute::
+
+ a = 1
+ t = 1.0
+
+ def onsite(site, pot):
+ return 4 * t + potential(site, pot)
+
+ lat = kwant.lattice.square(a)
+ syst = kwant.Builder()
+
+ syst[(lat(x, y) for x in range(L) for y in range(W))] = onsite
+ syst[lat.neighbors()] = -t
+
+.. jupyter-execute::
+ :hide-code:
+
+ #### Define and attach the leads. ####
+ lead = kwant.Builder(kwant.TranslationalSymmetry((-a, 0)))
+ lead[(lat(0, j) for j in range(W))] = 4 * t
+ lead[lat.neighbors()] = -t
+ syst.attach_lead(lead)
+ syst.attach_lead(lead.reversed())
+
+ syst = syst.finalized()
For each lattice point, the corresponding site is then passed as the
first argument to the function `onsite`. The values of any additional
@@ -180,9 +325,21 @@ of the lead -- this should be kept in mind.
Finally, we compute the transmission probability:
-.. literalinclude:: /code/include/quantum_well.py
- :start-after: #HIDDEN_BEGIN_sqvr
- :end-before: #HIDDEN_END_sqvr
+.. jupyter-execute::
+
+ def plot_conductance(syst, energy, welldepths):
+
+ # Compute conductance
+ data = []
+ for welldepth in welldepths:
+ smatrix = kwant.smatrix(syst, energy, params=dict(pot=-welldepth))
+ data.append(smatrix.transmission(1, 0))
+
+ pyplot.figure()
+ pyplot.plot(welldepths, data)
+ pyplot.xlabel("well depth [t]")
+ pyplot.ylabel("conductance [e^2/h]")
+ pyplot.show()
``kwant.smatrix`` allows us to specify a dictionary, `params`, that contains
the additional arguments required by the Hamiltonian matrix elements.
@@ -191,7 +348,11 @@ of the potential well by passing the potential value (remember above
we defined our `onsite` function that takes a parameter named `pot`).
We obtain the result:
-.. image:: /code/figure/quantum_well_result.*
+.. jupyter-execute::
+ :hide-code:
+
+ plot_conductance(syst, energy=0.2,
+ welldepths=[0.01 * i for i in range(100)])
Starting from no potential (well depth = 0), we observe the typical
oscillatory transmission behavior through resonances in the quantum well.
@@ -218,13 +379,44 @@ Nontrivial shapes
.. seealso::
The complete source code of this example can be found in
- :download:`ab_ring.py </code/download/ab_ring.py>`
+ :jupyter-download:script:`ab_ring`
+
+.. jupyter-kernel::
+ :id: ab_ring
+
+.. jupyter-execute::
+ :hide-code:
+
+ # 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
+
+.. jupyter-execute:: boilerplate.py
+ :hide-code:
Up to now, we only dealt with simple wire geometries. Now we turn to the case
of a more complex geometry, namely transport through a quantum ring
that is pierced by a magnetic flux :math:`\Phi`:
-.. image:: /code/figure/ab_ring_sketch.*
+.. image:: /figure/ab_ring_sketch.*
For a flux line, it is possible to choose a gauge such that a
charged particle acquires a phase :math:`e\Phi/h` whenever it
@@ -238,9 +430,14 @@ First, define a boolean function defining the desired shape, i.e. a function
that returns ``True`` whenever a point is inside the shape, and
``False`` otherwise:
-.. literalinclude:: /code/include/ab_ring.py
- :start-after: #HIDDEN_BEGIN_eusz
- :end-before: #HIDDEN_END_eusz
+.. jupyter-execute::
+
+ r1, r2 = 10, 20
+
+ def ring(pos):
+ (x, y) = pos
+ rsq = x ** 2 + y ** 2
+ return (r1 ** 2 < rsq < r2 ** 2)
Note that this function takes a real-space position as argument (not a
`~kwant.builder.Site`).
@@ -249,9 +446,16 @@ We can now simply add all of the lattice points inside this shape at
once, using the function `~kwant.lattice.Square.shape`
provided by the lattice:
-.. literalinclude:: /code/include/ab_ring.py
- :start-after: #HIDDEN_BEGIN_lcak
- :end-before: #HIDDEN_END_lcak
+.. jupyter-execute::
+
+ a = 1
+ t = 1.0
+
+ lat = kwant.lattice.square(a)
+ syst = kwant.Builder()
+
+ syst[lat.shape(ring, (0, r1 + 1))] = 4 * t
+ syst[lat.neighbors()] = -t
Here, ``lat.shape`` takes as a second parameter a (real-space) point that is
inside the desired shape. The hoppings can still be added using
@@ -264,9 +468,30 @@ along the branch cut in the lower arm of the ring. For this we select
all hoppings in x-direction that are of the form `(lat(1, j), lat(0, j))`
with ``j<0``:
-.. literalinclude:: /code/include/ab_ring.py
- :start-after: #HIDDEN_BEGIN_lvkt
- :end-before: #HIDDEN_END_lvkt
+.. jupyter-execute::
+
+ # In order to introduce a flux through the ring, we introduce a phase on
+ # the hoppings on the line cut through one of the arms. Since we want to
+ # change the flux without modifying the Builder instance repeatedly, we
+ # define the modified hoppings as a function that takes the flux as its
+ # parameter phi.
+ def 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
Here, `crosses_branchcut` is a boolean function that returns ``True`` for
the desired hoppings. We then use again a generator (this time with
@@ -279,9 +504,21 @@ by the parameter `phi`.
For the leads, we can also use the ``lat.shape``-functionality:
-.. literalinclude:: /code/include/ab_ring.py
- :start-after: #HIDDEN_BEGIN_qwgr
- :end-before: #HIDDEN_END_qwgr
+.. jupyter-execute::
+
+ #### Define the leads. ####
+ W = 10
+
+ 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
Here, the shape must be compatible with the translational symmetry
of the lead ``sym_lead``. In particular, this means that it should extend to
@@ -290,9 +527,12 @@ no restriction on ``x`` in ``lead_shape``) [#]_.
Attaching the leads is done as before:
-.. literalinclude:: /code/include/ab_ring.py
- :start-after: #HIDDEN_BEGIN_skbz
- :end-before: #HIDDEN_END_skbz
+.. jupyter-execute::
+ :hide-output:
+
+ #### Attach the leads ####
+ syst.attach_lead(lead)
+ syst.attach_lead(lead.reversed())
In fact, attaching leads seems not so simple any more for the current
structure with a scattering region very much different from the lead
@@ -305,16 +545,40 @@ back (going opposite to the direction of the translational vector)
until it intersects the scattering region. At this intersection,
the lead is attached:
-.. image:: /code/figure/ab_ring_sketch2.*
+.. image:: /figure/ab_ring_sketch2.*
After the lead has been attached, the system should look like this:
-.. image:: /code/figure/ab_ring_syst.*
+.. jupyter-execute::
+ :hide-code:
+
+ kwant.plot(syst);
The computation of the conductance goes in the same fashion as before.
Finally you should get the following result:
-.. image:: /code/figure/ab_ring_result.*
+
+.. jupyter-execute::
+ :hide-code:
+
+ 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, params=dict(phi=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()
+
+ # We should see a conductance that is periodic with the flux quantum
+ plot_conductance(syst.finalized(), energy=0.15,
+ fluxes=[0.01 * i * 3 * 2 * pi for i in range(100)])
where one can observe the conductance oscillations with the
period of one flux quantum.
@@ -330,7 +594,47 @@ period of one flux quantum.
becomes more apparent if we attach the leads a bit further away
from the central axis o the ring, as was done in this example:
- .. image:: /code/figure/ab_ring_note1.*
+ .. jupyter-kernel::
+ :id: ab_ring_note1
+
+ .. jupyter-execute::
+ :hide-code:
+
+ import kwant
+ from matplotlib import pyplot
+
+ .. jupyter-execute:: boilerplate.py
+ :hide-code:
+
+ .. jupyter-execute::
+ :hide-code:
+
+ a = 1
+ t = 1.0
+ W = 10
+ r1, r2 = 10, 20
+
+ lat = kwant.lattice.square()
+ syst = kwant.Builder()
+ def ring(pos):
+ (x, y) = pos
+ rsq = x**2 + y**2
+ return ( r1**2 < rsq < r2**2)
+ syst[lat.shape(ring, (0, 11))] = 4 * t
+ syst[lat.neighbors()] = -t
+ sym_lead0 = kwant.TranslationalSymmetry((-a, 0))
+ lead0 = kwant.Builder(sym_lead0)
+ def lead_shape(pos):
+ (x, y) = pos
+ return (-1 < x < 1) and ( 0.5 * W < y < 1.5 * W )
+ lead0[lat.shape(lead_shape, (0, W))] = 4 * t
+ lead0[lat.neighbors()] = -t
+ lead1 = lead0.reversed()
+ syst.attach_lead(lead0)
+ syst.attach_lead(lead1)
+
+ kwant.plot(syst);
+
- Per default, `~kwant.builder.Builder.attach_lead` attaches
the lead to the "outside" of the structure, by tracing the
@@ -345,7 +649,46 @@ period of one flux quantum.
starts the trace-back in the middle of the ring, resulting
in the lead being attached to the inner circle:
- .. image:: /code/figure/ab_ring_note2.*
+ .. jupyter-kernel::
+ :id: ab_ring_note2
+
+ .. jupyter-execute::
+ :hide-code:
+
+ import kwant
+ from matplotlib import pyplot
+
+ .. jupyter-execute:: boilerplate.py
+ :hide-code:
+
+ .. jupyter-execute::
+ :hide-code:
+
+ a = 1
+ t = 1.0
+ W = 10
+ r1, r2 = 10, 20
+
+ lat = kwant.lattice.square(a)
+ syst = kwant.Builder()
+ def ring(pos):
+ (x, y) = pos
+ rsq = x**2 + y**2
+ return ( r1**2 < rsq < r2**2)
+ syst[lat.shape(ring, (0, 11))] = 4 * t
+ syst[lat.neighbors()] = -t
+ sym_lead0 = kwant.TranslationalSymmetry((-a, 0))
+ lead0 = kwant.Builder(sym_lead0)
+ def lead_shape(pos):
+ (x, y) = pos
+ return (-1 < x < 1) and ( -W/2 < y < W/2 )
+ lead0[lat.shape(lead_shape, (0, 0))] = 4 * t
+ lead0[lat.neighbors()] = -t
+ lead1 = lead0.reversed()
+ syst.attach_lead(lead0)
+ syst.attach_lead(lead1, lat(0, 0))
+
+ kwant.plot(syst);
Note that here the lead is treated as if it would pass over
the other arm of the ring, without intersecting it.
diff --git a/doc/source/tutorial/superconductors.rst b/doc/source/tutorial/superconductors.rst
index 6113cb78a7be06bf57b7bf7f39b65095505b1f29..78f004d34820c3c4ef6219a855681ab542f0aa7a 100644
--- a/doc/source/tutorial/superconductors.rst
+++ b/doc/source/tutorial/superconductors.rst
@@ -3,7 +3,40 @@ Superconductors: orbital degrees of freedom, conservation laws and symmetries
.. seealso::
The complete source code of this example can be found in
- :download:`superconductor.py </code/download/superconductor.py>`
+ :jupyter-download:script:`superconductor`
+
+.. jupyter-kernel::
+ :id: superconductor
+
+.. jupyter-execute::
+ :hide-code:
+
+ # Tutorial 2.6. "Superconductors": orbitals, conservation laws and symmetries
+ # ===========================================================================
+ #
+ # Physics background
+ # ------------------
+ # - conductance of a NS-junction (Andreev reflection, superconducting gap)
+ #
+ # Kwant features highlighted
+ # --------------------------
+ # - Implementing electron and hole ("orbital") degrees of freedom
+ # using conservation laws.
+ # - Use of discrete symmetries to relate scattering states.
+
+ import kwant
+
+ from matplotlib import pyplot
+
+ import tinyarray
+ import numpy as np
+
+ tau_x = tinyarray.array([[0, 1], [1, 0]])
+ tau_y = tinyarray.array([[0, -1j], [1j, 0]])
+ tau_z = tinyarray.array([[1, 0], [0, -1]])
+
+.. jupyter-execute:: boilerplate.py
+ :hide-code:
This example deals with superconductivity on the level of the
Bogoliubov-de Gennes (BdG) equation. In this framework, the Hamiltonian
@@ -51,16 +84,48 @@ of freedom.
Let us consider a system that consists of a normal lead on the left,
a superconductor on the right, and a tunnel barrier in between:
-.. image:: /code/figure/superconductor_transport_sketch.*
+.. image:: /figure/superconductor_transport_sketch.*
We implement the BdG Hamiltonian in Kwant using a 2x2 matrix structure
for all Hamiltonian matrix elements, as we did
-previously in the :ref:`spin example <tutorial_spinorbit>`. We declare
-the square lattice and construct the scattering region with the following:
+previously in the :ref:`spin example <tutorial_spinorbit>`.
+We start by declaring some parameters that will be used in the following code:
+
+.. jupyter-execute::
+
+ a = 1
+ W, L = 10, 10
+ barrier = 1.5
+ barrierpos = (3, 4)
+ mu = 0.4
+ Delta = 0.1
+ Deltapos = 4
+ t = 1.0
+
+and we declare the square lattice and construct the scattering region with the following:
+
+.. jupyter-execute::
-.. literalinclude:: /code/include/superconductor.py
- :start-after: #HIDDEN_BEGIN_nbvn
- :end-before: #HIDDEN_END_nbvn
+ # Start with an empty tight-binding system. On each site, there
+ # are now electron and hole orbitals, so we must specify the
+ # number of orbitals per site. The orbital structure is the same
+ # as in the Hamiltonian.
+ lat = kwant.lattice.square(norbs=2)
+ syst = kwant.Builder()
+
+ #### Define the scattering region. ####
+ # The superconducting order parameter couples electron and hole orbitals
+ # on each site, and hence enters as an onsite potential.
+ # The pairing is only included beyond the point 'Deltapos' in the scattering region.
+ syst[(lat(x, y) for x in range(Deltapos) for y in range(W))] = (4 * t - mu) * tau_z
+ syst[(lat(x, y) for x in range(Deltapos, L) for y in range(W))] = (4 * t - mu) * tau_z + Delta * tau_x
+
+ # The tunnel barrier
+ syst[(lat(x, y) for x in range(barrierpos[0], barrierpos[1])
+ for y in range(W))] = (4 * t + barrier - mu) * tau_z
+
+ # Hoppings
+ syst[lat.neighbors()] = -t * tau_z
Note the new argument ``norbs`` to `~kwant.lattice.square`. This is
the number of orbitals per site in the discretized BdG Hamiltonian - of course,
@@ -80,9 +145,16 @@ of it). The next step towards computing conductance is to attach leads.
Let's attach two leads: a normal one to the left end, and a superconducting
one to the right end. Starting with the left lead, we have:
-.. literalinclude:: /code/include/superconductor.py
- :start-after: #HIDDEN_BEGIN_ttth
- :end-before: #HIDDEN_END_ttth
+.. jupyter-execute::
+
+ #### Define the leads. ####
+ # Left lead - normal, so the order parameter is zero.
+ sym_left = kwant.TranslationalSymmetry((-a, 0))
+ # Specify the conservation law used to treat electrons and holes separately.
+ # We only do this in the left lead, where the pairing is zero.
+ lead0 = kwant.Builder(sym_left, conservation_law=-tau_z, particle_hole=tau_y)
+ lead0[(lat(0, j) for j in range(W))] = (4 * t - mu) * tau_z
+ lead0[lat.neighbors()] = -t * tau_z
Note the two new new arguments in `~kwant.builder.Builder`, ``conservation_law``
and ``particle_hole``. For the purpose of computing conductance, ``conservation_law``
@@ -120,9 +192,19 @@ of ascending eigenvalues of the conservation law.
In order to move on with the conductance calculation, let's attach the second
lead to the right side of the scattering region:
-.. literalinclude:: /code/include/superconductor.py
- :start-after: #HIDDEN_BEGIN_zuuw
- :end-before: #HIDDEN_END_zuuw
+.. jupyter-execute::
+
+ # Right lead - superconducting, so the order parameter is included.
+ sym_right = kwant.TranslationalSymmetry((a, 0))
+ lead1 = kwant.Builder(sym_right)
+ lead1[(lat(0, j) for j in range(W))] = (4 * t - mu) * tau_z + Delta * tau_x
+ lead1[lat.neighbors()] = -t * tau_z
+
+ #### Attach the leads and finalize the system. ####
+ syst.attach_lead(lead0)
+ syst.attach_lead(lead1)
+
+ syst = syst.finalized()
The second (right) lead is superconducting, such that the electron and hole
blocks are coupled. Of course, this means that we can not separate them into
@@ -140,9 +222,23 @@ confused by the fact that it says ``transmission`` -- transmission
into the same lead is reflection), and reflection from electrons to holes
is ``smatrix.transmission((0, 1), (0, 0))``:
-.. literalinclude:: /code/include/superconductor.py
- :start-after: #HIDDEN_BEGIN_jbjt
- :end-before: #HIDDEN_END_jbjt
+.. jupyter-execute::
+
+ def plot_conductance(syst, energies):
+ # Compute conductance
+ data = []
+ for energy in energies:
+ smatrix = kwant.smatrix(syst, energy)
+ # Conductance is N - R_ee + R_he
+ data.append(smatrix.submatrix((0, 0), (0, 0)).shape[0] -
+ smatrix.transmission((0, 0), (0, 0)) +
+ smatrix.transmission((0, 1), (0, 0)))
+
+ pyplot.figure()
+ pyplot.plot(energies, data)
+ pyplot.xlabel("energy [t]")
+ pyplot.ylabel("conductance [e^2/h]")
+ pyplot.show()
Note that ``smatrix.submatrix((0, 0), (0, 0))`` returns the block concerning
reflection of electrons to electrons, and from its size we can extract the number of modes
@@ -150,7 +246,10 @@ reflection of electrons to electrons, and from its size we can extract the numbe
For the default parameters, we obtain the following conductance:
-.. image:: /code/figure/superconductor_transport_result.*
+.. jupyter-execute::
+ :hide-code:
+
+ plot_conductance(syst, energies=[0.002 * i for i in range(-10, 100)])
We a see a conductance that is proportional to the square of the tunneling
probability within the gap, and proportional to the tunneling probability
@@ -176,13 +275,34 @@ the electron and hole blocks. The exception is of course at zero energy, in whic
case particle-hole symmetry transforms between the electron and hole blocks, resulting
in a symmetric scattering matrix. We can check the symmetry explicitly with
-.. literalinclude:: /code/include/superconductor.py
- :start-after: #HIDDEN_BEGIN_pqmp
- :end-before: #HIDDEN_END_pqmp
+.. jupyter-execute::
+
+ def check_PHS(syst):
+ # Scattering matrix
+ s = kwant.smatrix(syst, energy=0)
+ # Electron to electron block
+ s_ee = s.submatrix((0,0), (0,0))
+ # Hole to hole block
+ s_hh = s.submatrix((0,1), (0,1))
+ print('s_ee: \n', np.round(s_ee, 3))
+ print('s_hh: \n', np.round(s_hh[::-1, ::-1], 3))
+ print('s_ee - s_hh^*: \n',
+ np.round(s_ee - s_hh[::-1, ::-1].conj(), 3), '\n')
+ # Electron to hole block
+ s_he = s.submatrix((0,1), (0,0))
+ # Hole to electron block
+ s_eh = s.submatrix((0,0), (0,1))
+ print('s_he: \n', np.round(s_he, 3))
+ print('s_eh: \n', np.round(s_eh[::-1, ::-1], 3))
+ print('s_he + s_eh^*: \n',
+ np.round(s_he + s_eh[::-1, ::-1].conj(), 3))
which yields the output
-.. literalinclude:: /code/figure/check_PHS_out.txt
+.. jupyter-execute::
+ :hide-code:
+
+ check_PHS(syst)
Note that :math:`\mathcal{P}` flips the sign of momentum, and for the parameters
we consider here, there are two electron and two hole modes active at zero energy.