doc/source/code/figure/faq.py.diff

-@@ -1,450 +1,479 @@
- # Frequently asked questions
- #
- # This script is a disorganized collection of code snippets.  As a whole, it is
- # not meant as an example of good programming practice.
- 
-+import _defs
- import kwant
- import numpy as np
- import tinyarray
- import matplotlib
- from matplotlib import pyplot as plt
- matplotlib.rcParams['figure.figsize'] = (3.5, 3.5)
- 
- 
-+def save_figure(file_name):
-+    for extension in ('pdf', 'png'):
-+        plt.savefig('.'.join((file_name, extension)),
-+                    dpi=_defs.dpi, bbox_inches='tight')
-+
-+
- #### What is a Site?
- 
- #HIDDEN_BEGIN_site
- a = 1
- lat = kwant.lattice.square(a)
- syst = kwant.Builder()
- 
- syst[lat(1, 0)] = 4
- syst[lat(1, 1)] = 4
- 
- kwant.plot(syst)
-+save_figure("faq_site")
- #HIDDEN_END_site
- 
- 
- #### What is a hopping?
- 
- a = 1
- lat = kwant.lattice.square(a)
- syst = kwant.Builder()
- 
- syst[lat(1, 0)] = 4
- syst[lat(1, 1)] = 4
- #HIDDEN_BEGIN_hopping
- syst[(lat(1, 0), lat(1, 1))] = 1j
- #HIDDEN_END_hopping
- 
- kwant.plot(syst)
-+save_figure("faq_hopping")
- 
- 
- #### What is a lattice?
- 
- #HIDDEN_BEGIN_lattice_monatomic
- # Two monatomic lattices
- primitive_vectors = [(1, 0), (0, 1)]
- lat_a = kwant.lattice.Monatomic(primitive_vectors, offset=(0, 0))
- lat_b = kwant.lattice.Monatomic(primitive_vectors, offset=(0.5, 0.5))
- # lat1 is equivalent to kwant.lattice.square()
- 
- syst = kwant.Builder()
- 
- syst[lat_a(0, 0)] = 4
- syst[lat_b(0, 0)] = 4
- 
- kwant.plot(syst)
-+save_figure("faq_lattice")
- #HIDDEN_END_lattice_monatomic
- 
- #HIDDEN_BEGIN_lattice_polyatomic
- # One polyatomic lattice containing two sublattices
- lat = kwant.lattice.Polyatomic([(1, 0), (0, 1)], [(0, 0), (0.5, 0.5)])
- sub_a, sub_b = lat.sublattices
- #HIDDEN_END_lattice_polyatomic
- 
- #### How to make a hole in a system?
- 
- #HIDDEN_BEGIN_hole
- # Define the lattice and the (empty) system
- a = 2
- lat = kwant.lattice.cubic(a)
- syst = kwant.Builder()
- 
- L = 10
- W = 10
- H = 2
- 
- # Add sites to the system in a cuboid
- 
- syst[(lat(i, j, k) for i in range(L) for j in range(W) for k in range(H))] = 4
- kwant.plot(syst)
-+save_figure("faq_hole1")
- 
- # Delete sites to create a hole
- 
- def in_hole(site):
-     x, y, z = site.pos / a - (L/2, W/2, H/2)  # position relative to centre
-     return abs(x) < L / 4 and abs(y) < W / 4
- 
- for site in filter(in_hole, list(syst.sites())):
-     del syst[site]
- 
- kwant.plot(syst)
-+save_figure("faq_hole2")
- #HIDDEN_END_hole
- 
- 
- #### How can we get access to the sites of our system?
- 
- builder = kwant.Builder()
- lat = kwant.lattice.square()
- builder[(lat(i, j) for i in range(3) for j in range(3))] = 4
- #HIDDEN_BEGIN_sites1
- # Before finalizing the system
- 
- sites = list(builder.sites())  # sites() doe *not* return a list
- #HIDDEN_END_sites1
- #HIDDEN_BEGIN_sites2
- # After finalizing the system
- syst = builder.finalized()
- sites = syst.sites  # syst.sites is an actual list
- #HIDDEN_END_sites2
- #HIDDEN_BEGIN_sites3
- i = syst.id_by_site[lat(0, 2)]  # we want the id of the site lat(0, 2)
- #HIDDEN_END_sites3
- 
- 
- #### How to plot a polyatomic lattice with different colors?
- 
- #HIDDEN_BEGIN_colors1
- lat = kwant.lattice.kagome()
- syst = kwant.Builder()
- 
- a, b, c = lat.sublattices  # The kagome lattice has 3 sublattices
- #HIDDEN_END_colors1
- 
- #HIDDEN_BEGIN_colors2
- # Plot sites from different families in different colors
- def family_color(site):
-     if site.family == a:
-         return 'red'
-     if site.family == b:
-         return 'green'
-     else:
-         return 'blue'
- 
- def plot_system(syst):
-     kwant.plot(syst, site_lw=0.1, site_color=family_color)
- 
- ## Add sites and hoppings.
- for i in range(4):
-     for j in range (4):
-         syst[a(i, j)] = 4
-         syst[b(i, j)] = 4
-         syst[c(i, j)] = 4
- 
- syst[lat.neighbors()] = -1
- 
- ## Plot the system.
- plot_system(syst)
-+save_figure("faq_colors")
- #HIDDEN_END_colors2
- 
--
- #### How to create all hoppings in a given direction using Hoppingkind?
- 
- # Monatomic lattice
- 
- #HIDDEN_BEGIN_direction1
- # Create hopping between neighbors with HoppingKind
- a = 1
- syst = kwant.Builder()
- lat = kwant.lattice.square(a)
- syst[ (lat(i, j) for i in range(5) for j in range(5)) ] = 4
- 
- syst[kwant.builder.HoppingKind((1, 0), lat)] = -1
- kwant.plot(syst)
-+save_figure("faq_direction1")
- #HIDDEN_END_direction1
- 
- # Polyatomic lattice
- 
- lat = kwant.lattice.kagome()
- syst = kwant.Builder()
- 
- a, b, c = lat.sublattices  # The kagome lattice has 3 sublattices
- 
- def family_color(site):
-     if site.family == a:
-         return 'red'
-     if site.family == b:
-         return 'green'
-     else:
-         return 'blue'
- 
- def plot_system(syst):
-     kwant.plot(syst, site_size=0.15, site_lw=0.05, site_color=family_color)
- 
- 
- for i in range(4):
-     for j in range (4):
-         syst[a(i, j)] = 4
-         syst[b(i, j)] = 4
-         syst[c(i, j)] = 4
- 
- 
- #HIDDEN_BEGIN_direction2
- # equivalent to syst[kwant.builder.HoppingKind((0, 1), b)] = -1
- syst[kwant.builder.HoppingKind((0, 1), b, b)] = -1
- #HIDDEN_END_direction2
- plot_system(syst)
-+save_figure("faq_direction2")
- # Delete the hoppings previously created
- del syst[kwant.builder.HoppingKind((0, 1), b, b)]
- #HIDDEN_BEGIN_direction3
- syst[kwant.builder.HoppingKind((0, 0), a, b)] = -1
- syst[kwant.builder.HoppingKind((0, 0), a, c)] = -1
- syst[kwant.builder.HoppingKind((0, 0), c, b)] = -1
- #HIDDEN_END_direction3
- plot_system(syst)
-+save_figure("faq_direction3")
- 
- 
- #### How to create the hoppings between adjacent sites?
- 
- # Monatomic lattice
- 
- #HIDDEN_BEGIN_adjacent1
- # Create hoppings with lat.neighbors()
- syst = kwant.Builder()
- lat = kwant.lattice.square()
- syst[(lat(i, j) for i in range(3) for j in range(3))] = 4
- 
- syst[lat.neighbors()] = -1  # Equivalent to lat.neighbors(1)
- kwant.plot(syst)
-+save_figure("faq_adjacent1")
- 
- del syst[lat.neighbors()]  # Delete all nearest-neighbor hoppings
- syst[lat.neighbors(2)] = -1
- 
- kwant.plot(syst)
-+save_figure("faq_adjacent2")
- #HIDDEN_END_adjacent1
- 
- # Polyatomic lattice
- 
- #HIDDEN_BEGIN_FAQ6
- 
- # Hoppings using .neighbors()
- #HIDDEN_BEGIN_adjacent2
- # Create the system
- lat = kwant.lattice.kagome()
- syst = kwant.Builder()
- a, b, c = lat.sublattices  # The kagome lattice has 3 sublattices
- 
- for i in range(4):
-     for j in range (4):
-         syst[a(i, j)] = 4  # red
-         syst[b(i, j)] = 4  # green
-         syst[c(i, j)] = 4  # blue
- 
- syst[lat.neighbors()] = -1
- #HIDDEN_END_adjacent2
- plot_system(syst)
-+save_figure("faq_adjacent3")
- del syst[lat.neighbors()]  # Delete the hoppings previously created
- #HIDDEN_BEGIN_adjacent3
- syst[a.neighbors()] = -1
- #HIDDEN_END_adjacent3
- plot_system(syst)
-+save_figure("faq_adjacent4")
- del syst[a.neighbors()]  # Delete the hoppings previously created
- 
- 
- syst[lat.neighbors(2)] = -1
- plot_system(syst)
- del syst[lat.neighbors(2)]
- 
- 
- #### How to create a lead with a lattice different from the scattering region?
- 
- # Plot sites from different families in different colors
- 
- def plot_system(syst):
- 
-     def family_color(site):
-         if site.family == subA:
-             return 'blue'
-         if site.family == subB:
-             return 'yellow'
-         else:
-             return 'green'
- 
-     kwant.plot(syst, site_lw=0.1, site_color=family_color)
- 
- 
- #HIDDEN_BEGIN_different_lattice1
- # Define the scattering Region
- L = 5
- W = 5
- 
- lat = kwant.lattice.honeycomb()
- subA, subB = lat.sublattices
- 
- syst = kwant.Builder()
- syst[(subA(i, j) for i in range(L) for j in range(W))] = 4
- syst[(subB(i, j) for i in range(L) for j in range(W))] = 4
- syst[lat.neighbors()] = -1
- #HIDDEN_END_different_lattice1
- plot_system(syst)
-+save_figure("faq_different_lattice1")
- 
- #HIDDEN_BEGIN_different_lattice2
- # Create a lead
- lat_lead = kwant.lattice.square()
- sym_lead1 = kwant.TranslationalSymmetry((0, 1))
- 
- lead1 = kwant.Builder(sym_lead1)
- lead1[(lat_lead(i, 0) for i in range(2, 7))] = 4
- lead1[lat_lead.neighbors()] = -1
- #HIDDEN_END_different_lattice2
- plot_system(lead1)
-+save_figure("faq_different_lattice2")
- 
- #HIDDEN_BEGIN_different_lattice3
- syst[(lat_lead(i, 5) for i in range(2, 7))] = 4
- syst[lat_lead.neighbors()] = -1
- 
- # Manually attach sites from graphene to square lattice
- syst[((lat_lead(i+2, 5), subB(i, 4)) for i in range(5))] = -1
- #HIDDEN_END_different_lattice3
- plot_system(syst)
-+save_figure("faq_different_lattice3")
- 
- #HIDDEN_BEGIN_different_lattice4
- syst.attach_lead(lead1)
- #HIDDEN_END_different_lattice4
- plot_system(syst)
-+save_figure("faq_different_lattice4")
- 
- 
- #### How to cut a finite system out of a system with translationnal symmetries?
- 
- #HIDDEN_BEGIN_fill1
- # Create 3d model.
- cubic = kwant.lattice.cubic()
- sym_3d = kwant.TranslationalSymmetry([1, 0, 0], [0, 1, 0], [0, 0, 1])
- model = kwant.Builder(sym_3d)
- model[cubic(0, 0, 0)] = 4
- model[cubic.neighbors()] = -1
- #HIDDEN_END_fill1
- 
- #HIDDEN_BEGIN_fill2
- # Build scattering region (white).
- def cuboid_shape(site):
-     x, y, z = abs(site.pos)
-     return x < 4 and y < 10 and z < 3
- 
- cuboid = kwant.Builder()
- cuboid.fill(model, cuboid_shape, (0, 0, 0));
- #HIDDEN_END_fill2
- kwant.plot(cuboid);
-+save_figure("faq_fill2")
- 
- #HIDDEN_BEGIN_fill3
- # Build electrode (black).
- def electrode_shape(site):
-     x, y, z = site.pos - (0, 5, 2)
-     return y**2 + z**2 < 2.3**2
- 
- electrode = kwant.Builder(kwant.TranslationalSymmetry([1, 0, 0]))
- electrode.fill(model, electrode_shape, (0, 5, 2))  # lead
- 
- # Scattering region
- cuboid.fill(electrode, lambda s: abs(s.pos[0]) < 7, (0, 5, 4))
- 
- cuboid.attach_lead(electrode)
- #HIDDEN_END_fill3
- kwant.plot(cuboid);
-+save_figure("faq_fill3")
- 
- 
- #### How does Kwant order the propagating modes of a lead?
- 
- #HIDDEN_BEGIN_pm
- lat = kwant.lattice.square()
- 
- lead = kwant.Builder(kwant.TranslationalSymmetry((-1, 0)))
- lead[(lat(0, i) for i in range(3))] = 4
- lead[lat.neighbors()] = -1
- 
- flead = lead.finalized()
- 
- E = 2.5
- prop_modes, _ = flead.modes(energy=E)
- #HIDDEN_END_pm
- 
- def plot_and_label_modes(lead, E):
-     # Plot the different modes
-     pmodes, _ = lead.modes(energy=E)
-     kwant.plotter.bands(lead, show=False)
-     for i, k in enumerate(pmodes.momenta):
-         plt.plot(k, E, 'ko')
-         plt.annotate(str(i), xy=(k, E), xytext=(-5, 8),
-                      textcoords='offset points',
-                      bbox=dict(boxstyle='round,pad=0.1',fc='white', alpha=0.7))
-     plt.plot([-3, 3], [E, E], 'r--')
-     plt.ylim(E-1, E+1)
-     plt.xlim(-2, 2)
-     plt.xlabel("momentum")
-     plt.ylabel("energy")
-     plt.show()
- 
- plot_and_label_modes(flead, E)
-+save_figure('faq_pm1')
- 
- # More involved example
- 
- s0 = np.eye(2)
- sz = np.array([[1, 0], [0, -1]])
- 
- lead2 = kwant.Builder(kwant.TranslationalSymmetry((-1, 0)))
- 
- lead2[(lat(0, i) for i in range(2))] = np.diag([1.8, -1])
- lead2[lat.neighbors()] = -1 * sz
- 
- flead2 = lead2.finalized()
- 
- plot_and_label_modes(flead2, 1)
-+save_figure('faq_pm2')
- 
- 
- #### How does Kwant order components of an individual wavefunction?
- def circle(R):
-     return lambda r: np.linalg.norm(r) < R
- 
- 
- def make_system(lat):
-     norbs = lat.norbs
-     syst = kwant.Builder()
-     syst[lat.shape(circle(3), (0, 0))] = 4 * np.eye(norbs)
-     syst[lat.neighbors()] = -1 * np.eye(norbs)
- 
-     lead = kwant.Builder(kwant.TranslationalSymmetry((-1, 0)))
-     lead[(lat(0, i) for i in range(-1, 2))] = 4 * np.eye(norbs)
-     lead[lat.neighbors()] = -1 * np.eye(norbs)
- 
-     syst.attach_lead(lead)
-     syst.attach_lead(lead.reversed())
- 
-     return syst.finalized()
- 
- 
- #HIDDEN_BEGIN_ord1
- lat = kwant.lattice.square(norbs=1)
- syst = make_system(lat)
- scattering_states = kwant.wave_function(syst, energy=1)
- wf = scattering_states(0)[0]  # scattering state from lead 0 incoming in mode 0
- 
- idx = syst.id_by_site[lat(0, 0)]  # look up index of site
- 
--print('wavefunction on lat(0, 0): ', wf[idx])
-+with open('faq_ord1.txt', 'w') as f:
-+    print('wavefunction on lat(0, 0): ', wf[idx], file=f)
- #HIDDEN_END_ord1
- 
- #HIDDEN_BEGIN_ord2
- lat = kwant.lattice.square(norbs=2)
- syst = make_system(lat)
- scattering_states = kwant.wave_function(syst, energy=1)
- wf = scattering_states(0)[0]  # scattering state from lead 0 incoming in mode 0
- 
- idx = syst.id_by_site[lat(0, 0)]  # look up index of site
- 
- # Group consecutive degrees of freedom from 'wf' together; these correspond
- # to degrees of freedom on the same site.
- wf = wf.reshape(-1, 2)
- 
--print('wavefunction on lat(0, 0): ', wf[idx])
-+with open('faq_ord2.txt', 'w') as f:
-+    print('wavefunction on lat(0, 0): ', wf[idx], file=f)
- #HIDDEN_END_ord2

