convert graphene example to jupyter-sphinx

doc/source/tutorial/graphene.rst

@@ -5,7 +5,7 @@ Beyond square lattices: graphene
 
 .. seealso::
     The complete source code of this example can be found in
-    :download:`graphene.py </code/download/graphene.py>`
+    :jupyter-download:script:`graphene`
 
 In the following example, we are going to calculate the
 conductance through a graphene quantum dot with a p-n junction
@@ -18,9 +18,40 @@ We begin by defining the honeycomb lattice of graphene. This is
 in principle already done in `kwant.lattice.honeycomb`, but we do it
 explicitly here to show how to define a new lattice:
 
-.. literalinclude:: /code/include/graphene.py
-    :start-after: #HIDDEN_BEGIN_hnla
-    :end-before: #HIDDEN_END_hnla
+.. jupyter-kernel::
+    :id: graphene
+
+.. jupyter-execute::
+    :hide-code:
+
+    # Tutorial 2.5. Beyond square lattices: graphene
+    # ==============================================
+    #
+    # Physics background
+    # ------------------
+    #  Transport through a graphene quantum dot with a pn-junction
+    #
+    # Kwant features highlighted
+    # --------------------------
+    #  - Application of all the aspects of tutorials 1-3 to a more complicated
+    #    lattice, namely graphene
+
+    from math import pi, sqrt, tanh
+
+    from matplotlib import pyplot
+
+    import kwant
+
+    # For computing eigenvalues
+    import scipy.sparse.linalg as sla
+
+    sin_30, cos_30 = (1 / 2, sqrt(3) / 2)
+
+.. jupyter-execute::
+
+    graphene = kwant.lattice.general([(1, 0), (sin_30, cos_30)],
+                                     [(0, 0), (0, 1 / sqrt(3))])
+    a, b = graphene.sublattices
 
 The first argument to the `~kwant.lattice.general` function is the list of
 primitive vectors of the lattice; the second one is the coordinates of basis
@@ -31,9 +62,30 @@ itself forms a regular lattice of the same type as well, and those
 In the next step we define the shape of the scattering region (circle again)
 and add all lattice points using the ``shape``-functionality:
 
-.. literalinclude:: /code/include/graphene.py
-    :start-after: #HIDDEN_BEGIN_shzy
-    :end-before: #HIDDEN_END_shzy
+.. jupyter-execute::
+    :hide-code:
+
+    r = 10
+    w = 2.0
+    pot = 0.1
+
+.. jupyter-execute::
+
+    #### Define the scattering region. ####
+    # circular scattering region
+    def circle(pos):
+        x, y = pos
+        return x ** 2 + y ** 2 < r ** 2
+
+    syst = kwant.Builder()
+
+    # w: width and pot: potential maximum of the p-n junction
+    def potential(site):
+        (x, y) = site.pos
+        d = y * cos_30 + x * sin_30
+        return pot * tanh(d / w)
+
+    syst[graphene.shape(circle, (0, 0))] = potential
 
 As you can see, this works exactly the same for any kind of lattice.
 We add the onsite energies using a function describing the p-n junction;
@@ -45,9 +97,11 @@ As a next step we add the hoppings, making use of
 `~kwant.builder.HoppingKind`. For illustration purposes we define
 the hoppings ourselves instead of using ``graphene.neighbors()``:
 
-.. literalinclude:: /code/include/graphene.py
-    :start-after: #HIDDEN_BEGIN_hsmc
-    :end-before: #HIDDEN_END_hsmc
+.. jupyter-execute::
+
+    # specify the hoppings of the graphene lattice in the
+    # format expected by builder.HoppingKind
+    hoppings = (((0, 0), a, b), ((0, 1), a, b), ((-1, 1), a, b))
 
 The nearest-neighbor model for graphene contains only
 hoppings between different basis atoms. For this type of
@@ -63,27 +117,50 @@ respect to the two primitive vectors ``[(1, 0), (sin_30, cos_30)]``.
 
 Adding the hoppings however still works the same way:
 
-.. literalinclude:: /code/include/graphene.py
-    :start-after: #HIDDEN_BEGIN_bfwb
-    :end-before: #HIDDEN_END_bfwb
+.. jupyter-execute::
+
+    syst[[kwant.builder.HoppingKind(*hopping) for hopping in hoppings]] = -1
 
 Modifying the scattering region is also possible as before. Let's
 do something crazy, and remove an atom in sublattice A
 (which removes also the hoppings from/to this site) as well
 as add an additional link:
 
-.. literalinclude:: /code/include/graphene.py
-    :start-after: #HIDDEN_BEGIN_efut
-    :end-before: #HIDDEN_END_efut
+.. jupyter-execute::
+
+    # Modify the scattering region
+    del syst[a(0, 0)]
+    syst[a(-2, 1), b(2, 2)] = -1
 
 Note again that the conversion from a tuple `(i,j)` to site
 is done by the sublattices `a` and `b`.
 
 The leads are defined almost as before:
 
-.. literalinclude:: /code/include/graphene.py
-    :start-after: #HIDDEN_BEGIN_aakh
-    :end-before: #HIDDEN_END_aakh
+.. jupyter-execute::
+
+    #### Define the leads. ####
+    # left lead
+    sym0 = kwant.TranslationalSymmetry(graphene.vec((-1, 0)))
+
+    def lead0_shape(pos):
+        x, y = pos
+        return (-0.4 * r < y < 0.4 * r)
+
+    lead0 = kwant.Builder(sym0)
+    lead0[graphene.shape(lead0_shape, (0, 0))] = -pot
+    lead0[[kwant.builder.HoppingKind(*hopping) for hopping in hoppings]] = -1
+
+    # The second lead, going to the top right
+    sym1 = kwant.TranslationalSymmetry(graphene.vec((0, 1)))
+
+    def lead1_shape(pos):
+        v = pos[1] * sin_30 - pos[0] * cos_30
+        return (-0.4 * r < v < 0.4 * r)
+
+    lead1 = kwant.Builder(sym1)
+    lead1[graphene.shape(lead1_shape, (0, 0))] = pot
+    lead1[[kwant.builder.HoppingKind(*hopping) for hopping in hoppings]] = -1
 
 Note the method `~kwant.lattice.Polyatomic.vec` used in calculating the
 parameter for `~kwant.lattice.TranslationalSymmetry`.  The latter expects a
