diff --git a/.gitignore b/.gitignore
index 64f23822bf30e9700f24c185f626343d379dd239..7357693774aea277987ca67c03ecd729581debb1 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 8fe4686dd04f742a21594ec6f885676af474769b..1b0a494f0927f1f781fb7bfa653caf02e24fb036 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 985a4ef95df64706afd9f2344c73f41fc4d99b88..8df4d59b00d3f81b876e8a3249cca7ddb5eab9c8 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 ded6e1426f8aeba1694994b425b39ac9454bb541..50a10d1920c488552360b3bfc388d661da86d88f 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 90d735445d34d33297a7a5532ff9b867892cb34b..8a2652072aea317f4a6e5090497c37b43afdef15 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 e33c3c2c86340a909b0959b42fdece6138e86100..7f1de4715780c3c04f7ccf668301027617fff545 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 7ac6dfe381a9289e7d0c1c7766222e0664e56a81..3d6c20249528f78846a946f1cd53d926f1db3c7c 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 4df4fd94bffa976cf028011428df2d4d1c19468d..7c8889dceec701572fde707e7c5aed3d5cf9f427 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 d78525ea67151e22aad3ef50ccbf20b0bf5f73e0..0000000000000000000000000000000000000000
--- 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 4871213278c445098a82cc72bd2c2cacc5cf8480..1352ed270df555a9cd1aa5c6ef51632590cd7337 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 d58b3836bd85bb6ec8fe06f930dfc41dedf6788a..9e778f10d77979805d3587eedf23d4447e49807c 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 bc7acb4dcf1b7bda417c3dc3458eb61655435e05..0b1f90cb6bd17d88a97b7e1924641fb640b15a63 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 74263c4f1ce1f523ba6c4205b4a5f99c8e92de41..82c7555424f6a0206c68cbce0235da3e6fa7ef73 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 8db900c966ddc680fd2bf0b5d0282830cc3b9695..0e0957092919800658a8ffae081cf47491ac4861 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 5a4bf03839c56c4732dd349e2fa628b76a99d02f..ad5007871331dd28bf8653ab5be63e102d58205b 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 eec8c54023f956381a24dc8cdbcd4795c7f7ac99..7b81c8668eadb48b0dbcceff7cfb0a04a7d78d64 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 0b9aa16dfb1d1ed1765e516369899c862fb08578..396ab255359afa8c53f7790a1c53b7b9fe2de178 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 a044351f55751cd6b833b4ad6851ebdbbd2cb380..ba32e81221f498d1b7057502b566632fbdf202df 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 06bf2198fb7a378d02170d9eff225c6bfcf4e0c5..156bdac66a114f713af8edc5d65d68da572df0bc 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 974ae08b427a7f8e2974b88f6ff869a931aea40c..0000000000000000000000000000000000000000
--- 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 2af92ec9c659812649911ea4586b615525cede85..0000000000000000000000000000000000000000
--- 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 5b223d3e96f0f474e7b684f2eae779c160698388..0000000000000000000000000000000000000000
--- 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 b192b18374981d45a3e98f23eeeb179fbe2c9295..0000000000000000000000000000000000000000
--- 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 03c04db3ec3da9c1a4782ce8dac7c214179d1cdc..0000000000000000000000000000000000000000
--- 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 3a00860b375fbdc558e141bfeae803c4fd19ccc7..0000000000000000000000000000000000000000
--- 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 512ca7e8d68f2554299d48548769cb017427c1fd..0000000000000000000000000000000000000000
--- 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 226893b45cd2f142992bbf262db7d172643628cb..0000000000000000000000000000000000000000
--- 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 89c7acfb793aecce1126a50fa45f31b360c2e053..0000000000000000000000000000000000000000
--- 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 e81e331b2bedf7de7241b857a50611a562d51316..0000000000000000000000000000000000000000
--- 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 65b2ee65115e9a7224886af3d632912eba972f6f..0000000000000000000000000000000000000000
--- 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 6b75a54e0df3d4b977ab1ec51a00960ef4f6f262..0000000000000000000000000000000000000000
--- 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 780052e9be77dbd6d2a52a9b6742a495d6481232..0000000000000000000000000000000000000000
--- 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 992a7dc37fb17309da2f9c2f51d2c68bcc90371d..0000000000000000000000000000000000000000
--- 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 f4c84b6c8984062aeb155f25f258ddd0345434e1..0000000000000000000000000000000000000000
--- 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 06caf0c4049280ceed1cc71aa7384730ff45e145..0000000000000000000000000000000000000000
--- 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()