doc/source/code/figure/graphene.py.diff

-@@ -1,179 +1,203 @@
- # Tutorial 2.5. Beyond square lattices: graphene
- # ==============================================
- #
- # Physics background
- # ------------------
- #  Transport through a graphene quantum dot with a pn-junction
- #
- # Kwant features highlighted
- # --------------------------
- #  - Application of all the aspects of tutorials 1-3 to a more complicated
- #    lattice, namely graphene
- 
-+import _defs
- from math import pi, sqrt, tanh
- 
- import kwant
- 
- # For computing eigenvalues
- import scipy.sparse.linalg as sla
- 
- # For plotting
- from matplotlib import pyplot
- 
- 
- # Define the graphene lattice
- sin_30, cos_30 = (1 / 2, sqrt(3) / 2)
- #HIDDEN_BEGIN_hnla
- graphene = kwant.lattice.general([(1, 0), (sin_30, cos_30)],
-                                  [(0, 0), (0, 1 / sqrt(3))])
- a, b = graphene.sublattices
- #HIDDEN_END_hnla
- 
- 
- #HIDDEN_BEGIN_shzy
- def make_system(r=10, w=2.0, pot=0.1):
- 
-     #### Define the scattering region. ####
-     # circular scattering region
-     def circle(pos):
-         x, y = pos
-         return x ** 2 + y ** 2 < r ** 2
- 
-     syst = kwant.Builder()
- 
-     # w: width and pot: potential maximum of the p-n junction
-     def potential(site):
-         (x, y) = site.pos
-         d = y * cos_30 + x * sin_30
-         return pot * tanh(d / w)
- 
-     syst[graphene.shape(circle, (0, 0))] = potential
- #HIDDEN_END_shzy
- 
-     # specify the hoppings of the graphene lattice in the
-     # format expected by builder.HoppingKind
- #HIDDEN_BEGIN_hsmc
-     hoppings = (((0, 0), a, b), ((0, 1), a, b), ((-1, 1), a, b))
- #HIDDEN_END_hsmc
- #HIDDEN_BEGIN_bfwb
-     syst[[kwant.builder.HoppingKind(*hopping) for hopping in hoppings]] = -1
- #HIDDEN_END_bfwb
- 
-     # Modify the scattering region
- #HIDDEN_BEGIN_efut
-     del syst[a(0, 0)]
-     syst[a(-2, 1), b(2, 2)] = -1
- #HIDDEN_END_efut
- 
-     #### Define the leads. ####
- #HIDDEN_BEGIN_aakh
-     # left lead
-     sym0 = kwant.TranslationalSymmetry(graphene.vec((-1, 0)))
- 
-     def lead0_shape(pos):
-         x, y = pos
-         return (-0.4 * r < y < 0.4 * r)
- 
-     lead0 = kwant.Builder(sym0)
-     lead0[graphene.shape(lead0_shape, (0, 0))] = -pot
-     lead0[[kwant.builder.HoppingKind(*hopping) for hopping in hoppings]] = -1
- 
-     # The second lead, going to the top right
-     sym1 = kwant.TranslationalSymmetry(graphene.vec((0, 1)))
- 
-     def lead1_shape(pos):
-         v = pos[1] * sin_30 - pos[0] * cos_30
-         return (-0.4 * r < v < 0.4 * r)
- 
-     lead1 = kwant.Builder(sym1)
-     lead1[graphene.shape(lead1_shape, (0, 0))] = pot
-     lead1[[kwant.builder.HoppingKind(*hopping) for hopping in hoppings]] = -1
- #HIDDEN_END_aakh
- 
- #HIDDEN_BEGIN_kmmw
-     return syst, [lead0, lead1]
- #HIDDEN_END_kmmw
- 
- 
- #HIDDEN_BEGIN_zydk
- def compute_evs(syst):
-     # Compute some eigenvalues of the closed system
-     sparse_mat = syst.hamiltonian_submatrix(sparse=True)
- 
-     evs = sla.eigs(sparse_mat, 2)[0]
-     print(evs.real)
- #HIDDEN_END_zydk
- 
- 
- def plot_conductance(syst, energies):
-     # Compute transmission as a function of energy
-     data = []
-     for energy in energies:
-         smatrix = kwant.smatrix(syst, energy)
-         data.append(smatrix.transmission(0, 1))
- 
--    pyplot.figure()
-+    fig = pyplot.figure()
-     pyplot.plot(energies, data)
--    pyplot.xlabel("energy [t]")
--    pyplot.ylabel("conductance [e^2/h]")
--    pyplot.show()
-+    pyplot.xlabel("energy [t]",
-+                 fontsize=_defs.mpl_label_size)
-+    pyplot.ylabel("conductance [e^2/h]",
-+                 fontsize=_defs.mpl_label_size)
-+    pyplot.setp(fig.get_axes()[0].get_xticklabels(),
-+               fontsize=_defs.mpl_tick_size)
-+    pyplot.setp(fig.get_axes()[0].get_yticklabels(),
-+               fontsize=_defs.mpl_tick_size)
-+    fig.set_size_inches(_defs.mpl_width_in, _defs.mpl_width_in * 3. / 4.)
-+    fig.subplots_adjust(left=0.15, right=0.95, top=0.95, bottom=0.15)
-+    for extension in ('pdf', 'png'):
-+        fig.savefig("graphene_result." + extension, dpi=_defs.dpi)
- 
- 
- def plot_bandstructure(flead, momenta):
-     bands = kwant.physics.Bands(flead)
-     energies = [bands(k) for k in momenta]
- 
--    pyplot.figure()
-+    fig = pyplot.figure()
-     pyplot.plot(momenta, energies)
--    pyplot.xlabel("momentum [(lattice constant)^-1]")
--    pyplot.ylabel("energy [t]")
--    pyplot.show()
-+    pyplot.xlabel("momentum [(lattice constant)^-1]",
-+                 fontsize=_defs.mpl_label_size)
-+    pyplot.ylabel("energy [t]",
-+                 fontsize=_defs.mpl_label_size)
-+    pyplot.setp(fig.get_axes()[0].get_xticklabels(),
-+               fontsize=_defs.mpl_tick_size)
-+    pyplot.setp(fig.get_axes()[0].get_yticklabels(),
-+               fontsize=_defs.mpl_tick_size)
-+    fig.set_size_inches(_defs.mpl_width_in, _defs.mpl_width_in * 3. / 4.)
-+    fig.subplots_adjust(left=0.15, right=0.95, top=0.95, bottom=0.15)
-+    for extension in ('pdf', 'png'):
-+        fig.savefig("graphene_bs." + extension, dpi=_defs.dpi)
- 
- 
- #HIDDEN The part of the following code block which begins with family_colors
- #HIDDEN is included verbatim in the tutorial text because nested code examples
- #HIDDEN are not supported.  Remember to update the tutorial text when you
- #HIDDEN modify this block.
- #HIDDEN_BEGIN_itkk
- def main():
-     pot = 0.1
-     syst, leads = make_system(pot=pot)
- 
-     # To highlight the two sublattices of graphene, we plot one with
-     # a filled, and the other one with an open circle:
-     def family_colors(site):
-         return 0 if site.family == a else 1
- 
--    # Plot the closed system without leads.
--    kwant.plot(syst, site_color=family_colors, site_lw=0.1, colorbar=False)
-+    size = (_defs.figwidth_in, _defs.figwidth_in)
-+    for extension in ('pdf', 'png'):
-+        kwant.plot(syst, site_color=family_colors, site_lw=0.1, colorbar=False,
-+                   file="graphene_syst1." + extension,
-+                   fig_size=size, dpi=_defs.dpi)
- #HIDDEN_END_itkk
- 
-     # Compute some eigenvalues.
- #HIDDEN_BEGIN_jmbi
-     compute_evs(syst.finalized())
- #HIDDEN_END_jmbi
- 
-     # Attach the leads to the system.
-     for lead in leads:
-         syst.attach_lead(lead)
- 
--    # Then, plot the system with leads.
--    kwant.plot(syst, site_color=family_colors, site_lw=0.1,
--               lead_site_lw=0, colorbar=False)
-+    size = (_defs.figwidth_in, 0.9 * _defs.figwidth_in)
-+    for extension in ('pdf', 'png'):
-+        kwant.plot(syst, site_color=family_colors, colorbar=False, site_lw=0.1,
-+                   file="graphene_syst2." + extension,
-+                   fig_size=size, dpi=_defs.dpi, lead_site_lw=0)
- 
-     # Finalize the system.
-     syst = syst.finalized()
- 
-     # Compute the band structure of lead 0.
-     momenta = [-pi + 0.02 * pi * i for i in range(101)]
-     plot_bandstructure(syst.leads[0], momenta)
- 
-     # Plot conductance.
-     energies = [-2 * pot + 4. / 50. * pot * i for i in range(51)]
-     plot_conductance(syst, energies)
- 
- 
- # Call the main function if the script gets executed (as opposed to imported).
- # See <http://docs.python.org/library/__main__.html>.
- if __name__ == '__main__':
-     main()

