From 001dafbb70c7902129b8c62e100f90396b995fab Mon Sep 17 00:00:00 2001
From: Christoph Groth <christoph.groth@cea.fr>
Date: Wed, 19 Jul 2017 17:49:01 +0200
Subject: [PATCH] improve the tutorial script and image buiding machinery
The basic idea remains the same, but should work much better now. The
difference is that images/example.py.diff is now a patch with complete
context and becomes the authoritative source for both the visible flavor
of an example as well as for its figure-generating variant. Both
tutorial/example.py and images/example.py are extracted from this file
by 'make html'. Thanks to the complete context the diffs are quite
readable and may be modified directly.
Alternatively, one may also modify the generated scripts. When
tutorial/example.py has been modified, it will be patched and saved as
images/example.py. The patching is done using the tool 'wiggle' that
works much better than 'patch'. If a conflict occurs, conflict markers
are added to the output file and its dated back to the dawn of time
(i.e. 1970) in order to mark the conflicts as not yet resolved. After
resolving, 'make html' is simply run again.
Upon a successful execution of the figure-generating script the diff
gets recreated automatically, 'generate-diffs' no longer exists.
---
.gitignore | 1 +
INSTALL.rst | 18 +
doc/Makefile | 63 ++-
doc/source/images/ab_ring.py.diff | 98 +++-
doc/source/images/band_structure.py.diff | 41 +-
doc/source/images/closed_system.py.diff | 101 +++-
doc/source/images/discretize.py.diff | 173 ++++++-
doc/source/images/faq.py.diff | 333 ++++++++++++-
doc/source/images/generate-diffs.sh | 16 -
doc/source/images/graphene.py.diff | 141 +++++-
.../images/kernel_polynomial_method.py.diff | 161 ++++++-
doc/source/images/magnetic_texture.py.diff | 169 ++++++-
doc/source/images/plot_graphene.py.diff | 70 ++-
doc/source/images/plot_qahe.py.diff | 63 ++-
doc/source/images/plot_zincblende.py.diff | 44 +-
doc/source/images/quantum_well.py.diff | 68 ++-
doc/source/images/quantum_wire.py.diff | 106 ++++-
doc/source/images/spin_orbit.py.diff | 84 +++-
doc/source/images/superconductor.py.diff | 102 +++-
doc/source/tutorial/README | 11 -
doc/source/tutorial/ab_ring.py | 127 -----
doc/source/tutorial/band_structure.py | 52 --
doc/source/tutorial/closed_system.py | 140 ------
doc/source/tutorial/discretize.py | 222 ---------
doc/source/tutorial/faq.py | 450 ------------------
doc/source/tutorial/graphene.py | 179 -------
.../tutorial/kernel_polynomial_method.py | 219 ---------
doc/source/tutorial/magnetic_texture.py | 245 ----------
doc/source/tutorial/plot_graphene.py | 112 -----
doc/source/tutorial/plot_qahe.py | 71 ---
doc/source/tutorial/plot_zincblende.py | 59 ---
doc/source/tutorial/quantum_well.py | 88 ----
doc/source/tutorial/quantum_wire.py | 121 -----
doc/source/tutorial/spin_orbit.py | 104 ----
doc/source/tutorial/superconductor.py | 132 -----
35 files changed, 1708 insertions(+), 2476 deletions(-)
delete mode 100755 doc/source/images/generate-diffs.sh
delete mode 100644 doc/source/tutorial/README
delete mode 100644 doc/source/tutorial/ab_ring.py
delete mode 100644 doc/source/tutorial/band_structure.py
delete mode 100644 doc/source/tutorial/closed_system.py
delete mode 100644 doc/source/tutorial/discretize.py
delete mode 100644 doc/source/tutorial/faq.py
delete mode 100644 doc/source/tutorial/graphene.py
delete mode 100644 doc/source/tutorial/kernel_polynomial_method.py
delete mode 100644 doc/source/tutorial/magnetic_texture.py
delete mode 100644 doc/source/tutorial/plot_graphene.py
delete mode 100644 doc/source/tutorial/plot_qahe.py
delete mode 100644 doc/source/tutorial/plot_zincblende.py
delete mode 100644 doc/source/tutorial/quantum_well.py
delete mode 100644 doc/source/tutorial/quantum_wire.py
delete mode 100644 doc/source/tutorial/spin_orbit.py
delete mode 100644 doc/source/tutorial/superconductor.py
diff --git a/.gitignore b/.gitignore
index 64f23822..73576937 100644
--- a/.gitignore
+++ b/.gitignore
@@ -8,6 +8,7 @@
/build
/dist
/doc/build
+/doc/source/tutorial/*.py
/doc/source/reference/generated/
/doc/source/images/*.png
/doc/source/images/*.pdf
diff --git a/INSTALL.rst b/INSTALL.rst
index 8fe4686d..1b0a494f 100644
--- a/INSTALL.rst
+++ b/INSTALL.rst
@@ -171,6 +171,24 @@ 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
+<http://neil.brown.name/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 985a4ef9..8df4d59b 100644
--- a/doc/Makefile
+++ b/doc/Makefile
@@ -35,6 +35,7 @@ GENERATEDPDF = $(patsubst %.svg,%.pdf,$(IMAGESOURCES))
# deleting the corresponding flag file.
SCRIPTS = $(patsubst source/images/%.diff,%,$(wildcard source/images/*.py.diff))
FLAGS = $(patsubst %.py,source/images/.%_flag,$(SCRIPTS))
+INCLUDES = $(patsubst %,source/tutorial/%,$(SCRIPTS))
.PHONY: help clean realclean html dirhtml pickle json htmlhelp qthelp latex changes linkcheck doctest
@@ -59,37 +60,38 @@ clean:
realclean: clean
-rm -f $(FLAGS)
+ -rm -f $(INCLUDES)
-rm -f $(patsubst %,source/images/%,$(SCRIPTS))
-rm -f $(patsubst %.py,source/images/%_*.png,$(SCRIPTS))
-rm -f $(patsubst %.py,source/images/%_*.pdf,$(SCRIPTS))
-html: $(FLAGS)
+html: $(FLAGS) $(INCLUDES)
$(SPHINXBUILD) -b html $(ALLSPHINXOPTS) $(BUILDDIR)/html
@echo
@echo "Build finished. The HTML pages are in $(BUILDDIR)/html."
-dirhtml: $(FLAGS)
+dirhtml: $(FLAGS) $(INCLUDES)
$(SPHINXBUILD) -b dirhtml $(ALLSPHINXOPTS) $(BUILDDIR)/dirhtml
@echo
@echo "Build finished. The HTML pages are in $(BUILDDIR)/dirhtml."
-pickle: $(FLAGS)
+pickle: $(FLAGS) $(INCLUDES)
$(SPHINXBUILD) -b pickle $(ALLSPHINXOPTS) $(BUILDDIR)/pickle
@echo
@echo "Build finished; now you can process the pickle files."
-json: $(FLAGS)
+json: $(FLAGS) $(INCLUDES)
$(SPHINXBUILD) -b json $(ALLSPHINXOPTS) $(BUILDDIR)/json
@echo
@echo "Build finished; now you can process the JSON files."
-htmlhelp: $(FLAGS)
+htmlhelp: $(FLAGS) $(INCLUDES)
$(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: $(FLAGS)
+qthelp: $(FLAGS) $(INCLUDES)
$(SPHINXBUILD) -b qthelp $(ALLSPHINXOPTS) $(BUILDDIR)/qthelp
@echo
@echo "Build finished; now you can run "qcollectiongenerator" with the" \
@@ -98,7 +100,7 @@ qthelp: $(FLAGS)
@echo "To view the help file:"
@echo "# assistant -collectionFile $(BUILDDIR)/qthelp/kwant.qhc"
-latex: $(GENERATEDPDF) $(FLAGS)
+latex: $(GENERATEDPDF) $(FLAGS) $(INCLUDES)
$(SPHINXBUILD) -b latex $(ALLSPHINXOPTS) $(BUILDDIR)/latex
@echo
@echo "Build finished; the LaTeX files are in $(BUILDDIR)/latex."
@@ -124,19 +126,48 @@ doctest:
%.pdf: %.svg
rsvg-convert -f pdf -o $@ $<
-# Make the image generation scripts by patching tutorial scipts.
+# Make tutorial scripts by extracting the (complete!) context of the "patches".
+# We make sure not to use 'wiggle' here.
.SECONDARY:
-%.py: %.py.diff
- @grep -v '^#HIDDEN' source/tutorial/$(notdir $@) >$@
- @patch $@ $<
-
-# The image generation scripts depend on their unpatched originals
+source/tutorial/%.py: source/images/%.py.diff
+ @sed -n '/^[- ]/ s/^.//p' <$< >$@
+ @touch -r $< $@
+
+# Make the image generation scripts by patching tutorial scripts. If the
+# tutorial scripts haven't been modified, don't patch but directly extract the
+# image generation scripts. This means that 'wiggle' is only needed when the
+# tutorial scripts have been modified.
+.SECONDARY:
+source/images/%.py: source/tutorial/%.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 image generation scripts also depend on the diffs.
define makedep
-source/images/$(1): source/tutorial/$(1)
+source/images/$(1): source/images/$(1).diff
endef
$(foreach name,$(SCRIPTS),$(eval $(call makedep,$(name))))
-# Generation of images
-.%_flag: %.py
+# Run an image 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 image
+# generation script is generated from it by patching.
+source/images/.%_flag: source/images/%.py
cd $(dir $<) && python3 $(notdir $<)
+ @if [ ! -f $<.diff -o $< -nt $<.diff ]; then \
+ wiggle --diff --lines source/tutorial/$(notdir $<) $< >$<.diff; \
+ touch -r $< $<.diff; \
+ fi
+ @rm -f $<.porig
@touch $@
diff --git a/doc/source/images/ab_ring.py.diff b/doc/source/images/ab_ring.py.diff
index ded6e142..50a10d19 100644
--- a/doc/source/images/ab_ring.py.diff
+++ b/doc/source/images/ab_ring.py.diff
@@ -1,6 +1,15 @@
---- original
-+++ modified
-@@ -12,6 +12,7 @@
+@@ -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
@@ -8,7 +17,79 @@
from cmath import exp
from math import pi
-@@ -81,6 +82,50 @@
+ 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
@@ -59,7 +140,9 @@
def plot_conductance(syst, energy, fluxes):
# compute conductance
-@@ -90,18 +135,31 @@
+ normalized_fluxes = [flux / (2 * pi) for flux in fluxes]
+ data = []
+ for flux in fluxes:
smatrix = kwant.smatrix(syst, energy, args=[flux])
data.append(smatrix.transmission(1, 0))
@@ -96,7 +179,9 @@
# Finalize the system.
syst = syst.finalized()
-@@ -111,6 +169,17 @@
+
+ # 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)])
@@ -114,3 +199,4 @@
# 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/images/band_structure.py.diff b/doc/source/images/band_structure.py.diff
index 90d73544..8a265207 100644
--- a/doc/source/images/band_structure.py.diff
+++ b/doc/source/images/band_structure.py.diff
@@ -1,6 +1,12 @@
---- original
-+++ modified
-@@ -9,6 +9,7 @@
+@@ -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.
@@ -8,8 +14,31 @@
import kwant
# For plotting.
-@@ -36,10 +37,19 @@
+ 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)
@@ -29,6 +58,10 @@
+ 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/images/closed_system.py.diff b/doc/source/images/closed_system.py.diff
index e33c3c2c..7f1de471 100644
--- a/doc/source/images/closed_system.py.diff
+++ b/doc/source/images/closed_system.py.diff
@@ -1,6 +1,14 @@
---- original
-+++ modified
-@@ -11,6 +11,7 @@
+@@ -1,140 +1,157 @@
+ # 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.
@@ -8,7 +16,62 @@
from cmath import exp
import numpy as np
import kwant
-@@ -68,24 +69,39 @@
+
+ # 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=0):
+ # The magnetic field is controlled by the parameter B
+ y = site1.pos[1]
+ return -t * exp(-1j * B * y)
+
+ syst[lat.shape(circle, (0, 0))] = 4 * t
+ # hoppings in x-direction
+ syst[kwant.builder.HoppingKind((1, 0), lat, lat)] = hopx
+ # hoppings in y-directions
+ syst[kwant.builder.HoppingKind((0, 1), lat, lat)] = -t
+
+ # It's a closed system for a change, so no leads
+ return syst
+ #HIDDEN_END_qlyd
+
+
+ #HIDDEN_BEGIN_yvri
+ def plot_spectrum(syst, Bfields):
+
+ # In the following, we compute the spectrum of the quantum dot
+ # using dense matrix methods. This works in this toy example, as
+ # the system is tiny. In a real example, one would want to use
+ # sparse matrix methods
+
+ energies = []
+ for B in Bfields:
+ # Obtain the Hamiltonian as a dense matrix
+ ham_mat = syst.hamiltonian_submatrix(args=[B], sparse=True)
+
+ # we only calculate the 15 lowest eigenvalues
+ ev = sla.eigsh(ham_mat, k=15, which='SM', return_eigenvectors=False)
energies.append(ev)
@@ -29,8 +92,10 @@
+ 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
+ #HIDDEN_BEGIN_wave
def plot_wave_function(syst):
+ size = (_defs.figwidth_in, _defs.figwidth_in)
+
@@ -46,15 +111,17 @@
+ 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):
+ size = (_defs.figwidth_in, _defs.figwidth_in)
+
# Calculate the wave functions in the system.
ham_mat = syst.hamiltonian_submatrix(sparse=True)
evecs = sla.eigsh(ham_mat, k=20, which='SM')[1]
-@@ -93,15 +109,16 @@
+
# Calculate and plot the local current of the 10th eigenmode.
J = kwant.operator.Current(syst)
current = J(evecs[:, 9])
@@ -64,6 +131,7 @@
+ syst, current, colorbar=False,
+ file="closed_system_current." + extension,
+ fig_size=size, dpi=_defs.dpi)
+ #HIDDEN_END_current
def main():
@@ -75,3 +143,26 @@
# 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/images/discretize.py.diff b/doc/source/images/discretize.py.diff
index 7ac6dfe3..3d6c2024 100644
--- a/doc/source/images/discretize.py.diff
+++ b/doc/source/images/discretize.py.diff
@@ -1,6 +1,12 @@
---- original
-+++ modified
-@@ -9,6 +9,7 @@
+@@ -1,222 +1,236 @@
+ # Tutorial 2.9. Processing continuum Hamiltonians with discretize
+ # ===============================================================
+ #
+ # Physics background
+ # ------------------
+ # - tight-binding approximation of continuous Hamiltonians
+ #
+ # Kwant features highlighted
# --------------------------
# - kwant.continuum.discretize
@@ -8,7 +14,11 @@
import kwant
import scipy.sparse.linalg
-@@ -20,10 +21,19 @@
+ import scipy.linalg
+ import numpy as np
+
+ # For plotting
+ import matplotlib as mpl
from matplotlib import pyplot as plt
@@ -21,34 +31,120 @@
+
+
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
-@@ -44,7 +54,7 @@
+ 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):
-@@ -91,7 +101,8 @@
+ #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_spacing=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)
-@@ -119,7 +130,7 @@
+ # 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')
@@ -57,7 +153,46 @@
def lattice_spacing():
-@@ -160,7 +171,7 @@
+ #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_spacing=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)
@@ -66,25 +201,45 @@
def substitutions():
-@@ -173,15 +184,18 @@
+ #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/images/faq.py.diff b/doc/source/images/faq.py.diff
index 4df4fd94..7c8889dc 100644
--- a/doc/source/images/faq.py.diff
+++ b/doc/source/images/faq.py.diff
@@ -1,6 +1,6 @@
---- original
-+++ modified
-@@ -3,6 +3,7 @@
+@@ -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.
@@ -8,7 +8,8 @@
import kwant
import numpy as np
import tinyarray
-@@ -11,6 +12,12 @@
+ import matplotlib
+ from matplotlib import pyplot as plt
matplotlib.rcParams['figure.figsize'] = (3.5, 3.5)
@@ -20,32 +21,72 @@
+
#### What is a Site?
+ #HIDDEN_BEGIN_site
a = 1
-@@ -21,6 +28,7 @@
+ 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?
-@@ -34,6 +42,7 @@
+
+ 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?
-@@ -50,6 +59,7 @@
+
+ #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)])
-@@ -70,6 +80,7 @@
+ 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)
@@ -53,47 +94,143 @@
# Delete sites to create a hole
-@@ -81,6 +92,7 @@
+ 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?
-@@ -127,7 +139,7 @@
+
+ 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?
-@@ -141,6 +153,7 @@
+ # 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
-@@ -171,12 +184,14 @@
+ 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?
-@@ -190,11 +205,13 @@
+
+ # 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)
@@ -104,66 +241,175 @@
kwant.plot(syst)
+save_figure("faq_adjacent2")
+ #HIDDEN_END_adjacent1
# Polyatomic lattice
-@@ -213,9 +230,11 @@
+ #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
-@@ -253,6 +272,7 @@
+ 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()
-@@ -262,6 +282,7 @@
+ 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
-@@ -269,9 +290,11 @@
+
# 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?
-@@ -291,6 +314,7 @@
+
+ #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):
-@@ -305,6 +329,7 @@
+ 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?
-@@ -337,6 +362,7 @@
+
+ #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)
@@ -171,7 +417,14 @@
# More involved example
-@@ -351,6 +377,7 @@
+ 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)
@@ -179,20 +432,52 @@
#### How does Kwant order components of an individual wavefunction?
-@@ -381,7 +408,8 @@
+ 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)
-@@ -394,4 +422,5 @@
+ 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/images/generate-diffs.sh b/doc/source/images/generate-diffs.sh
deleted file mode 100755
index d78525ea..00000000
--- a/doc/source/images/generate-diffs.sh
+++ /dev/null
@@ -1,16 +0,0 @@
-# !/bin/sh
-
-# This script regenerates the .diff files in this directory. It's these files
-# that are kept under vesion control instead of the scripts themselves.
-
-for f in [a-zA-Z]*.py; do
- # We use custom labels to suppress the time stamps which are unnecessary
- # here and would only lead to noise in version control.
- grep -v '#HIDDEN' ../tutorial/$f |
- diff -u --label original --label modified - $f >$f.diff_
- if cmp $f.diff_ $f.diff >/dev/null 2>&1; then
- rm $f.diff_
- else
- mv $f.diff_ $f.diff
- fi
-done
diff --git a/doc/source/images/graphene.py.diff b/doc/source/images/graphene.py.diff
index 48712132..1352ed27 100644
--- a/doc/source/images/graphene.py.diff
+++ b/doc/source/images/graphene.py.diff
@@ -1,6 +1,13 @@
---- original
-+++ modified
-@@ -10,6 +10,7 @@
+@@ -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
@@ -8,7 +15,102 @@
from math import pi, sqrt, tanh
import kwant
-@@ -96,22 +97,40 @@
+
+ # 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))
@@ -56,8 +158,17 @@
+ 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():
-@@ -123,8 +142,11 @@
+ 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
@@ -68,10 +179,14 @@
+ 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())
-@@ -133,9 +155,11 @@
+ #HIDDEN_END_jmbi
+
+ # Attach the leads to the system.
for lead in leads:
syst.attach_lead(lead)
@@ -86,3 +201,17 @@
# 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/images/kernel_polynomial_method.py.diff b/doc/source/images/kernel_polynomial_method.py.diff
index d58b3836..9e778f10 100644
--- a/doc/source/images/kernel_polynomial_method.py.diff
+++ b/doc/source/images/kernel_polynomial_method.py.diff
@@ -1,6 +1,14 @@
---- original
-+++ modified
-@@ -11,6 +11,9 @@
+@@ -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
@@ -10,7 +18,27 @@
# For plotting
from matplotlib import pyplot as plt
-@@ -36,13 +39,13 @@
+ #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
# Plot several density of states curves on the same axes.
@@ -26,7 +54,14 @@
plt.clf()
-@@ -57,10 +60,18 @@
+ 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(fsyst, axes, titles_to_data, file_name=None):
+ for ax, (title, ldos) in zip(axes, titles_to_data):
+ site_size = site_size_conversion(ldos) # convert LDoS to sizes
kwant.plot(fsyst, site_size=site_size, site_color=(0, 0, 1, 0.3), ax=ax)
ax.set_title(title)
ax.set(adjustable='box-forced', aspect='equal')
@@ -44,9 +79,21 @@
+
+
def simple_dos_example():
+ #HIDDEN_BEGIN_kpm1
fsyst = make_syst().finalized()
-@@ -74,18 +85,24 @@
+ 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)),
@@ -59,11 +106,14 @@
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)
@@ -71,10 +121,19 @@
+ 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):
-@@ -99,7 +116,9 @@
+ 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),
@@ -83,9 +142,13 @@
+ file_name='kpm_dos_acc'
+ )
+ #HIDDEN_BEGIN_acc2
spectrum.add_moments(100)
spectrum.add_vectors(5)
-@@ -109,7 +128,9 @@
+ #HIDDEN_END_acc2
+
+ increased_moments_dos = spectrum()
+
plot_dos([
('density', original_dos),
('higher number of moments', increased_moments_dos),
@@ -96,7 +159,36 @@
def operator_example(fsyst):
-@@ -139,7 +160,9 @@
+ #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)
+
+ _, axes = plt.subplots(1, 2, figsize=(12, 7))
plot_ldos(fsyst, axes,[
('energy = 0', zero_energy_ldos),
('energy = 1', finite_energy_ldos),
@@ -104,6 +196,57 @@
+ ],
+ file_name='kpm_ldos'
+ )
+ #HIDDEN_END_op4
def vector_factory_example(fsyst):
+ spectrum = kwant.kpm.SpectralDensity(fsyst)
+ #HIDDEN_BEGIN_fact1
+ # construct vectors with n random elements -1 or +1.
+ def binary_vectors(n):
+ return 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 main():
+ simple_dos_example()
+
+ fsyst = make_syst().finalized()
+
+ dos_integrating_example(fsyst)
+ increasing_accuracy_example(fsyst)
+ operator_example(fsyst)
+ ldos_example(fsyst)
+ vector_factory_example(fsyst)
+ bilinear_map_operator_example(fsyst)
+
+
+ # 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/images/magnetic_texture.py.diff b/doc/source/images/magnetic_texture.py.diff
index bc7acb4d..0b1f90cb 100644
--- a/doc/source/images/magnetic_texture.py.diff
+++ b/doc/source/images/magnetic_texture.py.diff
@@ -1,6 +1,14 @@
---- original
-+++ modified
-@@ -11,6 +11,9 @@
+@@ -1,245 +1,268 @@
+ # Tutorial 2.7. Spin textures
+ # ===========================
+ #
+ # Physics background
+ # ------------------
+ # - Spin textures
+ # - Skyrmions
+ #
+ # Kwant features highlighted
+ # --------------------------
# - operators
# - plotting vector fields
@@ -10,8 +18,73 @@
from math import sin, cos, tanh, pi
import itertools
import numpy as np
-@@ -79,7 +82,13 @@
+ 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):
@@ -25,7 +98,8 @@
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))
-@@ -88,28 +97,38 @@
+
+ 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)
@@ -52,11 +126,11 @@
+ 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__'):
@@ -72,34 +146,80 @@
def main():
-@@ -119,7 +138,7 @@
+ 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]
-@@ -144,7 +163,7 @@
+ 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)
-@@ -159,7 +178,7 @@
+ 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)
-@@ -173,7 +192,7 @@
+ 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),
@@ -107,8 +227,11 @@
+ ], fname='spin_current_comparison')
+ #HIDDEN_BEGIN_density_cut
def circle(site):
-@@ -183,7 +202,9 @@
+ 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)
@@ -116,10 +239,18 @@
+ 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
-@@ -197,8 +218,10 @@
+
+ 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)
@@ -129,10 +260,21 @@
+ 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)
-@@ -214,7 +237,7 @@
+
+ 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),
@@ -141,3 +283,4 @@
if __name__ == '__main__':
+ main()
diff --git a/doc/source/images/plot_graphene.py.diff b/doc/source/images/plot_graphene.py.diff
index 74263c4f..82c75554 100644
--- a/doc/source/images/plot_graphene.py.diff
+++ b/doc/source/images/plot_graphene.py.diff
@@ -1,6 +1,12 @@
---- original
-+++ modified
-@@ -9,6 +9,7 @@
+@@ -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
@@ -8,7 +14,14 @@
import kwant
from matplotlib import pyplot
-@@ -22,7 +23,7 @@
+ #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()
@@ -17,11 +30,16 @@
syst[lat.neighbors()] = t
syst.eradicate_dangling()
if tp:
-@@ -32,9 +33,11 @@
+ 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
@@ -31,8 +49,10 @@
+ fig_size=size, dpi=_defs.dpi)
# use color and linewidths to get a better plot
+ #HIDDEN_BEGIN_plotsyst2
def family_color(site):
-@@ -43,7 +46,11 @@
+ return 'black' if site.family == a else 'white'
+
def hopping_lw(site1, site2):
return 0.04 if site1.family == site2.family else 0.1
@@ -42,23 +62,38 @@
+ 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):
-@@ -58,7 +65,11 @@
+ 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):
-@@ -68,15 +79,22 @@
+ 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'
@@ -70,8 +105,10 @@
+ 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()
@@ -82,6 +119,23 @@
+ 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/images/plot_qahe.py.diff b/doc/source/images/plot_qahe.py.diff
index 8db900c9..0e095709 100644
--- a/doc/source/images/plot_qahe.py.diff
+++ b/doc/source/images/plot_qahe.py.diff
@@ -1,6 +1,14 @@
---- original
-+++ modified
-@@ -11,6 +11,7 @@
+@@ -1,71 +1,75 @@
+ # 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
@@ -8,8 +16,53 @@
import math
import matplotlib.pyplot
import kwant
-@@ -60,7 +61,10 @@
+
+
+ # 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_spacing=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)
@@ -17,6 +70,8 @@
+ 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/images/plot_zincblende.py.diff b/doc/source/images/plot_zincblende.py.diff
index 5a4bf038..ad500787 100644
--- a/doc/source/images/plot_zincblende.py.diff
+++ b/doc/source/images/plot_zincblende.py.diff
@@ -1,6 +1,12 @@
---- original
-+++ modified
-@@ -9,6 +9,7 @@
+@@ -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
@@ -8,8 +14,30 @@
import kwant
from matplotlib import pyplot
-@@ -33,7 +34,10 @@
+ #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)
@@ -17,10 +45,12 @@
+ 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)
-@@ -41,8 +45,12 @@
+
def family_colors(site):
return 'r' if site.family == a else 'g'
@@ -32,6 +62,10 @@
+ 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/images/quantum_well.py.diff b/doc/source/images/quantum_well.py.diff
index eec8c540..7b81c866 100644
--- a/doc/source/images/quantum_well.py.diff
+++ b/doc/source/images/quantum_well.py.diff
@@ -1,6 +1,12 @@
---- original
-+++ modified
-@@ -9,6 +9,7 @@
+@@ -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
@@ -8,7 +14,51 @@
import kwant
# For plotting
-@@ -55,19 +56,25 @@
+ 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=0):
+ 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, args=[-welldepth])
data.append(smatrix.transmission(1, 0))
@@ -30,6 +80,7 @@
+ 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():
@@ -41,3 +92,12 @@
# 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/images/quantum_wire.py.diff b/doc/source/images/quantum_wire.py.diff
index 0b9aa16d..396ab255 100644
--- a/doc/source/images/quantum_wire.py.diff
+++ b/doc/source/images/quantum_wire.py.diff
@@ -1,29 +1,124 @@
---- original
-+++ modified
-@@ -11,6 +11,7 @@
+@@ -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
-@@ -69,7 +70,10 @@
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()
-@@ -90,8 +94,13 @@
+ #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)
@@ -38,3 +133,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/images/spin_orbit.py.diff b/doc/source/images/spin_orbit.py.diff
index a044351f..ba32e812 100644
--- a/doc/source/images/spin_orbit.py.diff
+++ b/doc/source/images/spin_orbit.py.diff
@@ -1,6 +1,16 @@
---- original
-+++ modified
-@@ -13,6 +13,7 @@
+@@ -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
@@ -8,7 +18,65 @@
import kwant
# For plotting
-@@ -70,19 +71,24 @@
+ 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(a=1, 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(a)
+
+ 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
+ # hoppings in y-directions
+ syst[kwant.builder.HoppingKind((0, 1), lat, lat)] = \
+ -t * sigma_0 + 1j * alpha * sigma_x
+ #HIDDEN_END_uxrm
+
+ #### Define the left lead. ####
+ lead = kwant.Builder(kwant.TranslationalSymmetry((-a, 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
+ # hoppings in y-directions
+ lead[kwant.builder.HoppingKind((0, 1), lat, lat)] = \
+ -t * sigma_0 + 1j * alpha * sigma_x
+ #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))
@@ -40,3 +108,11 @@
# 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/images/superconductor.py.diff b/doc/source/images/superconductor.py.diff
index 06bf2198..156bdac6 100644
--- a/doc/source/images/superconductor.py.diff
+++ b/doc/source/images/superconductor.py.diff
@@ -1,6 +1,14 @@
---- original
-+++ modified
-@@ -11,6 +11,7 @@
+@@ -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.
@@ -8,7 +16,7 @@
import kwant
import tinyarray
-@@ -18,6 +19,7 @@
+ import numpy as np
# For plotting
from matplotlib import pyplot
@@ -16,10 +24,67 @@
tau_x = tinyarray.array([[0, 1], [1, 0]])
tau_y = tinyarray.array([[0, -1j], [1j, 0]])
-@@ -74,11 +76,19 @@
+ 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, phs=True):
+ # 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)
@@ -36,9 +101,28 @@
+ fig.savefig("superconductor_transport_result." + extension,
+ dpi=_defs.dpi)
+ #HIDDEN_BEGIN_pqmp
def check_PHS(syst):
# Scattering matrix
-@@ -103,14 +113,13 @@
+ 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)
@@ -56,3 +140,9 @@
# 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/tutorial/README b/doc/source/tutorial/README
deleted file mode 100644
index 974ae08b..00000000
--- a/doc/source/tutorial/README
+++ /dev/null
@@ -1,11 +0,0 @@
-Note for Kwant developers
--------------------------
-
-Say you have modified SCRIPT.py in this directory. Make sure that the tutorial
-image generation patches still apply by deleting doc/source/images/SCRIPT.py
-and running ``make doc/source/images/SCRIPT.py`` in doc. Now examine the newly
-created doc/source/images/SCRIPT.py. If you do not like the result or the
-patch did not apply, edit doc/source/images/SCRIPT.py until you like it. You
-can run `make html` during your edits to check things. Finally, even if you
-did not edit the script, execute generate-diffs.sh in doc/source/images. If
-the patches applied cleanly the diff files will usually stay unchanged.
diff --git a/doc/source/tutorial/ab_ring.py b/doc/source/tutorial/ab_ring.py
deleted file mode 100644
index 2af92ec9..00000000
--- a/doc/source/tutorial/ab_ring.py
+++ /dev/null
@@ -1,127 +0,0 @@
-# Tutorial 2.3.3. Nontrivial shapes
-# =================================
-#
-# Physics background
-# ------------------
-# Flux-dependent transmission through a quantum ring
-#
-# Kwant features highlighted
-# --------------------------
-# - More complex shapes with lattices
-# - Allows for discussion of subtleties of `attach_lead` (not in the
-# example, but in the tutorial main text)
-# - Modifcations of hoppings/sites after they have been added
-
-from cmath import exp
-from math import pi
-
-import kwant
-
-# For plotting
-from matplotlib import pyplot
-
-
-#HIDDEN_BEGIN_eusz
-def make_system(a=1, t=1.0, W=10, r1=10, r2=20):
- # Start with an empty tight-binding system and a single square lattice.
- # `a` is the lattice constant (by default set to 1 for simplicity).
-
- lat = kwant.lattice.square(a)
-
- syst = kwant.Builder()
-
- #### Define the scattering region. ####
- # Now, we aim for a more complex shape, namely a ring (or annulus)
- def ring(pos):
- (x, y) = pos
- rsq = x ** 2 + y ** 2
- return (r1 ** 2 < rsq < r2 ** 2)
-#HIDDEN_END_eusz
-
- # and add the corresponding lattice points using the `shape`-function
-#HIDDEN_BEGIN_lcak
- syst[lat.shape(ring, (0, r1 + 1))] = 4 * t
- syst[lat.neighbors()] = -t
-#HIDDEN_END_lcak
-
- # In order to introduce a flux through the ring, we introduce a phase on
- # the hoppings on the line cut through one of the arms. Since we want to
- # change the flux without modifying the Builder instance repeatedly, we
- # define the modified hoppings as a function that takes the flux as its
- # parameter phi.
-#HIDDEN_BEGIN_lvkt
- def hopping_phase(site1, site2, phi):
- return -t * exp(1j * phi)
-
- def crosses_branchcut(hop):
- ix0, iy0 = hop[0].tag
-
- # builder.HoppingKind with the argument (1, 0) below
- # returns hoppings ordered as ((i+1, j), (i, j))
- return iy0 < 0 and ix0 == 1 # ix1 == 0 then implied
-
- # Modify only those hopings in x-direction that cross the branch cut
- def hops_across_cut(syst):
- for hop in kwant.builder.HoppingKind((1, 0), lat, lat)(syst):
- if crosses_branchcut(hop):
- yield hop
- syst[hops_across_cut] = hopping_phase
-#HIDDEN_END_lvkt
-
- #### Define the leads. ####
- # left lead
-#HIDDEN_BEGIN_qwgr
- sym_lead = kwant.TranslationalSymmetry((-a, 0))
- lead = kwant.Builder(sym_lead)
-
- def lead_shape(pos):
- (x, y) = pos
- return (-W / 2 < y < W / 2)
-
- lead[lat.shape(lead_shape, (0, 0))] = 4 * t
- lead[lat.neighbors()] = -t
-#HIDDEN_END_qwgr
-
- #### Attach the leads and return the system. ####
-#HIDDEN_BEGIN_skbz
- syst.attach_lead(lead)
- syst.attach_lead(lead.reversed())
-#HIDDEN_END_skbz
-
- return syst
-
-
-def plot_conductance(syst, energy, fluxes):
- # compute conductance
-
- normalized_fluxes = [flux / (2 * pi) for flux in fluxes]
- data = []
- for flux in fluxes:
- smatrix = kwant.smatrix(syst, energy, args=[flux])
- data.append(smatrix.transmission(1, 0))
-
- pyplot.figure()
- pyplot.plot(normalized_fluxes, data)
- pyplot.xlabel("flux [flux quantum]")
- pyplot.ylabel("conductance [e^2/h]")
- pyplot.show()
-
-
-def main():
- syst = make_system()
-
- # Check that the system looks as intended.
- kwant.plot(syst)
-
- # Finalize the system.
- syst = syst.finalized()
-
- # We should see a conductance that is periodic with the flux quantum
- plot_conductance(syst, energy=0.15, fluxes=[0.01 * i * 3 * 2 * pi
- for i in range(100)])
-
-
-# Call the main function if the script gets executed (as opposed to imported).
-# See <http://docs.python.org/library/__main__.html>.
-if __name__ == '__main__':
- main()
diff --git a/doc/source/tutorial/band_structure.py b/doc/source/tutorial/band_structure.py
deleted file mode 100644
index 5b223d3e..00000000
--- a/doc/source/tutorial/band_structure.py
+++ /dev/null
@@ -1,52 +0,0 @@
-# Tutorial 2.4.1. Band structure calculations
-# ===========================================
-#
-# Physics background
-# ------------------
-# band structure of a simple quantum wire in tight-binding approximation
-#
-# Kwant features highlighted
-# --------------------------
-# - Computing the band structure of a finalized lead.
-
-import kwant
-
-# For plotting.
-from matplotlib import pyplot
-
-#HIDDEN_BEGIN_zxip
-def make_lead(a=1, t=1.0, W=10):
- # Start with an empty lead with a single square lattice
- lat = kwant.lattice.square(a)
-
- sym_lead = kwant.TranslationalSymmetry((-a, 0))
- lead = kwant.Builder(sym_lead)
-
- # build up one unit cell of the lead, and add the hoppings
- # to the next unit cell
- for j in range(W):
- lead[lat(0, j)] = 4 * t
-
- if j > 0:
- lead[lat(0, j), lat(0, j - 1)] = -t
-
- lead[lat(1, j), lat(0, j)] = -t
-
- return lead
-#HIDDEN_END_zxip
-
-
-#HIDDEN_BEGIN_pejz
-def main():
- lead = make_lead().finalized()
- kwant.plotter.bands(lead, show=False)
- pyplot.xlabel("momentum [(lattice constant)^-1]")
- pyplot.ylabel("energy [t]")
- pyplot.show()
-#HIDDEN_END_pejz
-
-
-# Call the main function if the script gets executed (as opposed to imported).
-# See <http://docs.python.org/library/__main__.html>.
-if __name__ == '__main__':
- main()
diff --git a/doc/source/tutorial/closed_system.py b/doc/source/tutorial/closed_system.py
deleted file mode 100644
index b192b183..00000000
--- a/doc/source/tutorial/closed_system.py
+++ /dev/null
@@ -1,140 +0,0 @@
-# Tutorial 2.4.2. Closed systems
-# ==============================
-#
-# Physics background
-# ------------------
-# Fock-darwin spectrum of a quantum dot (energy spectrum in
-# as a function of a magnetic field)
-#
-# Kwant features highlighted
-# --------------------------
-# - Use of `hamiltonian_submatrix` in order to obtain a Hamiltonian
-# matrix.
-
-from cmath import exp
-import numpy as np
-import kwant
-
-# For eigenvalue computation
-#HIDDEN_BEGIN_tibv
-import scipy.sparse.linalg as sla
-#HIDDEN_END_tibv
-
-# For plotting
-from matplotlib import pyplot
-
-
-def make_system(a=1, t=1.0, r=10):
- # Start with an empty tight-binding system and a single square lattice.
- # `a` is the lattice constant (by default set to 1 for simplicity).
-
-#HIDDEN_BEGIN_qlyd
- lat = kwant.lattice.square(a, 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=0):
- # The magnetic field is controlled by the parameter B
- y = site1.pos[1]
- return -t * exp(-1j * B * y)
-
- syst[lat.shape(circle, (0, 0))] = 4 * t
- # hoppings in x-direction
- syst[kwant.builder.HoppingKind((1, 0), lat, lat)] = hopx
- # hoppings in y-directions
- syst[kwant.builder.HoppingKind((0, 1), lat, lat)] = -t
-
- # It's a closed system for a change, so no leads
- return syst
-#HIDDEN_END_qlyd
-
-
-#HIDDEN_BEGIN_yvri
-def plot_spectrum(syst, Bfields):
-
- # In the following, we compute the spectrum of the quantum dot
- # using dense matrix methods. This works in this toy example, as
- # the system is tiny. In a real example, one would want to use
- # sparse matrix methods
-
- energies = []
- for B in Bfields:
- # Obtain the Hamiltonian as a dense matrix
- ham_mat = syst.hamiltonian_submatrix(args=[B], sparse=True)
-
- # we only calculate the 15 lowest eigenvalues
- ev = sla.eigsh(ham_mat, k=15, which='SM', return_eigenvectors=False)
-
- energies.append(ev)
-
- pyplot.figure()
- pyplot.plot(Bfields, energies)
- pyplot.xlabel("magnetic field [arbitrary units]")
- pyplot.ylabel("energy [t]")
- pyplot.show()
-#HIDDEN_END_yvri
-
-
-#HIDDEN_BEGIN_wave
-def plot_wave_function(syst):
- # Calculate the wave functions in the system.
- ham_mat = syst.hamiltonian_submatrix(sparse=True)
- evecs = sla.eigsh(ham_mat, k=20, which='SM')[1]
-
- # Plot the probability density of the 10th eigenmode.
- kwant.plotter.map(syst, np.abs(evecs[:, 9])**2,
- colorbar=False, oversampling=1)
-#HIDDEN_END_wave
-
-
-#HIDDEN_BEGIN_current
-def plot_current(syst):
- # Calculate the wave functions in the system.
- ham_mat = syst.hamiltonian_submatrix(sparse=True)
- evecs = sla.eigsh(ham_mat, k=20, which='SM')[1]
-
- # Calculate and plot the local current of the 10th eigenmode.
- J = kwant.operator.Current(syst)
- current = J(evecs[:, 9])
- kwant.plotter.current(syst, current, colorbar=False)
-#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/tutorial/discretize.py b/doc/source/tutorial/discretize.py
deleted file mode 100644
index 03c04db3..00000000
--- a/doc/source/tutorial/discretize.py
+++ /dev/null
@@ -1,222 +0,0 @@
-# Tutorial 2.9. Processing continuum Hamiltonians with discretize
-# ===============================================================
-#
-# Physics background
-# ------------------
-# - tight-binding approximation of continuous Hamiltonians
-#
-# Kwant features highlighted
-# --------------------------
-# - kwant.continuum.discretize
-
-
-import kwant
-import scipy.sparse.linalg
-import scipy.linalg
-import numpy as np
-
-# For plotting
-import matplotlib as mpl
-from matplotlib import pyplot as plt
-
-
-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)
-#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()
-
-
-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_spacing=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()
-#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()
-
-
-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_spacing=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()
-
-
-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])
-#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))
-#HIDDEN_END_subs_2
-
-
-def main():
-#HIDDEN_BEGIN_symbolic_discretization
- template = kwant.continuum.discretize('k_x * A(x) * k_x')
- print(template)
-#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/tutorial/faq.py b/doc/source/tutorial/faq.py
deleted file mode 100644
index 3a00860b..00000000
--- a/doc/source/tutorial/faq.py
+++ /dev/null
@@ -1,450 +0,0 @@
-# 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 kwant
-import numpy as np
-import tinyarray
-import matplotlib
-from matplotlib import pyplot as plt
-matplotlib.rcParams['figure.figsize'] = (3.5, 3.5)
-
-
-#### 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)
-#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)
-
-
-#### 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)
-#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)
-
-# 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)
-#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)
-#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)
-#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)
-# 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)
-
-
-#### 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)
-
-del syst[lat.neighbors()] # Delete all nearest-neighbor hoppings
-syst[lat.neighbors(2)] = -1
-
-kwant.plot(syst)
-#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)
-del syst[lat.neighbors()] # Delete the hoppings previously created
-#HIDDEN_BEGIN_adjacent3
-syst[a.neighbors()] = -1
-#HIDDEN_END_adjacent3
-plot_system(syst)
-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)
-
-#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)
-
-#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)
-
-#HIDDEN_BEGIN_different_lattice4
-syst.attach_lead(lead1)
-#HIDDEN_END_different_lattice4
-plot_system(syst)
-
-
-#### 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);
-
-#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);
-
-
-#### 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)
-
-# 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)
-
-
-#### 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])
-#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])
-#HIDDEN_END_ord2
diff --git a/doc/source/tutorial/graphene.py b/doc/source/tutorial/graphene.py
deleted file mode 100644
index 512ca7e8..00000000
--- a/doc/source/tutorial/graphene.py
+++ /dev/null
@@ -1,179 +0,0 @@
-# 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
-
-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()
- 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()
-
-
-#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)
-#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)
-
- # 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/tutorial/kernel_polynomial_method.py b/doc/source/tutorial/kernel_polynomial_method.py
deleted file mode 100644
index 226893b4..00000000
--- a/doc/source/tutorial/kernel_polynomial_method.py
+++ /dev/null
@@ -1,219 +0,0 @@
-# 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
-
-# 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
-
-
-# Plot several density of states curves on the same axes.
-def plot_dos(labels_to_data):
- 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("DoS [a.u.]")
- plt.show()
- 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(fsyst, axes, titles_to_data, file_name=None):
- for ax, (title, ldos) in zip(axes, titles_to_data):
- site_size = site_size_conversion(ldos) # convert LDoS to sizes
- kwant.plot(fsyst, site_size=site_size, site_color=(0, 0, 1, 0.3), ax=ax)
- ax.set_title(title)
- ax.set(adjustable='box-forced', aspect='equal')
- plt.show()
- plt.clf()
-
-
-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)),
- ])
-
-
-def dos_integrating_example(fsyst):
- spectrum = kwant.kpm.SpectralDensity(fsyst)
-
-#HIDDEN_BEGIN_int1
- 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))
-#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),
- ])
-
-#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),
- ])
-
-
-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)
-
- _, axes = plt.subplots(1, 2, figsize=(12, 7))
- plot_ldos(fsyst, axes,[
- ('energy = 0', zero_energy_ldos),
- ('energy = 1', finite_energy_ldos),
- ])
-#HIDDEN_END_op4
-
-
-def vector_factory_example(fsyst):
- spectrum = kwant.kpm.SpectralDensity(fsyst)
-#HIDDEN_BEGIN_fact1
- # construct vectors with n random elements -1 or +1.
- def binary_vectors(n):
- return 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 main():
- simple_dos_example()
-
- fsyst = make_syst().finalized()
-
- dos_integrating_example(fsyst)
- increasing_accuracy_example(fsyst)
- operator_example(fsyst)
- ldos_example(fsyst)
- vector_factory_example(fsyst)
- bilinear_map_operator_example(fsyst)
-
-
-# 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/tutorial/magnetic_texture.py b/doc/source/tutorial/magnetic_texture.py
deleted file mode 100644
index 89c7acfb..00000000
--- a/doc/source/tutorial/magnetic_texture.py
+++ /dev/null
@@ -1,245 +0,0 @@
-# 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)
-
-#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):
- 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))
- 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))
- 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()
-
-
-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))
-
-#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),
- ])
-
-#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),
- ])
-
-#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),
- ])
-
-
-#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)
-#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))
-#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),
- ])
-
-
-if __name__ == '__main__':
- main()
diff --git a/doc/source/tutorial/plot_graphene.py b/doc/source/tutorial/plot_graphene.py
deleted file mode 100644
index e81e331b..00000000
--- a/doc/source/tutorial/plot_graphene.py
+++ /dev/null
@@ -1,112 +0,0 @@
-# Tutorial 2.8.1. 2D example: graphene quantum dot
-# ================================================
-#
-# Physics background
-# ------------------
-# - graphene edge states
-#
-# Kwant features highlighted
-# --------------------------
-# - demonstrate different ways of plotting
-
-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.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
-
- # 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)
-#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')
-#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')
-#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)
-#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/tutorial/plot_qahe.py b/doc/source/tutorial/plot_qahe.py
deleted file mode 100644
index 65b2ee65..00000000
--- a/doc/source/tutorial/plot_qahe.py
+++ /dev/null
@@ -1,71 +0,0 @@
-# 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
-
-
-# 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_spacing=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)
-#HIDDEN_END_current
-
-
-if __name__ == '__main__':
- main()
diff --git a/doc/source/tutorial/plot_zincblende.py b/doc/source/tutorial/plot_zincblende.py
deleted file mode 100644
index 6b75a54e..00000000
--- a/doc/source/tutorial/plot_zincblende.py
+++ /dev/null
@@ -1,59 +0,0 @@
-# Tutorial 2.8.2. 3D example: zincblende structure
-# ================================================
-#
-# Physical background
-# -------------------
-# - 3D Bravais lattices
-#
-# Kwant features highlighted
-# --------------------------
-# - demonstrate different ways of plotting in 3D
-
-import kwant
-from matplotlib import pyplot
-
-#HIDDEN_BEGIN_zincblende1
-lat = kwant.lattice.general([(0, 0.5, 0.5), (0.5, 0, 0.5), (0.5, 0.5, 0)],
- [(0, 0, 0), (0.25, 0.25, 0.25)])
-a, b = lat.sublattices
-#HIDDEN_END_zincblende1
-
-#HIDDEN_BEGIN_zincblende2
-def make_cuboid(a=15, b=10, c=5):
- def cuboid_shape(pos):
- x, y, z = pos
- return 0 <= x < a and 0 <= y < b and 0 <= z < c
-
- syst = kwant.Builder()
- syst[lat.shape(cuboid_shape, (0, 0, 0))] = None
- syst[lat.neighbors()] = None
-
- return syst
-#HIDDEN_END_zincblende2
-
-
-def main():
- # the standard plotting style for 3D is mainly useful for
- # checking shapes:
-#HIDDEN_BEGIN_plot1
- syst = make_cuboid()
-
- kwant.plot(syst)
-#HIDDEN_END_plot1
-
- # visualize the crystal structure better for a very small system
-#HIDDEN_BEGIN_plot2
- syst = make_cuboid(a=1.5, b=1.5, c=1.5)
-
- def family_colors(site):
- return 'r' if site.family == a else 'g'
-
- kwant.plot(syst, site_size=0.18, site_lw=0.01, hop_lw=0.05,
- site_color=family_colors)
-#HIDDEN_END_plot2
-
-
-# Call the main function if the script gets executed (as opposed to imported).
-# See <http://docs.python.org/library/__main__.html>.
-if __name__ == '__main__':
- main()
diff --git a/doc/source/tutorial/quantum_well.py b/doc/source/tutorial/quantum_well.py
deleted file mode 100644
index 780052e9..00000000
--- a/doc/source/tutorial/quantum_well.py
+++ /dev/null
@@ -1,88 +0,0 @@
-# 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
-
-
-#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=0):
- 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, args=[-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()
-#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/tutorial/quantum_wire.py b/doc/source/tutorial/quantum_wire.py
deleted file mode 100644
index 992a7dc3..00000000
--- a/doc/source/tutorial/quantum_wire.py
+++ /dev/null
@@ -1,121 +0,0 @@
-# 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
-
-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)
-#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()
-pyplot.plot(energies, data)
-pyplot.xlabel("energy [t]")
-pyplot.ylabel("conductance [e^2/h]")
-pyplot.show()
-#HIDDEN_END_lliv
diff --git a/doc/source/tutorial/spin_orbit.py b/doc/source/tutorial/spin_orbit.py
deleted file mode 100644
index f4c84b6c..00000000
--- a/doc/source/tutorial/spin_orbit.py
+++ /dev/null
@@ -1,104 +0,0 @@
-# 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
-
-# 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(a=1, 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(a)
-
- 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
- # hoppings in y-directions
- syst[kwant.builder.HoppingKind((0, 1), lat, lat)] = \
- -t * sigma_0 + 1j * alpha * sigma_x
-#HIDDEN_END_uxrm
-
- #### Define the left lead. ####
- lead = kwant.Builder(kwant.TranslationalSymmetry((-a, 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
- # hoppings in y-directions
- lead[kwant.builder.HoppingKind((0, 1), lat, lat)] = \
- -t * sigma_0 + 1j * alpha * sigma_x
-#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()
- pyplot.plot(energies, data)
- pyplot.xlabel("energy [t]")
- pyplot.ylabel("conductance [e^2/h]")
- pyplot.show()
-
-
-def main():
- syst = make_system()
-
- # Check that the system looks as intended.
- kwant.plot(syst)
-
- # Finalize the system.
- syst = syst.finalized()
-
- # We should see 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/tutorial/superconductor.py b/doc/source/tutorial/superconductor.py
deleted file mode 100644
index 06caf0c4..00000000
--- a/doc/source/tutorial/superconductor.py
+++ /dev/null
@@ -1,132 +0,0 @@
-# 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
-
-import tinyarray
-import numpy as np
-
-# For plotting
-from matplotlib import pyplot
-
-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, phs=True):
- # 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()
- pyplot.plot(energies, data)
- pyplot.xlabel("energy [t]")
- pyplot.ylabel("conductance [e^2/h]")
- pyplot.show()
-
-#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)
-
- # 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()
--
GitLab