From d919806ba021c1394706a6b9127bee14ddd0c1be Mon Sep 17 00:00:00 2001
From: Rafal Skolasinski <r.j.skolasinski@gmail.com>
Date: Tue, 21 Mar 2017 12:29:05 +0100
Subject: [PATCH] add discretizer tutorial
---
.gitlab-ci.yml | 2 +
.../images/continuum_discretizer.py.diff | 93 +++++
doc/source/pre/whatsnew/1.3.rst | 8 +
doc/source/reference/index.rst | 1 +
doc/source/reference/kwant.continuum.rst | 21 ++
doc/source/tutorial/continuum_discretizer.py | 216 ++++++++++++
doc/source/tutorial/index.rst | 1 +
doc/source/tutorial/tutorial1.rst | 2 +
doc/source/tutorial/tutorial8.rst | 318 ++++++++++++++++++
kwant/continuum/_common.py | 6 +-
kwant/continuum/discretizer.py | 14 +-
kwant/continuum/tests/test_common.py | 10 +
12 files changed, 682 insertions(+), 10 deletions(-)
create mode 100644 doc/source/images/continuum_discretizer.py.diff
create mode 100644 doc/source/reference/kwant.continuum.rst
create mode 100644 doc/source/tutorial/continuum_discretizer.py
create mode 100644 doc/source/tutorial/tutorial8.rst
diff --git a/.gitlab-ci.yml b/.gitlab-ci.yml
index 2e9536a6..9d0e4c4d 100644
--- a/.gitlab-ci.yml
+++ b/.gitlab-ci.yml
@@ -57,6 +57,7 @@ check for dependencies installed:
build documentation:
stage: test
script:
+ - pip3 install sympy
- make -C doc realclean; make -C doc html SPHINXOPTS='-A website_deploy=True -n -W'
artifacts:
paths:
@@ -77,6 +78,7 @@ run tests:
check for broken links in doc:
stage: test
script:
+ - pip3 install sympy
- make -C doc linkcheck
allow_failure: true
diff --git a/doc/source/images/continuum_discretizer.py.diff b/doc/source/images/continuum_discretizer.py.diff
new file mode 100644
index 00000000..885ab797
--- /dev/null
+++ b/doc/source/images/continuum_discretizer.py.diff
@@ -0,0 +1,93 @@
+--- original
++++ modified
+@@ -9,6 +9,8 @@
+ # --------------------------
+ # - discretizer module from kwant.continuum
+
++import _defs
++from contextlib import redirect_stdout
+
+ import kwant
+ import scipy.sparse.linalg
+@@ -20,9 +22,19 @@
+ from matplotlib import pyplot as plt
+
+
++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 stadium_system(L=20, W=20):
+ hamiltonian = "k_x**2 + k_y**2 + V(x, y)"
+- template = kwant.continuum.discretize(hamiltonian, verbose=True)
++ with open('discretizer_verbose.txt', 'w') as f:
++ with redirect_stdout(f):
++ template = kwant.continuum.discretize(hamiltonian, verbose=True)
+
+ def stadium(site):
+ (x, y) = site.pos
+@@ -43,7 +55,7 @@
+ ham = syst.hamiltonian_submatrix(params=dict(V=potential), sparse=True)
+ evecs = scipy.sparse.linalg.eigsh(ham, k=10, which='SM')[1]
+ kwant.plotter.map(syst, abs(evecs[:, n])**2, show=False)
+- plt.show()
++ save_figure('discretizer_gs')
+
+
+ def qsh_system(a=20, L=2000, W=1000):
+@@ -90,7 +102,8 @@
+ plt.ylim(-0.05, 0.05)
+ plt.xlabel('momentum [1/A]')
+ plt.ylabel('energy [eV]')
+- plt.show()
++ save_figure('discretizer_qsh_band')
++
+ # get scattering wave functions at E=0
+ wf = kwant.wave_function(syst, energy=0, params=params)
+
+@@ -118,7 +131,7 @@
+
+ ax1.set_title('Probability density')
+ ax2.set_title('Spin density')
+- plt.show()
++ save_figure('discretizer_qsh_wf')
+
+
+ def lattice_spacing():
+@@ -155,7 +168,7 @@
+
+ plot(ax1, a=1)
+ plot(ax2, a=.25)
+- plt.show()
++ save_figure('discretizer_lattice_spacing')
+
+
+ def substitutions():
+@@ -168,14 +181,20 @@
+ sympify('k_x**2 * sz + alpha * k_x * sx', substitutions=subs),
+ )
+
+- print(e[0] == e[1] == e[2])
++ with open('discretizer_subs_1.txt', 'w') as f:
++ with redirect_stdout(f):
++ print(e[0] == e[1] == e[2])
+
+ subs = {'A': 'A(x) + B', 'V': 'V(x) + V_0', 'C': 5}
+- print(sympify('k_x * A * k_x + V + C', substitutions=subs))
++ with open('discretizer_subs_2.txt', 'w') as f:
++ with redirect_stdout(f):
++ print(sympify('k_x * A * k_x + V + C', substitutions=subs))
+
+
+ def main():
+- template = kwant.continuum.discretize('k_x * A(x) * k_x', verbose=True)
++ with open('discretizer_intro_verbose.txt', 'w') as f:
++ with redirect_stdout(f):
++ kwant.continuum.discretize('k_x * A(x) * k_x', verbose=True)
+
+ syst = stadium_system()
+ plot_eigenstate(syst)
diff --git a/doc/source/pre/whatsnew/1.3.rst b/doc/source/pre/whatsnew/1.3.rst
index 325c581a..3bdf5259 100644
--- a/doc/source/pre/whatsnew/1.3.rst
+++ b/doc/source/pre/whatsnew/1.3.rst
@@ -133,3 +133,11 @@ Reference implementation of the Kernel Polynomial Method
The kernel polynomial method is now implemented within Kwant to obtain the
density of states or, more generally, the spectral density of a given operator
acting on a system or Hamiltonian.
+
+Tools for coninuum Hamiltonians
+-------------------------------
+The `~kwant.continuum` sub-package is a collection of tools for working with
+continuum models and for discretizing them into tight-binding models. It aims
+at providing a handy interface to convert symbolic Hamiltonian into a builder
+with N-D translational symmetry that can be use to calculate TB band structures
+or construct systems with different/lower symmetry.
diff --git a/doc/source/reference/index.rst b/doc/source/reference/index.rst
index 479356d6..cd9be611 100644
--- a/doc/source/reference/index.rst
+++ b/doc/source/reference/index.rst
@@ -38,4 +38,5 @@ The following modules provide functionality for special applications.
kwant.digest
kwant.rmt
kwant.kpm
+ kwant.continuum
kwant.wraparound
diff --git a/doc/source/reference/kwant.continuum.rst b/doc/source/reference/kwant.continuum.rst
new file mode 100644
index 00000000..d6799438
--- /dev/null
+++ b/doc/source/reference/kwant.continuum.rst
@@ -0,0 +1,21 @@
+:mod:`kwant.continuum` -- Tools for continuum systems
+=====================================================
+
+.. module:: kwant.continuum
+
+Discretizer
+-----------
+.. autosummary::
+ :toctree: generated/
+
+ discretize
+ discretize_symbolic
+ build_discretized
+
+Symbolic helpers
+----------------
+.. autosummary::
+ :toctree: generated/
+
+ sympify
+ lambdify
diff --git a/doc/source/tutorial/continuum_discretizer.py b/doc/source/tutorial/continuum_discretizer.py
new file mode 100644
index 00000000..70ef0589
--- /dev/null
+++ b/doc/source/tutorial/continuum_discretizer.py
@@ -0,0 +1,216 @@
+# Tutorial 2.9. Processing continuum Hamiltonians with discretizer
+# ================================================================
+#
+# Physics background
+# ------------------
+# - tight-binding approximation of continous Hamiltonians
+#
+# Kwant features highlighted
+# --------------------------
+# - discretizer module from kwant.continuum
+
+
+import kwant
+import scipy.sparse.linalg
+import scipy.linalg
+import numpy as np
+
+# For plotting
+import matplotlib as mpl
+from matplotlib import pyplot as plt
+
+
+def stadium_system(L=20, W=20):
+#HIDDEN_BEGIN_template
+ hamiltonian = "k_x**2 + k_y**2 + V(x, y)"
+ template = kwant.continuum.discretize(hamiltonian, verbose=True)
+#HIDDEN_END_template
+
+#HIDDEN_BEGIN_fill
+ def stadium(site):
+ (x, y) = site.pos
+ x = max(abs(x) - 20, 0)
+ return x**2 + y**2 < 30**2
+
+ syst = kwant.Builder()
+ syst.fill(template, stadium, (0, 0));
+ syst = syst.finalized()
+#HIDDEN_END_fill
+ return syst
+
+
+#HIDDEN_BEGIN_plot_eigenstate
+def plot_eigenstate(syst, n=2, Vx=.0003, Vy=.0005):
+
+ def potential(x, y):
+ return Vx * x + Vy * y
+
+ ham = syst.hamiltonian_submatrix(params=dict(V=potential), sparse=True)
+ evecs = scipy.sparse.linalg.eigsh(ham, k=10, which='SM')[1]
+ kwant.plotter.map(syst, abs(evecs[:, n])**2, show=False)
+#HIDDEN_END_plot_eigenstate
+ plt.show()
+
+
+def qsh_system(a=20, L=2000, W=1000):
+#HIDDEN_BEGIN_define_qsh
+ hamiltonian = """
+ + C * identity(4) + M * kron(sigma_0, sigma_z)
+ - B * (k_x**2 + k_y**2) * kron(sigma_0, sigma_z)
+ - D * (k_x**2 + k_y**2) * kron(sigma_0, sigma_0)
+ + A * k_x * kron(sigma_z, sigma_x)
+ - A * k_y * kron(sigma_0, sigma_y)
+ """
+
+ template = kwant.continuum.discretize(hamiltonian, grid_spacing=a)
+#HIDDEN_END_define_qsh
+
+#HIDDEN_BEGIN_define_qsh_build
+ def shape(site):
+ (x, y) = site.pos
+ return (0 <= y < W and 0 <= x < L)
+
+ def lead_shape(site):
+ (x, y) = site.pos
+ return (0 <= y < W)
+
+ syst = kwant.Builder()
+ syst.fill(template, shape, (0, 0))
+
+ lead = kwant.Builder(kwant.TranslationalSymmetry([-a, 0]))
+ lead.fill(template, lead_shape, (0, 0))
+
+ syst.attach_lead(lead)
+ syst.attach_lead(lead.reversed())
+
+ syst = syst.finalized()
+#HIDDEN_END_define_qsh_build
+ return syst
+
+
+def analyze_qsh():
+ params = dict(A=3.65, B=-68.6, D=-51.1, M=-0.01, C=0)
+ syst = qsh_system()
+
+#HIDDEN_BEGIN_plot_qsh_band
+ kwant.plotter.bands(syst.leads[0], params=params,
+ momenta=np.linspace(-0.3, 0.3, 201), show=False)
+#HIDDEN_END_plot_qsh_band
+
+ plt.grid()
+ plt.xlim(-.3, 0.3)
+ plt.ylim(-0.05, 0.05)
+ plt.xlabel('momentum [1/A]')
+ plt.ylabel('energy [eV]')
+ plt.show()
+#HIDDEN_BEGIN_plot_qsh_wf
+ # get scattering wave functions at E=0
+ wf = kwant.wave_function(syst, energy=0, params=params)
+
+ # prepare density operators
+ sigma_z = np.array([[1, 0], [0, -1]])
+ prob_density = kwant.operator.Density(syst, np.eye(4))
+ spin_density = kwant.operator.Density(syst, np.kron(sigma_z, np.eye(2)))
+
+ # calculate expectation values and plot them
+ wf_sqr = sum(prob_density(psi) for psi in wf(0)) # states from left lead
+ rho_sz = sum(spin_density(psi) for psi in wf(0)) # states from left lead
+
+ fig, (ax1, ax2) = plt.subplots(1, 2, sharey=True, figsize=(16, 4))
+ kwant.plotter.map(syst, wf_sqr, ax=ax1)
+ kwant.plotter.map(syst, rho_sz, ax=ax2)
+#HIDDEN_END_plot_qsh_wf
+ ax = ax1
+ im = [obj for obj in ax.get_children()
+ if isinstance(obj, mpl.image.AxesImage)][0]
+ fig.colorbar(im, ax=ax)
+
+ ax = ax2
+ im = [obj for obj in ax.get_children()
+ if isinstance(obj, mpl.image.AxesImage)][0]
+ fig.colorbar(im, ax=ax)
+
+ ax1.set_title('Probability density')
+ ax2.set_title('Spin density')
+ plt.show()
+
+
+def lattice_spacing():
+#HIDDEN_BEGIN_ls_def
+ hamiltonian = "k_x**2 * identity(2) + alpha * k_x * sigma_y"
+ params = dict(alpha=.5)
+#HIDDEN_END_ls_def
+
+ def plot(ax, a=1):
+#HIDDEN_BEGIN_ls_hk_cont
+ h_k = kwant.continuum.lambdify(hamiltonian, substitutions=params)
+ k_cont = np.linspace(-4, 4, 201)
+ e_cont = [scipy.linalg.eigvalsh(h_k(k_x=ki)) for ki in k_cont]
+#HIDDEN_END_ls_hk_cont
+
+#HIDDEN_BEGIN_ls_hk_tb
+ lead = kwant.continuum.discretize(hamiltonian, grid_spacing=a)
+ band = kwant.physics.Bands(lead.finalized(), params=params)
+
+ k_tb = np.linspace(-np.pi/a, np.pi/a, 201)
+ e_tb = np.array([band(a*ki) for ki in k_tb])
+#HIDDEN_END_ls_hk_tb
+
+ ax.plot(k_cont, e_cont, 'r-')
+ ax.plot(k_tb, e_tb, 'k-')
+
+ ax.plot([], [], 'r-', label='continuum')
+ ax.plot([], [], 'k-', label='tight-binding')
+
+ ax.set_xlim(-4, 4)
+ ax.set_ylim(-1, 14)
+ ax.set_title('a={}'.format(a))
+
+ ax.set_xlabel('momentum [a.u.]')
+ ax.set_ylabel('energy [a.u.]')
+ ax.grid()
+ ax.legend()
+
+ _, (ax1, ax2) = plt.subplots(1, 2, sharey=True, figsize=(12, 4))
+
+ plot(ax1, a=1)
+ plot(ax2, a=.25)
+ plt.show()
+
+
+def substitutions():
+#HIDDEN_BEGIN_subs_1
+ sympify = kwant.continuum.sympify
+ subs = {'sx': [[0, 1], [1, 0]], 'sz': [[1, 0], [0, -1]]}
+
+ e = (
+ sympify('[[k_x**2, alpha * k_x], [k_x * alpha, -k_x**2]]'),
+ sympify('k_x**2 * sigma_z + alpha * k_x * sigma_x'),
+ sympify('k_x**2 * sz + alpha * k_x * sx', substitutions=subs),
+ )
+
+ print(e[0] == e[1] == e[2])
+#HIDDEN_END_subs_1
+
+#HIDDEN_BEGIN_subs_2
+ subs = {'A': 'A(x) + B', 'V': 'V(x) + V_0', 'C': 5}
+ print(sympify('k_x * A * k_x + V + C', substitutions=subs))
+#HIDDEN_END_subs_2
+
+
+def main():
+#HIDDEN_BEGIN_symbolic_discretization
+ template = kwant.continuum.discretize('k_x * A(x) * k_x', verbose=True)
+#HIDDEN_END_symbolic_discretization
+
+ syst = stadium_system()
+ plot_eigenstate(syst)
+
+ analyze_qsh()
+ lattice_spacing()
+ substitutions()
+
+# Call the main function if the script gets executed (as opposed to imported).
+# See <http://docs.python.org/library/__main__.html>.
+if __name__ == '__main__':
+ main()
diff --git a/doc/source/tutorial/index.rst b/doc/source/tutorial/index.rst
index dcd702f4..1729a231 100644
--- a/doc/source/tutorial/index.rst
+++ b/doc/source/tutorial/index.rst
@@ -10,3 +10,4 @@ Tutorial: learning Kwant through examples
tutorial5
tutorial6
tutorial7
+ tutorial8
diff --git a/doc/source/tutorial/tutorial1.rst b/doc/source/tutorial/tutorial1.rst
index b4698c0b..6aac3667 100644
--- a/doc/source/tutorial/tutorial1.rst
+++ b/doc/source/tutorial/tutorial1.rst
@@ -1,6 +1,8 @@
First steps: setting up a simple system and computing conductance
-----------------------------------------------------------------
+.. _tutorial_discretization_schrodinger:
+
Discretization of a Schrödinger Hamiltonian
...........................................
diff --git a/doc/source/tutorial/tutorial8.rst b/doc/source/tutorial/tutorial8.rst
new file mode 100644
index 00000000..216f24b1
--- /dev/null
+++ b/doc/source/tutorial/tutorial8.rst
@@ -0,0 +1,318 @@
+Discretizing continuous Hamiltonians
+------------------------------------
+
+Introduction
+............
+
+In ":ref:`tutorial_discretization_schrodinger`" we have learnt that Kwant works
+with tight-binding Hamiltonians. Often, however, one will start with a
+continuum model and will subsequently need to discretize it to arrive at a
+tight-binding model.
+Although discretizing a Hamiltonian is usually a simple
+process, it is tedious and repetitive. The situation is further exacerbated
+when one introduces additional on-site degrees of freedom, and tracking all
+the necessary terms becomes a chore.
+The `~kwant.continuum` sub-package aims to be a solution to this problem.
+It is a collection of tools for working with
+continuum models and for discretizing them into tight-binding models.
+
+.. seealso::
+ The complete source code of this tutorial can be found in
+ :download:`tutorial/continuum_discretizer.py <../../../tutorial/continuum_discretizer.py>`
+
+
+.. _tutorial_discretizer_introduction:
+
+Discretizing by hand
+....................
+
+As an example, let us consider the following continuum Schrödinger equation
+for a semiconducting heterostructure (using the effective mass approximation):
+
+.. math::
+
+ \left( k_x \frac{\hbar^2}{2 m(x)} k_x \right) \psi(x) = E \, \psi(x).
+
+Replacing the momenta by their corresponding differential operators
+
+.. math::
+ k_\alpha = -i \partial_\alpha,
+
+for :math:`\alpha = x, y` or :math:`z`, and discretizing on a regular lattice of
+points with spacing :math:`a`, we obtain the tight-binding model
+
+.. math::
+
+ H = - \frac{1}{a^2} \sum_i A\left(x+\frac{a}{2}\right)
+ \big(\ket{i}\bra{i+1} + h.c.\big)
+ + \frac{1}{a^2} \sum_i
+ \left( A\left(x+\frac{a}{2}\right) + A\left(x-\frac{a}{2}\right)\right)
+ \ket{i} \bra{i},
+
+with :math:`A(x) = \frac{\hbar^2}{2 m(x)}`.
+
+Using `~kwant.continuum.discretize` to obtain a template
+........................................................
+
+The function `kwant.continuum.discretize` takes a symbolic Hamiltonian and
+turns it into a `~kwant.builder.Builder` instance with appropriate spatial
+symmetry that serves as a template.
+(We will see how to use the template to build systems with a particular
+shape later).
+
+.. literalinclude:: continuum_discretizer.py
+ :start-after: #HIDDEN_BEGIN_symbolic_discretization
+ :end-before: #HIDDEN_END_symbolic_discretization
+
+It is worth noting that ``discretize`` treats ``k_x`` and ``x`` as
+non-commuting operators, and so their order is preserved during the
+discretization process.
+
+Setting the ``verbose`` parameter to ``True`` prints extra information about the
+onsite and hopping functions assigned to the ``Builder`` produced
+by ``discretize``:
+
+.. literalinclude:: ../images/discretizer_intro_verbose.txt
+
+.. specialnote:: Technical details
+
+ - ``kwant.continuum`` uses ``sympy`` internally to handle symbolic
+ expressions. Strings are converted using `kwant.continuum.sympify`,
+ which essentially applies some Kwant-specific rules (such as treating
+ ``k_x`` and ``x`` as non-commutative) before calling ``sympy.sympify``
+
+ - The builder returned by ``discretize`` will have an N-D
+ translational symmetry, where ``N`` is the number of dimensions that were
+ discretized. This is the case, even if there are expressions in the input
+ (e.g. ``V(x, y)``) which in principle *may not* have this symmetry. When
+ using the returned builder directly, or when using it as a template to
+ construct systems with different/lower symmetry, it is important to
+ ensure that any functional parameters passed to the system respect the
+ symmetry of the system. Kwant provides no consistency check for this.
+
+ - The discretization process consists of taking input
+ :math:`H(k_x, k_y, k_z)`, multiplying it from the right by
+ :math:`\psi(x, y, z)` and iteratively applying a second-order accurate
+ central derivative approximation for every
+ :math:`k_\alpha=-i\partial_\alpha`:
+
+ .. math::
+ \partial_\alpha \psi(\alpha) =
+ \frac{1}{a} \left( \psi\left(\alpha + \frac{a}{2}\right)
+ -\psi\left(\alpha - \frac{a}{2}\right)\right).
+
+ This process is done separately for every summand in Hamiltonian.
+ Once all symbols denoting operators are applied internal algorithm is
+ calculating ``gcd`` for hoppings coming from each summand in order to
+ find best possible approximation. Please see source code for details.
+
+ - Instead of using ``discretize`` one can use
+ `~kwant.continuum.discretize_symbolic` to obtain symbolic output.
+ When working interactively in `Jupyter notebooks <https://jupyter.org/>`_
+ it can be useful to use this to see a symbolic representation of
+ the discretized Hamiltonian. This works best when combined with ``sympy``
+ `Pretty Printing <http://docs.sympy.org/latest/tutorial/printing.html#setting-up-pretty-printing>`_.
+
+ - The symbolic result of discretization obtained with
+ ``discretize_symbolic`` can be converted into a
+ builder using `~kwant.continuum.build_discretized`.
+ This can be useful if one wants to alter the tight-binding Hamiltonian
+ before building the system.
+
+
+Building a Kwant system from the template
+.........................................
+
+Let us now use the output of ``discretize`` as a template to
+build a system and plot some of its energy eigenstate. For this example the
+Hamiltonian will be
+
+.. math::
+
+ H = k_x^2 + k_y^2 + V(x, y),
+
+where :math:`V(x, y)` is some arbitrary potential.
+
+First, use ``discretize`` to obtain a
+builder that we will use as a template:
+
+.. literalinclude:: continuum_discretizer.py
+ :start-after: #HIDDEN_BEGIN_template
+ :end-before: #HIDDEN_END_template
+
+We now use this system with the `~kwant.builder.Builder.fill`
+method of `~kwant.builder.Builder` to construct the system we
+want to investigate:
+
+.. literalinclude:: continuum_discretizer.py
+ :start-after: #HIDDEN_BEGIN_fill
+ :end-before: #HIDDEN_END_fill
+
+After finalizing this system, we can plot one of the system's
+energy eigenstates:
+
+.. literalinclude:: continuum_discretizer.py
+ :start-after: #HIDDEN_BEGIN_plot_eigenstate
+ :end-before: #HIDDEN_END_plot_eigenstate
+
+.. image:: ../images/discretizer_gs.*
+
+Note in the above that we provided the function ``V`` to
+``syst.hamiltonian_submatrix`` using ``params=dict(V=potential)``, rather than
+via ``args``.
+
+In addition, the function passed as ``V`` expects two input parameters ``x``
+and ``y``, the same as in the initial continuum Hamiltonian.
+
+
+Models with more structure: Bernevig-Hughes-Zhang
+.................................................
+
+When working with multi-band systems, like the Bernevig-Hughes-Zhang (BHZ)
+model [1]_ [2]_, one can provide matrix input to `~kwant.continuum.discretize`
+using ``identity`` and ``kron``. For example, the definition of the BHZ model can be
+written succinctly as:
+
+.. literalinclude:: continuum_discretizer.py
+ :start-after: #HIDDEN_BEGIN_define_qsh
+ :end-before: #HIDDEN_END_define_qsh
+
+We can then make a ribbon out of this template system:
+
+.. literalinclude:: continuum_discretizer.py
+ :start-after: #HIDDEN_BEGIN_define_qsh_build
+ :end-before: #HIDDEN_END_define_qsh_build
+
+and plot its dispersion using `kwant.plotter.bands`:
+
+.. literalinclude:: continuum_discretizer.py
+ :start-after: #HIDDEN_BEGIN_plot_qsh_band
+ :end-before: #HIDDEN_END_plot_qsh_band
+
+.. image:: ../images/discretizer_qsh_band.*
+
+In the above we see the edge states of the quantum spin Hall effect, which
+we can visualize using `kwant.plotter.map`:
+
+.. literalinclude:: continuum_discretizer.py
+ :start-after: #HIDDEN_BEGIN_plot_qsh_wf
+ :end-before: #HIDDEN_END_plot_qsh_wf
+
+.. image:: ../images/discretizer_qsh_wf.*
+
+
+Limitations of discretization
+.............................
+
+It is important to remember that the discretization of a continuum
+model is an *approximation* that is only valid in the low-energy
+limit. For example, the quadratic continuum Hamiltonian
+
+.. math::
+
+ H_\textrm{continuous}(k_x) = \frac{\hbar^2}{2m}k_x^2
+
+
+and its discretized approximation
+
+.. math::
+
+ H_\textrm{tight-binding}(k_x) = 2t \big(1 - \cos(k_x a)\big),
+
+
+where :math:`t=\frac{\hbar^2}{2ma^2}`, are only valid in the limit
+:math:`E \lt t`. The grid spacing :math:`a` must be chosen according
+to how high in energy you need your tight-binding model to be valid.
+
+It is possible to set :math:`a` through the ``grid_spacing`` parameter
+to `~kwant.continuum.discretize`, as we will illustrate in the following
+example. Let us start from the continuum Hamiltonian
+
+.. math::
+
+ H(k) = k_x^2 \mathbb{1}_{2\times2} + α k_x \sigma_y.
+
+We start by defining this model as a string and setting the value of the
+:math:`α` parameter:
+
+.. literalinclude:: continuum_discretizer.py
+ :start-after: #HIDDEN_BEGIN_ls_def
+ :end-before: #HIDDEN_END_ls_def
+
+Now we can use `kwant.continuum.lambdify` to obtain a function that computes
+:math:`H(k)`:
+
+.. literalinclude:: continuum_discretizer.py
+ :start-after: #HIDDEN_BEGIN_ls_hk_cont
+ :end-before: #HIDDEN_END_ls_hk_cont
+
+We can also construct a discretized approximation using
+`kwant.continuum.discretize`, in a similar manner to previous examples:
+
+.. literalinclude:: continuum_discretizer.py
+ :start-after: #HIDDEN_BEGIN_ls_hk_tb
+ :end-before: #HIDDEN_END_ls_hk_tb
+
+Below we can see the continuum and tight-binding dispersions for two
+different values of the discretization grid spacing :math:`a`:
+
+.. image:: ../images/discretizer_lattice_spacing.*
+
+
+We clearly see that the smaller grid spacing is, the better we approximate
+the original continuous dispersion. It is also worth remembering that the
+Brillouin zone also scales with grid spacing: :math:`[-\frac{\pi}{a},
+\frac{\pi}{a}]`.
+
+
+Advanced topics
+...............
+
+The input to `kwant.continuum.discretize` and `kwant.continuum.lambdify` can be
+not only a ``string``, as we saw above, but also a ``sympy`` expression or
+a ``sympy`` matrix.
+This functionality will probably be mostly useful to people who
+are already experienced with ``sympy``.
+
+
+It is possible to use ``identity`` (for identity matrix), ``kron`` (for Kronecker product), as well as Pauli matrices ``sigma_0``,
+``sigma_x``, ``sigma_y``, ``sigma_z`` in the input to
+`~kwant.continuum.lambdify` and `~kwant.continuum.discretize`, in order to simplify
+expressions involving matrices. Matrices can also be provided explicitly using
+square ``[]`` brackets. For example, all following expressions are equivalent:
+
+.. literalinclude:: continuum_discretizer.py
+ :start-after: #HIDDEN_BEGIN_subs_1
+ :end-before: #HIDDEN_END_subs_1
+
+.. literalinclude:: ../images/discretizer_subs_1.txt
+
+We can use the ``substitutions`` keyword parameter to substitute expressions
+and numerical values:
+
+.. literalinclude:: continuum_discretizer.py
+ :start-after: #HIDDEN_BEGIN_subs_2
+ :end-before: #HIDDEN_END_subs_2
+
+.. literalinclude:: ../images/discretizer_subs_2.txt
+
+Symbolic expressions obtained in this way can be directly passed to all
+``discretizer`` functions.
+
+.. specialnote:: Technical details
+
+ Because of the way that ``sympy`` handles commutation relations all symbols
+ representing position and momentum operators are set to be non commutative.
+ This means that the order of momentum and position operators in the input
+ expression is preserved. Note that it is not possible to define individual
+ commutation relations within ``sympy``, even expressions such :math:`x k_y x`
+ will not be simplified, even though mathematically :math:`[x, k_y] = 0`.
+
+
+.. rubric:: References
+
+.. [1] `Science, 314, 1757 (2006)
+ <https://arxiv.org/abs/cond-mat/0611399>`_.
+
+.. [2] `Phys. Rev. B 82, 045122 (2010)
+ <https://arxiv.org/abs/1005.1682>`_.
diff --git a/kwant/continuum/_common.py b/kwant/continuum/_common.py
index 0eabf86a..f614f2bb 100644
--- a/kwant/continuum/_common.py
+++ b/kwant/continuum/_common.py
@@ -30,7 +30,7 @@ pauli = (msigma(0), msigma(1), msigma(2), msigma(3))
_clash = _clash.copy()
_clash.update({s.name: s for s in momentum_operators})
_clash.update({s.name: s for s in position_operators})
-_clash.update({'kron': TensorProduct, 'eye': sympy.eye})
+_clash.update({'kron': TensorProduct, 'eye': sympy.eye, 'identity': sympy.eye})
_clash.update({'sigma_{}'.format(c): p for c, p in zip('0xyz', pauli)})
@@ -49,8 +49,8 @@ def lambdify(hamiltonian, substitutions=None):
----------
hamiltonian : sympy.Expr or sympy.Matrix, or string
Symbolic representation of a continous Hamiltonian. When providing
- a sympy expression, it is recommended to use operators
- defined in `~kwant.continuum.momentum_operators` in order to provide
+ a sympy expression, it is recommended to use ``momentum_operators``
+ defined in `~kwant.continuum` in order to provide
proper commutation properties. If a string is provided it will be
converted to a sympy expression using `kwant.continuum.sympify`.
substitutions : dict, defaults to empty
diff --git a/kwant/continuum/discretizer.py b/kwant/continuum/discretizer.py
index 5e3eab13..c10c6e07 100644
--- a/kwant/continuum/discretizer.py
+++ b/kwant/continuum/discretizer.py
@@ -46,8 +46,8 @@ def discretize(hamiltonian, discrete_coordinates=None, *, grid_spacing=1,
----------
hamiltonian : sympy.Expr or sympy.Matrix, or string
Symbolic representation of a continous Hamiltonian. When providing
- a sympy expression, it is recommended to use operators
- defined in `~kwant.continuum.momentum_operators` in order to provide
+ a sympy expression, it is recommended to use ``momentum_operators``
+ defined in `~kwant.continuum` in order to provide
proper commutation properties. If a string is provided it will be
converted to a sympy expression using `kwant.continuum.sympify`.
discrete_coordinates : sequence of strings, or ``None`` (default)
@@ -89,8 +89,8 @@ def discretize_symbolic(hamiltonian, discrete_coordinates=None, *,
----------
hamiltonian : sympy.Expr or sympy.Matrix, or string
Symbolic representation of a continous Hamiltonian. When providing
- a sympy expression, it is recommended to use operators
- defined in `~kwant.continuum.momentum_operators` in order to provide
+ a sympy expression, it is recommended to use ``momentum_operators``
+ defined in `~kwant.continuum` in order to provide
proper commutation properties. If a string is provided it will be
converted to a sympy expression using `kwant.continuum.sympify`.
discrete_coordinates : sequence of strings, or ``None`` (default)
@@ -522,8 +522,8 @@ def _value_function(expr, discrete_coordinates, grid_spacing, onsite,
"""
expr = expr.subs({sympy.Symbol('a'): grid_spacing})
- return_string, map_func_calls, const_symbols, _cache = \
- _return_string(expr, discrete_coordinates=discrete_coordinates)
+ return_string, map_func_calls, const_symbols, _cache = _return_string(
+ expr, discrete_coordinates=discrete_coordinates)
# first check if value function needs to read coordinates
atoms_names = {s.name for s in expr.atoms(sympy.Symbol)}
@@ -554,7 +554,7 @@ def _value_function(expr, discrete_coordinates, grid_spacing, onsite,
lines.append(return_string)
- separator = '\n' + 4 * ' '
+ separator = '\n '
# 'site_string' is tightly coupled to the symbols used in '_assign_symbol'
site_string = 'site' if onsite else 'site1, site2'
if required_kwargs:
diff --git a/kwant/continuum/tests/test_common.py b/kwant/continuum/tests/test_common.py
index 1ab80022..313f3193 100644
--- a/kwant/continuum/tests/test_common.py
+++ b/kwant/continuum/tests/test_common.py
@@ -13,6 +13,7 @@ import pytest
import tinyarray as ta
from sympy.physics.matrices import msigma
+from sympy.physics.quantum import TensorProduct
import sympy
from kwant.continuum._common import position_operators, momentum_operators
@@ -33,6 +34,11 @@ kx, ky, kz = momentum_operators
('[[k_x*A(x)*k_x, B(x, y)*k_x], [k_x*B(x, y), C*k_y**2]]',
sympy.Matrix([[kx*com_A(x_op)*kx, com_B(x_op, y_op)*kx],
[kx*com_B(x_op, y_op), com_C*ky**2]])),
+ ('kron(sigma_x, sigma_y)', TensorProduct(msigma(1), msigma(2))),
+ ('identity(2)', sympy.eye(2)),
+ ('eye(2)', sympy.eye(2)),
+ ('1 * sigma_x + 2 * sigma_y + 3 * sigma_z',
+ msigma(1) + 2 * msigma(2) + 3 * msigma(3))
])
def test_sympify(input_expr, output_expr):
assert sympify(input_expr) == output_expr
@@ -60,8 +66,12 @@ def test_sympify_substitutions(input_expr, output_expr, subs):
@pytest.mark.parametrize('input_expr, output_expr, subs', [
('A + k_x**2 * eye(2)',
+ kx**2 * sympy.eye(2) + msigma(2),
+ {'A': "[[0, -1j], [1j, 0]]"}),
+ ('A + k_x**2 * identity(2)',
kx**2 * sympy.eye(2) + msigma(2),
{'A': "[[0, -1j], [1j, 0]]"})
+
])
def test_sympify_mix_symbol_and_matrx(input_expr, output_expr, subs):
assert sympify(input_expr, substitutions=subs) == output_expr
--
GitLab