doc/source/code/figure/kernel_polynomial_method.py.diff

-@@ -1,341 +1,367 @@
- # Tutorial 2.8. Calculating spectral density with the Kernel Polynomial Method
- # ============================================================================
- #
- # Physics background
- # ------------------
- #  - Chebyshev polynomials, random trace approximation, spectral densities.
- #
- # Kwant features highlighted
- # --------------------------
- #  - kpm module,kwant operators.
- 
- import scipy
- 
-+import _defs
-+from contextlib import redirect_stdout
-+
- # For plotting
- from matplotlib import pyplot as plt
- 
- #HIDDEN_BEGIN_sys1
- # necessary imports
- import kwant
- import numpy as np
- 
- 
- # define the system
- def make_syst(r=30, t=-1, a=1):
-     syst = kwant.Builder()
-     lat = kwant.lattice.honeycomb(a, norbs=1)
- 
-     def circle(pos):
-         x, y = pos
-         return x ** 2 + y ** 2 < r ** 2
- 
-     syst[lat.shape(circle, (0, 0))] = 0.
-     syst[lat.neighbors()] = t
-     syst.eradicate_dangling()
- 
-     return syst
- #HIDDEN_END_sys1
- 
- #HIDDEN_BEGIN_sys2
- # define a Haldane system
- def make_syst_topo(r=30, a=1, t=1, t2=0.5):
-     syst = kwant.Builder()
-     lat = kwant.lattice.honeycomb(a, norbs=1, name=['a', 'b'])
- 
-     def circle(pos):
-         x, y = pos
-         return x ** 2 + y ** 2 < r ** 2
- 
-     syst[lat.shape(circle, (0, 0))] = 0.
-     syst[lat.neighbors()] = t
-     # add second neighbours hoppings
-     syst[lat.a.neighbors()] = 1j * t2
-     syst[lat.b.neighbors()] = -1j * t2
-     syst.eradicate_dangling()
- 
-     return lat, syst.finalized()
- #HIDDEN_END_sys2
- 
- 
- #HIDDEN_BEGIN_sys3
- # define the system
- def make_syst_staggered(r=30, t=-1, a=1, m=0.1):
-     syst = kwant.Builder()
-     lat = kwant.lattice.honeycomb(a, norbs=1)
- 
-     def circle(pos):
-         x, y = pos
-         return x ** 2 + y ** 2 < r ** 2
- 
-     syst[lat.a.shape(circle, (0, 0))] = m
-     syst[lat.b.shape(circle, (0, 0))] = -m
-     syst[lat.neighbors()] = t
-     syst.eradicate_dangling()
- 
-     return syst
- #HIDDEN_END_sys3
- 
- # Plot several density of states curves on the same axes.
--def plot_dos(labels_to_data):
-+def plot_dos(labels_to_data, file_name=None, ylabel="DoS [a.u.]"):
-     plt.figure(figsize=(5,4))
-     for label, (x, y) in labels_to_data:
-         plt.plot(x, y.real, label=label, linewidth=2)
-     plt.legend(loc=2, framealpha=0.5)
-     plt.xlabel("energy [t]")
-     plt.ylabel(ylabel)
--    plt.show()
-+    save_figure(file_name)
-     plt.clf()
- 
- 
- # Plot fill density of states plus curves on the same axes.
--def plot_dos_and_curves(dos labels_to_data):
-+def plot_dos_and_curves(dos, labels_to_data, file_name=None, ylabel="DoS [a.u.]"):
-     plt.figure(figsize=(5,4))
-     plt.fill_between(dos[0], dos[1], label="DoS [a.u.]",
-                      alpha=0.5, color='gray')
-     for label, (x, y) in labels_to_data:
-         plt.plot(x, y, label=label, linewidth=2)
-     plt.legend(loc=2, framealpha=0.5)
-     plt.xlabel("energy [t]")
-     plt.ylabel(ylabel)
--    plt.show()
-+    save_figure(file_name)
-     plt.clf()
- 
- def site_size_conversion(densities):
-     return 3 * np.abs(densities) / max(densities)
- 
- 
- # Plot several local density of states maps in different subplots
- def plot_ldos(syst, densities, file_name=None):
-     fig, axes = plt.subplots(1, len(densities), figsize=(7*len(densities), 7))
-     for ax, (title, rho) in zip(axes, densities):
-         kwant.plotter.density(syst, rho.real, ax=ax)
-         ax.set_title(title)
-         ax.set(adjustable='box', aspect='equal')
--    plt.show()
-+    save_figure(file_name)
-     plt.clf()
- 
-+def save_figure(file_name):
-+    if not file_name:
-+        return
-+    for extension in ('pdf', 'png'):
-+        plt.savefig('.'.join((file_name,extension)),
-+                    dpi=_defs.dpi, bbox_inches='tight')
-+
-+
- def simple_dos_example():
- #HIDDEN_BEGIN_kpm1
-     fsyst = make_syst().finalized()
- 
-     spectrum = kwant.kpm.SpectralDensity(fsyst, rng=0)
- #HIDDEN_END_kpm1
- 
- #HIDDEN_BEGIN_kpm2
-     energies, densities = spectrum()
- #HIDDEN_END_kpm2
- 
- #HIDDEN_BEGIN_kpm3
-     energy_subset = np.linspace(0, 2)
-     density_subset = spectrum(energy_subset)
- #HIDDEN_END_kpm3
- 
-     plot_dos([
-         ('densities', (energies, densities)),
-         ('density subset', (energy_subset, density_subset)),
--    ])
-+     ],
-+     file_name='kpm_dos'
-+    )
- 
- 
- def dos_integrating_example(fsyst):
-     spectrum = kwant.kpm.SpectralDensity(fsyst, rng=0)
- 
- #HIDDEN_BEGIN_int1
--    print('identity resolution:', spectrum.integrate())
-+    with open('kpm_normalization.txt', 'w') as f:
-+        with redirect_stdout(f):
-+            print('identity resolution:', spectrum.integrate())
- #HIDDEN_END_int1
- 
- #HIDDEN_BEGIN_int2
-     # Fermi energy 0.1 and temperature 0.2
-     fermi = lambda E: 1 / (np.exp((E - 0.1) / 0.2) + 1)
- 
--    print('number of filled states:', spectrum.integrate(fermi))
-+    with open('kpm_total_states.txt', 'w') as f:
-+        with redirect_stdout(f):
-+            print('number of filled states:', spectrum.integrate(fermi))
- #HIDDEN_END_int2
- 
- 
- def increasing_accuracy_example(fsyst):
-     spectrum = kwant.kpm.SpectralDensity(fsyst, rng=0)
-     original_dos = spectrum()  # get unaltered DoS
- 
- #HIDDEN_BEGIN_acc1
-     spectrum.add_moments(energy_resolution=0.03)
- #HIDDEN_END_acc1
- 
-     increased_resolution_dos = spectrum()
- 
-     plot_dos([
-         ('density', original_dos),
-         ('higher energy resolution', increased_resolution_dos),
--    ])
-+     ],
-+     file_name='kpm_dos_acc'
-+    )
- 
- #HIDDEN_BEGIN_acc2
-     spectrum.add_moments(100)
-     spectrum.add_vectors(5)
- #HIDDEN_END_acc2
- 
-     increased_moments_dos = spectrum()
- 
-     plot_dos([
-         ('density', original_dos),
-         ('higher number of moments', increased_moments_dos),
--    ])
-+     ],
-+     file_name='kpm_dos_r'
-+    )
- 
- 
- def operator_example(fsyst):
- #HIDDEN_BEGIN_op1
-     # identity matrix
-     matrix_op = scipy.sparse.eye(len(fsyst.sites))
-     matrix_spectrum = kwant.kpm.SpectralDensity(fsyst, operator=matrix_op, rng=0)
- #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, rng=0)
- #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, rng=0)
- #HIDDEN_END_op3
- 
- #HIDDEN_BEGIN_op4
-     zero_energy_ldos = local_dos(energy=0)
-     finite_energy_ldos = local_dos(energy=1)
-     plot_ldos(fsyst, [
-         ('energy = 0', zero_energy_ldos),
-         ('energy = 1', finite_energy_ldos)
--    ])
-+     ],
-+     file_name='kpm_ldos'
-+    )
- #HIDDEN_END_op4
- 
- 
- def ldos_sites_example():
-     fsyst = make_syst_staggered().finalized()
- #HIDDEN_BEGIN_op5
-     # find 'A' and 'B' sites in the unit cell at the center of the disk
-     center_tag = np.array([0, 0])
-     where = lambda s: s.tag == center_tag
-     # make local vectors
-     vector_factory = kwant.kpm.LocalVectors(fsyst, where)
- #HIDDEN_END_op5
- 
- #HIDDEN_BEGIN_op6
-     # 'num_vectors' can be unspecified when using 'LocalVectors'
-     local_dos = kwant.kpm.SpectralDensity(fsyst, num_vectors=None,
-                                           vector_factory=vector_factory,
-                                           mean=False,
-                                           rng=0)
-     energies, densities = local_dos()
-     plot_dos([
-         ('A sublattice', (energies, densities[:, 0])),
-         ('B sublattice', (energies, densities[:, 1])),
--    ])
-+     ],
-+     file_name='kpm_ldos_sites'
-+    )
- #HIDDEN_END_op6
- 
- 
- def vector_factory_example(fsyst):
-     spectrum = kwant.kpm.SpectralDensity(fsyst, rng=0)
- #HIDDEN_BEGIN_fact1
-     # construct a generator of vectors with n random elements -1 or +1.
-     n = fsyst.hamiltonian_submatrix(sparse=True).shape[0]
-     def binary_vectors():
-         while True:
-            yield np.rint(np.random.random_sample(n)) * 2 - 1
- 
-     custom_factory = kwant.kpm.SpectralDensity(fsyst,
-                                                vector_factory=binary_vectors(),
-                                                rng=0)
- #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, rng=0)
-     rho_alt_spectrum = kwant.kpm.SpectralDensity(fsyst, operator=rho_alt, rng=0)
- #HIDDEN_END_blm
- 
-     plot_dos([
-         ('kwant.operator.Density', rho_spectrum()),
-         ('bilinear operator', rho_alt_spectrum()),
-     ])
- 
- def conductivity_example():
- #HIDDEN_BEGIN_cond
-     # construct the Haldane model
-     lat, fsyst = make_syst_topo()
-     # find 'A' and 'B' sites in the unit cell at the center of the disk
-     where = lambda s: np.linalg.norm(s.pos) < 1
- 
-     # component 'xx'
-     s_factory = kwant.kpm.LocalVectors(fsyst, where)
-     cond_xx = kwant.kpm.conductivity(fsyst, alpha='x', beta='x', mean=True,
-                                      num_vectors=None, vector_factory=s_factory,
-                                      rng=0)
-     # component 'xy'
-     s_factory = kwant.kpm.LocalVectors(fsyst, where)
-     cond_xy = kwant.kpm.conductivity(fsyst, alpha='x', beta='y', mean=True,
-                                      num_vectors=None, vector_factory=s_factory,
-                                      rng=0)
- 
-     energies = cond_xx.energies
-     cond_array_xx = np.array([cond_xx(e, temperature=0.01) for e in energies])
-     cond_array_xy = np.array([cond_xy(e, temperature=0.01) for e in energies])
- 
-     # area of the unit cell per site
-     area_per_site = np.abs(np.cross(*lat.prim_vecs)) / len(lat.sublattices)
-     cond_array_xx /= area_per_site
-     cond_array_xy /= area_per_site
- #HIDDEN_END_cond
-     # ldos
-     s_factory = kwant.kpm.LocalVectors(fsyst, where)
-     spectrum = kwant.kpm.SpectralDensity(fsyst, num_vectors=None,
-                                          vector_factory=s_factory,
-                                          rng=0)
- 
-     plot_dos_and_curves(
-     (spectrum.energies, spectrum.densities * 8),
-     [
-         (r'Longitudinal conductivity $\sigma_{xx} / 4$',
-          (energies, cond_array_xx / 4)),
-         (r'Hall conductivity $\sigma_{xy}$',
-          (energies, cond_array_xy))],
-         ylabel=r'$\sigma [e^2/h]$',
-         file_name='kpm_cond'
-     )
- 
- 
- def main():
-     simple_dos_example()
- 
-     fsyst = make_syst().finalized()
- 
-     dos_integrating_example(fsyst)
-     increasing_accuracy_example(fsyst)
-     operator_example(fsyst)
-     ldos_example(fsyst)
-     ldos_sites_example()
-     vector_factory_example(fsyst)
-     bilinear_map_operator_example(fsyst)
-     conductivity_example()
- 
- 
- # Call the main function if the script gets executed (as opposed to imported).
- # See <http://docs.python.org/library/__main__.html>.
- if __name__ == '__main__':
-     main()

