From 72064255c6de0796d3c556328650a37abfecfc51 Mon Sep 17 00:00:00 2001
From: Joseph Weston <joseph.weston08@gmail.com>
Date: Tue, 21 Mar 2017 14:44:35 +0100
Subject: [PATCH] add operator tutorial and tutorial script
---
doc/source/images/magnetic_texture.py.diff | 143 ++++++++
doc/source/tutorial/magnetic_texture.py | 245 ++++++++++++++
doc/source/tutorial/tutorial6.rst | 367 +++++++++++++++++++++
3 files changed, 755 insertions(+)
create mode 100644 doc/source/images/magnetic_texture.py.diff
create mode 100644 doc/source/tutorial/magnetic_texture.py
create mode 100644 doc/source/tutorial/tutorial6.rst
diff --git a/doc/source/images/magnetic_texture.py.diff b/doc/source/images/magnetic_texture.py.diff
new file mode 100644
index 00000000..bc7acb4d
--- /dev/null
+++ b/doc/source/images/magnetic_texture.py.diff
@@ -0,0 +1,143 @@
+--- original
++++ modified
+@@ -11,6 +11,9 @@
+ # - operators
+ # - plotting vector fields
+
++import _defs
++from contextlib import redirect_stdout
++
+ from math import sin, cos, tanh, pi
+ import itertools
+ import numpy as np
+@@ -79,7 +82,13 @@
+ return 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))
+@@ -88,28 +97,38 @@
+ m_i = np.reshape(m_i, x.shape + (3,))
+ m_i = np.rollaxis(m_i, 2, 0)
+
+- fig, ax = plt.subplots(1, 1)
++ fig, ax = plt.subplots(1, 1, figsize=(9, 7))
+ im = ax.quiver(x, y, *m_i, pivot='mid', scale=75)
+ fig.colorbar(im)
+- plt.show()
++ if fname:
++ save_plot(fname)
++ else:
++ plt.show()
+
+
+-def plot_densities(syst, densities):
+- fig, axes = plt.subplots(1, len(densities))
++def plot_densities(syst, densities, fname=None):
++ fig, axes = plt.subplots(1, len(densities), figsize=(7*len(densities), 7))
+ for ax, (title, rho) in zip(axes, densities):
+ kwant.plotter.map(syst, rho, ax=ax, a=4)
+ ax.set_title(title)
+- plt.show()
++ if fname:
++ save_plot(fname)
++ else:
++ plt.show()
++
+
+
+-def plot_currents(syst, currents):
+- fig, axes = plt.subplots(1, len(currents))
++def plot_currents(syst, currents, fname=None):
++ fig, axes = plt.subplots(1, len(currents), figsize=(7*len(currents), 7))
+ if not hasattr(axes, '__len__'):
+ axes = (axes,)
+ for ax, (title, current) in zip(axes, currents):
+ kwant.plotter.current(syst, current, ax=ax, colorbar=False)
+ ax.set_title(title)
+- plt.show()
++ if fname:
++ save_plot(fname)
++ else:
++ plt.show()
+
+
+ def main():
+@@ -119,7 +138,7 @@
+ wf = kwant.wave_function(syst, energy=-1, params=params)
+ psi = wf(0)[0]
+
+- plot_vector_field(syst, dict(r0=20, delta=10))
++ plot_vector_field(syst, dict(r0=20, delta=10), fname='mag_field_direction')
+
+ # even (odd) indices correspond to spin up (down)
+ up, down = psi[::2], psi[1::2]
+@@ -144,7 +163,7 @@
+ ('$σ_0$', density),
+ ('$σ_z$', spin_z),
+ ('$σ_y$', spin_y),
+- ])
++ ], fname='spin_densities')
+
+ J_0 = kwant.operator.Current(syst)
+ J_z = kwant.operator.Current(syst, sigma_z)
+@@ -159,7 +178,7 @@
+ ('$J_{σ_0}$', current),
+ ('$J_{σ_z}$', spin_z_current),
+ ('$J_{σ_y}$', spin_y_current),
+- ])
++ ], fname='spin_currents')
+
+ def following_m_i(site, r0, delta):
+ m_i = field_direction(site.pos, r0, delta)
+@@ -173,7 +192,7 @@
+ plot_currents(syst, [
+ (r'$J_{\mathbf{m}_i}$', m_current),
+ ('$J_{σ_z}$', spin_z_current),
+- ])
++ ], fname='spin_current_comparison')
+
+
+ def circle(site):
+@@ -183,7 +202,9 @@
+
+ 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)
+
+ def left_cut(site_to, site_from):
+ return site_from.pos[0] <= -39 and site_to.pos[0] > -39
+@@ -197,8 +218,10 @@
+ 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))
+
+ J_m = kwant.operator.Current(syst, following_m_i)
+ J_z = kwant.operator.Current(syst, sigma_z)
+@@ -214,7 +237,7 @@
+ 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__':
diff --git a/doc/source/tutorial/magnetic_texture.py b/doc/source/tutorial/magnetic_texture.py
new file mode 100644
index 00000000..89c7acfb
--- /dev/null
+++ b/doc/source/tutorial/magnetic_texture.py
@@ -0,0 +1,245 @@
+# 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/tutorial6.rst b/doc/source/tutorial/tutorial6.rst
new file mode 100644
index 00000000..c2c314ac
--- /dev/null
+++ b/doc/source/tutorial/tutorial6.rst
@@ -0,0 +1,367 @@
+Computing local quantities: densities and currents
+==================================================
+In the previous tutorials we have mainly concentrated on calculating *global*
+properties such as conductance and band structures. Often, however, insight can
+be gained from calculating *locally-defined* quantities, that is, quantities
+defined over individual sites or hoppings in your system. In the
+:ref:`closed-systems` tutorial we saw how we could visualize the density
+associated with the eigenstates of a system using `kwant.plotter.map`.
+
+In this tutorial we will see how we can calculate more general quantities than
+simple densities by studying spin transport in a system with a magnetic
+texture.
+
+.. seealso::
+ The complete source code of this example can be found in
+ :download:`tutorial/magnetic_texture.py <../../../tutorial/magnetic_texture.py>`
+
+
+Introduction
+------------
+Our starting point will be the following spinful tight-binding model on
+a square lattice:
+
+.. math::
+ H = - \sum_{⟨ij⟩}\sum_{α} |iα⟩⟨jα|
+ + J \sum_{i}\sum_{αβ} \mathbf{m}_i⋅ \mathbf{σ}_{αβ} |iα⟩⟨iβ|,
+
+where latin indices run over sites, and greek indices run over spin. We can
+identify the first term as a nearest-neighbor hopping between like-spins, and
+the second as a term that couples spins on the same site. The second term acts
+like a magnetic field of strength :math:`J` that varies from site to site and
+that, on site :math:`i`, points in the direction of the unit vector
+:math:`\mathbf{m}_i`. :math:`\mathbf{σ}_{αβ}` is a vector of Pauli matrices.
+We shall take the following form for :math:`\mathbf{m}_i`:
+
+.. math::
+ \mathbf{m}_i &=&\ \left(
+ \frac{x_i}{x_i^2 + y_i^2} \sin θ_i,\
+ \frac{y_i}{x_i^2 + y_i^2} \sin θ_i,\
+ \cos θ_i \right)^T,
+ \\
+ θ_i &=&\ \frac{π}{2} \tanh \frac{r_i - r_0}{δ},
+
+where :math:`x_i` and :math:`y_i` are the :math:`x` and :math:`y` coordinates
+of site :math:`i`, and :math:`r_i = \sqrt{x_i^2 + y_i^2}`.
+
+To define this model in Kwant we start as usual by defining functions that
+depend on the model parameters:
+
+.. literalinclude:: magnetic_texture.py
+ :start-after: #HIDDEN_BEGIN_model
+ :end-before: #HIDDEN_END_model
+
+and define our system as a square shape on a square lattice with two orbitals
+per site, with leads attached on the left and right:
+
+.. literalinclude:: magnetic_texture.py
+ :start-after: #HIDDEN_BEGIN_syst
+ :end-before: #HIDDEN_END_syst
+
+Below is a plot of a projection of :math:`\mathbf{m}_i` onto the x-y plane
+inside the scattering region. The z component is shown by the color scale:
+
+.. image:: ../images/mag_field_direction.*
+
+We will now be interested in analyzing the form of the scattering states
+that originate from the left lead:
+
+.. literalinclude:: magnetic_texture.py
+ :start-after: #HIDDEN_BEGIN_wavefunction
+ :end-before: #HIDDEN_END_wavefunction
+
+Local densities
+---------------
+If we were simulating a spinless system with only a single degree of freedom,
+then calculating the density on each site would be as simple as calculating the
+absolute square of the wavefunction like::
+
+ density = np.abs(psi)**2
+
+When there are multiple degrees of freedom per site, however, one has to be
+more careful. In the present case with two (spin) degrees of freedom per site
+one could calculate the per-site density like:
+
+.. literalinclude:: magnetic_texture.py
+ :start-after: #HIDDEN_BEGIN_ldos
+ :end-before: #HIDDEN_END_ldos
+
+With more than one degree of freedom per site we have more freedom as to what
+local quantities we can meaningfully compute. For example, we may wish to
+calculate the local z-projected spin density. We could calculate
+this in the following way:
+
+.. literalinclude:: magnetic_texture.py
+ :start-after: #HIDDEN_BEGIN_lsdz
+ :end-before: #HIDDEN_END_lsdz
+
+If we wanted instead to calculate the local y-projected spin density, we would
+need to use an even more complicated expression:
+
+.. literalinclude:: magnetic_texture.py
+ :start-after: #HIDDEN_BEGIN_lsdy
+ :end-before: #HIDDEN_END_lsdy
+
+The `kwant.operator` module aims to alleviate somewhat this tedious
+book-keeping by providing a simple interface for defining operators that act on
+wavefunctions. To calculate the above quantities we would use the
+`~kwant.operator.Density` operator like so:
+
+.. literalinclude:: magnetic_texture.py
+ :start-after: #HIDDEN_BEGIN_lden
+ :end-before: #HIDDEN_END_lden
+
+`~kwant.operator.Density` takes a `~kwant.system.System` as its first parameter
+as well as (optionally) a square matrix that defines the quantity that you wish
+to calculate per site. When an instance of a `~kwant.operator.Density` is then
+evaluated with a wavefunction, the quantity
+
+.. math:: Ï_i = \mathbf{ψ}^†_i \mathbf{M} \mathbf{ψ}_i
+
+is calculated for each site :math:`i`, where :math:`\mathbf{ψ}_{i}` is a vector
+consisting of the wavefunction components on that site and :math:`\mathbf{M}`
+is the square matrix referred to previously.
+
+Below we can see colorplots of the above-calculated quantities. The array that
+is returned by evaluating a `~kwant.operator.Density` can be used directly with
+`kwant.plotter.map`:
+
+.. image:: ../images/spin_densities.*
+
+
+.. specialnote:: Technical Details
+
+ Although we refer loosely to "densities" and "operators" above, a
+ `~kwant.operator.Density` actually represents a *collection* of linear
+ operators. This can be made clear by rewriting the above definition
+ of :math:`Ï_i` in the following way:
+
+ .. math::
+ Ï_i = \sum_{αβ} ψ^*_{α} \mathcal{M}_{iαβ} ψ_{β}
+
+ where greek indices run over the degrees of freedom in the Hilbert space of
+ the scattering region and latin indices run over sites. We can this
+ identify :math:`\mathcal{M}_{iαβ}` as the components of a rank-3 tensor and can
+ represent them as a "vector of matrices":
+
+ .. math::
+ \mathcal{M} = \left[
+ \left(\begin{matrix}
+ \mathbf{M} & 0 & … \\
+ 0 & 0 & … \\
+ ⋮ & ⋮ & ⋱
+ \end{matrix}\right)
+ ,\
+ \left(\begin{matrix}
+ 0 & 0 & … \\
+ 0 & \mathbf{M} & … \\
+ ⋮ & ⋮ & ⋱
+ \end{matrix}\right)
+ , … \right]
+
+ where :math:`\mathbf{M}` is defined as in the main text, and the :math:`0`
+ are zero matrices of the same shape as :math:`\mathbf{M}`.
+
+
+Local currents
+--------------
+`kwant.operator` also has a class `~kwant.operator.Current` for calculating
+local currents, analogously to the local "densities" described above. If
+one has defined a density via a matrix :math:`\mathbf{M}` and the above
+equation, then one can define a local current flowing from site :math:`b`
+to site :math:`a`:
+
+.. math:: J_{ab} = i \left(
+ \mathbf{ψ}^†_b (\mathbf{H}_{ab})^†\mathbf{M} \mathbf{ψ}_a
+ - \mathbf{ψ}^†_a \mathbf{M} \mathbf{H}_{ab} \mathbf{ψ}_b
+ \right),
+
+where :math:`\mathbf{H}_{ab}` is the hopping matrix from site :math:`b` to site
+:math:`a`. For example, to calculate the local current and
+spin current:
+
+.. literalinclude:: magnetic_texture.py
+ :start-after: #HIDDEN_BEGIN_current
+ :end-before: #HIDDEN_END_current
+
+Evaluating a `~kwant.operator.Current` operator on a wavefunction returns a
+1D array of values that can be directly used with `kwant.plotter.current`:
+
+.. image:: ../images/spin_currents.*
+
+.. note::
+
+ Evaluating a `~kwant.operator.Current` operator on a wavefunction
+ returns a 1D array of the same length as the number of hoppings in the
+ system, ordered in the same way as the edges in the system's graph.
+
+.. specialnote:: Technical Details
+
+ Similarly to how we saw in the previous section that `~kwant.operator.Density`
+ can be thought of as a collection of operators, `~kwant.operator.Current`
+ can be defined in a similar way. Starting from the definition of a "density":
+
+ .. math:: Ï_a = \sum_{αβ} ψ^*_{α} \mathcal{M}_{aαβ} ψ_{β},
+
+ we can define *currents* :math:`J_{ab}` via the continuity equation:
+
+ .. math:: \frac{∂Ï_a}{∂t} - \sum_{b} J_{ab} = 0
+
+ where the sum runs over sites :math:`b` neigboring site :math:`a`.
+ Plugging in the definition for :math:`Ï_a`, along with the Schrödinger
+ equation and the assumption that :math:`\mathcal{M}` is time independent,
+ gives:
+
+ .. math:: J_{ab} = \sum_{αβ}
+ ψ^*_α \left(i \sum_{γ}
+ \mathcal{H}^*_{abγα} \mathcal{M}_{aγβ}
+ - \mathcal{M}_{aαγ} \mathcal{H}_{abγβ}
+ \right) ψ_β,
+
+ where latin indices run over sites and greek indices run over the Hilbert
+ space degrees of freedom, and
+
+ .. math:: \mathcal{H}_{ab} = \left(\begin{matrix}
+ ⋱ & ⋮ & ⋮ & ⋮ & ⋰ \\
+ ⋯ & ⋱ & 0 & \mathbf{H}_{ab} & ⋯ \\
+ ⋯ & 0 & ⋱ & 0 & ⋯ \\
+ ⋯ & 0 & 0 & ⋱ & ⋯ \\
+ ⋰ & ⋮ & ⋮ & ⋮ & ⋱
+ \end{matrix}\right).
+
+ i.e. :math:`\mathcal{H}_{ab}` is a matrix that is zero everywhere
+ except on elements connecting *from* site :math:`b` *to* site :math:`a`,
+ where it is equal to the hopping matrix :math:`\mathbf{H}_{ab}` between
+ these two sites.
+
+ This allows us to identify the rank-4 quantity
+
+ .. math:: \mathcal{J}_{abαβ} = i \sum_{γ}
+ \mathcal{H}^*_{abγα} \mathcal{M}_{aγβ}
+ - \mathcal{M}_{aαγ} \mathcal{H}_{abγβ}
+
+ as the local current between connected sites.
+
+ The diagonal part of this quantity, :math:`\mathcal{J}_{aa}`,
+ represents the extent to which the density defined by :math:`\mathcal{M}_a`
+ is not conserved on site :math:`a`. It can be calculated using
+ `~kwant.operator.Source`, rather than `~kwant.operator.Current`, which
+ only computes the off-diagonal part.
+
+
+Spatially varying operators
+---------------------------
+The above examples are reasonably simple in the sense that the book-keeping
+required to manually calculate the various densities and currents is still
+manageable. Now we shall look at the case where we wish to calculate some
+projected spin currents, but where the spin projection axis varies from place
+to place. More specifically, we want to visualize the spin current along the
+direction of :math:`\mathbf{m}_i`, which changes continuously over the whole
+scattering region.
+
+Doing this is as simple as passing a *function* when instantiating
+the `~kwant.operator.Current`, instead of a constant matrix:
+
+.. literalinclude:: magnetic_texture.py
+ :start-after: #HIDDEN_BEGIN_following
+ :end-before: #HIDDEN_END_following
+
+The function must take a `~kwant.builder.Site` as its first parameter,
+and may optionally take other parameters (i.e. it must have the same
+signature as a Hamiltonian onsite function), and must return the square
+matrix that defines the operator we wish to calculate.
+
+.. note::
+
+ In the above example we had to pass the extra parameters needed by the
+ ``following_operator`` function via the ``param`` keyword argument. In
+ general you must pass all the parameters needed by the Hamiltonian via
+ ``params`` (as you would when calling `~kwant.solvers.default.smatrix` or
+ `~kwant.solvers.default.wave_function`). In the previous examples,
+ however, we used the fact that the system hoppings do not depend on any
+ parameters (these are the only Hamiltonian elements required to calculate
+ currents) to avoid passing the system parameters for the sake of brevity.
+
+Using this we can see that the spin current is essentially oriented along
+the direction of :math:`m_i` in the present regime where the onsite term
+in the Hamiltonian is dominant:
+
+.. image:: ../images/spin_current_comparison.*
+
+.. note:: Although this example used exclusively `~kwant.operator.Current`,
+ you can do the same with `~kwant.operator.Density`.
+
+
+Defining operators over parts of a system
+-----------------------------------------
+
+Another useful feature of `kwant.operator` is the ability to calculate
+operators over selected parts of a system. For example, we may wish to
+calculate the total density of states in a certain part
+of the system, or the current flowing through a cut in the system.
+We can do this selection when creating the operator by using the
+keyword parameter ``where``.
+
+Density of states in a circle
+*****************************
+
+To calculate the density of states inside a circle of radius
+20 we can simply do:
+
+.. literalinclude:: magnetic_texture.py
+ :start-after: #HIDDEN_BEGIN_density_cut
+ :end-before: #HIDDEN_END_density_cut
+
+.. literalinclude:: ../images/circle_dos.txt
+
+note that we also provide ``sum=True``, which means that evaluating the
+operator on a wavefunction will produce a single scalar. This is semantically
+equivalent to providing ``sum=False`` (the default) and running ``numpy.sum``
+on the output.
+
+Current flowing through a cut
+*****************************
+
+Below we calculate the probability current and z-projected spin current near
+the interfaces with the left and right leads.
+
+.. literalinclude:: magnetic_texture.py
+ :start-after: #HIDDEN_BEGIN_current_cut
+ :end-before: #HIDDEN_END_current_cut
+
+.. literalinclude:: ../images/current_cut.txt
+
+We see that the probability current is conserved across the scattering region,
+but the z-projected spin current is not due to the fact that the Hamiltonian
+does not commute with :math:`σ_z` everywhere in the scattering region.
+
+.. note:: ``where`` can also be provided as a sequence of `~kwant.builder.Site`
+ or a sequence of hoppings (i.e. pairs of `~kwant.builder.Site`),
+ rather than a function.
+
+
+Advanced Topics
+---------------
+
+Using ``bind`` for speed
+************************
+In most of the above examples we only used each operator *once* after creating
+it. Often one will want to evaluate an operator with many different
+wavefunctions, for example with all scattering wavefunctions at a certain
+energy, but with the *same set of parameters*. In such cases it is best to tell
+the operator to pre-compute the onsite matrices and any necessary Hamiltonian
+elements using the given set of parameters, so that this work is not duplicated
+every time the operator is evaluated.
+
+This can be achieved with `~kwant.operator.Current.bind`:
+
+.. warning:: Take care that you do not use an operator that was bound to a
+ particular set of parameters with wavefunctions calculated with a
+ *different* set of parameters. This will almost certainly give
+ incorrect results.
+
+.. literalinclude:: magnetic_texture.py
+ :start-after: #HIDDEN_BEGIN_bind
+ :end-before: #HIDDEN_END_bind
+
+.. image:: ../images/bound_current.*
--
GitLab