@@ -98,28 +175,63 @@ in a square lattice -- they follow the non-orthogonal primitive vectors defined
 in the beginning.
 
 Later, we will compute some eigenvalues of the closed scattering region without
-leads. This is why we postpone attaching the leads to the system. Instead,
-we return the scattering region and the leads separately.
+leads. This is why we postpone attaching the leads to the system.
 
-.. literalinclude:: /code/include/graphene.py
-    :start-after: #HIDDEN_BEGIN_kmmw
-    :end-before: #HIDDEN_END_kmmw
 
 The computation of some eigenvalues of the closed system is done
 in the following piece of code:
 
-.. literalinclude:: /code/include/graphene.py
-    :start-after: #HIDDEN_BEGIN_zydk
-    :end-before: #HIDDEN_END_zydk
+.. jupyter-execute::
+
+    def compute_evs(syst):
+        # Compute some eigenvalues of the closed system
+        sparse_mat = syst.hamiltonian_submatrix(sparse=True)
+
+        evs = sla.eigs(sparse_mat, 2)[0]
+        print(evs.real)
 
 The code for computing the band structure and the conductance is identical
 to the previous examples, and needs not be further explained here.
 
-Finally, in the `main` function we make and plot the system:
+Finally, we plot the system:
+
+.. jupyter-execute::
+    :hide-code:
+
+    def plot_conductance(syst, energies):
+        # Compute transmission as a function of energy
+        data = []
+        for energy in energies:
+            smatrix = kwant.smatrix(syst, energy)
+            data.append(smatrix.transmission(0, 1))
+
+        pyplot.figure()
+        pyplot.plot(energies, data)
+        pyplot.xlabel("energy [t]")
+        pyplot.ylabel("conductance [e^2/h]")
+        pyplot.show()
+
 
-.. literalinclude:: /code/include/graphene.py
-    :start-after: #HIDDEN_BEGIN_itkk
-    :end-before: #HIDDEN_END_itkk
+    def plot_bandstructure(flead, momenta):
+        bands = kwant.physics.Bands(flead)
+        energies = [bands(k) for k in momenta]
+
+        pyplot.figure()
+        pyplot.plot(momenta, energies)
+        pyplot.xlabel("momentum [(lattice constant)^-1]")
+        pyplot.ylabel("energy [t]")
+        pyplot.show()
+
+
+.. jupyter-execute::
+
+    # To highlight the two sublattices of graphene, we plot one with
+    # a filled, and the other one with an open circle:
+    def family_colors(site):
+        return 0 if site.family == a else 1
+
+    # Plot the closed system without leads.
+    kwant.plot(syst, site_color=family_colors, site_lw=0.1, colorbar=False);
 
 We customize the plotting: we set the `site_colors` argument of
 `~kwant.plotter.plot` to a function which returns 0 for
@@ -132,27 +244,43 @@ The function `~kwant.plotter.plot` shows these values using a color scale
 (grayscale by default). The symbol `size` is specified in points, and is
 independent on the overall figure size.
 
-Plotting the closed system gives this result:
-
-.. image:: /code/figure/graphene_syst1.*
 
 Computing the eigenvalues of largest magnitude,
 
-.. literalinclude:: /code/include/graphene.py
-    :start-after: #HIDDEN_BEGIN_jmbi
-    :end-before: #HIDDEN_END_jmbi
+.. jupyter-execute::
+
+    compute_evs(syst.finalized())
 
-should yield two eigenvalues equal to ``[ 3.07869311,
+yields two eigenvalues equal to ``[ 3.07869311,
 -3.06233144]``.
 
-The remaining code of `main` attaches the leads to the system and plots it
+The remaining code attaches the leads to the system and plots it
 again:
 
-.. image:: /code/figure/graphene_syst2.*
+.. jupyter-execute::
+
+    # Attach the leads to the system.
+    for lead in [lead0, lead1]:
+        syst.attach_lead(lead)
+
+    # Then, plot the system with leads.
+    kwant.plot(syst, site_color=family_colors, site_lw=0.1,
+               lead_site_lw=0, colorbar=False);
+
+We then finalize the system:
 
-It computes the band structure of one of lead 0:
+.. jupyter-execute::
 
-.. image:: /code/figure/graphene_bs.*
+    syst = syst.finalized()
+
+and compute the band structure of one of lead 0:
+
+.. jupyter-execute::
+
+
+    # Compute the band structure of lead 0.
+    momenta = [-pi + 0.02 * pi * i for i in range(101)]
+    plot_bandstructure(syst.leads[0], momenta)
 
 showing all the features of a zigzag lead, including the flat
 edge state bands (note that the band structure is not symmetric around
@@ -160,7 +288,11 @@ zero energy, due to a potential in the leads).
 
 Finally the transmission through the system is computed,
 
-.. image:: /code/figure/graphene_result.*
+.. jupyter-execute::
+
+    # Plot conductance.
+    energies = [-2 * pot + 4. / 50. * pot * i for i in range(51)]
+    plot_conductance(syst, energies)
 
 showing the typical resonance-like transmission probability through
 an open quantum dot