doc/source/code/figure/magnetic_texture.py.diff

-@@ -1,245 +1,268 @@
- # Tutorial 2.7. Spin textures
- # ===========================
- #
- # Physics background
- # ------------------
- #  - Spin textures
- #  - Skyrmions
- #
- # Kwant features highlighted
- # --------------------------
- #  - operators
- #  - plotting vector fields
- 
-+import _defs
-+from contextlib import redirect_stdout
-+
- from math import sin, cos, tanh, pi
- import itertools
- import numpy as np
- import tinyarray as ta
- import matplotlib.pyplot as plt
- 
- import kwant
- 
- sigma_0 = ta.array([[1, 0], [0, 1]])
- sigma_x = ta.array([[0, 1], [1, 0]])
- sigma_y = ta.array([[0, -1j], [1j, 0]])
- sigma_z = ta.array([[1, 0], [0, -1]])
- 
- # vector of Pauli matrices σ_αiβ where greek
- # letters denote spinor indices
- sigma = np.rollaxis(np.array([sigma_x, sigma_y, sigma_z]), 1)
- 
- #HIDDEN_BEGIN_model
- def field_direction(pos, r0, delta):
-     x, y = pos
-     r = np.linalg.norm(pos)
-     r_tilde = (r - r0) / delta
-     theta = (tanh(r_tilde) - 1) * (pi / 2)
- 
-     if r == 0:
-         m_i = [0, 0, -1]
-     else:
-         m_i = [
-             (x / r) * sin(theta),
-             (y / r) * sin(theta),
-             cos(theta),
-         ]
- 
-     return np.array(m_i)
- 
- 
- def scattering_onsite(site, r0, delta, J):
-     m_i = field_direction(site.pos, r0, delta)
-     return J * np.dot(m_i, sigma)
- 
- 
- def lead_onsite(site, J):
-     return J * sigma_z
- #HIDDEN_END_model
- 
- 
- #HIDDEN_BEGIN_syst
- lat = kwant.lattice.square(norbs=2)
- 
- def make_system(L=80):
- 
-     syst = kwant.Builder()
- 
-     def square(pos):
-         return all(-L/2 < p < L/2 for p in pos)
- 
-     syst[lat.shape(square, (0, 0))] = scattering_onsite
-     syst[lat.neighbors()] = -sigma_0
- 
-     lead = kwant.Builder(kwant.TranslationalSymmetry((-1, 0)),
-                          conservation_law=-sigma_z)
- 
-     lead[lat.shape(square, (0, 0))] = lead_onsite
-     lead[lat.neighbors()] = -sigma_0
- 
-     syst.attach_lead(lead)
-     syst.attach_lead(lead.reversed())
- 
-     return syst
- #HIDDEN_END_syst
- 
- 
--def plot_vector_field(syst, params):
-+def save_plot(fname):
-+    for extension in ('.pdf', '.png'):
-+        plt.savefig(fname + extension,
-+                    dpi=_defs.dpi, bbox_inches='tight')
-+
-+
-+def plot_vector_field(syst, params, fname=None):
-     xmin, ymin = min(s.tag for s in syst.sites)
-     xmax, ymax = max(s.tag for s in syst.sites)
-     x, y = np.meshgrid(np.arange(xmin, xmax+1), np.arange(ymin, ymax+1))
- 
-     m_i = [field_direction(p, **params) for p in zip(x.flat, y.flat)]
-     m_i = np.reshape(m_i, x.shape + (3,))
-     m_i = np.rollaxis(m_i, 2, 0)
- 
--    fig, ax = plt.subplots(1, 1)
-+    fig, ax = plt.subplots(1, 1, figsize=(9, 7))
-     im = ax.quiver(x, y, *m_i, pivot='mid', scale=75)
-     fig.colorbar(im)
--    plt.show()
-+    if fname:
-+        save_plot(fname)
-+    else:
-+        plt.show()
- 
- 
--def plot_densities(syst, densities):
--    fig, axes = plt.subplots(1, len(densities))
-+def plot_densities(syst, densities, fname=None):
-+    fig, axes = plt.subplots(1, len(densities), figsize=(7*len(densities), 7))
-     for ax, (title, rho) in zip(axes, densities):
-         kwant.plotter.density(syst, rho, ax=ax)
-         ax.set_title(title)
--    plt.show()
-+    if fname:
-+        save_plot(fname)
-+    else:
-+        plt.show()
- 
- 
--def plot_currents(syst, currents):
--    fig, axes = plt.subplots(1, len(currents))
-+
-+def plot_currents(syst, currents, fname=None):
-+    fig, axes = plt.subplots(1, len(currents), figsize=(7*len(currents), 7))
-     if not hasattr(axes, '__len__'):
-         axes = (axes,)
-     for ax, (title, current) in zip(axes, currents):
-         kwant.plotter.current(syst, current, ax=ax, colorbar=False)
-         ax.set_title(title)
--    plt.show()
-+    if fname:
-+        save_plot(fname)
-+    else:
-+        plt.show()
- 
- 
- def main():
-     syst = make_system().finalized()
- 
- #HIDDEN_BEGIN_wavefunction
-     params = dict(r0=20, delta=10, J=1)
-     wf = kwant.wave_function(syst, energy=-1, params=params)
-     psi = wf(0)[0]
- #HIDDEN_END_wavefunction
- 
--    plot_vector_field(syst, dict(r0=20, delta=10))
-+    plot_vector_field(syst, dict(r0=20, delta=10), fname='mag_field_direction')
- 
- #HIDDEN_BEGIN_ldos
-     # even (odd) indices correspond to spin up (down)
-     up, down = psi[::2], psi[1::2]
-     density = np.abs(up)**2 + np.abs(down)**2
- #HIDDEN_END_ldos
- 
- #HIDDEN_BEGIN_lsdz
-     # spin down components have a minus sign
-     spin_z = np.abs(up)**2 - np.abs(down)**2
- #HIDDEN_END_lsdz
- 
- #HIDDEN_BEGIN_lsdy
-     # spin down components have a minus sign
-     spin_y = 1j * (down.conjugate() * up - up.conjugate() * down)
- #HIDDEN_END_lsdy
- 
- #HIDDEN_BEGIN_lden
-     rho = kwant.operator.Density(syst)
-     rho_sz = kwant.operator.Density(syst, sigma_z)
-     rho_sy = kwant.operator.Density(syst, sigma_y)
- 
-     # calculate the expectation values of the operators with 'psi'
-     density = rho(psi)
-     spin_z = rho_sz(psi)
-     spin_y = rho_sy(psi)
- #HIDDEN_END_lden
- 
-     plot_densities(syst, [
-         ('$σ_0$', density),
-         ('$σ_z$', spin_z),
-         ('$σ_y$', spin_y),
--    ])
-+    ], fname='spin_densities')
- 
- #HIDDEN_BEGIN_current
-     J_0 = kwant.operator.Current(syst)
-     J_z = kwant.operator.Current(syst, sigma_z)
-     J_y = kwant.operator.Current(syst, sigma_y)
- 
-     # calculate the expectation values of the operators with 'psi'
-     current = J_0(psi)
-     spin_z_current = J_z(psi)
-     spin_y_current = J_y(psi)
- #HIDDEN_END_current
- 
-     plot_currents(syst, [
-         ('$J_{σ_0}$', current),
-         ('$J_{σ_z}$', spin_z_current),
-         ('$J_{σ_y}$', spin_y_current),
--    ])
-+    ], fname='spin_currents')
- 
- #HIDDEN_BEGIN_following
-     def following_m_i(site, r0, delta):
-         m_i = field_direction(site.pos, r0, delta)
-         return np.dot(m_i, sigma)
- 
-     J_m = kwant.operator.Current(syst, following_m_i)
- 
-     # evaluate the operator
-     m_current = J_m(psi, params=dict(r0=25, delta=10))
- #HIDDEN_END_following
- 
-     plot_currents(syst, [
-         (r'$J_{\mathbf{m}_i}$', m_current),
-         ('$J_{σ_z}$', spin_z_current),
--    ])
-+    ], fname='spin_current_comparison')
- 
- 
- #HIDDEN_BEGIN_density_cut
-     def circle(site):
-         return np.linalg.norm(site.pos) < 20
- 
-     rho_circle = kwant.operator.Density(syst, where=circle, sum=True)
- 
-     all_states = np.vstack((wf(0), wf(1)))
-     dos_in_circle = sum(rho_circle(p) for p in all_states) / (2 * pi)
--    print('density of states in circle:', dos_in_circle)
-+    with open('circle_dos.txt', 'w') as f:
-+        with redirect_stdout(f):
-+            print('density of states in circle:', dos_in_circle)
- #HIDDEN_END_density_cut
- 
- #HIDDEN_BEGIN_current_cut
-     def left_cut(site_to, site_from):
-         return site_from.pos[0] <= -39 and site_to.pos[0] > -39
- 
-     def right_cut(site_to, site_from):
-         return site_from.pos[0] < 39 and site_to.pos[0] >= 39
- 
-     J_left = kwant.operator.Current(syst, where=left_cut, sum=True)
-     J_right = kwant.operator.Current(syst, where=right_cut, sum=True)
- 
-     Jz_left = kwant.operator.Current(syst, sigma_z, where=left_cut, sum=True)
-     Jz_right = kwant.operator.Current(syst, sigma_z, where=right_cut, sum=True)
- 
--    print('J_left:', J_left(psi), ' J_right:', J_right(psi))
--    print('Jz_left:', Jz_left(psi), ' Jz_right:', Jz_right(psi))
-+    with open('current_cut.txt', 'w') as f:
-+        with redirect_stdout(f):
-+            print('J_left:', J_left(psi), ' J_right:', J_right(psi))
-+            print('Jz_left:', Jz_left(psi), ' Jz_right:', Jz_right(psi))
- #HIDDEN_END_current_cut
- 
- #HIDDEN_BEGIN_bind
-     J_m = kwant.operator.Current(syst, following_m_i)
-     J_z = kwant.operator.Current(syst, sigma_z)
- 
-     J_m_bound = J_m.bind(params=dict(r0=25, delta=10, J=1))
-     J_z_bound = J_z.bind(params=dict(r0=25, delta=10, J=1))
- 
-     # Sum current local from all scattering states on the left at energy=-1
-     wf_left = wf(0)
-     J_m_left = sum(J_m_bound(p) for p in wf_left)
-     J_z_left = sum(J_z_bound(p) for p in wf_left)
- #HIDDEN_END_bind
- 
-     plot_currents(syst, [
-         (r'$J_{\mathbf{m}_i}$ (from left)', J_m_left),
-         (r'$J_{σ_z}$ (from left)', J_z_left),
--    ])
-+    ], fname='bound_current')
- 
- 
- if __name__ == '__main__':
-     main()

doc/source/code/figure/plot_graphene.py.diff

-@@ -1,112 +1,130 @@
- # Tutorial 2.8.1. 2D example: graphene quantum dot
- # ================================================
- #
- # Physics background
- # ------------------
- #  - graphene edge states
- #
- # Kwant features highlighted
- # --------------------------
- #  - demonstrate different ways of plotting
- 
-+import _defs
- import kwant
- from matplotlib import pyplot
- 
- #HIDDEN_BEGIN_makesyst
- lat = kwant.lattice.honeycomb()
- a, b = lat.sublattices
- 
- def make_system(r=8, t=-1, tp=-0.1):
- 
-     def circle(pos):
-         x, y = pos
-         return x**2 + y**2 < r**2
- 
-     syst = kwant.Builder()
--    syst[lat.shape(circle, (0, 0))] = 0
-+    syst[lat.shape(circle, (0,0))] = 0
-     syst[lat.neighbors()] = t
-     syst.eradicate_dangling()
-     if tp:
-         syst[lat.neighbors(2)] = tp
- 
-     return syst
- #HIDDEN_END_makesyst
- 
- 
- #HIDDEN_BEGIN_plotsyst1
- def plot_system(syst):
--    kwant.plot(syst)
- #HIDDEN_END_plotsyst1
--    # the standard plot is ok, but not very intelligible. One can do
--    # better by playing wioth colors and linewidths
-+    # standard plot - not very intelligible for this particular situation
-+    size = (_defs.figwidth_in, _defs.figwidth_in)
-+    for extension in ('pdf', 'png'):
-+        kwant.plot(syst, file="plot_graphene_syst1." + extension,
-+                   fig_size=size, dpi=_defs.dpi)
- 
-     # use color and linewidths to get a better plot
- #HIDDEN_BEGIN_plotsyst2
-     def family_color(site):
-         return 'black' if site.family == a else 'white'
- 
-     def hopping_lw(site1, site2):
-         return 0.04 if site1.family == site2.family else 0.1
- 
--    kwant.plot(syst, site_lw=0.1, site_color=family_color, hop_lw=hopping_lw)
-+    size = (_defs.figwidth_in, _defs.figwidth_in)
-+    for extension in ('pdf', 'png'):
-+        kwant.plot(syst, site_lw=0.1, site_color=family_color,
-+                   hop_lw=hopping_lw, file="plot_graphene_syst2." + extension,
-+                   fig_size=size, dpi=_defs.dpi)
- #HIDDEN_END_plotsyst2
- 
- 
- #HIDDEN_BEGIN_plotdata1
- def plot_data(syst, n):
-     import scipy.linalg as la
- 
-     syst = syst.finalized()
-     ham = syst.hamiltonian_submatrix()
-     evecs = la.eigh(ham)[1]
- 
-     wf = abs(evecs[:, n])**2
- #HIDDEN_END_plotdata1
- 
-     # the usual - works great in general, looks just a bit crufty for
-     # small systems
- 
- #HIDDEN_BEGIN_plotdata2
--    kwant.plotter.map(syst, wf, oversampling=10, cmap='gist_heat_r')
-+    size = (_defs.figwidth_in, _defs.figwidth_in)
-+    for extension in ('pdf', 'png'):
-+        kwant.plotter.map(syst, wf, oversampling=10, cmap='gist_heat_r',
-+                          file="plot_graphene_data1." + extension,
-+                          fig_size=size, dpi=_defs.dpi)
- #HIDDEN_END_plotdata2
- 
-     # use two different sort of triangles to cleanly fill the space
- #HIDDEN_BEGIN_plotdata3
-     def family_shape(i):
-         site = syst.sites[i]
-         return ('p', 3, 180) if site.family == a else ('p', 3, 0)
- 
-     def family_color(i):
-         return 'black' if syst.sites[i].family == a else 'white'
- 
--    kwant.plot(syst, site_color=wf, site_symbol=family_shape,
--               site_size=0.5, hop_lw=0, cmap='gist_heat_r')
-+    size = (_defs.figwidth_in, _defs.figwidth_in)
-+    for extension in ('pdf', 'png'):
-+        kwant.plot(syst, site_color=wf, site_symbol=family_shape,
-+                   site_size=0.5, hop_lw=0, cmap='gist_heat_r',
-+                   file="plot_graphene_data2." + extension,
-+                   fig_size=size, dpi=_defs.dpi)
- #HIDDEN_END_plotdata3
- 
-     # plot by changing the symbols itself
- #HIDDEN_BEGIN_plotdata4
-     def site_size(i):
-         return 3 * wf[i] / wf.max()
- 
--    kwant.plot(syst, site_size=site_size, site_color=(0, 0, 1, 0.3),
--               hop_lw=0.1)
-+    size = (_defs.figwidth_in, _defs.figwidth_in)
-+    for extension in ('pdf', 'png'):
-+        kwant.plot(syst, site_size=site_size, site_color=(0,0,1,0.3),
-+                   hop_lw=0.1, file="plot_graphene_data3." + extension,
-+                   fig_size=size, dpi=_defs.dpi)
- #HIDDEN_END_plotdata4
- 
- 
- def main():
-     # plot the graphene system in different styles
-     syst = make_system()
- 
-     plot_system(syst)
- 
-     # compute a wavefunction (number 225) and plot it in different
-     # styles
-     syst = make_system(tp=0)
- 
-     plot_data(syst, 225)
- 
- 
- # Call the main function if the script gets executed (as opposed to imported).
- # See <http://docs.python.org/library/__main__.html>.
- if __name__ == '__main__':
-     main()

doc/source/code/figure/plot_qahe.py.diff

-@@ -1,72 +1,76 @@
- # Comprehensive example: quantum anomalous Hall effect
- # ====================================================
- #
- # Physics background
- # ------------------
- # + Quantum anomalous Hall effect
- #
- # Features highlighted
- # --------------------
- # + Use of `kwant.continuum` to discretize a continuum Hamiltonian
- # + Use of `kwant.operator` to compute local current
- # + Use of `kwant.plotter.current` to plot local current
- 
-+import _defs
- import math
- import matplotlib.pyplot
- import kwant
- import kwant.continuum
- 
- 
- # 2 band model exhibiting quantum anomalous Hall effect
- #HIDDEN_BEGIN_model
- def make_model(a):
-     ham = ("alpha * (k_x * sigma_x - k_y * sigma_y)"
-            "+ (m + beta * kk) * sigma_z"
-            "+ (gamma * kk + U) * sigma_0")
-     subs = {"kk": "k_x**2 + k_y**2"}
-     return kwant.continuum.discretize(ham, locals=subs, grid=a)
- #HIDDEN_END_model
- 
- 
- def make_system(model, L):
-     def lead_shape(site):
-         x, y = site.pos / L
-         return abs(y) < 0.5
- 
-     # QPC shape: a rectangle with 2 gaussians
-     # etched out of the top and bottom edge.
-     def central_shape(site):
-         x, y = site.pos / L
-         return abs(x) < 3/5 and abs(y) < 0.5 - 0.4 * math.exp(-40 * x**2)
- 
-     lead = kwant.Builder(kwant.TranslationalSymmetry(
-         model.lattice.vec((-1, 0))))
-     lead.fill(model, lead_shape, (0, 0))
- 
-     syst = kwant.Builder()
-     syst.fill(model, central_shape, (0, 0))
-     syst.attach_lead(lead)
-     syst.attach_lead(lead.reversed())
- 
-     return syst.finalized()
- 
- 
- def main():
-     # Set up our model and system, and define the model parameters.
-     params = dict(alpha=0.365, beta=0.686, gamma=0.512, m=-0.01, U=0)
-     model = make_model(1)
-     syst = make_system(model, 70)
-     kwant.plot(syst)
- 
-     # Calculate the scattering states at energy 'm' coming from the left
-     # lead, and the associated particle current.
-     psi = kwant.wave_function(syst, energy=params['m'], params=params)(0)
- #HIDDEN_BEGIN_current
-     J = kwant.operator.Current(syst).bind(params=params)
-     current = sum(J(p) for p in psi)
--    kwant.plotter.current(syst, current)
-+    for extension in ('pdf', 'png'):
-+        kwant.plotter.current(syst, current,
-+                              file="plot_qahe_current." + extension,
-+                              dpi=_defs.dpi)
- #HIDDEN_END_current
- 
- 
- if __name__ == '__main__':
-     main()

doc/source/code/figure/plot_zincblende.py.diff

-@@ -1,59 +1,67 @@
- # Tutorial 2.8.2. 3D example: zincblende structure
- # ================================================
- #
- # Physical background
- # -------------------
- #  - 3D Bravais lattices
- #
- # Kwant features highlighted
- # --------------------------
- #  - demonstrate different ways of plotting in 3D
- 
-+import _defs
- import kwant
- from matplotlib import pyplot
- 
- #HIDDEN_BEGIN_zincblende1
- lat = kwant.lattice.general([(0, 0.5, 0.5), (0.5, 0, 0.5), (0.5, 0.5, 0)],
-                             [(0, 0, 0), (0.25, 0.25, 0.25)])
- a, b = lat.sublattices
- #HIDDEN_END_zincblende1
- 
- #HIDDEN_BEGIN_zincblende2
- def make_cuboid(a=15, b=10, c=5):
-     def cuboid_shape(pos):
-         x, y, z = pos
-         return 0 <= x < a and 0 <= y < b and 0 <= z < c
- 
-     syst = kwant.Builder()
-     syst[lat.shape(cuboid_shape, (0, 0, 0))] = None
-     syst[lat.neighbors()] = None
- 
-     return syst
- #HIDDEN_END_zincblende2
- 
- 
- def main():
-     # the standard plotting style for 3D is mainly useful for
-     # checking shapes:
- #HIDDEN_BEGIN_plot1
-     syst = make_cuboid()
- 
--    kwant.plot(syst)
-+    size = (_defs.figwidth_in, _defs.figwidth_in)
-+    for extension in ('pdf', 'png'):
-+        kwant.plot(syst, file="plot_zincblende_syst1." + extension,
-+                   fig_size=size, dpi=_defs.dpi)
- #HIDDEN_END_plot1
- 
-     # visualize the crystal structure better for a very small system
- #HIDDEN_BEGIN_plot2
-     syst = make_cuboid(a=1.5, b=1.5, c=1.5)
- 
-     def family_colors(site):
-         return 'r' if site.family == a else 'g'
- 
--    kwant.plot(syst, site_size=0.18, site_lw=0.01, hop_lw=0.05,
--               site_color=family_colors)
-+    size = (_defs.figwidth_in, _defs.figwidth_in)
-+    for extension in ('pdf', 'png'):
-+        kwant.plot(syst, site_size=0.18, site_lw=0.01, hop_lw=0.05,
-+                   site_color=family_colors,
-+                   file="plot_zincblende_syst2." + extension,
-+                   fig_size=size, dpi=_defs.dpi)
- #HIDDEN_END_plot2
- 
- 
- # Call the main function if the script gets executed (as opposed to imported).
- # See <http://docs.python.org/library/__main__.html>.
- if __name__ == '__main__':
-     main()

doc/source/code/figure/quantum_well.py.diff

-@@ -1,88 +1,95 @@
- # Tutorial 2.3.2. Spatially dependent values through functions
- # ============================================================
- #
- # Physics background
- # ------------------
- #  transmission through a quantum well
- #
- # Kwant features highlighted
- # --------------------------
- #  - Functions as values in Builder
- 
-+import _defs
- import kwant
- 
- # For plotting
- from matplotlib import pyplot
- 
- 
- #HIDDEN_BEGIN_ehso
- def make_system(a=1, t=1.0, W=10, L=30, L_well=10):
-     # Start with an empty tight-binding system and a single square lattice.
-     # `a` is the lattice constant (by default set to 1 for simplicity.
-     lat = kwant.lattice.square(a)
- 
-     syst = kwant.Builder()
- 
-     #### Define the scattering region. ####
-     # Potential profile
-     def potential(site, pot):
-         (x, y) = site.pos
-         if (L - L_well) / 2 < x < (L + L_well) / 2:
-             return pot
-         else:
-             return 0
- #HIDDEN_END_ehso
- 
- #HIDDEN_BEGIN_coid
-     def onsite(site, pot):
-         return 4 * t + potential(site, pot)
- 
-     syst[(lat(x, y) for x in range(L) for y in range(W))] = onsite
-     syst[lat.neighbors()] = -t
- #HIDDEN_END_coid
- 
-     #### Define and attach the leads. ####
-     lead = kwant.Builder(kwant.TranslationalSymmetry((-a, 0)))
-     lead[(lat(0, j) for j in range(W))] = 4 * t
-     lead[lat.neighbors()] = -t
-     syst.attach_lead(lead)
-     syst.attach_lead(lead.reversed())
- 
-     return syst
- 
- 
- def plot_conductance(syst, energy, welldepths):
- #HIDDEN_BEGIN_sqvr
- 
-     # Compute conductance
-     data = []
-     for welldepth in welldepths:
-         smatrix = kwant.smatrix(syst, energy, params=dict(pot=-welldepth))
-         data.append(smatrix.transmission(1, 0))
- 
--    pyplot.figure()
-+    fig = pyplot.figure()
-     pyplot.plot(welldepths, data)
--    pyplot.xlabel("well depth [t]")
--    pyplot.ylabel("conductance [e^2/h]")
--    pyplot.show()
-+    pyplot.xlabel("well depth [t]",
-+                 fontsize=_defs.mpl_label_size)
-+    pyplot.ylabel("conductance [e^2/h]",
-+                 fontsize=_defs.mpl_label_size)
-+    pyplot.setp(fig.get_axes()[0].get_xticklabels(),
-+               fontsize=_defs.mpl_tick_size)
-+    pyplot.setp(fig.get_axes()[0].get_yticklabels(),
-+               fontsize=_defs.mpl_tick_size)
-+    fig.set_size_inches(_defs.mpl_width_in, _defs.mpl_width_in * 3. / 4.)
-+    fig.subplots_adjust(left=0.15, right=0.95, top=0.95, bottom=0.15)
-+    for extension in ('pdf', 'png'):
-+        fig.savefig("quantum_well_result." + extension, dpi=_defs.dpi)
- #HIDDEN_END_sqvr
- 
- 
- def main():
-     syst = make_system()
- 
--    # Check that the system looks as intended.
--    kwant.plot(syst)
--
-     # Finalize the system.
-     syst = syst.finalized()
- 
-     # We should see conductance steps.
-     plot_conductance(syst, energy=0.2,
-                      welldepths=[0.01 * i for i in range(100)])
- 
- 
- # Call the main function if the script gets executed (as opposed to imported).
- # See <http://docs.python.org/library/__main__.html>.
- if __name__ == '__main__':
-     main()

doc/source/code/figure/quantum_wire.py.diff

-@@ -1,121 +1,130 @@
- # Tutorial 2.2.2. Transport through a quantum wire
- # ================================================
- #
- # Physics background
- # ------------------
- #  Conductance of a quantum wire; subbands
- #
- # Kwant features highlighted
- # --------------------------
- #  - Builder for setting up transport systems easily
- #  - Making scattering region and leads
- #  - Using the simple sparse solver for computing Landauer conductance
- 
-+import _defs
- from matplotlib import pyplot
- #HIDDEN_BEGIN_dwhx
- import kwant
- #HIDDEN_END_dwhx
- 
- # First, define the tight-binding system
- 
- #HIDDEN_BEGIN_goiq
- syst = kwant.Builder()
- #HIDDEN_END_goiq
- 
- # Here, we are only working with square lattices
- #HIDDEN_BEGIN_suwo
- a = 1
- lat = kwant.lattice.square(a)
- #HIDDEN_END_suwo
- 
- #HIDDEN_BEGIN_zfvr
- t = 1.0
- W = 10
- L = 30
- 
- # Define the scattering region
- 
- for i in range(L):
-     for j in range(W):
-         # On-site Hamiltonian
-         syst[lat(i, j)] = 4 * t
- 
-         # Hopping in y-direction
-         if j > 0:
-             syst[lat(i, j), lat(i, j - 1)] = -t
- 
-         # Hopping in x-direction
-         if i > 0:
-             syst[lat(i, j), lat(i - 1, j)] = -t
- #HIDDEN_END_zfvr
- 
- # Then, define and attach the leads:
- 
- # First the lead to the left
- # (Note: TranslationalSymmetry takes a real-space vector)
- #HIDDEN_BEGIN_xcmc
- sym_left_lead = kwant.TranslationalSymmetry((-a, 0))
- left_lead = kwant.Builder(sym_left_lead)
- #HIDDEN_END_xcmc
- 
- #HIDDEN_BEGIN_ndez
- for j in range(W):
-     left_lead[lat(0, j)] = 4 * t
-     if j > 0:
-         left_lead[lat(0, j), lat(0, j - 1)] = -t
-     left_lead[lat(1, j), lat(0, j)] = -t
- #HIDDEN_END_ndez
- 
- #HIDDEN_BEGIN_fskr
- syst.attach_lead(left_lead)
- #HIDDEN_END_fskr
- 
- # Then the lead to the right
- #HIDDEN_BEGIN_xhqc
- sym_right_lead = kwant.TranslationalSymmetry((a, 0))
- right_lead = kwant.Builder(sym_right_lead)
- 
- for j in range(W):
-     right_lead[lat(0, j)] = 4 * t
-     if j > 0:
-         right_lead[lat(0, j), lat(0, j - 1)] = -t
-     right_lead[lat(1, j), lat(0, j)] = -t
- 
- syst.attach_lead(right_lead)
- #HIDDEN_END_xhqc
- 
- # Plot it, to make sure it's OK
- #HIDDEN_BEGIN_wsgh
--kwant.plot(syst)
-+size = (_defs.figwidth_in, 0.3 * _defs.figwidth_in)
-+for extension in ('pdf', 'png'):
-+    kwant.plot(syst, file="quantum_wire_syst." + extension,
-+               fig_size=size, dpi=_defs.dpi)
- #HIDDEN_END_wsgh
- 
- # Finalize the system
- #HIDDEN_BEGIN_dngj
- syst = syst.finalized()
- #HIDDEN_END_dngj
- 
- # Now that we have the system, we can compute conductance
- #HIDDEN_BEGIN_buzn
- energies = []
- data = []
- for ie in range(100):
-     energy = ie * 0.01
- 
-     # compute the scattering matrix at a given energy
-     smatrix = kwant.smatrix(syst, energy)
- 
-     # compute the transmission probability from lead 0 to
-     # lead 1
-     energies.append(energy)
-     data.append(smatrix.transmission(1, 0))
- #HIDDEN_END_buzn
- 
- # Use matplotlib to write output
- # We should see conductance steps
- #HIDDEN_BEGIN_lliv
--pyplot.figure()
-+fig = pyplot.figure()
- pyplot.plot(energies, data)
--pyplot.xlabel("energy [t]")
--pyplot.ylabel("conductance [e^2/h]")
--pyplot.show()
-+pyplot.xlabel("energy [t]", fontsize=_defs.mpl_label_size)
-+pyplot.ylabel("conductance [e^2/h]", fontsize=_defs.mpl_label_size)
-+pyplot.setp(fig.get_axes()[0].get_xticklabels(), fontsize=_defs.mpl_tick_size)
-+pyplot.setp(fig.get_axes()[0].get_yticklabels(), fontsize=_defs.mpl_tick_size)
-+fig.set_size_inches(_defs.mpl_width_in, _defs.mpl_width_in * 3. / 4.)
-+fig.subplots_adjust(left=0.15, right=0.95, top=0.95, bottom=0.15)
-+for extension in ('pdf', 'png'):
-+    fig.savefig("quantum_wire_result." + extension, dpi=_defs.dpi)
- #HIDDEN_END_lliv

doc/source/code/figure/spin_orbit.py.diff

-@@ -1,104 +1,110 @@
- # Tutorial 2.3.1. Matrix structure of on-site and hopping elements
- # ================================================================
- #
- # Physics background
- # ------------------
- #  Gaps in quantum wires with spin-orbit coupling and Zeeman splititng,
- #  as theoretically predicted in
- #   http://prl.aps.org/abstract/PRL/v90/i25/e256601
- #  and (supposedly) experimentally oberved in
- #   http://www.nature.com/nphys/journal/v6/n5/abs/nphys1626.html
- #
- # Kwant features highlighted
- # --------------------------
- #  - Numpy matrices as values in Builder
- 
-+import _defs
- import kwant
- 
- # For plotting
- from matplotlib import pyplot
- 
- # For matrix support
- #HIDDEN_BEGIN_xumz
- import tinyarray
- #HIDDEN_END_xumz
- 
- # define Pauli-matrices for convenience
- #HIDDEN_BEGIN_hwbt
- sigma_0 = tinyarray.array([[1, 0], [0, 1]])
- sigma_x = tinyarray.array([[0, 1], [1, 0]])
- sigma_y = tinyarray.array([[0, -1j], [1j, 0]])
- sigma_z = tinyarray.array([[1, 0], [0, -1]])
- #HIDDEN_END_hwbt
- 
- 
- def make_system(t=1.0, alpha=0.5, e_z=0.08, W=10, L=30):
-     # Start with an empty tight-binding system and a single square lattice.
-     # `a` is the lattice constant (by default set to 1 for simplicity).
-     lat = kwant.lattice.square()
- 
-     syst = kwant.Builder()
- 
-     #### Define the scattering region. ####
- #HIDDEN_BEGIN_uxrm
-     syst[(lat(x, y) for x in range(L) for y in range(W))] = \
-         4 * t * sigma_0 + e_z * sigma_z
-     # hoppings in x-direction
-     syst[kwant.builder.HoppingKind((1, 0), lat, lat)] = \
-         -t * sigma_0 + 1j * alpha * sigma_y / 2
-     # hoppings in y-directions
-     syst[kwant.builder.HoppingKind((0, 1), lat, lat)] = \
-         -t * sigma_0 - 1j * alpha * sigma_x / 2
- #HIDDEN_END_uxrm
- 
-     #### Define the left lead. ####
-     lead = kwant.Builder(kwant.TranslationalSymmetry((-1, 0)))
- 
- #HIDDEN_BEGIN_yliu
-     lead[(lat(0, j) for j in range(W))] = 4 * t * sigma_0 + e_z * sigma_z
-     # hoppings in x-direction
-     lead[kwant.builder.HoppingKind((1, 0), lat, lat)] = \
-         -t * sigma_0 + 1j * alpha * sigma_y / 2
-     # hoppings in y-directions
-     lead[kwant.builder.HoppingKind((0, 1), lat, lat)] = \
-         -t * sigma_0 - 1j * alpha * sigma_x / 2
- #HIDDEN_END_yliu
- 
-     #### Attach the leads and return the finalized system. ####
-     syst.attach_lead(lead)
-     syst.attach_lead(lead.reversed())
- 
-     return syst
- 
- 
- def plot_conductance(syst, energies):
-     # Compute conductance
-     data = []
-     for energy in energies:
-         smatrix = kwant.smatrix(syst, energy)
-         data.append(smatrix.transmission(1, 0))
- 
--    pyplot.figure()
-+    fig = pyplot.figure()
-     pyplot.plot(energies, data)
--    pyplot.xlabel("energy [t]")
--    pyplot.ylabel("conductance [e^2/h]")
--    pyplot.show()
-+    pyplot.xlabel("energy [t]", fontsize=_defs.mpl_label_size)
-+    pyplot.ylabel("conductance [e^2/h]",
-+                 fontsize=_defs.mpl_label_size)
-+    pyplot.setp(fig.get_axes()[0].get_xticklabels(),
-+               fontsize=_defs.mpl_tick_size)
-+    pyplot.setp(fig.get_axes()[0].get_yticklabels(),
-+               fontsize=_defs.mpl_tick_size)
-+    fig.set_size_inches(_defs.mpl_width_in, _defs.mpl_width_in * 3. / 4.)
-+    fig.subplots_adjust(left=0.15, right=0.95, top=0.95, bottom=0.15)
-+    for extension in ('pdf', 'png'):
-+        fig.savefig("spin_orbit_result." + extension, dpi=_defs.dpi)
- 
- 
- def main():
-     syst = make_system()
- 
--    # Check that the system looks as intended.
--    kwant.plot(syst)
--
-     # Finalize the system.
-     syst = syst.finalized()
- 
-     # We should see non-monotonic conductance steps.
-     plot_conductance(syst, energies=[0.01 * i - 0.3 for i in range(100)])
- 
- 
- # Call the main function if the script gets executed (as opposed to imported).
- # See <http://docs.python.org/library/__main__.html>.
- if __name__ == '__main__':
-     main()

doc/source/code/figure/superconductor.py.diff

-@@ -1,132 +1,141 @@
- # Tutorial 2.6. "Superconductors": orbitals, conservation laws and symmetries
- # ===========================================================================
- #
- # Physics background
- # ------------------
- # - conductance of a NS-junction (Andreev reflection, superconducting gap)
- #
- # Kwant features highlighted
- # --------------------------
- # - Implementing electron and hole ("orbital") degrees of freedom
- #   using conservation laws.
- # - Use of discrete symmetries to relate scattering states.
- 
-+import _defs
- import kwant
- 
- import tinyarray
- import numpy as np
- 
- # For plotting
- from matplotlib import pyplot
-+from contextlib import redirect_stdout
- 
- tau_x = tinyarray.array([[0, 1], [1, 0]])
- tau_y = tinyarray.array([[0, -1j], [1j, 0]])
- tau_z = tinyarray.array([[1, 0], [0, -1]])
- 
- #HIDDEN_BEGIN_nbvn
- def make_system(a=1, W=10, L=10, barrier=1.5, barrierpos=(3, 4),
-                 mu=0.4, Delta=0.1, Deltapos=4, t=1.0):
-     # Start with an empty tight-binding system. On each site, there
-     # are now electron and hole orbitals, so we must specify the
-     # number of orbitals per site. The orbital structure is the same
-     # as in the Hamiltonian.
-     lat = kwant.lattice.square(norbs=2)
-     syst = kwant.Builder()
- 
-     #### Define the scattering region. ####
-     # The superconducting order parameter couples electron and hole orbitals
-     # on each site, and hence enters as an onsite potential.
-     # The pairing is only included beyond the point 'Deltapos' in the scattering region.
-     syst[(lat(x, y) for x in range(Deltapos) for y in range(W))] = (4 * t - mu) * tau_z
-     syst[(lat(x, y) for x in range(Deltapos, L) for y in range(W))] = (4 * t - mu) * tau_z + Delta * tau_x
- 
-     # The tunnel barrier
-     syst[(lat(x, y) for x in range(barrierpos[0], barrierpos[1])
-          for y in range(W))] = (4 * t + barrier - mu) * tau_z
- 
-     # Hoppings
-     syst[lat.neighbors()] = -t * tau_z
- #HIDDEN_END_nbvn
- #HIDDEN_BEGIN_ttth
-     #### Define the leads. ####
-     # Left lead - normal, so the order parameter is zero.
-     sym_left = kwant.TranslationalSymmetry((-a, 0))
-     # Specify the conservation law used to treat electrons and holes separately.
-     # We only do this in the left lead, where the pairing is zero.
-     lead0 = kwant.Builder(sym_left, conservation_law=-tau_z, particle_hole=tau_y)
-     lead0[(lat(0, j) for j in range(W))] = (4 * t - mu) * tau_z
-     lead0[lat.neighbors()] = -t * tau_z
- #HIDDEN_END_ttth
- #HIDDEN_BEGIN_zuuw
-     # Right lead - superconducting, so the order parameter is included.
-     sym_right = kwant.TranslationalSymmetry((a, 0))
-     lead1 = kwant.Builder(sym_right)
-     lead1[(lat(0, j) for j in range(W))] = (4 * t - mu) * tau_z + Delta * tau_x
-     lead1[lat.neighbors()] = -t * tau_z
- 
-     #### Attach the leads and return the system. ####
-     syst.attach_lead(lead0)
-     syst.attach_lead(lead1)
- 
-     return syst
- #HIDDEN_END_zuuw
- 
- #HIDDEN_BEGIN_jbjt
- def plot_conductance(syst, energies):
-     # Compute conductance
-     data = []
-     for energy in energies:
-         smatrix = kwant.smatrix(syst, energy)
-         # Conductance is N - R_ee + R_he
-         data.append(smatrix.submatrix((0, 0), (0, 0)).shape[0] -
-                     smatrix.transmission((0, 0), (0, 0)) +
-                     smatrix.transmission((0, 1), (0, 0)))
- #HIDDEN_END_jbjt
--    pyplot.figure()
-+    fig = pyplot.figure()
-     pyplot.plot(energies, data)
-     pyplot.xlabel("energy [t]")
-     pyplot.ylabel("conductance [e^2/h]")
--    pyplot.show()
-+    pyplot.setp(fig.get_axes()[0].get_xticklabels(),
-+                fontsize=_defs.mpl_tick_size)
-+    pyplot.setp(fig.get_axes()[0].get_yticklabels(),
-+                fontsize=_defs.mpl_tick_size)
-+    fig.set_size_inches(_defs.mpl_width_in, _defs.mpl_width_in * 3. / 4.)
-+    fig.subplots_adjust(left=0.15, right=0.95, top=0.95, bottom=0.15)
-+    for extension in ('pdf', 'png'):
-+        fig.savefig("superconductor_transport_result." + extension,
-+                    dpi=_defs.dpi)
- 
- #HIDDEN_BEGIN_pqmp
- def check_PHS(syst):
-     # Scattering matrix
-     s = kwant.smatrix(syst, energy=0)
-     # Electron to electron block
-     s_ee = s.submatrix((0,0), (0,0))
-     # Hole to hole block
-     s_hh = s.submatrix((0,1), (0,1))
-     print('s_ee: \n', np.round(s_ee, 3))
-     print('s_hh: \n', np.round(s_hh[::-1, ::-1], 3))
-     print('s_ee - s_hh^*: \n',
-           np.round(s_ee - s_hh[::-1, ::-1].conj(), 3), '\n')
-     # Electron to hole block
-     s_he = s.submatrix((0,1), (0,0))
-     # Hole to electron block
-     s_eh = s.submatrix((0,0), (0,1))
-     print('s_he: \n', np.round(s_he, 3))
-     print('s_eh: \n', np.round(s_eh[::-1, ::-1], 3))
-     print('s_he + s_eh^*: \n',
-           np.round(s_he + s_eh[::-1, ::-1].conj(), 3))
- #HIDDEN_END_pqmp
- 
- def main():
-     syst = make_system(W=10)
- 
--    # Check that the system looks as intended.
--    kwant.plot(syst)
--
-     # Finalize the system.
-     syst = syst.finalized()
- 
-     # Check particle-hole symmetry of the scattering matrix
--    check_PHS(syst)
-+    with open('check_PHS_out.txt', 'w') as f:
-+        with redirect_stdout(f):
-+            check_PHS(syst)
- 
-     # Compute and plot the conductance
-     plot_conductance(syst, energies=[0.002 * i for i in range(-10, 100)])
- 
- 
- # Call the main function if the script gets executed (as opposed to imported).
- # See <http://docs.python.org/library/__main__.html>.
- if __name__ == '__main__':
-     main()

doc/source/code/include/quantum_wire_revisited.py

-# Tutorial 2.2.3. Building the same system with less code
-# =======================================================
-#
-# Physics background
-# ------------------
-#  Conductance of a quantum wire; subbands
-#
-# Kwant features highlighted
-# --------------------------
-#  - Using iterables and builder.HoppingKind for making systems
-#  - introducing `reversed()` for the leads
-#
-# Note: Does the same as tutorial1a.py, but using other features of Kwant.
-
-#HIDDEN_BEGIN_xkzy
-import kwant
-
-# For plotting
-from matplotlib import pyplot
-
-
-def make_system(a=1, t=1.0, W=10, L=30):
-    # Start with an empty tight-binding system and a single square lattice.
-    # `a` is the lattice constant (by default set to 1 for simplicity.
-    lat = kwant.lattice.square(a)
-
-    syst = kwant.Builder()
-#HIDDEN_END_xkzy
-
-    #### Define the scattering region. ####
-#HIDDEN_BEGIN_vvjt
-    syst[(lat(x, y) for x in range(L) for y in range(W))] = 4 * t
-#HIDDEN_END_vvjt
-#HIDDEN_BEGIN_nooi
-    syst[lat.neighbors()] = -t
-#HIDDEN_END_nooi
-
-    #### Define and attach the leads. ####
-    # Construct the left lead.
-#HIDDEN_BEGIN_iepx
-    lead = kwant.Builder(kwant.TranslationalSymmetry((-a, 0)))
-    lead[(lat(0, j) for j in range(W))] = 4 * t
-    lead[lat.neighbors()] = -t
-#HIDDEN_END_iepx
-
-    # Attach the left lead and its reversed copy.
-#HIDDEN_BEGIN_yxot
-    syst.attach_lead(lead)
-    syst.attach_lead(lead.reversed())
-
-    return syst
-#HIDDEN_END_yxot
-
-
-#HIDDEN_BEGIN_ayuk
-def plot_conductance(syst, energies):
-    # Compute conductance
-    data = []
-    for energy in energies:
-        smatrix = kwant.smatrix(syst, energy)
-        data.append(smatrix.transmission(1, 0))
-
-    pyplot.figure()
-    pyplot.plot(energies, data)
-    pyplot.xlabel("energy [t]")
-    pyplot.ylabel("conductance [e^2/h]")
-    pyplot.show()
-#HIDDEN_END_ayuk
-
-
-#HIDDEN_BEGIN_cjel
-def main():
-    syst = make_system()
-
-    # Check that the system looks as intended.
-    kwant.plot(syst)
-
-    # Finalize the system.
-    syst = syst.finalized()
-
-    # We should see conductance steps.
-    plot_conductance(syst, energies=[0.01 * i for i in range(100)])
-#HIDDEN_END_cjel
-
-
-# Call the main function if the script gets executed (as opposed to imported).
-# See <http://docs.python.org/library/__main__.html>.
-#HIDDEN_BEGIN_ypbj
-if __name__ == '__main__':
-    main()
-#HIDDEN_END_ypbj

doc/source/conf.py

@@ -15,8 +15,15 @@ import sys
 import os
 import string
 from distutils.util import get_platform
-sys.path.insert(0, "../../build/lib.{0}-{1}.{2}".format(
-        get_platform(), *sys.version_info[:2]))
+
+package_path = os.path.abspath(
+    "../../build/lib.{0}-{1}.{2}"
+    .format(get_platform(), *sys.version_info[:2]))
+
+# Insert into sys.path so that we can import kwant here
+sys.path.insert(0, package_path)
+# Insert into PYTHONPATH so that jupyter-sphinx will pick it up
+os.environ['PYTHONPATH'] = ':'.join((package_path, os.environ.get('PYTHONPATH','')))
 
 import kwant
 import kwant.qsymm
@@ -30,8 +37,9 @@ import kwant.continuum  # sphinx gets confused with lazy loading
 sys.path.insert(0, os.path.abspath('../sphinxext'))
 
 extensions = ['sphinx.ext.autodoc', 'sphinx.ext.autosummary',
-              'sphinx.ext.todo', 'sphinx.ext.mathjax', 'numpydoc',
-              'kwantdoc', 'sphinx.ext.linkcode']
+              'sphinx.ext.todo', 'sphinx.ext.mathjax', 'sphinx.ext.napoleon',
+              'sphinx.ext.linkcode', 'jupyter_sphinx', 'sphinx_togglebutton',
+              'sphinxcontrib.rsvgconverter']
 
 # Add any paths that contain templates here, relative to this directory.
 templates_path = ['../templates']
@@ -99,14 +107,20 @@ pygments_style = 'sphinx'
 # Do not show all class members automatically in the class documentation
 numpydoc_show_class_members = False
 
+# Jupyter Sphinx config
+
+jupyter_sphinx_thebelab_config = {
+    "binderOptions": {
+        "repo": "kwant-project/binder",
+        "ref": "master",
+    }
+}
+
 # -- Options for HTML output ---------------------------------------------------
 
 # The theme to use for HTML and HTML Help pages.  Major themes that come with
 # Sphinx are currently 'default' and 'sphinxdoc'.
-html_theme = 'kwantdoctheme'
-html_theme_path = ['..']
-html_theme_options = {'collapsiblesidebar': True}
-html_style = 'kwant.css'
+html_theme = 'sphinx_book_theme'
 
 # Theme options are theme-specific and customize the look and feel of a theme
 # further.  For a list of options available for each theme, see the
@@ -125,7 +139,7 @@ html_style = 'kwant.css'
 
 # The name of an image file (relative to this directory) to place at the top
 # of the sidebar.
-#html_logo = None
+html_logo = "_static/kwant_logo.png"
 
 # The name of an image file (within the static path) to use as favicon of the
 # docs.  This file should be a Windows icon file (.ico) being 16x16 or 32x32
@@ -175,13 +189,9 @@ html_domain_indices = False
 # If nonempty, this is the file name suffix for HTML files (e.g. ".xhtml").
 #html_file_suffix = ''
 
-# Output file base name for HTML help builder.
-htmlhelp_basename = 'kwantdoc'
-
 
 # -- Options for LaTeX output --------------------------------------------------
 
-# http://thread.gmane.org/gmane.comp.python.sphinx.devel/4220/focus=4238
 latex_elements = {'papersize': 'a4paper',
                   'release': '',
                   'releasename': '',
@@ -218,7 +228,7 @@ latex_documents = [
 
 latex_engine = 'xelatex'
 
-latex_use_xindy = False
+latex_use_xindy = False  # Xindy not installable in CI environment
 
 # The name of an image file (relative to this directory) to place at the top of
 # the title page.
@@ -239,8 +249,9 @@ latex_domain_indices = False
 autosummary_generate = True
 
 autoclass_content = "both"
-autodoc_default_flags = ['show-inheritance']
-
+autodoc_default_options = {
+    'show-inheritance': True,
+}
 
 # -- Teach Sphinx to document bound methods like functions ---------------------
 import types
@@ -253,7 +264,7 @@ class BoundMethodDocumenter(autodoc.FunctionDocumenter):
     @classmethod
     def can_document_member(cls, member, membername, isattr, parent):
         # Return True iff `member` is a bound method.  Taken from
-        # <http://stackoverflow.com/a/1260881>.
+        # <https://stackoverflow.com/a/1260881>.
         return (isinstance(member, types.MethodType) and
                 member.__self__ is not None and
                 not issubclass(member.__self__.__class__, type) and
@@ -286,10 +297,6 @@ nitpick_ignore = [('py:class', 'Warning'), ('py:class', 'Exception'),
                   ('py:class', 'kwant.builder._FinalizedBuilderMixin')]
 
 
-# Use custom MathJax CDN, as cdn.mathjax.org will soon shut down
-mathjax_path = "https://cdnjs.cloudflare.com/ajax/libs/mathjax/2.7.0/MathJax.js?config=TeX-AMS-MML_HTMLorMML"
-
-
 # -- Make Sphinx insert source code links --------------------------------------
 def linkcode_resolve(domain, info):
     def find_source():

doc/source/code/figure/ab_ring_sketch.svg → doc/source/figure/ab_ring_sketch.svg

@@ -3,7 +3,7 @@
 
 <svg
    xmlns:dc="http://purl.org/dc/elements/1.1/"
-   xmlns:cc="http://creativecommons.org/ns#"
+   xmlns:cc="https://creativecommons.org/ns#"
    xmlns:rdf="http://www.w3.org/1999/02/22-rdf-syntax-ns#"
    xmlns:svg="http://www.w3.org/2000/svg"
    xmlns="http://www.w3.org/2000/svg"

doc/source/code/figure/ab_ring_sketch2.svg → doc/source/figure/ab_ring_sketch2.svg

@@ -3,7 +3,7 @@
 
 <svg
    xmlns:dc="http://purl.org/dc/elements/1.1/"
-   xmlns:cc="http://creativecommons.org/ns#"
+   xmlns:cc="https://creativecommons.org/ns#"
    xmlns:rdf="http://www.w3.org/1999/02/22-rdf-syntax-ns#"
    xmlns:svg="http://www.w3.org/2000/svg"
    xmlns="http://www.w3.org/2000/svg"

doc/source/code/figure/superconductor_transport_sketch.svg → doc/source/figure/superconductor_transport_sketch.svg

@@ -3,7 +3,7 @@
 
 <svg
    xmlns:dc="http://purl.org/dc/elements/1.1/"
-   xmlns:cc="http://creativecommons.org/ns#"
+   xmlns:cc="https://creativecommons.org/ns#"
    xmlns:rdf="http://www.w3.org/1999/02/22-rdf-syntax-ns#"
    xmlns:svg="http://www.w3.org/2000/svg"
    xmlns="http://www.w3.org/2000/svg"

doc/source/pre/whatsnew/0.2.rst

@@ -21,7 +21,7 @@ New MUMPS-based solver
 ----------------------
 The code for sparse matrix solvers has been reorganized and a new solver has
 been added next to `kwant.solvers.sparse`: `kwant.solvers.mumps`.  The new
-solver uses the `MUMPS <http://graal.ens-lyon.fr/MUMPS/>`_ software package and
+solver uses the `MUMPS <https://graal.ens-lyon.fr/MUMPS/>`_ software package and
 is much (typically several times) faster than the UMFPACK-based old solver.
 In addition, MUMPS uses considerably less memory for a given system while at
 the same time it is able to take advantage of more than 2 GiB of RAM.

doc/source/pre/whatsnew/1.0.rst

@@ -3,8 +3,8 @@ What's new in Kwant 1.0
 
 This article explains the new features in Kwant 1.0 compared to Kwant 0.2.
 Kwant 1.0 was released on 9 September 2013.  Please consult the `full list of
-changes in Kwant <https://git.kwant-project.org/kwant/log/?h=v1.0.5>`_ for all
-the changes up to the most recent bugfix release.
+changes in Kwant <https://gitlab.kwant-project.org/kwant/kwant/-/commits/v1.0.5>`_
+for all the changes up to the most recent bugfix release.
 
 
 Lattice and shape improvements

doc/source/pre/whatsnew/1.1.rst

@@ -4,7 +4,7 @@ What's new in Kwant 1.1
 This article explains the user-visible changes in Kwant 1.1.0, released on 21
 October 2015.  See also the `full list of changes up to the most recent bugfix
 release of the 1.1 series
-<https://gitlab.kwant-project.org/kwant/kwant/compare/v1.1.0...latest-1.1>`_.
+<https://gitlab.kwant-project.org/kwant/kwant/-/compare/v1.1.0...latest-1.1>`_.
 
 Harmonize `~kwant.physics.Bands` with `~kwant.physics.modes`
 ------------------------------------------------------------

doc/source/pre/whatsnew/1.2.rst

@@ -4,7 +4,7 @@ What's new in Kwant 1.2
 This article explains the user-visible changes in Kwant 1.2.2, released on 9
 December 2015.  See also the `full list of changes up to the most recent bugfix
 release of the 1.2 series
-<https://gitlab.kwant-project.org/kwant/kwant/compare/v1.2.2...latest-1.2>`_.
+<https://gitlab.kwant-project.org/kwant/kwant/-/compare/v1.2.2...latest-1.2>`_.
 
 Kwant 1.2 is identical to Kwant 1.1 except that it has been updated to run on
 Python 3.4 and above.  Bugfix releases for the 1.1 and 1.2 series will mirror