Skip to content
Snippets Groups Projects

Compare revisions

Changes are shown as if the source revision was being merged into the target revision. Learn more about comparing revisions.

Source

Select target project
No results found

Target

Select target project
  • kwant/kwant
  • jbweston/kwant
  • anton-akhmerov/kwant
  • cwg/kwant
  • Mathieu/kwant
  • slavoutich/kwant
  • pacome/kwant
  • behrmann/kwant
  • michaelwimmer/kwant
  • albeercik/kwant
  • eunjongkim/kwant
  • basnijholt/kwant
  • r-j-skolasinski/kwant
  • sahmed95/kwant
  • pablopiskunow/kwant
  • mare/kwant
  • dvarjas/kwant
  • Paul/kwant
  • bbuijtendorp/kwant
  • tkloss/kwant
  • torosdahl/kwant
  • kel85uk/kwant
  • kpoyhonen/kwant
  • Fromeworld/kwant
  • quaeritis/kwant
  • marwahaha/kwant
  • fernandodfufrpe/kwant
  • oly/kwant
  • jiamingh/kwant
  • mehdi2369/kwant
  • ValFadeev/kwant
  • Kostas/kwant
  • chelseabaptiste03/kwant
33 results
Show changes
Hide whitespace changes
Inline Side-by-side
Showing
with 3803 additions and 798 deletions
doc/source/tutorial/first_steps.rst
View file @ 1d8db7bb
......@@ -13,7 +13,7 @@ Schrödinger equation
.. math::
H = \frac{-\hbar^2}{2m}(\partial_x^2 + \partial_y^2) + V(y)
with a hard-wall confinement :math:`V(y)` in y-direction.
with a hard-wall confinement :math:`V(y)` in the y-direction.
To be able to implement the quantum wire with Kwant, the continuous Hamiltonian
:math:`H` has to be discretized thus turning it into a tight-binding
......@@ -60,25 +60,58 @@ simplicity, we choose to work in such units that :math:`t = a = 1`.
Transport through a quantum wire
................................
.. seealso::
You can execute the code examples live in your browser by
activating thebelab:
.. thebe-button:: Activate Thebelab
.. seealso::
The complete source code of this example can be found in
:download:`quantum_wire.py </code/download/quantum_wire.py>`
:jupyter-download:script:`quantum_wire`
In order to use Kwant, we need to import it:
.. literalinclude:: /code/include/quantum_wire.py
:start-after: #HIDDEN_BEGIN_dwhx
:end-before: #HIDDEN_END_dwhx
.. jupyter-kernel::
:id: quantum_wire
.. jupyter-execute::
:hide-code:
# 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
# For plotting
from matplotlib import pyplot
.. jupyter-execute:: boilerplate.py
:hide-code:
.. jupyter-execute::
import kwant
Enabling Kwant is as easy as this [#]_ !
The first step is now the definition of the system with scattering region and
The first step is now to define the system with scattering region and
leads. For this we make use of the `~kwant.builder.Builder` type that allows to
define a system in a convenient way. We need to create an instance of it:
.. literalinclude:: /code/include/quantum_wire.py
:start-after: #HIDDEN_BEGIN_goiq
:end-before: #HIDDEN_END_goiq
.. jupyter-execute::
# First define the tight-binding system
syst = kwant.Builder()
Observe that we just accessed `~kwant.builder.Builder` by the name
``kwant.Builder``. We could have just as well written
......@@ -92,16 +125,19 @@ Apart from `~kwant.builder.Builder` we also need to specify
what kind of sites we want to add to the system. Here we work with
a square lattice. For simplicity, we set the lattice constant to unity:
.. literalinclude:: /code/include/quantum_wire.py
:start-after: #HIDDEN_BEGIN_suwo
:end-before: #HIDDEN_END_suwo
.. jupyter-execute::
a = 1
lat = kwant.lattice.square(a, norbs=1)
Since we work with a square lattice, we label the points with two
integer coordinates `(i, j)`. `~kwant.builder.Builder` then
allows us to add matrix elements corresponding to lattice points:
``syst[lat(i, j)] = ...`` sets the on-site energy for the point `(i, j)`,
and ``syst[lat(i1, j1), lat(i2, j2)] = ...`` the hopping matrix element
**from** point `(i2, j2)` **to** point `(i1, j1)`.
**from** point `(i2, j2)` **to** point `(i1, j1)`. In specifying ``norbs=1``
in the definition of the lattice we tell Kwant that there is 1 degree
of freedom per lattice site.
Note that we need to specify sites for `~kwant.builder.Builder`
in the form ``lat(i, j)``. The lattice object `lat` does the
......@@ -111,9 +147,25 @@ needed in Builder (more about that in the technical details below).
We now build a rectangular scattering region that is `W`
lattice points wide and `L` lattice points long:
.. literalinclude:: /code/include/quantum_wire.py
:start-after: #HIDDEN_BEGIN_zfvr
:end-before: #HIDDEN_END_zfvr
.. jupyter-execute::
t = 1.0
W, L = 10, 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
Observe how the above code corresponds directly to the terms of the discretized
Hamiltonian:
......@@ -140,9 +192,14 @@ Next, we define the leads. Leads are also constructed using
`~kwant.builder.Builder`, but in this case, the
system must have a translational symmetry:
.. literalinclude:: /code/include/quantum_wire.py
:start-after: #HIDDEN_BEGIN_xcmc
:end-before: #HIDDEN_END_xcmc
.. jupyter-execute::
# Then, define and attach the leads:
# First the lead to the left
# (Note: TranslationalSymmetry takes a real-space vector)
sym_left_lead = kwant.TranslationalSymmetry((-a, 0))
left_lead = kwant.Builder(sym_left_lead)
Here, the `~kwant.builder.Builder` takes a
`~kwant.lattice.TranslationalSymmetry` as the optional parameter. Note that the
......@@ -155,17 +212,22 @@ as the hoppings inside one unit cell and to the next unit cell of the lead.
For a square lattice, and a lead in y-direction the unit cell is
simply a vertical line of points:
.. literalinclude:: /code/include/quantum_wire.py
:start-after: #HIDDEN_BEGIN_ndez
:end-before: #HIDDEN_END_ndez
.. jupyter-execute::
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
Note that here it doesn't matter if you add the hoppings to the next or the
previous unit cell -- the translational symmetry takes care of that. The
isolated, infinite is attached at the correct position using
.. literalinclude:: /code/include/quantum_wire.py
:start-after: #HIDDEN_BEGIN_fskr
:end-before: #HIDDEN_END_fskr
.. jupyter-execute::
:hide-output:
syst.attach_lead(left_lead)
This call returns the lead number which will be used to refer to the lead when
computing transmissions (further down in this tutorial). More details about
......@@ -175,9 +237,21 @@ We also want to add a lead on the right side. The only difference to
the left lead is that the vector of the translational
symmetry must point to the right, the remaining code is the same:
.. literalinclude:: /code/include/quantum_wire.py
:start-after: #HIDDEN_BEGIN_xhqc
:end-before: #HIDDEN_END_xhqc
.. jupyter-execute::
:hide-output:
# Then the lead to the right
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)
Note that here we added points with x-coordinate 0, just as for the left lead.
You might object that the right lead should be placed `L`
......@@ -187,13 +261,9 @@ you do not need to worry about that.
Now we have finished building our system! We plot it, to make sure we didn't
make any mistakes:
.. literalinclude:: /code/include/quantum_wire.py
:start-after: #HIDDEN_BEGIN_wsgh
:end-before: #HIDDEN_END_wsgh
This should bring up this picture:
.. jupyter-execute::
.. image:: /code/figure/quantum_wire_syst.*
kwant.plot(syst);
The system is represented in the usual way for tight-binding systems:
dots represent the lattice points `(i, j)`, and for every
......@@ -203,16 +273,29 @@ fading color.
In order to use our system for a transport calculation, we need to finalize it
.. literalinclude:: /code/include/quantum_wire.py
:start-after: #HIDDEN_BEGIN_dngj
:end-before: #HIDDEN_END_dngj
.. jupyter-execute::
# Finalize the system
syst = syst.finalized()
Having successfully created a system, we now can immediately start to compute
its conductance as a function of energy:
.. literalinclude:: /code/include/quantum_wire.py
:start-after: #HIDDEN_BEGIN_buzn
:end-before: #HIDDEN_END_buzn
.. jupyter-execute::
# Now that we have the system, we can compute conductance
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))
We use ``kwant.smatrix`` which is a short name for
`kwant.solvers.default.smatrix` of the default solver module
......@@ -227,13 +310,15 @@ Finally we can use ``matplotlib`` to make a plot of the computed data
(although writing to file and using an external viewer such as
gnuplot or xmgrace is just as viable)
.. literalinclude:: /code/include/quantum_wire.py
:start-after: #HIDDEN_BEGIN_lliv
:end-before: #HIDDEN_END_lliv
.. jupyter-execute::
This should yield the result
.. image:: /code/figure/quantum_wire_result.*
# Use matplotlib to write output
# We should see conductance steps
pyplot.figure()
pyplot.plot(energies, data)
pyplot.xlabel("energy [t]")
pyplot.ylabel("conductance [e^2/h]")
pyplot.show()
We see a conductance quantized in units of :math:`e^2/h`,
increasing in steps as the energy is increased. The
......@@ -273,13 +358,13 @@ subbands that increases with energy.
**sites** as indices. Sites themselves have a certain type, and
belong to a **site family**. A site family is also used to convert
something that represents a site (like a tuple) into a
proper `~kwant.builder.Site` object that can be used with
proper `~kwant.system.Site` object that can be used with
`~kwant.builder.Builder`.
In the above example, `lat` is the site family. ``lat(i, j)``
then translates the description of a lattice site in terms of two
integer indices (which is the natural way to do here) into
a proper `~kwant.builder.Site` object.
a proper `~kwant.system.Site` object.
The concept of site families and sites allows `~kwant.builder.Builder`
to mix arbitrary lattices and site families
......@@ -333,7 +418,12 @@ Building the same system with less code
.. seealso::
The complete source code of this example can be found in
:download:`quantum_wire_revisited.py </code/download/quantum_wire_revisited.py>`
:jupyter-download:script:`quantum_wire_revisited`
.. jupyter-kernel::
:id: quantum_wire_revisited
Kwant allows for more than one way to build a system. The reason is that
`~kwant.builder.Builder` is essentially just a container that can be filled in
......@@ -347,19 +437,57 @@ the code into separate entities. In this example we therefore also aim at more
structure.
We begin the program collecting all imports in the beginning of the
file and put the build-up of the system into a separate function
`make_system`:
file and defining the a square lattice and empty scattering region.
.. jupyter-execute::
:hide-code:
# 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.
.. jupyter-execute::
import kwant
# For plotting
from matplotlib import pyplot
.. jupyter-execute:: boilerplate.py
:hide-code:
.. jupyter-execute::
a = 1
t = 1.0
W, L = 10, 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).
# Each lattice site has 1 degree of freedom, hence norbs=1.
lat = kwant.lattice.square(a, norbs=1)
syst = kwant.Builder()
.. literalinclude:: /code/include/quantum_wire_revisited.py
:start-after: #HIDDEN_BEGIN_xkzy
:end-before: #HIDDEN_END_xkzy
Previously, the scattering region was build using two ``for``-loops.
Instead, we now write:
.. literalinclude:: /code/include/quantum_wire_revisited.py
:start-after: #HIDDEN_BEGIN_vvjt
:end-before: #HIDDEN_END_vvjt
.. jupyter-execute::
syst[(lat(x, y) for x in range(L) for y in range(W))] = 4 * t
Here, all lattice points are added at once in the first line. The
construct ``((i, j) for i in range(L) for j in range(W))`` is a
......@@ -375,9 +503,9 @@ hoppings. In this case, an iterable like for the lattice
points becomes a bit cumbersome, and we use instead another
feature of Kwant:
.. literalinclude:: /code/include/quantum_wire_revisited.py
:start-after: #HIDDEN_BEGIN_nooi
:end-before: #HIDDEN_END_nooi
.. jupyter-execute::
syst[lat.neighbors()] = -t
In regular lattices, hoppings form large groups such that hoppings within a
group can be transformed into one another by lattice translations. In order to
......@@ -397,49 +525,48 @@ values for all the nth-nearest neighbors at once, one can similarly use
The left lead is constructed in an analogous way:
.. literalinclude:: /code/include/quantum_wire_revisited.py
:start-after: #HIDDEN_BEGIN_iepx
:end-before: #HIDDEN_END_iepx
.. jupyter-execute::
lead = kwant.Builder(kwant.TranslationalSymmetry((-a, 0)))
lead[(lat(0, j) for j in range(W))] = 4 * t
lead[lat.neighbors()] = -t
The previous example duplicated almost identical code for the left and
the right lead. The only difference was the direction of the translational
symmetry vector. Here, we only construct the left lead, and use the method
`~kwant.builder.Builder.reversed` of `~kwant.builder.Builder` to obtain a copy
of a lead pointing in the opposite direction. Both leads are attached as
before and the finished system returned:
before:
.. literalinclude:: /code/include/quantum_wire_revisited.py
:start-after: #HIDDEN_BEGIN_yxot
:end-before: #HIDDEN_END_yxot
.. jupyter-execute::
:hide-output:
The remainder of the script has been organized into two functions. One for the
plotting of the conductance.
syst.attach_lead(lead)
syst.attach_lead(lead.reversed())
.. literalinclude:: /code/include/quantum_wire_revisited.py
:start-after: #HIDDEN_BEGIN_ayuk
:end-before: #HIDDEN_END_ayuk
The remainder of the script proceeds identically. We first finalize the system:
And one ``main`` function.
.. jupyter-execute::
.. literalinclude:: /code/include/quantum_wire_revisited.py
:start-after: #HIDDEN_BEGIN_cjel
:end-before: #HIDDEN_END_cjel
syst = syst.finalized()
Finally, we use the following standard Python construct [#]_ to execute
``main`` if the program is used as a script (i.e. executed as
``python quantum_wire_revisited.py``):
and then calculate the transmission and plot:
.. literalinclude:: /code/include/quantum_wire_revisited.py
:start-after: #HIDDEN_BEGIN_ypbj
:end-before: #HIDDEN_END_ypbj
.. jupyter-execute::
If the example, however, is imported inside Python using ``import
quantum_wire_revisted as qw``, ``main`` is not executed automatically.
Instead, you can execute it manually using ``qw.main()``. On the other
hand, you also have access to the other functions, ``make_system`` and
``plot_conductance``, and can thus play with the parameters.
energies = []
data = []
for ie in range(100):
energy = ie * 0.01
smatrix = kwant.smatrix(syst, energy)
energies.append(energy)
data.append(smatrix.transmission(1, 0))
The result of the example should be identical to the previous one.
pyplot.figure()
pyplot.plot(energies, data)
pyplot.xlabel("energy [t]")
pyplot.ylabel("conductance [e^2/h]")
pyplot.show()
.. specialnote:: Technical details
......@@ -453,11 +580,9 @@ The result of the example should be identical to the previous one.
For technical reasons it is not possible to add several points
using a tuple of sites. Hence it is worth noting
a subtle detail in
a subtle detail in:
.. literalinclude:: /code/include/quantum_wire_revisited.py
:start-after: #HIDDEN_BEGIN_vvjt
:end-before: #HIDDEN_END_vvjt
>>> syst[(lat(x, y) for x in range(L) for y in range(W))] = 4 * t
Note that ``(lat(x, y) for x in range(L) for y in range(W))`` is not
a tuple, but a generator.
......@@ -512,6 +637,132 @@ The result of the example should be identical to the previous one.
it would be impossible to distinguish whether one would like to add two
separate sites, or one hopping.
Tips for organizing your simulation scripts
...........................................
.. seealso::
The complete source code of this example can be found in
:jupyter-download:script:`quantum_wire_organized`
.. jupyter-kernel::
:id: quantum_wire_organized
.. jupyter-execute::
:hide-code:
# Tutorial 2.2.4. Organizing a simulation script
# ==============================================
#
# Physics background
# ------------------
# Conductance of a quantum wire; subbands
#
# Note: Does the same as quantum_write_revisited.py, but features
# better code organization
The above two examples illustrate some of the core features of Kwant, however
the code was presented in a style which is good for exposition, but which is
bad for making your code understandable and reusable. In this example we will
lay out some best practices for writing your own simulation scripts.
In the above examples we constructed a single Kwant system, using global variables
for parameters such as the lattice constant and the length and width of the system.
Instead, it is preferable to create a *function* that you can call, and which will
return a Kwant ``Builder``:
.. jupyter-execute::
from matplotlib import pyplot
import kwant
.. jupyter-execute:: boilerplate.py
:hide-code:
.. jupyter-execute::
def make_system(L, W, a=1, t=1.0):
lat = kwant.lattice.square(a, norbs=1)
syst = kwant.Builder()
syst[(lat(i, j) for i in range(L) for j in range(W))] = 4 * t
syst[lat.neighbors()] = -t
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
By encapsulating system creation within ``make_system`` we *document* our code
by telling readers that *this* is how we create a system, and that creating a system
depends on *these* parameters (the length and width of the system, in this case, as well
as the lattice constant and the value for the hopping parameter). By defining a function
we also ensure that we can consistently create different systems (e.g. of different sizes)
of the same type (rectangular slab).
We similarly encapsulate the part of the script that does computation and plotting into
a function ``plot_conductance``:
.. jupyter-execute::
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()
And the ``main`` function that glues together the components that we previously defined:
.. jupyter-execute::
def main():
syst = make_system(W=10, L=30)
# Check that the system looks as intended.
kwant.plot(syst)
# Finalize the system.
fsyst = syst.finalized()
# We should see conductance steps.
plot_conductance(fsyst, energies=[0.01 * i for i in range(100)])
Finally, we use the following standard Python construct [#]_ to execute
``main`` if the program is used as a script (i.e. executed as
``python quantum_wire_organized.py``):
.. jupyter-execute::
# 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()
If the example, however, is imported inside Python using ``import
quantum_wire_organized as qw``, ``main`` is not executed automatically.
Instead, you can execute it manually using ``qw.main()``. On the other
hand, you also have access to the other functions, ``make_system`` and
``plot_conductance``, and can thus play with the parameters.
The result of this example should be identical to the previous one.
.. rubric:: Footnotes
.. [#] https://docs.python.org/3/library/__main__.html
doc/source/tutorial/graphene.rst
View file @ 1d8db7bb
......@@ -3,9 +3,15 @@
Beyond square lattices: graphene
--------------------------------
.. seealso::
You can execute the code examples live in your browser by
activating thebelab:
.. thebe-button:: Activate Thebelab
.. 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 +24,44 @@ 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:: boilerplate.py
:hide-code:
.. jupyter-execute::
graphene = kwant.lattice.general([(1, 0), (sin_30, cos_30)],
[(0, 0), (0, 1 / sqrt(3))],
norbs=1)
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 +72,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 +107,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 +127,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 +185,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 +254,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::
should yield two eigenvalues equal to ``[ 3.07869311,
compute_evs(syst.finalized())
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);
It computes the band structure of one of lead 0:
We then finalize the system:
.. image:: /code/figure/graphene_bs.*
.. jupyter-execute::
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 +298,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
......
doc/source/tutorial/index.rst
View file @ 1d8db7bb
......@@ -12,4 +12,5 @@ Tutorial: learning Kwant through examples
plotting
kpm
discretize
magnetic_field
faq
doc/source/tutorial/kpm.rst
View file @ 1d8db7bb
......@@ -15,10 +15,38 @@ expansion of the density of states. It can also be used to calculate the
spectral density of arbitrary operators. Kwant has an implementation of the
KPM method `kwant.kpm`, that is based on the algorithms presented in Ref. [1]_.
.. seealso::
You can execute the code examples live in your browser by
activating thebelab:
.. thebe-button:: Activate Thebelab
.. seealso::
The complete source code of this example can be found in
:download:`kernel_polynomial_method.py </code/download/kernel_polynomial_method.py>`
:jupyter-download:script:`kernel_polynomial_method`
.. jupyter-kernel::
:id: kernel_polynomial_method
.. jupyter-execute::
:hide-code:
# 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
from matplotlib import pyplot
.. jupyter-execute:: boilerplate.py
:hide-code:
Introduction
************
......@@ -103,35 +131,104 @@ to obtain the (global) density of states of a graphene disk.
We start by importing kwant and defining our system.
.. literalinclude:: /code/include/kernel_polynomial_method.py
:start-after: #HIDDEN_BEGIN_sys1
:end-before: #HIDDEN_END_sys1
.. jupyter-execute::
# 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
.. jupyter-execute::
:hide-code:
## common plotting routines ##
# Plot several density of states curves on the same axes.
def plot_dos(labels_to_data):
pyplot.figure()
for label, (x, y) in labels_to_data:
pyplot.plot(x, y.real, label=label, linewidth=2)
pyplot.legend(loc=2, framealpha=0.5)
pyplot.xlabel("energy [t]")
pyplot.ylabel("DoS [a.u.]")
pyplot.show()
# Plot fill density of states plus curves on the same axes.
def plot_dos_and_curves(dos, labels_to_data):
pyplot.figure()
pyplot.fill_between(dos[0], dos[1], label="DoS [a.u.]",
alpha=0.5, color='gray')
for label, (x, y) in labels_to_data:
pyplot.plot(x, y, label=label, linewidth=2)
pyplot.legend(loc=2, framealpha=0.5)
pyplot.xlabel("energy [t]")
pyplot.ylabel("$σ [e^2/h]$")
pyplot.show()
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):
fig, axes = pyplot.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')
pyplot.show()
After making a system we can then create a `~kwant.kpm.SpectralDensity`
object that represents the density of states for this system.
.. literalinclude:: /code/include/kernel_polynomial_method.py
:start-after: #HIDDEN_BEGIN_kpm1
:end-before: #HIDDEN_END_kpm1
.. jupyter-execute::
fsyst = make_syst().finalized()
spectrum = kwant.kpm.SpectralDensity(fsyst, rng=0)
The `~kwant.kpm.SpectralDensity` can then be called like a function to obtain a
sequence of energies in the spectrum of the Hamiltonian, and the corresponding
density of states at these energies.
.. literalinclude:: /code/include/kernel_polynomial_method.py
:start-after: #HIDDEN_BEGIN_kpm2
:end-before: #HIDDEN_END_kpm2
.. jupyter-execute::
energies, densities = spectrum()
When called with no arguments, an optimal set of energies is chosen (these are
not evenly distributed over the spectrum, see Ref. [1]_ for details), however
it is also possible to provide an explicit sequence of energies at which to
evaluate the density of states.
.. literalinclude:: /code/include/kernel_polynomial_method.py
:start-after: #HIDDEN_BEGIN_kpm3
:end-before: #HIDDEN_END_kpm3
.. jupyter-execute::
energy_subset = np.linspace(0, 2)
density_subset = spectrum(energy_subset)
.. image:: /code/figure/kpm_dos.*
.. jupyter-execute::
:hide-code:
plot_dos([
('densities', (energies, densities)),
('density subset', (energy_subset, density_subset)),
])
In addition to being called like functions, `~kwant.kpm.SpectralDensity`
objects also have a method `~kwant.kpm.SpectralDensity.integrate` which can be
......@@ -139,22 +236,21 @@ used to integrate the density of states against some distribution function over
the whole spectrum. If no distribution function is specified, then the uniform
distribution is used:
.. literalinclude:: /code/include/kernel_polynomial_method.py
:start-after: #HIDDEN_BEGIN_int1
:end-before: #HIDDEN_END_int1
.. jupyter-execute::
.. literalinclude:: /code/figure/kpm_normalization.txt
print('identity resolution:', spectrum.integrate())
We see that the integral of the density of states is normalized to the total
number of available states in the system. If we wish to calculate, say, the
number of states populated in equilibrium, then we should integrate with
respect to a Fermi-Dirac distribution:
.. literalinclude:: /code/include/kernel_polynomial_method.py
:start-after: #HIDDEN_BEGIN_int2
:end-before: #HIDDEN_END_int2
.. jupyter-execute::
# Fermi energy 0.1 and temperature 0.2
fermi = lambda E: 1 / (np.exp((E - 0.1) / 0.2) + 1)
.. literalinclude:: /code/figure/kpm_total_states.txt
print('number of filled states:', spectrum.integrate(fermi))
.. specialnote:: Stability and performance: spectral bounds
......@@ -194,23 +290,52 @@ In the following example, we compute the local density of states at the center
of the graphene disk, and we add a staggering potential between the two
sublattices.
.. literalinclude:: /code/include/kernel_polynomial_method.py
:start-after: #HIDDEN_BEGIN_sys3
:end-before: #HIDDEN_END_sys3
.. jupyter-execute::
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
Next, we choose one site of each sublattice "A" and "B",
.. literalinclude:: /code/include/kernel_polynomial_method.py
:start-after: #HIDDEN_BEGIN_op5
:end-before: #HIDDEN_END_op5
.. jupyter-execute::
fsyst_staggered = make_syst_staggered().finalized()
# 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_staggered, where)
and plot their respective local density of states.
.. literalinclude:: /code/include/kernel_polynomial_method.py
:start-after: #HIDDEN_BEGIN_op6
:end-before: #HIDDEN_END_op6
.. jupyter-execute::
.. image:: /code/figure/kpm_ldos_sites.*
# 'num_vectors' can be unspecified when using 'LocalVectors'
local_dos = kwant.kpm.SpectralDensity(fsyst_staggered, num_vectors=None,
vector_factory=vector_factory,
mean=False,
rng=0)
energies, densities = local_dos()
.. jupyter-execute::
:hide-code:
plot_dos([
('A sublattice', (energies, densities[:, 0])),
('B sublattice', (energies, densities[:, 1])),
])
Note that there is no noise comming from the random vectors.
......@@ -225,9 +350,15 @@ exactly is changed.
The simplest way to obtain a more accurate solution is to use the
``add_moments`` method:
.. literalinclude:: /code/include/kernel_polynomial_method.py
:start-after: #HIDDEN_BEGIN_acc1
:end-before: #HIDDEN_END_acc1
.. jupyter-execute::
:hide-code:
spectrum = kwant.kpm.SpectralDensity(fsyst, rng=0)
original_dos = spectrum()
.. jupyter-execute::
spectrum.add_moments(energy_resolution=0.03)
This will update the number of calculated moments and also the default
number of sampling points such that the maximum distance between successive
......@@ -237,12 +368,19 @@ Alternatively, you can directly increase the number of moments
with ``add_moments``, or the number of random vectors with ``add_vectors``.
The added vectors will be generated with the ``vector_factory``.
.. literalinclude:: /code/include/kernel_polynomial_method.py
:start-after: #HIDDEN_BEGIN_acc2
:end-before: #HIDDEN_END_acc2
.. jupyter-execute::
spectrum.add_moments(100)
spectrum.add_vectors(5)
.. image:: /code/figure/kpm_dos_r.*
.. jupyter-execute::
:hide-code:
increased_moments_dos = spectrum()
plot_dos([
('density', original_dos),
('higher number of moments', increased_moments_dos),
])
.. _operator_spectral_density:
......@@ -264,17 +402,19 @@ the spectral density of the given operator that is calculated.
If an explicit matrix is provided then it must have the same
shape as the system Hamiltonian.
.. literalinclude:: /code/include/kernel_polynomial_method.py
:start-after: #HIDDEN_BEGIN_op1
:end-before: #HIDDEN_END_op1
.. jupyter-execute::
# identity matrix
matrix_op = scipy.sparse.eye(len(fsyst.sites))
matrix_spectrum = kwant.kpm.SpectralDensity(fsyst, operator=matrix_op, rng=0)
Or, to do the same calculation using `kwant.operator.Density`:
.. literalinclude:: /code/include/kernel_polynomial_method.py
:start-after: #HIDDEN_BEGIN_op2
:end-before: #HIDDEN_END_op2
.. jupyter-execute::
# '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)
Spectral density with random vectors
^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
......@@ -283,9 +423,11 @@ Using operators from `kwant.operator` allows us to calculate quantities
such as the *local* density of states by telling the operator not to
sum over all the sites of the system:
.. literalinclude:: /code/include/kernel_polynomial_method.py
:start-after: #HIDDEN_BEGIN_op3
:end-before: #HIDDEN_END_op3
.. jupyter-execute::
# '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)
`~kwant.kpm.SpectralDensity` will properly handle this vector output,
and will average the local density obtained with random vectors.
......@@ -294,12 +436,18 @@ The accuracy of this approximation depends on the number of random vectors used.
This allows us to plot an approximate local density of states at different
points in the spectrum:
.. literalinclude:: /code/include/kernel_polynomial_method.py
:start-after: #HIDDEN_BEGIN_op4
:end-before: #HIDDEN_END_op4
.. jupyter-execute::
zero_energy_ldos = local_dos(energy=0)
finite_energy_ldos = local_dos(energy=1)
.. image:: /code/figure/kpm_ldos.*
.. jupyter-execute::
:hide-code:
plot_ldos(fsyst, [
('energy = 0', zero_energy_ldos),
('energy = 1', finite_energy_ldos)
])
Calculating Kubo conductivity
*****************************
......@@ -325,17 +473,55 @@ with nearest neigbours hoppings. To turn it into a topological
insulator we add imaginary second nearest neigbours hoppings, where
the sign changes for each sublattice.
.. literalinclude:: /code/include/kernel_polynomial_method.py
:start-after: #HIDDEN_BEGIN_sys2
:end-before: #HIDDEN_END_sys2
.. jupyter-execute::
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()
To calculate the bulk conductivity, we will select sites in the unit cell
in the middle of the sample, and create a vector factory that outputs local
vectors
.. literalinclude:: /code/include/kernel_polynomial_method.py
:start-after: #HIDDEN_BEGIN_cond
:end-before: #HIDDEN_END_cond
.. jupyter-execute::
# construct the Haldane model
lat, fsyst_topo = 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_topo, where)
cond_xx = kwant.kpm.conductivity(fsyst_topo, alpha='x', beta='x', mean=True,
num_vectors=None, vector_factory=s_factory,
rng=0)
# component 'xy'
s_factory = kwant.kpm.LocalVectors(fsyst_topo, where)
cond_xy = kwant.kpm.conductivity(fsyst_topo, 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
Note that the Kubo conductivity must be normalized with the area covered
by the vectors. In this case, each local vector represents a site, and
......@@ -345,8 +531,23 @@ value of the conductivity over large parts of the system. In this
case, the area that normalizes the result, is the area covered by
the random vectors.
.. image:: /code/figure/kpm_cond.*
.. jupyter-execute::
:hide-code:
s_factory = kwant.kpm.LocalVectors(fsyst_topo, where)
spectrum = kwant.kpm.SpectralDensity(fsyst_topo, 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.real / 4)),
(r'Hall conductivity $\sigma_{xy}$',
(energies, cond_array_xy.real))],
)
.. _advanced_topics:
......@@ -370,9 +571,17 @@ To change how the random vectors are generated, you need only specify a
function that takes the dimension of the Hilbert space as a single parameter,
and which returns a vector in that Hilbert space:
.. literalinclude:: /code/include/kernel_polynomial_method.py
:start-after: #HIDDEN_BEGIN_fact1
:end-before: #HIDDEN_END_fact1
.. jupyter-execute::
# 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)
Aditionally, a `~kwant.kpm.LocalVectors` generator is also available, that
returns local vectors that correspond to the sites passed. Note that
......@@ -407,9 +616,16 @@ first argument, and linear in its second argument. Below, we compare two
methods for computing the local density of states, one using
`kwant.operator.Density`, and the other using a custom function.
.. literalinclude:: /code/include/kernel_polynomial_method.py
:start-after: #HIDDEN_BEGIN_blm
:end-before: #HIDDEN_END_blm
.. jupyter-execute::
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)
__ operator_spectral_density_
......
doc/source/tutorial/magnetic_field.rst 0 → 100644
View file @ 1d8db7bb
Adding magnetic field
---------------------
Computing Landau levels in a harmonic oscillator basis
......................................................
.. seealso::
You can execute the code examples live in your browser by
activating thebelab:
.. thebe-button:: Activate Thebelab
.. seealso::
The complete source code of this example can be found in
:jupyter-download:script:`landau-levels`
.. jupyter-kernel::
:id: landau-levels
When electrons move in an external magnetic field, their motion perpendicular
to the field direction is quantized into discrete Landau levels. Kwant implements
an efficient scheme for computing the Landau levels of arbitrary continuum
Hamiltonians. The general scheme revolves around rewriting the Hamiltonian in terms
of a basis of harmonic oscillator states [#]_, and is commonly illustrated in textbooks
for quadratic Hamiltonians.
To demonstrate the general scheme, let us consider a magnetic field oriented along
the :math:`z` direction :math:`\vec{B} = (0, 0, B)`, such that electron motion
in the :math:`xy` plane is Landau quantized. The magnetic field enters the Hamiltonian
through the kinetic momentum
.. math:: \hbar \vec{k} = - i \hbar \nabla + e\vec{A}(x, y).
In the symmetric gauge :math:`\vec{A}(x, y) = (B/2)[-y, x, 0]`, we introduce ladder
operators with the substitution
.. math::
k_x = \frac{1}{\sqrt{2} l_B} (a + a^\dagger), \quad \quad
k_y = \frac{i}{\sqrt{2} l_B} (a - a^\dagger),
with the magnetic length :math:`l_B = \sqrt{\hbar/eB}`. The ladder operators obey the
commutation relation
.. math:: [a, a^\dagger] = 1,
and define a quantum harmonic oscillator. We can thus write any electron continuum
Hamiltonian in terms of :math:`a` and :math:`a^\dagger`. Such a Hamiltonian has a
simple matrix representation in the eigenbasis of the number operator :math:`a^\dagger a`.
The eigenstates satisfy :math:`a^\dagger a | n \rangle = n | n \rangle` with the integer
Landau level index :math:`n \geq 0`, and in coordinate representation are proportional to
.. math::
\psi_n (x, y) = \left( \frac{\partial}{ \partial w} - \frac{w^*}{4 l_B^2} \right)
w^n e^{-|w|^2/4l_B^2},
with :math:`w = x + i y`. The matrix elements of the ladder operators are
.. math::
\langle n | a | m \rangle = \sqrt{m}~\delta_{n, m-1}, \quad \quad
\langle n | a^\dagger | m \rangle = \sqrt{m + 1}~\delta_{n, m+1}.
Truncating the basis to the first :math:`N` Landau levels allows us to approximate
the Hamiltonian as a discrete, finite matrix.
We can now formulate the algorithm that Kwant uses to compute Landau levels.
1. We take a generic continuum Hamiltonian, written in terms of the kinetic
momentum :math:`\vec{k}`. The Hamiltonian must be translationally
invariant along the directions perpendicular to the field direction.
2. We substitute the momenta perpendicular to the magnetic field with the ladder
operators :math:`a` and :math:`a^\dagger`.
3. We construct a `kwant.builder.Builder` using a special lattice which includes
the Landau level index :math:`n` as a degree of freedom on each site. The directions
normal to the field direction are not included in the builder, because they are
encoded in the Landau level index.
This procedure is automated with `kwant.continuum.discretize_landau`.
As an example, let us take the Bernevig-Hughes-Zhang model that we first considered in the
discretizer tutorial ":ref:`discretize-bhz-model`":
.. math::
C + M σ_0 \otimes σ_z + F(k_x^2 + k_y^2) σ_0 \otimes σ_z + D(k_x^2 + k_y^2) σ_0 \otimes σ_0
+ A k_x σ_z \otimes σ_x + A k_y σ_0 \otimes σ_y.
We can discretize this Hamiltonian in a basis of Landau levels as follows
.. jupyter-execute::
import numpy as np
import scipy.linalg
from matplotlib import pyplot
import kwant
import kwant.continuum
.. jupyter-execute:: boilerplate.py
:hide-code:
.. jupyter-execute::
hamiltonian = """
+ C * identity(4) + M * kron(sigma_0, sigma_z)
- F * (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)
"""
syst = kwant.continuum.discretize_landau(hamiltonian, N=10)
syst = syst.finalized()
We can then plot the spectrum of the system as a function of magnetic field, and
observe a crossing of Landau levels at finite magnetic field near zero energy,
characteristic of a quantum spin Hall insulator with band inversion.
.. jupyter-execute::
params = dict(A=3.645, F =-68.6, D=-51.2, M=-0.01, C=0)
b_values = np.linspace(0.0001, 0.0004, 200)
fig = kwant.plotter.spectrum(syst, ('B', b_values), params=params, show=False)
pyplot.ylim(-0.1, 0.2);
Comparing with tight-binding
============================
In the limit where fewer than one quantum of flux is threaded through a plaquette of
the discretization lattice we can compare the discretization in Landau levels with
a discretization in realspace.
.. jupyter-execute::
lat = kwant.lattice.square(norbs=1)
syst = kwant.Builder(kwant.TranslationalSymmetry((-1, 0)))
def peierls(to_site, from_site, B):
y = from_site.tag[1]
return -1 * np.exp(-1j * B * y)
syst[(lat(0, j) for j in range(-19, 20))] = 4
syst[lat.neighbors()] = -1
syst[kwant.HoppingKind((1, 0), lat)] = peierls
syst = syst.finalized()
landau_syst = kwant.continuum.discretize_landau("k_x**2 + k_y**2", N=5)
landau_syst = landau_syst.finalized()
Here we plot the dispersion relation for the discretized ribbon and compare it
with the Landau levels shown as dashed lines.
.. jupyter-execute::
fig, ax = pyplot.subplots(1, 1)
ax.set_xlabel("momentum")
ax.set_ylabel("energy")
ax.set_ylim(0, 1)
params = dict(B=0.1)
kwant.plotter.bands(syst, ax=ax, params=params)
h = landau_syst.hamiltonian_submatrix(params=params)
for ev in scipy.linalg.eigvalsh(h):
ax.axhline(ev, linestyle='--')
The dispersion and the Landau levels diverge with increasing energy, because the real space
discretization of the ribbon gives a worse approximation to the dispersion at higher energies.
Discretizing 3D models
======================
Although the preceding examples have only included the plane perpendicular to the
magnetic field, the Landau level quantization also works if the direction
parallel to the field is included. In fact, although the system must be
translationally invariant in the plane perpendicular to the field, the system may
be arbitrary along the parallel direction. For example, it is therefore possible to
model a heterostructure and/or set up a scattering problem along the field direction.
Let's say that we wish to to model a heterostructure with a varying potential
:math:`V` along the direction of a magnetic field, :math:`z`, that includes
Zeeman splitting and Rashba spin-orbit coupling:
.. math::
\frac{\hbar^2}{2m}\sigma_0(k_x^2 + k_y^2 + k_z^2)
+ V(z)\sigma_0
+ \frac{\mu_B B}{2}\sigma_z
+ \hbar\alpha(\sigma_x k_y - \sigma_y k_x).
We can discretize this Hamiltonian in a basis of Landau levels as before:
.. jupyter-execute::
continuum_hamiltonian = """
(k_x**2 + k_y**2 + k_z**2) * sigma_0
+ V(z) * sigma_0
+ mu * B * sigma_z / 2
+ alpha * (sigma_x * k_y - sigma_y * k_x)
"""
template = kwant.continuum.discretize_landau(continuum_hamiltonian, N=10)
This creates a system with a single translational symmetry, along
the :math:`z` direction, which we can use as a template
to construct our heterostructure:
.. jupyter-execute::
def hetero_structure(site):
z, = site.pos
return 0 <= z < 10
def hetero_potential(z):
if z < 2:
return 0
elif z < 7:
return 0.5
else:
return 0.7
syst = kwant.Builder()
syst.fill(template, hetero_structure, (0,))
syst = syst.finalized()
params = dict(
B=0.5,
mu=0.2,
alpha=0.4,
V=hetero_potential,
)
syst.hamiltonian_submatrix(params=params);
.. rubric:: Footnotes
.. [#] `Wikipedia <https://en.wikipedia.org/wiki/Landau_quantization>`_ has
a nice introduction to Landau quantization.
doc/source/tutorial/operators.rst
View file @ 1d8db7bb
......@@ -12,9 +12,53 @@ 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:`magnetic_texture.py </code/download/magnetic_texture.py>`
You can execute the code examples live in your browser by
activating thebelab:
.. thebe-button:: Activate Thebelab
.. seealso::
The complete source code of this example can be found in
:jupyter-download:script:`magnetic_texture`
.. jupyter-kernel::
:id: magnetic_texture
.. jupyter-execute::
:hide-code:
# 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)
.. jupyter-execute:: boilerplate.py
:hide-code:
Introduction
------------
......@@ -47,28 +91,119 @@ 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:: /code/include/magnetic_texture.py
:start-after: #HIDDEN_BEGIN_model
:end-before: #HIDDEN_END_model
.. jupyter-execute::
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
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:: /code/include/magnetic_texture.py
:start-after: #HIDDEN_BEGIN_syst
:end-before: #HIDDEN_END_syst
.. jupyter-execute::
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
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:: /code/figure/mag_field_direction.*
.. jupyter-execute::
:hide-code:
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), figsize=(13, 10))
for ax, (title, rho) in zip(axes, densities):
kwant.plotter.density(syst, rho, ax=ax)
ax.set_title(title)
plt.show()
def plot_currents(syst, currents):
fig, axes = plt.subplots(1, len(currents), figsize=(13, 10))
if not hasattr(axes, '__len__'):
axes = (axes,)
for ax, (title, current) in zip(axes, currents):
kwant.plotter.current(syst, current, ax=ax, colorbar=False,
fig_size=(13, 10))
ax.set_title(title)
plt.show()
.. jupyter-execute::
:hide-code:
syst = make_system().finalized()
.. jupyter-execute::
:hide-code:
plot_vector_field(syst, dict(r0=20, delta=10))
We will now be interested in analyzing the form of the scattering states
that originate from the left lead:
.. literalinclude:: /code/include/magnetic_texture.py
:start-after: #HIDDEN_BEGIN_wavefunction
:end-before: #HIDDEN_END_wavefunction
.. jupyter-execute::
params = dict(r0=20, delta=10, J=1)
wf = kwant.wave_function(syst, energy=-1, params=params)
psi = wf(0)[0]
Local densities
---------------
......@@ -82,34 +217,45 @@ 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:: /code/include/magnetic_texture.py
:start-after: #HIDDEN_BEGIN_ldos
:end-before: #HIDDEN_END_ldos
.. jupyter-execute::
# even (odd) indices correspond to spin up (down)
up, down = psi[::2], psi[1::2]
density = np.abs(up)**2 + np.abs(down)**2
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:: /code/include/magnetic_texture.py
:start-after: #HIDDEN_BEGIN_lsdz
:end-before: #HIDDEN_END_lsdz
.. jupyter-execute::
# spin down components have a minus sign
spin_z = np.abs(up)**2 - np.abs(down)**2
If we wanted instead to calculate the local y-projected spin density, we would
need to use an even more complicated expression:
.. literalinclude:: /code/include/magnetic_texture.py
:start-after: #HIDDEN_BEGIN_lsdy
:end-before: #HIDDEN_END_lsdy
.. jupyter-execute::
# spin down components have a minus sign
spin_y = 1j * (down.conjugate() * up - up.conjugate() * down)
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:: /code/include/magnetic_texture.py
:start-after: #HIDDEN_BEGIN_lden
:end-before: #HIDDEN_END_lden
.. jupyter-execute::
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)
`~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
......@@ -126,8 +272,14 @@ 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.density`:
.. image:: /code/figure/spin_densities.*
.. jupyter-execute::
:hide-code:
plot_densities(syst, [
('$σ_0$', density),
('$σ_z$', spin_z),
('$σ_y$', spin_y),
])
.. specialnote:: Technical Details
......@@ -180,14 +332,27 @@ 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:: /code/include/magnetic_texture.py
:start-after: #HIDDEN_BEGIN_current
:end-before: #HIDDEN_END_current
.. jupyter-execute::
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)
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:: /code/figure/spin_currents.*
.. jupyter-execute::
plot_currents(syst, [
('$J_{σ_0}$', current),
('$J_{σ_z}$', spin_z_current),
('$J_{σ_y}$', spin_y_current),
])
.. note::
......@@ -262,11 +427,18 @@ scattering region.
Doing this is as simple as passing a *function* when instantiating
the `~kwant.operator.Current`, instead of a constant matrix:
.. literalinclude:: /code/include/magnetic_texture.py
:start-after: #HIDDEN_BEGIN_following
:end-before: #HIDDEN_END_following
.. jupyter-execute::
def following_m_i(site, r0, delta):
m_i = field_direction(site.pos, r0, delta)
return np.dot(m_i, sigma)
The function must take a `~kwant.builder.Site` as its first parameter,
J_m = kwant.operator.Current(syst, following_m_i)
# evaluate the operator
m_current = J_m(psi, params=dict(r0=25, delta=10))
The function must take a `~kwant.system.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.
......@@ -274,7 +446,7 @@ 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
``following_operator`` function via the ``params`` 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,
......@@ -286,7 +458,13 @@ 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:: /code/figure/spin_current_comparison.*
.. jupyter-execute::
:hide-code:
plot_currents(syst, [
(r'$J_{\mathbf{m}_i}$', m_current),
('$J_{σ_z}$', spin_z_current),
])
.. note:: Although this example used exclusively `~kwant.operator.Current`,
you can do the same with `~kwant.operator.Density`.
......@@ -308,11 +486,16 @@ Density of states in a circle
To calculate the density of states inside a circle of radius
20 we can simply do:
.. literalinclude:: /code/include/magnetic_texture.py
:start-after: #HIDDEN_BEGIN_density_cut
:end-before: #HIDDEN_END_density_cut
.. jupyter-execute::
def circle(site):
return np.linalg.norm(site.pos) < 20
.. literalinclude:: /code/figure/circle_dos.txt
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)
note that we also provide ``sum=True``, which means that evaluating the
operator on a wavefunction will produce a single scalar. This is semantically
......@@ -325,18 +508,29 @@ 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:: /code/include/magnetic_texture.py
:start-after: #HIDDEN_BEGIN_current_cut
:end-before: #HIDDEN_END_current_cut
.. jupyter-execute::
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
.. literalinclude:: /code/figure/current_cut.txt
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))
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`),
.. note:: ``where`` can also be provided as a sequence of `~kwant.system.Site`
or a sequence of hoppings (i.e. pairs of `~kwant.system.Site`),
rather than a function.
......@@ -360,8 +554,23 @@ This can be achieved with `~kwant.operator.Current.bind`:
*different* set of parameters. This will almost certainly give
incorrect results.
.. literalinclude:: /code/include/magnetic_texture.py
:start-after: #HIDDEN_BEGIN_bind
:end-before: #HIDDEN_END_bind
.. jupyter-execute::
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)
.. jupyter-execute::
:hide-code:
.. image:: /code/figure/bound_current.*
plot_currents(syst, [
(r'$J_{\mathbf{m}_i}$ (from left)', J_m_left),
(r'$J_{σ_z}$ (from left)', J_z_left),
])
doc/source/tutorial/plotting.rst
View file @ 1d8db7bb
......@@ -10,17 +10,66 @@ these options can be used to achieve various very different objectives.
2D example: graphene quantum dot
................................
.. seealso::
You can execute the code examples live in your browser by
activating thebelab:
.. thebe-button:: Activate Thebelab
.. seealso::
The complete source code of this example can be found in
:download:`plot_graphene.py </code/download/plot_graphene.py>`
:jupyter-download:script:`plot_graphene`
.. jupyter-kernel::
:id: plot_graphene
.. jupyter-execute::
:hide-code:
# Tutorial 2.8.1. 2D example: graphene quantum dot
# ================================================
#
# Physics background
# ------------------
# - graphene edge states
#
# Kwant features highlighted
# --------------------------
# - demonstrate different ways of plotting
import warnings
warnings.simplefilter("ignore")
from matplotlib import pyplot
import kwant
.. jupyter-execute:: boilerplate.py
:hide-code:
We begin by first considering a circular graphene quantum dot (similar to what
has been used in parts of the tutorial :ref:`tutorial-graphene`.) In contrast
to previous examples, we will also use hoppings beyond next-nearest neighbors:
.. literalinclude:: /code/include/plot_graphene.py
:start-after: #HIDDEN_BEGIN_makesyst
:end-before: #HIDDEN_END_makesyst
.. jupyter-execute::
lat = kwant.lattice.honeycomb(norbs=1)
a, b = lat.sublattices
def make_system(r=8, t=-1, tp=-0.1):
def circle(pos):
x, y = pos
return x**2 + y**2 < r**2
syst = kwant.Builder()
syst[lat.shape(circle, (0, 0))] = 0
syst[lat.neighbors()] = t
syst.eradicate_dangling()
if tp:
syst[lat.neighbors(2)] = tp
return syst
Note that adding hoppings hoppings to the `n`-th nearest neighbors can be
simply done by passing `n` as an argument to
......@@ -29,16 +78,14 @@ simply done by passing `n` as an argument to
out of the shape. It is necessary to do so *before* adding the
next-nearest-neighbor hopping [#]_.
Of course, the system can be plotted simply with default settings:
.. literalinclude:: /code/include/plot_graphene.py
:start-after: #HIDDEN_BEGIN_plotsyst1
:end-before: #HIDDEN_END_plotsyst1
However, due to the richer structure of the lattice, this results in a rather
Of course, the system can be plotted simply with default settings,
however, due to the richer structure of the lattice, this results in a rather
busy plot:
.. image:: /code/figure/plot_graphene_syst1.*
.. jupyter-execute::
syst = make_system()
kwant.plot(syst);
A much clearer plot can be obtained by using different colors for both
sublattices, and by having different line widths for different hoppings. This
......@@ -47,15 +94,19 @@ can be achieved by passing a function to the arguments of
must be a function taking one site as argument, for hoppings a function taking
the start end end site of hopping as arguments:
.. literalinclude:: /code/include/plot_graphene.py
:start-after: #HIDDEN_BEGIN_plotsyst2
:end-before: #HIDDEN_END_plotsyst2
.. jupyter-execute::
Note that since we are using an unfinalized Builder, a ``site`` is really an
instance of `~kwant.builder.Site`. With these adjustments we arrive at a plot
that carries the same information, but is much easier to interpret:
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
.. image:: /code/figure/plot_graphene_syst2.*
kwant.plot(syst, site_lw=0.1, site_color=family_color, hop_lw=hopping_lw);
Note that since we are using an unfinalized Builder, a ``site`` is really an
instance of `~kwant.system.Site`. With these adjustments we arrive at a plot
that carries the same information, but is much easier to interpret.
Apart from plotting the *system* itself, `~kwant.plotter.plot` can also be used
to plot *data* living on the system.
......@@ -67,23 +118,27 @@ nearest-neighbor hopping. Computing the wave functions is done in the usual
way (note that for a large-scale system, one would probably want to use sparse
linear algebra):
.. literalinclude:: /code/include/plot_graphene.py
:start-after: #HIDDEN_BEGIN_plotdata1
:end-before: #HIDDEN_END_plotdata1
In most cases, to plot the wave function probability, one wouldn't use
`~kwant.plotter.plot`, but rather `~kwant.plotter.map`. Here, we plot the
`n`-th wave function using it:
.. jupyter-execute::
.. literalinclude:: /code/include/plot_graphene.py
:start-after: #HIDDEN_BEGIN_plotdata2
:end-before: #HIDDEN_END_plotdata2
import scipy.linalg as la
syst = make_system(tp=0).finalized()
ham = syst.hamiltonian_submatrix()
evecs = la.eigh(ham)[1]
wf = abs(evecs[:, 225])**2
In most cases, to plot the wave function probability, one wouldn't use
`~kwant.plotter.plot`, but rather `~kwant.plotter.map`. Here, we plot the
`n`-th wave function using it.
This results in a standard pseudocolor plot, showing in this case (``n=225``) a
graphene edge state, i.e. a wave function mostly localized at the zigzag edges
of the quantum dot.
.. image:: /code/figure/plot_graphene_data1.*
.. jupyter-execute::
kwant.plotter.map(syst, wf, oversampling=10, cmap='gist_heat_r');
However although in general preferable, `~kwant.plotter.map` has a few
deficiencies for this small system: For example, there are a few distortions at
......@@ -99,9 +154,17 @@ the symbol shape depending on the sublattice. With a triangle pointing up and
down on the respective sublattice, the symbols used by plot fill the space
completely:
.. literalinclude:: /code/include/plot_graphene.py
:start-after: #HIDDEN_BEGIN_plotdata3
:end-before: #HIDDEN_END_plotdata3
.. jupyter-execute::
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');
Note that with ``hop_lw=0`` we deactivate plotting the hoppings (that would not
serve any purpose here). Moreover, ``site_size=0.5`` guarantees that the two
......@@ -112,33 +175,25 @@ this is interpreted as the radius of the inner circle).
Finally, note that since we are dealing with a finalized system now, a site `i`
is represented by an integer. In order to obtain the original
`~kwant.builder.Site`, ``syst.sites[i]`` can be used.
With this we arrive at
.. image:: /code/figure/plot_graphene_data2.*
with the same information as `~kwant.plotter.map`, but with a cleaner look.
`~kwant.system.Site`, ``syst.sites[i]`` can be used.
The way how data is presented of course influences what features of the data
are best visible in a given plot. With `~kwant.plotter.plot` one can easily go
beyond pseudocolor-like plots. For example, we can represent the wave function
probability using the symbols itself:
.. literalinclude:: /code/include/plot_graphene.py
:start-after: #HIDDEN_BEGIN_plotdata4
:end-before: #HIDDEN_END_plotdata4
.. jupyter-execute::
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);
Here, we choose the symbol size proportional to the wave function probability,
while the site color is transparent to also allow for overlapping symbols to be
visible. The hoppings are also plotted in order to show the underlying lattice.
With this, we arrive at
.. image:: /code/figure/plot_graphene_data3.*
which shows the edge state nature of the wave function most clearly.
.. rubric:: Footnotes
.. [#] A dangling site is defined as having only one hopping connecting it to
......@@ -150,7 +205,31 @@ which shows the edge state nature of the wave function most clearly.
.. seealso::
The complete source code of this example can be found in
:download:`plot_zincblende.py </code/download/plot_zincblende.py>`
:jupyter-download:script:`plot_zincblende`
.. jupyter-kernel::
:id: plot_zincblende
.. jupyter-execute::
:hide-code:
# Tutorial 2.8.2. 3D example: zincblende structure
# ================================================
#
# Physical background
# -------------------
# - 3D Bravais lattices
#
# Kwant features highlighted
# --------------------------
# - demonstrate different ways of plotting in 3D
from matplotlib import pyplot
import kwant
.. jupyter-execute:: boilerplate.py
:hide-code:
Zincblende is a very common crystal structure of semiconductors. It is a
face-centered cubic crystal with two inequivalent atoms in the unit cell
......@@ -159,18 +238,30 @@ structure, but two equivalent atoms per unit cell).
It is very easily generated in Kwant with `kwant.lattice.general`:
.. literalinclude:: /code/include/plot_zincblende.py
:start-after: #HIDDEN_BEGIN_zincblende1
:end-before: #HIDDEN_END_zincblende1
.. jupyter-execute::
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)],
norbs=1)
a, b = lat.sublattices
Note how we keep references to the two different sublattices for later use.
A three-dimensional structure is created as easily as in two dimensions, by
using the `~kwant.lattice.PolyatomicLattice.shape`-functionality:
.. literalinclude:: /code/include/plot_zincblende.py
:start-after: #HIDDEN_BEGIN_zincblende2
:end-before: #HIDDEN_END_zincblende2
.. jupyter-execute::
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
We restrict ourselves here to a simple cuboid, and do not bother to add real
values for onsite and hopping energies, but only the placeholder ``None`` (in a
......@@ -179,13 +270,11 @@ real calculation, several atomic orbitals would have to be considered).
`~kwant.plotter.plot` can plot 3D systems just as easily as its two-dimensional
counterparts:
.. literalinclude:: /code/include/plot_zincblende.py
:start-after: #HIDDEN_BEGIN_plot1
:end-before: #HIDDEN_END_plot1
.. jupyter-execute::
resulting in
syst = make_cuboid()
.. image:: /code/figure/plot_zincblende_syst1.*
kwant.plot(syst);
You might notice that the standard options for plotting are quite different in
3D than in 2D. For example, by default hoppings are not printed, but sites are
......@@ -207,14 +296,18 @@ Also for 3D it is possible to customize the plot. For example, we
can explicitly plot the hoppings as lines, and color sites differently
depending on the sublattice:
.. literalinclude:: /code/include/plot_zincblende.py
:start-after: #HIDDEN_BEGIN_plot2
:end-before: #HIDDEN_END_plot2
.. jupyter-execute::
which results in a 3D plot that allows to interactively (when plotted
in a window) explore the crystal structure:
syst = make_cuboid(a=1.5, b=1.5, c=1.5)
.. image:: /code/figure/plot_zincblende_syst2.*
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);
which results in a 3D plot that allows to interactively (when plotted
in a window) explore the crystal structure.
Hence, a few lines of code using Kwant allow to explore all the different
crystal lattices out there!
......
doc/source/tutorial/spectrum.rst
View file @ 1d8db7bb
......@@ -4,9 +4,40 @@ Beyond transport: Band structure and closed systems
Band structure calculations
...........................
.. seealso::
You can execute the code examples live in your browser by
activating thebelab:
.. thebe-button:: Activate Thebelab
.. seealso::
The complete source code of this example can be found in
:download:`band_structure.py </code/download/band_structure.py>`
:jupyter-download:script:`band_structure`
.. jupyter-kernel::
:id: band_structure
.. jupyter-execute::
:hide-code:
# Tutorial 2.4.1. Band structure calculations
# ===========================================
#
# Physics background
# ------------------
# band structure of a simple quantum wire in tight-binding approximation
#
# Kwant features highlighted
# --------------------------
# - Computing the band structure of a finalized lead.
import kwant
# For plotting
from matplotlib import pyplot
.. jupyter-execute:: boilerplate.py
:hide-code:
When doing transport simulations, one also often needs to know the band
structure of the leads, i.e. the energies of the propagating plane waves in the
......@@ -19,9 +50,26 @@ tight-binding wire.
Computing band structures in Kwant is easy. Just define a lead in the
usual way:
.. literalinclude:: /code/include/band_structure.py
:start-after: #HIDDEN_BEGIN_zxip
:end-before: #HIDDEN_END_zxip
.. jupyter-execute::
def make_lead(a=1, t=1.0, W=10):
# Start with an empty lead with a single square lattice
lat = kwant.lattice.square(a, norbs=1)
sym_lead = kwant.TranslationalSymmetry((-a, 0))
lead = kwant.Builder(sym_lead)
# build up one unit cell of the lead, and add the hoppings
# to the next unit cell
for j in range(W):
lead[lat(0, j)] = 4 * t
if j > 0:
lead[lat(0, j), lat(0, j - 1)] = -t
lead[lat(1, j), lat(0, j)] = -t
return lead
"Usual way" means defining a translational symmetry vector, as well
as one unit cell of the lead, and the hoppings to neighboring
......@@ -42,13 +90,24 @@ do something more profound with the dispersion relation these energies may be
calculated directly using `~kwant.physics.Bands`. For now we just plot the
bandstructure:
.. literalinclude:: /code/include/band_structure.py
:start-after: #HIDDEN_BEGIN_pejz
:end-before: #HIDDEN_END_pejz
.. jupyter-execute::
def main():
lead = make_lead().finalized()
kwant.plotter.bands(lead, show=False)
pyplot.xlabel("momentum [(lattice constant)^-1]")
pyplot.ylabel("energy [t]")
pyplot.show()
This gives the result:
.. image:: /code/figure/band_structure_result.*
.. jupyter-execute::
:hide-code:
# 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()
where we observe the cosine-like dispersion of the square lattice. Close
to ``k=0`` this agrees well with the quadratic dispersion this tight-binding
......@@ -61,7 +120,34 @@ Closed systems
.. seealso::
The complete source code of this example can be found in
:download:`closed_system.py </code/download/closed_system.py>`
:jupyter-download:script:`closed_system`
.. jupyter-kernel::
:id: closed_system
.. jupyter-execute::
:hide-code:
# Tutorial 2.4.2. Closed systems
# ==============================
#
# Physics background
# ------------------
# Fock-darwin spectrum of a quantum dot (energy spectrum in
# as a function of a magnetic field)
#
# Kwant features highlighted
# --------------------------
# - Use of `hamiltonian_submatrix` in order to obtain a Hamiltonian
# matrix.
from cmath import exp
import numpy as np
from matplotlib import pyplot
import kwant
.. jupyter-execute:: boilerplate.py
:hide-code:
Although Kwant is (currently) mainly aimed towards transport problems, it
can also easily be used to compute properties of closed systems -- after
......@@ -74,16 +160,49 @@ To compute the eigenenergies and eigenstates, we will make use of the sparse
linear algebra functionality of `SciPy <https://www.scipy.org>`_, which
interfaces the ARPACK package:
.. literalinclude:: /code/include/closed_system.py
:start-after: #HIDDEN_BEGIN_tibv
:end-before: #HIDDEN_END_tibv
.. jupyter-execute::
# For eigenvalue computation
import scipy.sparse.linalg as sla
We set up the system using the `shape`-function as in
:ref:`tutorial-abring`, but do not add any leads:
.. literalinclude:: /code/include/closed_system.py
:start-after: #HIDDEN_BEGIN_qlyd
:end-before: #HIDDEN_END_qlyd
.. jupyter-execute::
:hide-code:
a = 1
t = 1.0
r = 10
.. jupyter-execute::
def make_system(a=1, t=1.0, r=10):
lat = kwant.lattice.square(a, norbs=1)
syst = kwant.Builder()
# Define the quantum dot
def circle(pos):
(x, y) = pos
rsq = x ** 2 + y ** 2
return rsq < r ** 2
def hopx(site1, site2, B):
# The magnetic field is controlled by the parameter B
y = site1.pos[1]
return -t * exp(-1j * B * y)
syst[lat.shape(circle, (0, 0))] = 4 * t
# hoppings in x-direction
syst[kwant.builder.HoppingKind((1, 0), lat, lat)] = hopx
# hoppings in y-directions
syst[kwant.builder.HoppingKind((0, 1), lat, lat)] = -t
# It's a closed system for a change, so no leads
return syst
We add the magnetic field using a function and a global variable as we
did in the two previous tutorial. (Here, the gauge is chosen such that
......@@ -93,14 +212,40 @@ The spectrum can be obtained by diagonalizing the Hamiltonian of the
system, which in turn can be obtained from the finalized
system using `~kwant.system.System.hamiltonian_submatrix`:
.. literalinclude:: /code/include/closed_system.py
:start-after: #HIDDEN_BEGIN_yvri
:end-before: #HIDDEN_END_yvri
.. jupyter-execute::
def plot_spectrum(syst, Bfields):
energies = []
for B in Bfields:
# Obtain the Hamiltonian as a sparse matrix
ham_mat = syst.hamiltonian_submatrix(params=dict(B=B), sparse=True)
# we only calculate the 15 lowest eigenvalues
ev = sla.eigsh(ham_mat.tocsc(), k=15, sigma=0,
return_eigenvectors=False)
energies.append(ev)
pyplot.figure()
pyplot.plot(Bfields, energies)
pyplot.xlabel("magnetic field [arbitrary units]")
pyplot.ylabel("energy [t]")
pyplot.show()
Note that we use sparse linear algebra to efficiently calculate only a
few lowest eigenvalues. Finally, we obtain the result:
.. image:: /code/figure/closed_system_result.*
.. jupyter-execute::
:hide-code:
syst = make_system()
syst = syst.finalized()
# We should observe energy levels that flow towards Landau
# level energies with increasing magnetic field.
plot_spectrum(syst, [iB * 0.002 for iB in range(100)])
At zero magnetic field several energy levels are degenerate (since our
quantum dot is rather symmetric). These degeneracies are split
......@@ -110,11 +255,34 @@ Landau level energies at higher magnetic fields [#]_.
The eigenvectors are obtained very similarly, and can be plotted directly
using `~kwant.plotter.map`:
.. literalinclude:: /code/include/closed_system.py
:start-after: #HIDDEN_BEGIN_wave
:end-before: #HIDDEN_END_wave
.. jupyter-execute::
:hide-code:
def sorted_eigs(ev):
evals, evecs = ev
evals, evecs = map(np.array, zip(*sorted(zip(evals, evecs.transpose()))))
return evals, evecs.transpose()
.. jupyter-execute::
def plot_wave_function(syst, B=0.001):
# Calculate the wave functions in the system.
ham_mat = syst.hamiltonian_submatrix(sparse=True, params=dict(B=B))
evals, evecs = sorted_eigs(sla.eigsh(ham_mat.tocsc(), k=20, sigma=0))
# Plot the probability density of the 10th eigenmode.
kwant.plotter.map(syst, np.abs(evecs[:, 9])**2,
colorbar=False, oversampling=1)
.. image:: /code/figure/closed_system_eigenvector.*
.. jupyter-execute::
:hide-code:
syst = make_system(r=30)
# Plot an eigenmode of a circular dot. Here we create a larger system for
# better spatial resolution.
syst = make_system(r=30).finalized()
plot_wave_function(syst);
The last two arguments to `~kwant.plotter.map` are optional. The first prevents
a colorbar from appearing. The second, ``oversampling=1``, makes the image look
......@@ -127,11 +295,22 @@ that they can have *non-zero* local current. We can calculate the local
current due to a state by using `kwant.operator.Current` and plotting
it using `kwant.plotter.current`:
.. literalinclude:: /code/include/closed_system.py
:start-after: #HIDDEN_BEGIN_current
:end-before: #HIDDEN_END_current
.. jupyter-execute::
def plot_current(syst, B=0.001):
# Calculate the wave functions in the system.
ham_mat = syst.hamiltonian_submatrix(sparse=True, params=dict(B=B))
evals, evecs = sorted_eigs(sla.eigsh(ham_mat.tocsc(), k=20, sigma=0))
# Calculate and plot the local current of the 10th eigenmode.
J = kwant.operator.Current(syst)
current = J(evecs[:, 9], params=dict(B=B))
kwant.plotter.current(syst, current, colorbar=False)
.. jupyter-execute::
:hide-code:
.. image:: /code/figure/closed_system_current.*
plot_current(syst);
.. specialnote:: Technical details
......
doc/source/tutorial/spin_potential_shape.rst
View file @ 1d8db7bb
......@@ -9,9 +9,18 @@ very simple examples of the previous section.
Matrix structure of on-site and hopping elements
................................................
.. seealso::
You can execute the code examples live in your browser by
activating thebelab:
.. thebe-button:: Activate Thebelab
.. seealso::
The complete source code of this example can be found in
:download:`spin_orbit.py </code/download/spin_orbit.py>`
:jupyter-download:script:`spin_orbit`
.. jupyter-kernel::
:id: spin_orbit
We begin by extending the simple 2DEG-Hamiltonian by a Rashba spin-orbit
coupling and a Zeeman splitting due to an external magnetic field:
......@@ -42,25 +51,82 @@ use matrices in our program, we import the Tinyarray package. (`NumPy
<http://www.numpy.org/>`_ would work as well, but Tinyarray is much faster
for small arrays.)
.. literalinclude:: /code/include/spin_orbit.py
:start-after: #HIDDEN_BEGIN_xumz
:end-before: #HIDDEN_END_xumz
.. jupyter-execute::
:hide-code:
# Tutorial 2.3.1. Matrix structure of on-site and hopping elements
# ================================================================
#
# Physics background
# ------------------
# Gaps in quantum wires with spin-orbit coupling and Zeeman splititng,
# as theoretically predicted in
# http://prl.aps.org/abstract/PRL/v90/i25/e256601
# and (supposedly) experimentally oberved in
# http://www.nature.com/nphys/journal/v6/n5/abs/nphys1626.html
#
# Kwant features highlighted
# --------------------------
# - Numpy matrices as values in Builder
import kwant
# For plotting
from matplotlib import pyplot
.. jupyter-execute:: boilerplate.py
:hide-code:
.. jupyter-execute::
# For matrix support
import tinyarray
For convenience, we define the Pauli-matrices first (with :math:`\sigma_0` the
unit matrix):
.. literalinclude:: /code/include/spin_orbit.py
:start-after: #HIDDEN_BEGIN_hwbt
:end-before: #HIDDEN_END_hwbt
.. jupyter-execute::
# define Pauli-matrices for convenience
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]])
and we also define some other parameters useful for constructing our system:
.. jupyter-execute::
t = 1.0
alpha = 0.5
e_z = 0.08
W, L = 10, 30
Previously, we used numbers as the values of our matrix elements.
However, `~kwant.builder.Builder` also accepts matrices as values, and
we can simply write:
.. literalinclude:: /code/include/spin_orbit.py
:start-after: #HIDDEN_BEGIN_uxrm
:end-before: #HIDDEN_END_uxrm
.. jupyter-execute::
:hide-code:
lat = kwant.lattice.square(norbs=2)
syst = kwant.Builder()
.. jupyter-execute::
#### Define the scattering region. ####
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
Note that we specify ``norbs=2`` when creating the lattice, as each site
has 2 degrees of freedom associated with it, giving us 2x2 matrices as
onsite/hopping terms.
Note that the Zeeman energy adds to the onsite term, whereas the Rashba
spin-orbit term adds to the hoppings (due to the derivative operator).
Furthermore, the hoppings in x and y-direction have a different matrix
......@@ -85,14 +151,49 @@ when specifying `(1, 0)` it is not necessary to specify `(-1, 0)`),
The leads also allow for a matrix structure,
.. literalinclude:: /code/include/spin_orbit.py
:start-after: #HIDDEN_BEGIN_yliu
:end-before: #HIDDEN_END_yliu
.. jupyter-execute::
:hide-code:
#### Define the left lead. ####
lead = kwant.Builder(kwant.TranslationalSymmetry((-1, 0)))
.. jupyter-execute::
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
.. jupyter-execute::
:hide-code:
#### Attach the leads and finalize the system. ####
syst.attach_lead(lead)
syst.attach_lead(lead.reversed())
syst = syst.finalized()
The remainder of the code is unchanged, and as a result we should obtain
the following, clearly non-monotonic conductance steps:
.. image:: /code/figure/spin_orbit_result.*
.. jupyter-execute::
:hide-code:
# Compute conductance
energies=[0.01 * i - 0.3 for i in range(100)]
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()
.. specialnote:: Technical details
......@@ -129,7 +230,32 @@ Spatially dependent values through functions
.. seealso::
The complete source code of this example can be found in
:download:`quantum_well.py </code/download/quantum_well.py>`
:jupyter-download:script:`quantum_well`
.. jupyter-kernel::
:id: quantum_well
.. jupyter-execute::
:hide-code:
# Tutorial 2.3.2. Spatially dependent values through functions
# ============================================================
#
# Physics background
# ------------------
# transmission through a quantum well
#
# Kwant features highlighted
# --------------------------
# - Functions as values in Builder
import kwant
# For plotting
from matplotlib import pyplot
.. jupyter-execute:: boilerplate.py
:hide-code:
Up to now, all examples had position-independent matrix-elements
(and thus translational invariance along the wire, which
......@@ -148,22 +274,50 @@ changing the potential then implies the need to build up the system again.
Instead, we use a python *function* to define the onsite energies. We
define the potential profile of a quantum well as:
.. literalinclude:: /code/include/quantum_well.py
:start-after: #HIDDEN_BEGIN_ehso
:end-before: #HIDDEN_END_ehso
.. jupyter-execute::
W, L, L_well = 10, 30, 10
This function takes two arguments: the first of type `~kwant.builder.Site`,
def potential(site, pot):
(x, y) = site.pos
if (L - L_well) / 2 < x < (L + L_well) / 2:
return pot
else:
return 0
This function takes two arguments: the first of type `~kwant.system.Site`,
from which you can get the real-space coordinates using ``site.pos``, and the
value of the potential as the second. Note that in `potential` we can access
variables of the surrounding function: `L` and `L_well` are taken from the
namespace of `make_system`.
variables `L` and `L_well` that are defined globally.
Kwant now allows us to pass a function as a value to
`~kwant.builder.Builder`:
.. literalinclude:: /code/include/quantum_well.py
:start-after: #HIDDEN_BEGIN_coid
:end-before: #HIDDEN_END_coid
.. jupyter-execute::
a = 1
t = 1.0
def onsite(site, pot):
return 4 * t + potential(site, pot)
lat = kwant.lattice.square(a, norbs=1)
syst = kwant.Builder()
syst[(lat(x, y) for x in range(L) for y in range(W))] = onsite
syst[lat.neighbors()] = -t
.. jupyter-execute::
:hide-code:
#### 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())
syst = syst.finalized()
For each lattice point, the corresponding site is then passed as the
first argument to the function `onsite`. The values of any additional
......@@ -180,9 +334,21 @@ of the lead -- this should be kept in mind.
Finally, we compute the transmission probability:
.. literalinclude:: /code/include/quantum_well.py
:start-after: #HIDDEN_BEGIN_sqvr
:end-before: #HIDDEN_END_sqvr
.. jupyter-execute::
def plot_conductance(syst, energy, welldepths):
# Compute conductance
data = []
for welldepth in welldepths:
smatrix = kwant.smatrix(syst, energy, params=dict(pot=-welldepth))
data.append(smatrix.transmission(1, 0))
pyplot.figure()
pyplot.plot(welldepths, data)
pyplot.xlabel("well depth [t]")
pyplot.ylabel("conductance [e^2/h]")
pyplot.show()
``kwant.smatrix`` allows us to specify a dictionary, `params`, that contains
the additional arguments required by the Hamiltonian matrix elements.
......@@ -191,7 +357,11 @@ of the potential well by passing the potential value (remember above
we defined our `onsite` function that takes a parameter named `pot`).
We obtain the result:
.. image:: /code/figure/quantum_well_result.*
.. jupyter-execute::
:hide-code:
plot_conductance(syst, energy=0.2,
welldepths=[0.01 * i for i in range(100)])
Starting from no potential (well depth = 0), we observe the typical
oscillatory transmission behavior through resonances in the quantum well.
......@@ -205,10 +375,10 @@ oscillatory transmission behavior through resonances in the quantum well.
.. specialnote:: Technical details
- Functions can also be used for hoppings. In this case, they take
two `~kwant.builder.Site`'s as arguments and then an arbitrary number
two `~kwant.system.Site`'s as arguments and then an arbitrary number
of additional arguments.
- Apart from the real-space position `pos`, `~kwant.builder.Site` has also an
- Apart from the real-space position `pos`, `~kwant.system.Site` has also an
attribute `tag` containing the lattice indices of the site.
.. _tutorial-abring:
......@@ -218,13 +388,44 @@ Nontrivial shapes
.. seealso::
The complete source code of this example can be found in
:download:`ab_ring.py </code/download/ab_ring.py>`
:jupyter-download:script:`ab_ring`
.. jupyter-kernel::
:id: ab_ring
.. jupyter-execute::
:hide-code:
# Tutorial 2.3.3. Nontrivial shapes
# =================================
#
# Physics background
# ------------------
# Flux-dependent transmission through a quantum ring
#
# Kwant features highlighted
# --------------------------
# - More complex shapes with lattices
# - Allows for discussion of subtleties of `attach_lead` (not in the
# example, but in the tutorial main text)
# - Modifcations of hoppings/sites after they have been added
from cmath import exp
from math import pi
import kwant
# For plotting
from matplotlib import pyplot
.. jupyter-execute:: boilerplate.py
:hide-code:
Up to now, we only dealt with simple wire geometries. Now we turn to the case
of a more complex geometry, namely transport through a quantum ring
that is pierced by a magnetic flux :math:`\Phi`:
.. image:: /code/figure/ab_ring_sketch.*
.. image:: /figure/ab_ring_sketch.*
For a flux line, it is possible to choose a gauge such that a
charged particle acquires a phase :math:`e\Phi/h` whenever it
......@@ -238,20 +439,32 @@ First, define a boolean function defining the desired shape, i.e. a function
that returns ``True`` whenever a point is inside the shape, and
``False`` otherwise:
.. literalinclude:: /code/include/ab_ring.py
:start-after: #HIDDEN_BEGIN_eusz
:end-before: #HIDDEN_END_eusz
.. jupyter-execute::
r1, r2 = 10, 20
def ring(pos):
(x, y) = pos
rsq = x ** 2 + y ** 2
return (r1 ** 2 < rsq < r2 ** 2)
Note that this function takes a real-space position as argument (not a
`~kwant.builder.Site`).
`~kwant.system.Site`).
We can now simply add all of the lattice points inside this shape at
once, using the function `~kwant.lattice.Square.shape`
provided by the lattice:
.. literalinclude:: /code/include/ab_ring.py
:start-after: #HIDDEN_BEGIN_lcak
:end-before: #HIDDEN_END_lcak
.. jupyter-execute::
a = 1
t = 1.0
lat = kwant.lattice.square(a, norbs=1)
syst = kwant.Builder()
syst[lat.shape(ring, (0, r1 + 1))] = 4 * t
syst[lat.neighbors()] = -t
Here, ``lat.shape`` takes as a second parameter a (real-space) point that is
inside the desired shape. The hoppings can still be added using
......@@ -264,9 +477,30 @@ along the branch cut in the lower arm of the ring. For this we select
all hoppings in x-direction that are of the form `(lat(1, j), lat(0, j))`
with ``j<0``:
.. literalinclude:: /code/include/ab_ring.py
:start-after: #HIDDEN_BEGIN_lvkt
:end-before: #HIDDEN_END_lvkt
.. jupyter-execute::
# In order to introduce a flux through the ring, we introduce a phase on
# the hoppings on the line cut through one of the arms. Since we want to
# change the flux without modifying the Builder instance repeatedly, we
# define the modified hoppings as a function that takes the flux as its
# parameter phi.
def hopping_phase(site1, site2, phi):
return -t * exp(1j * phi)
def crosses_branchcut(hop):
ix0, iy0 = hop[0].tag
# builder.HoppingKind with the argument (1, 0) below
# returns hoppings ordered as ((i+1, j), (i, j))
return iy0 < 0 and ix0 == 1 # ix1 == 0 then implied
# Modify only those hopings in x-direction that cross the branch cut
def hops_across_cut(syst):
for hop in kwant.builder.HoppingKind((1, 0), lat, lat)(syst):
if crosses_branchcut(hop):
yield hop
syst[hops_across_cut] = hopping_phase
Here, `crosses_branchcut` is a boolean function that returns ``True`` for
the desired hoppings. We then use again a generator (this time with
......@@ -279,9 +513,21 @@ by the parameter `phi`.
For the leads, we can also use the ``lat.shape``-functionality:
.. literalinclude:: /code/include/ab_ring.py
:start-after: #HIDDEN_BEGIN_qwgr
:end-before: #HIDDEN_END_qwgr
.. jupyter-execute::
#### Define the leads. ####
W = 10
sym_lead = kwant.TranslationalSymmetry((-a, 0))
lead = kwant.Builder(sym_lead)
def lead_shape(pos):
(x, y) = pos
return (-W / 2 < y < W / 2)
lead[lat.shape(lead_shape, (0, 0))] = 4 * t
lead[lat.neighbors()] = -t
Here, the shape must be compatible with the translational symmetry
of the lead ``sym_lead``. In particular, this means that it should extend to
......@@ -290,9 +536,12 @@ no restriction on ``x`` in ``lead_shape``) [#]_.
Attaching the leads is done as before:
.. literalinclude:: /code/include/ab_ring.py
:start-after: #HIDDEN_BEGIN_skbz
:end-before: #HIDDEN_END_skbz
.. jupyter-execute::
:hide-output:
#### Attach the leads ####
syst.attach_lead(lead)
syst.attach_lead(lead.reversed())
In fact, attaching leads seems not so simple any more for the current
structure with a scattering region very much different from the lead
......@@ -305,16 +554,40 @@ back (going opposite to the direction of the translational vector)
until it intersects the scattering region. At this intersection,
the lead is attached:
.. image:: /code/figure/ab_ring_sketch2.*
.. image:: /figure/ab_ring_sketch2.*
After the lead has been attached, the system should look like this:
.. image:: /code/figure/ab_ring_syst.*
.. jupyter-execute::
:hide-code:
kwant.plot(syst);
The computation of the conductance goes in the same fashion as before.
Finally you should get the following result:
.. image:: /code/figure/ab_ring_result.*
.. jupyter-execute::
:hide-code:
def plot_conductance(syst, energy, fluxes):
# compute conductance
normalized_fluxes = [flux / (2 * pi) for flux in fluxes]
data = []
for flux in fluxes:
smatrix = kwant.smatrix(syst, energy, params=dict(phi=flux))
data.append(smatrix.transmission(1, 0))
pyplot.figure()
pyplot.plot(normalized_fluxes, data)
pyplot.xlabel("flux [flux quantum]")
pyplot.ylabel("conductance [e^2/h]")
pyplot.show()
# We should see a conductance that is periodic with the flux quantum
plot_conductance(syst.finalized(), energy=0.15,
fluxes=[0.01 * i * 3 * 2 * pi for i in range(100)])
where one can observe the conductance oscillations with the
period of one flux quantum.
......@@ -330,7 +603,47 @@ period of one flux quantum.
becomes more apparent if we attach the leads a bit further away
from the central axis o the ring, as was done in this example:
.. image:: /code/figure/ab_ring_note1.*
.. jupyter-kernel::
:id: ab_ring_note1
.. jupyter-execute::
:hide-code:
import kwant
from matplotlib import pyplot
.. jupyter-execute:: boilerplate.py
:hide-code:
.. jupyter-execute::
:hide-code:
a = 1
t = 1.0
W = 10
r1, r2 = 10, 20
lat = kwant.lattice.square(norbs=1)
syst = kwant.Builder()
def ring(pos):
(x, y) = pos
rsq = x**2 + y**2
return ( r1**2 < rsq < r2**2)
syst[lat.shape(ring, (0, 11))] = 4 * t
syst[lat.neighbors()] = -t
sym_lead0 = kwant.TranslationalSymmetry((-a, 0))
lead0 = kwant.Builder(sym_lead0)
def lead_shape(pos):
(x, y) = pos
return (-1 < x < 1) and ( 0.5 * W < y < 1.5 * W )
lead0[lat.shape(lead_shape, (0, W))] = 4 * t
lead0[lat.neighbors()] = -t
lead1 = lead0.reversed()
syst.attach_lead(lead0)
syst.attach_lead(lead1)
kwant.plot(syst);
- Per default, `~kwant.builder.Builder.attach_lead` attaches
the lead to the "outside" of the structure, by tracing the
......@@ -345,7 +658,46 @@ period of one flux quantum.
starts the trace-back in the middle of the ring, resulting
in the lead being attached to the inner circle:
.. image:: /code/figure/ab_ring_note2.*
.. jupyter-kernel::
:id: ab_ring_note2
.. jupyter-execute::
:hide-code:
import kwant
from matplotlib import pyplot
.. jupyter-execute:: boilerplate.py
:hide-code:
.. jupyter-execute::
:hide-code:
a = 1
t = 1.0
W = 10
r1, r2 = 10, 20
lat = kwant.lattice.square(a, norbs=1)
syst = kwant.Builder()
def ring(pos):
(x, y) = pos
rsq = x**2 + y**2
return ( r1**2 < rsq < r2**2)
syst[lat.shape(ring, (0, 11))] = 4 * t
syst[lat.neighbors()] = -t
sym_lead0 = kwant.TranslationalSymmetry((-a, 0))
lead0 = kwant.Builder(sym_lead0)
def lead_shape(pos):
(x, y) = pos
return (-1 < x < 1) and ( -W/2 < y < W/2 )
lead0[lat.shape(lead_shape, (0, 0))] = 4 * t
lead0[lat.neighbors()] = -t
lead1 = lead0.reversed()
syst.attach_lead(lead0)
syst.attach_lead(lead1, lat(0, 0))
kwant.plot(syst);
Note that here the lead is treated as if it would pass over
the other arm of the ring, without intersecting it.
......
doc/source/tutorial/superconductors.rst
View file @ 1d8db7bb
Superconductors: orbital degrees of freedom, conservation laws and symmetries
-----------------------------------------------------------------------------
.. seealso::
You can execute the code examples live in your browser by
activating thebelab:
.. thebe-button:: Activate Thebelab
.. seealso::
The complete source code of this example can be found in
:download:`superconductor.py </code/download/superconductor.py>`
:jupyter-download:script:`superconductor`
.. jupyter-kernel::
:id: superconductor
.. jupyter-execute::
:hide-code:
# Tutorial 2.6. "Superconductors": orbitals, conservation laws and symmetries
# ===========================================================================
#
# Physics background
# ------------------
# - conductance of a NS-junction (Andreev reflection, superconducting gap)
#
# Kwant features highlighted
# --------------------------
# - Implementing electron and hole ("orbital") degrees of freedom
# using conservation laws.
# - Use of discrete symmetries to relate scattering states.
import kwant
from matplotlib import pyplot
import tinyarray
import numpy as np
tau_x = tinyarray.array([[0, 1], [1, 0]])
tau_y = tinyarray.array([[0, -1j], [1j, 0]])
tau_z = tinyarray.array([[1, 0], [0, -1]])
.. jupyter-execute:: boilerplate.py
:hide-code:
This example deals with superconductivity on the level of the
Bogoliubov-de Gennes (BdG) equation. In this framework, the Hamiltonian
......@@ -51,18 +90,50 @@ of freedom.
Let us consider a system that consists of a normal lead on the left,
a superconductor on the right, and a tunnel barrier in between:
.. image:: /code/figure/superconductor_transport_sketch.*
.. image:: /figure/superconductor_transport_sketch.*
We implement the BdG Hamiltonian in Kwant using a 2x2 matrix structure
for all Hamiltonian matrix elements, as we did
previously in the :ref:`spin example <tutorial_spinorbit>`. We declare
the square lattice and construct the scattering region with the following:
.. literalinclude:: /code/include/superconductor.py
:start-after: #HIDDEN_BEGIN_nbvn
:end-before: #HIDDEN_END_nbvn
Note the new argument ``norbs`` to `~kwant.lattice.square`. This is
previously in the :ref:`spin example <tutorial_spinorbit>`.
We start by declaring some parameters that will be used in the following code:
.. jupyter-execute::
a = 1
W, L = 10, 10
barrier = 1.5
barrierpos = (3, 4)
mu = 0.4
Delta = 0.1
Deltapos = 4
t = 1.0
and we declare the square lattice and construct the scattering region with the following:
.. jupyter-execute::
# 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
Note the argument ``norbs`` to `~kwant.lattice.square`. This is
the number of orbitals per site in the discretized BdG Hamiltonian - of course,
``norbs = 2``, since each site has one electron orbital and one hole orbital.
It is necessary to specify ``norbs`` here, such that we may later separate the
......@@ -80,9 +151,16 @@ of it). The next step towards computing conductance is to attach leads.
Let's attach two leads: a normal one to the left end, and a superconducting
one to the right end. Starting with the left lead, we have:
.. literalinclude:: /code/include/superconductor.py
:start-after: #HIDDEN_BEGIN_ttth
:end-before: #HIDDEN_END_ttth
.. jupyter-execute::
#### 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
Note the two new new arguments in `~kwant.builder.Builder`, ``conservation_law``
and ``particle_hole``. For the purpose of computing conductance, ``conservation_law``
......@@ -120,9 +198,19 @@ of ascending eigenvalues of the conservation law.
In order to move on with the conductance calculation, let's attach the second
lead to the right side of the scattering region:
.. literalinclude:: /code/include/superconductor.py
:start-after: #HIDDEN_BEGIN_zuuw
:end-before: #HIDDEN_END_zuuw
.. jupyter-execute::
# 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 finalize the system. ####
syst.attach_lead(lead0)
syst.attach_lead(lead1)
syst = syst.finalized()
The second (right) lead is superconducting, such that the electron and hole
blocks are coupled. Of course, this means that we can not separate them into
......@@ -140,9 +228,23 @@ confused by the fact that it says ``transmission`` -- transmission
into the same lead is reflection), and reflection from electrons to holes
is ``smatrix.transmission((0, 1), (0, 0))``:
.. literalinclude:: /code/include/superconductor.py
:start-after: #HIDDEN_BEGIN_jbjt
:end-before: #HIDDEN_END_jbjt
.. jupyter-execute::
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)))
pyplot.figure()
pyplot.plot(energies, data)
pyplot.xlabel("energy [t]")
pyplot.ylabel("conductance [e^2/h]")
pyplot.show()
Note that ``smatrix.submatrix((0, 0), (0, 0))`` returns the block concerning
reflection of electrons to electrons, and from its size we can extract the number of modes
......@@ -150,7 +252,10 @@ reflection of electrons to electrons, and from its size we can extract the numbe
For the default parameters, we obtain the following conductance:
.. image:: /code/figure/superconductor_transport_result.*
.. jupyter-execute::
:hide-code:
plot_conductance(syst, energies=[0.002 * i for i in range(-10, 100)])
We a see a conductance that is proportional to the square of the tunneling
probability within the gap, and proportional to the tunneling probability
......@@ -176,13 +281,34 @@ the electron and hole blocks. The exception is of course at zero energy, in whic
case particle-hole symmetry transforms between the electron and hole blocks, resulting
in a symmetric scattering matrix. We can check the symmetry explicitly with
.. literalinclude:: /code/include/superconductor.py
:start-after: #HIDDEN_BEGIN_pqmp
:end-before: #HIDDEN_END_pqmp
.. jupyter-execute::
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))
which yields the output
.. literalinclude:: /code/figure/check_PHS_out.txt
.. jupyter-execute::
:hide-code:
check_PHS(syst)
Note that :math:`\mathcal{P}` flips the sign of momentum, and for the parameters
we consider here, there are two electron and two hole modes active at zero energy.
......
docker/Dockerfile.conda
View file @ 1d8db7bb
......@@ -2,6 +2,7 @@ FROM conda/miniconda3:latest
ENV LANG C.UTF-8
ENV LC_ALL C.UTF-8
ENV DEBIAN_FRONTEND noninteractive
RUN apt-get update && apt-get install -y --no-install-recommends \
# all the hard non-Python dependencies
......@@ -9,7 +10,7 @@ RUN apt-get update && apt-get install -y --no-install-recommends \
# Additional tools for running CI
file rsync openssh-client \
# Documentation building
inkscape texlive-full zip \
inkscape texlive-full zip librsvg2-bin \
&& apt-get clean && \
rm -rf /var/lib/apt/lists/*
......@@ -18,3 +19,5 @@ COPY kwant-latest.yml kwant-stable.yml kwant-stable-no-extras.yml /
RUN conda env create -qf kwant-stable.yml
RUN conda env create -qf kwant-stable-no-extras.yml
RUN conda env create -qf kwant-latest.yml
RUN /usr/local/envs/kwant-latest/bin/python -m ipykernel install --user --name kwant-latest
docker/Dockerfile.debian
View file @ 1d8db7bb
......@@ -2,11 +2,12 @@ FROM debian:latest
ENV LANG C.UTF-8
ENV LC_ALL C.UTF-8
ENV DEBIAN_FRONTEND noninteractive
RUN apt-get update && apt-get install -y --no-install-recommends \
gnupg dirmngr apt-transport-https ca-certificates curl software-properties-common
RUN echo "deb http://downloads.kwant-project.org/debian/ stretch-backports main" >> /etc/apt/sources.list && \
RUN echo "deb http://downloads.kwant-project.org/debian/ stable main" >> /etc/apt/sources.list && \
apt-key adv --no-tty --keyserver pool.sks-keyservers.net --recv-key C3F147F5980F3535 && \
apt-get update && apt-get install -y --no-install-recommends \
# all the hard non-Python dependencies
......@@ -20,6 +21,10 @@ RUN echo "deb http://downloads.kwant-project.org/debian/ stretch-backports main"
&& apt-get clean && \
rm -rf /var/lib/apt/lists/*
### install optional dependencies not available from the Debian repositories
RUN pip3 install \
qsymm==1.2.6
### install build and testing dependencies
RUN pip3 install \
cython \
......
docker/Dockerfile.ubuntu
View file @ 1d8db7bb
FROM ubuntu:16.04
FROM ubuntu:18.04
ENV LANG C.UTF-8
ENV LC_ALL C.UTF-8
ENV DEBIAN_FRONTEND noninteractive
RUN apt-get update && apt-get install -y --no-install-recommends \
gnupg dirmngr apt-transport-https software-properties-common
......@@ -18,6 +19,10 @@ RUN apt-add-repository -s ppa:kwant-project/ppa && \
&& apt-get clean && \
rm -rf /var/lib/apt/lists/*
### install optional dependencies not available from the Debian repositories
RUN pip3 install \
qsymm==1.2.6
### install build and testing dependencies
RUN pip3 install \
cython \
......
docker/kwant-latest.yml
View file @ 1d8db7bb
......@@ -8,6 +8,7 @@ dependencies:
- tinyarray
- sympy
- matplotlib
- qsymm
# Linear algebra libraraies
- mumps
- blas #=1.1 openblas
......@@ -23,6 +24,9 @@ dependencies:
- pytest-flakes
- pytest-pep8
# Documentation building
- sphinx=1.7.4 # later versions seem to have problems
- sphinx
- numpydoc
- requests
- jupyter_sphinx
- pip:
- sphinxcontrib-svg2pdfconverter
docker/kwant-stable-no-extras.yml
View file @ 1d8db7bb
......@@ -2,9 +2,9 @@ name: kwant-stable-no-extras
channels:
- conda-forge
dependencies:
- python=3.5
- numpy=1.11.0
- scipy=0.17 # No minor version because conda-forge has no 0.17.0
- python=3.6
- numpy=1.13.1
- scipy=0.19.1
- tinyarray=1.2.0
# Linear algebra libraraies
- blas #=1.1 openblas
......@@ -19,7 +19,3 @@ dependencies:
- pytest-cov
- pytest-flakes
- pytest-pep8
# Documentation building
- sphinx=1.7.4 # later versions seem to have problems
- numpydoc
- requests
docker/kwant-stable.yml
View file @ 1d8db7bb
......@@ -2,12 +2,13 @@ name: kwant-stable
channels:
- conda-forge
dependencies:
- python=3.5
- numpy=1.11.0
- scipy=0.17 # No minor version because conda-forge has no 0.17.0
- python=3.6
- numpy=1.13.3
- scipy=0.19.1
- tinyarray=1.2.0
- sympy=0.7.6
- matplotlib=1.5.1
- sympy=1.1.1
- matplotlib=2.1.1
- qsymm=1.2.6
# Linear algebra libraraies
- mumps
- blas #=1.1 openblas
......@@ -22,7 +23,3 @@ dependencies:
- pytest-cov
- pytest-flakes
- pytest-pep8
# Documentation building
- sphinx=1.7.4 # later versions seem to have problems
- numpydoc
- requests
kwant/__init__.py
View file @ 1d8db7bb
......@@ -29,22 +29,29 @@ except ImportError:
raise
from ._common import KwantDeprecationWarning, UserCodeError
__all__.extend(['KwantDeprecationWarning', 'UserCodeError'])
for module in ['system', 'builder', 'lattice', 'solvers', 'digest', 'rmt',
'operator', 'kpm', 'wraparound']:
exec('from . import {0}'.format(module))
__all__.append(module)
# Pre-import most submodules. (We make exceptions for submodules that have
# special dependencies or where that would take too long.)
from . import system
from . import builder
from . import lattice
from . import solvers
from . import digest
from . import rmt
from . import operator
from . import kpm
from . import wraparound
__all__.extend(['system', 'builder', 'lattice', 'solvers', 'digest', 'rmt',
'operator', 'kpm', 'wraparound'])
# Make selected functionality available directly in the root namespace.
available = [('builder', ['Builder', 'HoppingKind']),
('lattice', ['TranslationalSymmetry']),
('solvers.default',
['smatrix', 'greens_function', 'ldos', 'wave_function'])]
for module, names in available:
exec('from .{0} import {1}'.format(module, ', '.join(names)))
__all__.extend(names)
from .builder import Builder, HoppingKind
__all__.extend(['Builder', 'HoppingKind'])
from .lattice import TranslationalSymmetry
__all__.extend(['TranslationalSymmetry'])
from .solvers.default import smatrix, greens_function, ldos, wave_function
__all__.extend(['smatrix', 'greens_function', 'ldos', 'wave_function'])
# Importing plotter might not work, but this does not have to be a problem --
# only no plotting will be available.
......
kwant/_common.py
View file @ 1d8db7bb
......@@ -13,6 +13,7 @@ import inspect
import warnings
import importlib
import functools
import collections
from contextlib import contextmanager
__all__ = ['KwantDeprecationWarning', 'UserCodeError']
......@@ -191,3 +192,27 @@ class lazy_import:
package = sys.modules[self.__package]
setattr(package, self.__module, mod)
return getattr(mod, name)
def _hashable(obj):
return isinstance(obj, collections.abc.Hashable)
def memoize(f):
"""Decorator to memoize a function that works even with unhashable args.
This decorator will even work with functions whose args are not hashable.
The cache key is made up by the hashable arguments and the ids of the
non-hashable args. It is up to the user to make sure that non-hashable
args do not change during the lifetime of the decorator.
This decorator will keep reevaluating functions that return None.
"""
def lookup(*args):
key = tuple(arg if _hashable(arg) else id(arg) for arg in args)
result = cache.get(key)
if result is None:
cache[key] = result = f(*args)
return result
cache = {}
return lookup
kwant/_system.pyx
View file @ 1d8db7bb
# Copyright 2011-2013 Kwant authors.
# Copyright 2011-2019 Kwant authors.
#
# This file is part of Kwant. It is subject to the license terms in the file
# LICENSE.rst found in the top-level directory of this distribution and at
......@@ -12,12 +12,16 @@ import numpy as np
from scipy import sparse as sp
from itertools import chain
import types
import bisect
from .graph.core cimport CGraph, gintArraySlice
from .graph.defs cimport gint
from .graph.defs import gint_dtype
from ._common import deprecate_args
### Non-vectorized methods
msg = ('Hopping from site {0} to site {1} does not match the '
'dimensions of onsite Hamiltonians of these sites.')
......@@ -352,3 +356,383 @@ def hamiltonian_submatrix(self, args=(), to_sites=None, from_sites=None,
mat = func(ham, args, params, self.graph, diag, from_sites,
n_by_to_site, to_norb, to_off, from_norb, from_off)
return (mat, to_norb, from_norb) if return_norb else mat
### Vectorized methods
_shape_error_msg = (
"The following hoppings have matrix elements of incompatible shape "
"with other matrix elements: {}"
)
@cython.boundscheck(False)
def _vectorized_make_sparse(subgraphs, hams, long [:] norbs, long [:] orb_offsets,
long [:] site_offsets, shape=None):
ndata = sum(h.shape[0] * h.shape[1] * h.shape[2] for h in hams)
cdef long [:] rows_view, cols_view
cdef complex [:] data_view
rows = np.empty((ndata,), dtype=long)
cols = np.empty((ndata,), dtype=long)
data = np.empty((ndata,), dtype=complex)
rows_view = rows
cols_view = cols
data_view = data
cdef long i, j, k, m, N, M, P, to_off, from_off,\
ta, fa, to_norbs, from_norbs
cdef long [:] ts_offs, fs_offs
cdef complex [:, :, :] h
m = 0
# This outer loop zip() is pure Python, but that's ok, as it
# has very few entries and the inner loops are fully vectorized
for ((ta, fa), (ts_offs, fs_offs)), h in zip(subgraphs, hams):
N = h.shape[0]
M = h.shape[1]
P = h.shape[2]
if norbs[ta] != M or norbs[fa] != P:
to_sites = site_offsets[ta] + np.array(ts_offs)
from_sites = site_offsets[fa] + np.array(fs_offs)
hops = np.array([to_sites, from_sites]).transpose()
raise ValueError(_shape_error_msg.format(hops))
for i in range(N):
to_off = orb_offsets[ta] + norbs[ta] * ts_offs[i]
from_off = orb_offsets[fa] + norbs[fa] * fs_offs[i]
for j in range(M):
for k in range(P):
rows_view[m] = to_off + j
cols_view[m] = from_off + k
data_view[m] = h[i, j, k]
m += 1
if shape is None:
shape = (orb_offsets[-1], orb_offsets[-1])
return sp.coo_matrix((data, (rows, cols)), shape=shape)
@cython.boundscheck(False)
def _vectorized_make_dense(subgraphs, hams, long [:] norbs, long [:] orb_offsets,
long [:] site_offsets, shape=None):
if shape is None:
shape = (orb_offsets[-1], orb_offsets[-1])
mat = np.zeros(shape, dtype=complex)
cdef complex [:, :] mat_view
mat_view = mat
cdef long i, j, k, N, M, P, to_off, from_off,\
ta, fa, to_norbs, from_norbs
cdef long [:] ts_offs, fs_offs
cdef complex [:, :, :] h
# This outer loop zip() is pure Python, but that's ok, as it
# has very few entries and the inner loops are fully vectorized
for ((ta, fa), (ts_offs, fs_offs)), h in zip(subgraphs, hams):
N = h.shape[0]
M = h.shape[1]
P = h.shape[2]
if norbs[ta] != M or norbs[fa] != P:
to_sites = site_offsets[ta] + np.array(ts_offs)
from_sites = site_offsets[fa] + np.array(fs_offs)
hops = np.array([to_sites, from_sites]).transpose()
raise ValueError(_shape_error_msg.format(hops))
for i in range(N):
to_off = orb_offsets[ta] + norbs[ta] * ts_offs[i]
from_off = orb_offsets[fa] + norbs[fa] * fs_offs[i]
for j in range(M):
for k in range(P):
mat_view[to_off + j, from_off + k] = h[i, j, k]
return mat
def _count_norbs(syst, site_offsets, hams, args=(), params=None):
"""Return the norbs and orbital offset of each site family in 'syst'
Parameters
----------
syst : `kwant.system.System`
site_offsets : sequence of int
Contains the index of the first site for each site array
in 'syst.site_arrays'.
hams : sequence of arrays or 'None'
The Hamiltonian for each term in 'syst.terms'. 'None'
if the corresponding term has not been evaluated.
args, params : positional and keyword arguments to the system Hamiltonian
"""
site_ranges = syst.site_ranges
if site_ranges is None:
# NOTE: Remove this branch once we can rely on
# site families storing the norb information.
site_arrays = syst.site_arrays
family_norbs = {s.family: None for s in site_arrays}
# Compute the norbs per site family using already evaluated
# Hamiltonian terms where possible
for ham, t in zip(hams, syst.terms):
(to_w, from_w), _ = syst.subgraphs[t.subgraph]
if ham is not None:
family_norbs[site_arrays[to_w].family] = ham.shape[1]
family_norbs[site_arrays[from_w].family] = ham.shape[2]
# Evaluate Hamiltonian terms where we do not already have them
for n, t in enumerate(syst.terms):
(to_w, from_w), _ = syst.subgraphs[t.subgraph]
to_fam = site_arrays[to_w].family
from_fam = site_arrays[from_w].family
if family_norbs[to_fam] is None or family_norbs[from_fam] is None:
ham = syst.hamiltonian_term(n, args=args, params=params)
family_norbs[to_fam] = ham.shape[1]
family_norbs[from_fam] = ham.shape[2]
# This case is technically allowed by the format (some sites present
# but no hamiltonian terms that touch them) but is very unlikely.
if any(norbs is None for norbs in family_norbs.values()):
raise ValueError("Cannot determine the number of orbitals for "
"some site families.")
orb_offset = 0
site_ranges = []
for start, site_array in zip(site_offsets, syst.site_arrays):
norbs = family_norbs[site_array.family]
site_ranges.append((start, norbs, orb_offset))
orb_offset += len(site_array) * norbs
site_ranges.append((site_offsets[-1], 0, orb_offset))
site_ranges = np.array(site_ranges)
_, norbs, orb_offsets = site_ranges.transpose()
# The final (extra) element in orb_offsets is the total number of orbitals
assert len(orb_offsets) == len(syst.site_arrays) + 1
return norbs, orb_offsets
def _expand_norbs(compressed_norbs, site_offsets):
"Return norbs per site, given norbs per site array."
norbs = np.empty(site_offsets[-1], int)
for start, stop, norb in zip(site_offsets, site_offsets[1:],
compressed_norbs):
norbs[start:stop] = norb
return norbs
def _reverse_subgraph(subgraph):
(to_sa, from_sa), (to_off, from_off) = subgraph
return ((from_sa, to_sa), (from_off, to_off))
@deprecate_args
@cython.binding(True)
def vectorized_hamiltonian_submatrix(self, args=(), sparse=False,
return_norb=False, *, params=None):
"""Return The system Hamiltonian.
Parameters
----------
args : tuple, defaults to empty
Positional arguments to pass to ``hamiltonian_term``. Mutually
exclusive with 'params'.
sparse : bool
Whether to return a sparse or a dense matrix. Defaults to ``False``.
return_norb : bool
Whether to return arrays of numbers of orbitals. Defaults to ``False``.
params : dict, optional
Dictionary of parameter names and their values. Mutually exclusive
with 'args'.
Returns
-------
hamiltonian_part : numpy.ndarray or scipy.sparse.coo_matrix
The Hamiltonian of the system.
norb : array of integers
Numbers of orbitals on each site. Only returned when ``return_norb``
is true.
Notes
-----
Providing positional arguments via 'args' is deprecated,
instead, provide named parameters as a dictionary via 'params'.
"""
site_offsets = np.cumsum([0] + [len(arr) for arr in self.site_arrays])
subgraphs = [
self.subgraphs[t.subgraph]
for t in self.terms
]
# Add Hermitian conjugates
subgraphs += [
_reverse_subgraph(self.subgraphs[t.subgraph])
for t in self.terms
if t.hermitian
]
hams = [
self.hamiltonian_term(n, args=args, params=params)
for n, _ in enumerate(self.terms)
]
# Add Hermitian conjugates
hams += [
ham.conjugate().transpose(0, 2, 1)
for ham, t in zip(hams, self.terms)
if t.hermitian
]
norbs, orb_offsets = _count_norbs(self, site_offsets, hams,
args=args, params=params)
func = _vectorized_make_sparse if sparse else _vectorized_make_dense
mat = func(subgraphs, hams, norbs, orb_offsets, site_offsets)
if return_norb:
return (mat, _expand_norbs(norbs, site_offsets))
else:
return mat
@deprecate_args
@cython.binding(True)
def vectorized_cell_hamiltonian(self, args=(), sparse=False, *, params=None):
"""Hamiltonian of a single cell of the infinite system.
Providing positional arguments via 'args' is deprecated,
instead, provide named parameters as a dictionary via 'params'.
"""
site_offsets = np.cumsum([0] + [len(arr) for arr in self.site_arrays])
# Site array where next cell starts
next_cell = bisect.bisect(site_offsets, self.cell_size) - 1
def inside_cell(term):
(to_which, from_which), _= self.subgraphs[term.subgraph]
return to_which < next_cell and from_which < next_cell
cell_terms = [
n for n, t in enumerate(self.terms)
if inside_cell(t)
]
subgraphs = [
self.subgraphs[self.terms[n].subgraph]
for n in cell_terms
]
# Add Hermitian conjugates
subgraphs += [
_reverse_subgraph(self.subgraphs[self.terms[n].subgraph])
for n in cell_terms
if self.terms[n].hermitian
]
hams = [
self.hamiltonian_term(n, args=args, params=params)
for n in cell_terms
]
# Add Hermitian conjugates
hams += [
ham.conjugate().transpose(0, 2, 1)
for ham, n in zip(hams, cell_terms)
if self.terms[n].hermitian
]
# _count_norbs requires passing hamiltonians for all terms, or 'None'
# if it has not been evaluated
all_hams = [None] * len(self.terms)
for n, ham in zip(cell_terms, hams):
all_hams[n] = ham
norbs, orb_offsets = _count_norbs(self, site_offsets, all_hams,
args=args, params=params)
shape = (orb_offsets[next_cell], orb_offsets[next_cell])
func = _vectorized_make_sparse if sparse else _vectorized_make_dense
mat = func(subgraphs, hams, norbs, orb_offsets, site_offsets, shape=shape)
return mat
@deprecate_args
@cython.binding(True)
def vectorized_inter_cell_hopping(self, args=(), sparse=False, *, params=None):
"""Hopping Hamiltonian between two cells of the infinite system.
Providing positional arguments via 'args' is deprecated,
instead, provide named parameters as a dictionary via 'params'.
"""
# Take advantage of the fact that there are distinct entries in
# onsite_terms for sites inside the unit cell, and distinct entries
# in onsite_terms for hoppings between the unit cell sites and
# interface sites. This way we only need to check the first entry
# in onsite/hopping_terms
site_offsets = np.cumsum([0] + [len(arr) for arr in self.site_arrays])
# Site array where next cell starts
next_cell = bisect.bisect(site_offsets, self.cell_size) - 1
inter_cell_hopping_terms = [
n for n, t in enumerate(self.terms)
# *from* site is in interface
if (self.subgraphs[t.subgraph][0][1] >= next_cell
and self.subgraphs[t.subgraph][0][0] < next_cell)
]
reversed_inter_cell_hopping_terms = [
n for n, t in enumerate(self.terms)
# *to* site is in interface
if (self.subgraphs[t.subgraph][0][0] >= next_cell
and self.subgraphs[t.subgraph][0][1] < next_cell)
]
# Evaluate inter-cell hoppings only
inter_cell_hams = [
self.hamiltonian_term(n, args=args, params=params)
for n in inter_cell_hopping_terms
]
reversed_inter_cell_hams = [
self.hamiltonian_term(n, args=args, params=params)
.conjugate().transpose(0, 2, 1)
for n in reversed_inter_cell_hopping_terms
]
hams = inter_cell_hams + reversed_inter_cell_hams
subgraphs = [
self.subgraphs[self.terms[n].subgraph]
for n in inter_cell_hopping_terms
]
subgraphs += [
_reverse_subgraph(self.subgraphs[self.terms[n].subgraph])
for n in reversed_inter_cell_hopping_terms
]
# All the 'from' sites are in the previous domain, but to build a
# matrix we need to get the equivalent sites in the fundamental domain.
# We set the 'from' site array to the one from the fundamental domain.
subgraphs = [
((to_sa, from_sa - next_cell), (to_off, from_off))
for (to_sa, from_sa), (to_off, from_off) in subgraphs
]
# _count_norbs requires passing hamiltonians for all terms, or 'None'
# if it has not been evaluated
all_hams = [None] * len(self.terms)
for n, ham in zip(inter_cell_hopping_terms, inter_cell_hams):
all_hams[n] = ham
for n, ham in zip(reversed_inter_cell_hopping_terms,
reversed_inter_cell_hams):
# Transpose to get back correct shape wrt. original term
all_hams[n] = ham.transpose(0, 2, 1)
norbs, orb_offsets = _count_norbs(self, site_offsets, all_hams,
args=args, params=params)
shape = (orb_offsets[next_cell],
orb_offsets[len(self.site_arrays) - next_cell])
func = _vectorized_make_sparse if sparse else _vectorized_make_dense
mat = func(subgraphs, hams, norbs, orb_offsets, site_offsets, shape=shape)
return mat
kwant/builder.py
View file @ 1d8db7bb
# Copyright 2011-2016 Kwant authors.
# Copyright 2011-2019 Kwant authors.
#
# This file is part of Kwant. It is subject to the license terms in the file
# LICENSE.rst found in the top-level directory of this distribution and at
......@@ -13,189 +13,29 @@ import collections
import copy
from functools import total_ordering, wraps, update_wrapper
from itertools import islice, chain
import textwrap
import bisect
import numbers
import inspect
import tinyarray as ta
import numpy as np
from scipy import sparse
from . import system, graph, KwantDeprecationWarning, UserCodeError
from .system import Site, SiteArray, SiteFamily
from .linalg import lll
from .operator import Density
from .physics import DiscreteSymmetry
from .physics import DiscreteSymmetry, magnetic_gauge
from ._common import (ensure_isinstance, get_parameters, reraise_warnings,
interleave, deprecate_args)
interleave, deprecate_args, memoize)
__all__ = ['Builder', 'Site', 'SiteFamily', 'SimpleSiteFamily', 'Symmetry',
'HoppingKind', 'Lead', 'BuilderLead', 'SelfEnergyLead', 'ModesLead']
__all__ = ['Builder', 'SimpleSiteFamily', 'Symmetry', 'HoppingKind', 'Lead',
'BuilderLead', 'SelfEnergyLead', 'ModesLead', 'add_peierls_phase']
################ Sites and site families
class Site(tuple):
"""A site, member of a `SiteFamily`.
Sites are the vertices of the graph which describes the tight binding
system in a `Builder`.
A site is uniquely identified by its family and its tag.
Parameters
----------
family : an instance of `SiteFamily`
The 'type' of the site.
tag : a hashable python object
The unique identifier of the site within the site family, typically a
vector of integers.
Raises
------
ValueError
If `tag` is not a proper tag for `family`.
Notes
-----
For convenience, ``family(*tag)`` can be used instead of ``Site(family,
tag)`` to create a site.
The parameters of the constructor (see above) are stored as instance
variables under the same names. Given a site ``site``, common things to
query are thus ``site.family``, ``site.tag``, and ``site.pos``.
"""
__slots__ = ()
family = property(operator.itemgetter(0),
doc="The site family to which the site belongs.")
tag = property(operator.itemgetter(1), doc="The tag of the site.")
def __new__(cls, family, tag, _i_know_what_i_do=False):
if _i_know_what_i_do:
return tuple.__new__(cls, (family, tag))
try:
tag = family.normalize_tag(tag)
except (TypeError, ValueError) as e:
msg = 'Tag {0} is not allowed for site family {1}: {2}'
raise type(e)(msg.format(repr(tag), repr(family), e.args[0]))
return tuple.__new__(cls, (family, tag))
def __repr__(self):
return 'Site({0}, {1})'.format(repr(self.family), repr(self.tag))
def __str__(self):
sf = self.family
return '<Site {0} of {1}>'.format(self.tag, sf.name if sf.name else sf)
def __getnewargs__(self):
return (self.family, self.tag, True)
@property
def pos(self):
"""Real space position of the site.
This relies on ``family`` having a ``pos`` method (see `SiteFamily`).
"""
return self.family.pos(self.tag)
################ Site families
@total_ordering
class SiteFamily(metaclass=abc.ABCMeta):
"""Abstract base class for site families.
Site families are the 'type' of `Site` objects. Within a family, individual
sites are uniquely identified by tags. Valid tags must be hashable Python
objects, further details are up to the family.
Site families must be immutable and fully defined by their initial
arguments. They must inherit from this abstract base class and call its
__init__ function providing it with two arguments: a canonical
representation and a name. The canonical representation will be returned as
the objects representation and must uniquely identify the site family
instance. The name is a string used to distinguish otherwise identical site
families. It may be empty. ``norbs`` defines the number of orbitals
on sites associated with this site family; it may be `None`, in which case
the number of orbitals is not specified.
All site families must define the method `normalize_tag` which brings a tag
to the standard format for this site family.
Site families that are intended for use with plotting should also provide a
method `pos(tag)`, which returns a vector with real-space coordinates of the
site belonging to this family with a given tag.
If the ``norbs`` of a site family are provided, and sites of this family
are used to populate a `~kwant.builder.Builder`, then the associated
Hamiltonian values must have the correct shape. That is, if a site family
has ``norbs = 2``, then any on-site terms for sites belonging to this
family should be 2x2 matrices. Similarly, any hoppings to/from sites
belonging to this family must have a matrix structure where there are two
rows/columns. This condition applies equally to Hamiltonian values that
are given by functions. If this condition is not satisfied, an error will
be raised.
"""
def __init__(self, canonical_repr, name, norbs):
self.canonical_repr = canonical_repr
self.hash = hash(canonical_repr)
self.name = name
if norbs is not None:
if int(norbs) != norbs or norbs <= 0:
raise ValueError('The norbs parameter must be an integer > 0.')
norbs = int(norbs)
self.norbs = norbs
def __repr__(self):
return self.canonical_repr
def __str__(self):
if self.name:
msg = '<{0} site family {1}{2}>'
else:
msg = '<unnamed {0} site family{2}>'
orbs = ' with {0} orbitals'.format(self.norbs) if self.norbs else ''
return msg.format(self.__class__.__name__, self.name, orbs)
def __hash__(self):
return self.hash
def __eq__(self, other):
try:
return self.canonical_repr == other.canonical_repr
except AttributeError:
return False
def __ne__(self, other):
try:
return self.canonical_repr != other.canonical_repr
except AttributeError:
return True
def __lt__(self, other):
# If this raises an AttributeError, we were trying
# to compare it to something non-comparable anyway.
return self.canonical_repr < other.canonical_repr
@abc.abstractmethod
def normalize_tag(self, tag):
"""Return a normalized version of the tag.
Raises TypeError or ValueError if the tag is not acceptable.
"""
pass
def __call__(self, *tag):
"""
A convenience function.
This function allows to write fam(1, 2) instead of Site(fam, (1, 2)).
"""
# Catch a likely and difficult to find mistake.
if tag and isinstance(tag[0], tuple):
raise ValueError('Use site_family(1, 2) instead of '
'site_family((1, 2))!')
return Site(self, tag)
class SimpleSiteFamily(SiteFamily):
"""A site family used as an example and for testing.
......@@ -301,19 +141,44 @@ class Symmetry(metaclass=abc.ABCMeta):
def which(self, site):
"""Calculate the domain of the site.
Return the group element whose action on a certain site from the
fundamental domain will result in the given ``site``.
Parameters
----------
site : `~kwant.system.Site` or `~kwant.system.SiteArray`
Returns
-------
group_element : tuple or sequence of tuples
A single tuple if ``site`` is a Site, or a sequence of tuples if
``site`` is a SiteArray. The group element(s) whose action
on a certain site(s) from the fundamental domain will result
in the given ``site``.
"""
pass
@abc.abstractmethod
def act(self, element, a, b=None):
"""Act with a symmetry group element on a site or hopping."""
"""Act with symmetry group element(s) on site(s) or hopping(s).
Parameters
----------
element : tuple or sequence of tuples
Group element(s) with which to act on the provided site(s)
or hopping(s)
a, b : `~kwant.system.Site` or `~kwant.system.SiteArray`
If Site then ``element`` is a single tuple, if SiteArray then
``element`` is a sequence of tuples. If only ``a`` is provided then
``element`` acts on the site(s) of ``a``. If ``b`` is also provided
then ``element`` acts on the hopping(s) ``(a, b)``.
"""
pass
def to_fd(self, a, b=None):
"""Map a site or hopping to the fundamental domain.
Parameters
----------
a, b : `~kwant.system.Site` or `~kwant.system.SiteArray`
If ``b`` is None, return a site equivalent to ``a`` within the
fundamental domain. Otherwise, return a hopping equivalent to ``(a,
b)`` but where the first element belongs to the fundamental domain.
......@@ -323,11 +188,30 @@ class Symmetry(metaclass=abc.ABCMeta):
return self.act(-self.which(a), a, b)
def in_fd(self, site):
"""Tell whether ``site`` lies within the fundamental domain."""
for d in self.which(site):
if d != 0:
return False
return True
"""Tell whether ``site`` lies within the fundamental domain.
Parameters
----------
site : `~kwant.system.Site` or `~kwant.system.SiteArray`
Returns
-------
in_fd : bool or sequence of bool
single bool if ``site`` is a Site, or a sequence of
bool if ``site`` is a SiteArray. In the latter case
we return whether each site in the SiteArray is in
the fundamental domain.
"""
if isinstance(site, Site):
for d in self.which(site):
if d != 0:
return False
return True
elif isinstance(site, SiteArray):
which = self.which(site)
return np.logical_and.reduce(which != 0, axis=1)
else:
raise TypeError("'site' must be a Site or SiteArray")
@abc.abstractmethod
def subgroup(self, *generators):
......@@ -406,8 +290,8 @@ class HoppingKind(tuple):
----------
delta : Sequence of integers
The sequence is interpreted as a vector with integer elements.
family_a : `~kwant.builder.SiteFamily`
family_b : `~kwant.builder.SiteFamily` or ``None`` (default)
family_a : `~kwant.system.SiteFamily`
family_b : `~kwant.system.SiteFamily` or ``None`` (default)
The default value means: use the same family as `family_a`.
Notes
......@@ -564,7 +448,7 @@ class BuilderLead(Lead):
The tight-binding system of a lead. It has to possess appropriate
symmetry, and it may not contain hoppings between further than
neighboring images of the fundamental domain.
interface : sequence of `Site` instances
interface : sequence of `~kwant.system.Site` instances
Sequence of sites in the scattering region to which the lead is
attached.
......@@ -572,9 +456,9 @@ class BuilderLead(Lead):
----------
builder : `Builder`
The tight-binding system of a lead.
interface : list of `Site` instances
interface : list of `~kwant.system.Site` instances
A sorted list of interface sites.
padding : list of `Site` instances
padding : list of `~kwant.system.Site` instances
A sorted list of sites that originate from the lead, have the same
onsite Hamiltonian, and are connected by the same hoppings as the lead
sites.
......@@ -601,7 +485,10 @@ class BuilderLead(Lead):
The order of interface sites is kept during finalization.
"""
return InfiniteSystem(self.builder, self.interface)
if self.builder.vectorize:
return InfiniteVectorizedSystem(self.builder, self.interface)
else:
return InfiniteSystem(self.builder, self.interface)
def _ensure_signature(func):
......@@ -632,7 +519,7 @@ class SelfEnergyLead(Lead):
selfenergy_func : function
Has the same signature as `selfenergy` (without the ``self``
parameter) and returns the self energy matrix for the interface sites.
interface : sequence of `Site` instances
interface : sequence of `~kwant.system.Site` instances
parameters : sequence of strings
The parameters on which the lead depends.
"""
......@@ -663,7 +550,7 @@ class ModesLead(Lead):
and returns the modes of the lead as a tuple of
`~kwant.physics.PropagatingModes` and `~kwant.physics.StabilizedModes`.
interface :
sequence of `Site` instances
sequence of `~kwant.system.Site` instances
parameters : sequence of strings
The parameters on which the lead depends.
"""
......@@ -791,10 +678,11 @@ class Builder:
This is one of the central types in Kwant. It is used to construct tight
binding systems in a flexible way.
The nodes of the graph are `Site` instances. The edges, i.e. the hoppings,
are pairs (2-tuples) of sites. Each node and each edge has a value
associated with it. The values associated with nodes are interpreted as
on-site Hamiltonians, the ones associated with edges as hopping integrals.
The nodes of the graph are `~kwant.system.Site` instances. The edges,
i.e. the hoppings, are pairs (2-tuples) of sites. Each node and each
edge has a value associated with it. The values associated with nodes
are interpreted as on-site Hamiltonians, the ones associated with edges
as hopping integrals.
To make the graph accessible in a way that is natural within the Python
language it is exposed as a *mapping* (much like a built-in Python
......@@ -821,6 +709,14 @@ class Builder:
chiral : 2D array, dictionary, function or `None`
The unitary part of the onsite chiral symmetry operator.
Same format as that of `conservation_law`.
vectorize : bool, default: False
If True then the finalized Builder will evaluate its Hamiltonian in
a vectorized manner. This requires that all value functions take
`~kwant.system.SiteArray` as first/second parameter rather than
`~kwant.system.Site`, and return an array of values. The returned
array must have the same length as the provided SiteArray(s), and
may contain either scalars (i.e. a vector of values) or matrices
(i.e. a 3d array of values).
Notes
-----
......@@ -899,7 +795,7 @@ class Builder:
"""
def __init__(self, symmetry=None, *, conservation_law=None, time_reversal=None,
particle_hole=None, chiral=None):
particle_hole=None, chiral=None, vectorize=False):
if symmetry is None:
symmetry = NoSymmetry()
else:
......@@ -909,6 +805,7 @@ class Builder:
self.time_reversal = time_reversal
self.particle_hole = particle_hole
self.chiral = chiral
self.vectorize = vectorize
self.leads = []
self.H = {}
......@@ -992,6 +889,7 @@ class Builder:
result.time_reversal = self.time_reversal
result.particle_hole = self.particle_hole
result.chiral = self.chiral
result.vectorize = self.vectorize
result.leads = self.leads
result.H = self.H
return result
......@@ -1446,7 +1344,7 @@ class Builder:
A boolean function of site returning whether the site should be
added to the target builder or not. The shape must be compatible
with the symmetry of the target builder.
start : `Site` instance or iterable thereof or iterable of numbers
start : `~kwant.system.Site` or iterable thereof or iterable of numbers
The site(s) at which the the flood-fill starts. If start is an
iterable of numbers, the starting site will be
``template.closest(start)``.
......@@ -1456,7 +1354,7 @@ class Builder:
Returns
-------
added_sites : list of `Site` objects that were added to the system.
added_sites : list of `~kwant.system.Site` that were added to the system.
Raises
------
......@@ -1608,7 +1506,7 @@ class Builder:
----------
lead_builder : `Builder` with 1D translational symmetry
Builder of the lead which has to be attached.
origin : `Site`
origin : `~kwant.system.Site`
The site which should belong to a domain where the lead should
begin. It is used to attach a lead inside the system, e.g. to an
inner radius of a ring.
......@@ -1618,7 +1516,7 @@ class Builder:
Returns
-------
added_sites : list of `Site` objects that were added to the system.
added_sites : list of `~kwant.system.Site`
Raises
------
......@@ -1671,6 +1569,13 @@ class Builder:
if add_cells < 0 or int(add_cells) != add_cells:
raise ValueError('add_cells must be an integer >= 0.')
if self.vectorize != lead_builder.vectorize:
raise ValueError(
"Sites of the lead were added to the scattering "
"region, but only one of these systems is vectorized. "
"Vectorize both systems to allow attaching leads."
)
sym = lead_builder.symmetry
H = lead_builder.H
......@@ -1691,6 +1596,7 @@ class Builder:
time_reversal=lead_builder.time_reversal,
particle_hole=lead_builder.particle_hole,
chiral=lead_builder.chiral,
vectorize=lead_builder.vectorize,
)
with reraise_warnings():
new_lead.fill(lead_builder, lambda site: True,
......@@ -1787,9 +1693,15 @@ class Builder:
`Symmetry` can be finalized.
"""
if self.symmetry.num_directions == 0:
return FiniteSystem(self)
if self.vectorize:
return FiniteVectorizedSystem(self)
else:
return FiniteSystem(self)
elif self.symmetry.num_directions == 1:
return InfiniteSystem(self)
if self.vectorize:
return InfiniteVectorizedSystem(self)
else:
return InfiniteSystem(self)
else:
raise ValueError('Currently, only builders without or with a 1D '
'translational symmetry can be finalized.')
......@@ -1807,6 +1719,173 @@ class Builder:
inter_cell_hopping = cell_hamiltonian = precalculated = \
_require_system
def _add_peierls_phase(syst, peierls_parameter):
@memoize
def wrap_hopping(hop):
def f(*args):
a, b, *args, phi = args
return hop(a, b, *args) * phi(a, b)
params = ('_site0', '_site1')
if callable(hop):
params += get_parameters(hop)[2:] # cut off both site parameters
params += (peierls_parameter,)
f.__signature__ = inspect.Signature(
inspect.Parameter(name, inspect.Parameter.POSITIONAL_ONLY)
for name in params)
return f
const_hoppings = {}
def const_hopping(a, b, phi):
return const_hoppings[a, b] * phi(a, b)
# fix the signature
params = ('_site0', '_site1', peierls_parameter)
const_hopping.__signature__ = inspect.Signature(
inspect.Parameter(name, inspect.Parameter.POSITIONAL_ONLY)
for name in params)
ret = Builder(
syst.symmetry,
# No time reversal, as magnetic fields break it
conservation_law=syst.conservation_law,
particle_hole=syst.particle_hole,
chiral=syst.chiral,
)
for i, lead in enumerate(syst.leads):
ret.leads.append(
BuilderLead(
_add_peierls_phase(
lead.builder,
peierls_parameter + '_lead{}'.format(i),
),
lead.interface,
lead.padding,
)
)
for site, val in syst.site_value_pairs():
ret[site] = val
for hop, val in syst.hopping_value_pairs():
if callable(val):
ret[hop] = wrap_hopping(val)
else:
const_hoppings[hop] = val
ret[hop] = const_hopping
return ret
def add_peierls_phase(syst, peierls_parameter='phi', fix_gauge=True):
"""Add a Peierls phase parameter to a Builder.
Parameters
----------
syst : kwant.Builder
The system to which to add a Peierls phase.
peierls_parameter : str
The name of the Peierls phase parameter to value functions.
fix_gauge : bool (default: True)
If True, fix the gauge of the system and return it also
Returns
-------
syst : kwant.builder.Finitesystem or kwant.builder.InfiniteSystem
A system where all hoppings are given by value functions
that have an additional parameter 'phi' that should be
the Peierls phase returned by `kwant.physics.magnetic_gauge`.
Any leads have similar parameters named 'phi_lead0',
'phi_lead1' etc.
gauge : callable
Only returned if 'fix_gauge' is True. Called with magnetic
field(s), and returns a parameter dictionary that can be
passed to the system as 'params' (see the example below).
Examples
--------
>>> import numpy as np
>>> import kwant
>>>
>>> syst = make_system()
>>> lead = make_lead()
>>> syst.attach_lead(lead)
>>> syst.attach_lead(lead.reversed())
>>>
>>> fsyst, phase = kwant.builder.add_peierls_phase(syst)
>>>
>>> def B_syst(pos):
>>> return np.exp(-np.sum(pos * pos))
>>>
>>> kwant.smatrix(fsyst, parameters=dict(t=-1, **phase(B_syst, 0, 0)))
"""
def wrap_gauge(gauge):
phase_names = [peierls_parameter]
phase_names += [peierls_parameter + '_lead{}'.format(i)
for i in range(len(syst.leads))]
doc = textwrap.dedent("""\
Return the Peierls phase for a magnetic field configuration.
Parameters
----------
syst_field : scalar, vector or callable
The magnetic field to apply to the scattering region.
If callable, takes a position and returns the
magnetic field at that position. Can be a scalar if
the system is 1D or 2D, otherwise must be a vector.
Magnetic field is expressed in units :math:`φ₀ / l²`,
where :math:`φ₀` is the magnetic flux quantum and
:math:`l` is the unit of length.
*lead_fields : scalar, vector or callable
The magnetic fields to apply to each of the leads, in
the same format as 'syst_field'. In addition, if a callable
is provided, then the magnetic field must have the symmetry
of the associated lead.
tol : float, default: 1E-8
The tolerance to which to calculate the flux through each
hopping loop in the system.
average : bool, default: False
If True, estimate the magnetic flux through each hopping loop
in the system by evaluating the magnetic field at a single
position inside the loop and multiplying it by the area of the
loop. If False, then ``scipy.integrate.quad`` is used to integrate
the magnetic field. This parameter is only used when 'syst_field'
or 'lead_fields' are callable.
Returns
-------
phases : dict of callables
To be passed directly as the parameters of a system wrapped with
'kwant.builder.add_peierls_phase'.
""")
@wraps(gauge)
def f(*args, **kwargs):
phases = gauge(*args, **kwargs)
return dict(zip(phase_names, phases))
f.__doc__ = doc
return f
if syst.vectorize:
raise TypeError("'add_peierls_phase' does not work with "
"vectorized Builders")
ret = _add_peierls_phase(syst, peierls_parameter).finalized()
if fix_gauge:
return ret, wrap_gauge(magnetic_gauge(ret))
else:
return ret
################ Finalized systems
......@@ -1946,6 +2025,135 @@ class _FinalizedBuilderMixin:
return DiscreteSymmetry(projectors, *(evaluate(symm) for symm in
self._symmetries))
def pos(self, i):
return self.sites[i].pos
class _VectorizedFinalizedBuilderMixin(_FinalizedBuilderMixin):
"""Common functionality for all vectorized finalized builders
Attributes
----------
_term_values : sequence
Each item is either an array of values (if 'param_names is None')
or a vectorized value function (in which case 'param_names' is a
sequence of strings: the extra parameters to the value function).
_onsite_term_by_site_id : sequence of int
Maps site (integer) to the term that evaluates the associated
onsite matrix element.
_hopping_term_by_edge_id : sequence of int
Maps edge id in the graph to the term that evaluates the associated
hopping matrix element. If negative, then the associated hopping is
the Hermitian conjugate of another hopping; if the number stored is
't' (< 0) then the associated hopping is stored in term '-t - 1'.
"""
def hamiltonian(self, i, j, *args, params=None):
warnings.warn(
"Querying individual matrix elements with 'hamiltonian' is "
"deprecated, and now takes O(log N) time where N is the number "
"of matrix elements per hamiltonian term. Query many matrix "
"elements at once using 'hamiltonian_term'.",
KwantDeprecationWarning
)
site_offsets = np.cumsum([0] + [len(s) for s in self.site_arrays])
if i == j:
which_term = self._onsite_term_by_site_id[i]
(w, _), (off, _) = self.subgraphs[self.terms[which_term].subgraph]
i_off = i - site_offsets[w]
selector = np.searchsorted(off, i_off)
onsite = self.hamiltonian_term(
which_term, [selector], args, params=params)
return onsite[0]
else:
edge_id = self.graph.first_edge_id(i, j)
which_term = self._hopping_term_by_edge_id[edge_id]
herm_conj = which_term < 0
if herm_conj:
which_term = -which_term - 1
term = self.terms[which_term]
# To search for hopping (i, j) in hopping_terms, we need
# to treat hopping_terms as a record array of integer pairs
dtype = np.dtype([('f0', int), ('f1', int)])
# For hermitian conjugate terms search through the
# *other* term's hoppings because those are sorted.
(to_w, from_w), hoppings = self.subgraphs[term.subgraph]
hoppings = hoppings.transpose() # to get array-of-pairs
if herm_conj:
hop = (j - site_offsets[to_w], i - site_offsets[from_w])
else:
hop = (i - site_offsets[to_w], j - site_offsets[from_w])
hop = np.array(hop, dtype=dtype)
hoppings = hoppings.view(dtype).reshape(-1)
selector = np.recarray.searchsorted(hoppings, hop)
h = self.hamiltonian_term(
which_term, [selector], args, params=params)
h = h[0]
if herm_conj:
h = h.conjugate().transpose()
return h
def hamiltonian_term(self, term_number, selector=slice(None),
args=(), params=None):
if args and params:
raise TypeError("'args' and 'params' are mutually exclusive.")
term = self.terms[term_number]
val = self._term_values[term_number]
if not callable(val):
return val[selector]
# Construct site arrays to pass to the vectorized value function.
subgraph = self.subgraphs[self.terms[term_number].subgraph]
(to_which, from_which), (to_off, from_off) = subgraph
is_onsite = to_off is from_off
to_off = to_off[selector]
from_off = from_off[selector]
assert len(to_off) == len(from_off)
to_family = self.site_arrays[to_which].family
to_tags = self.site_arrays[to_which].tags
to_site_array = SiteArray(to_family, to_tags[to_off])
if not is_onsite:
from_family = self.site_arrays[from_which].family
from_tags = self.site_arrays[from_which].tags
from_site_array = SiteArray(from_family, from_tags[from_off])
# Construct args from params
if params:
# There was a problem extracting parameter names from the value
# function (probably an illegal signature) and we are using
# keyword parameters.
if self._term_errors[term_number] is not None:
raise self._term_errors[term_number]
try:
args = [params[p] for p in term.parameters]
except KeyError:
missing = [p for p in term.parameters if p not in params]
msg = ('System is missing required arguments: ',
', '.join(map('"{}"'.format, missing)))
raise TypeError(''.join(msg))
if is_onsite:
try:
ham = val(to_site_array, *args)
except Exception as exc:
_raise_user_error(exc, val)
else:
try:
ham = val(
*self.symmetry.to_fd(
to_site_array,
from_site_array),
*args)
except Exception as exc:
_raise_user_error(exc, val)
ham = _normalize_term_value(ham, len(to_site_array))
return ham
# The same (value, parameters) pair will be used for many sites/hoppings,
# so we cache it to avoid wasting extra memory.
......@@ -1978,6 +2186,64 @@ def _value_params_pair_cache(nstrip):
return get
def _make_graph(H, id_by_site):
g = graph.Graph()
g.num_nodes = len(id_by_site) # Some sites could not appear in any edge.
for tail, hvhv in H.items():
for head in islice(hvhv, 2, None, 2):
if tail == head:
continue
g.add_edge(id_by_site[tail], id_by_site[head])
return g.compressed()
def _finalize_leads(leads, id_by_site):
#### Connect leads.
finalized_leads = []
lead_interfaces = []
lead_paddings = []
for lead_nr, lead in enumerate(leads):
try:
with warnings.catch_warnings(record=True) as ws:
warnings.simplefilter("always")
# The following line is the whole "payload" of the entire
# try-block.
finalized_leads.append(lead.finalized())
except ValueError as e:
# Re-raise the exception with an additional message.
msg = 'Problem finalizing lead {0}:'.format(lead_nr)
e.args = (' '.join((msg,) + e.args),)
raise
else:
for w in ws:
# Re-raise any warnings with an additional message and the
# proper stacklevel.
w = w.message
msg = 'When finalizing lead {0}:'.format(lead_nr)
warnings.warn(w.__class__(' '.join((msg,) + w.args)),
stacklevel=3)
try:
interface = [id_by_site[isite] for isite in lead.interface]
except KeyError as e:
msg = ("Lead {0} is attached to a site that does not "
"belong to the scattering region:\n {1}")
raise ValueError(msg.format(lead_nr, e.args[0]))
lead_interfaces.append(np.array(interface))
padding = getattr(lead, 'padding', [])
# Some padding sites might have been removed after the lead was
# attached. Unlike in the case of the interface, this is not a
# problem.
finalized_padding = [
id_by_site[isite] for isite in padding if isite in id_by_site
]
lead_paddings.append(np.array(finalized_padding))
return finalized_leads, lead_interfaces, lead_paddings
class FiniteSystem(_FinalizedBuilderMixin, system.FiniteSystem):
"""Finalized `Builder` with leads.
......@@ -1986,11 +2252,11 @@ class FiniteSystem(_FinalizedBuilderMixin, system.FiniteSystem):
Attributes
----------
sites : sequence
``sites[i]`` is the `~kwant.builder.Site` instance that corresponds
``sites[i]`` is the `~kwant.system.Site` instance that corresponds
to the integer-labeled site ``i`` of the low-level system. The sites
are ordered first by their family and then by their tag.
id_by_site : dict
The inverse of ``sites``; maps high-level `~kwant.builder.Site`
id_by_site : mapping
The inverse of ``sites``; maps high-level `~kwant.system.Site`
instances to their integer label.
Satisfies ``id_by_site[sites[i]] == i``.
"""
......@@ -2004,57 +2270,10 @@ class FiniteSystem(_FinalizedBuilderMixin, system.FiniteSystem):
for site_id, site in enumerate(sites):
id_by_site[site] = site_id
#### Make graph.
g = graph.Graph()
g.num_nodes = len(sites) # Some sites could not appear in any edge.
for tail, hvhv in builder.H.items():
for head in islice(hvhv, 2, None, 2):
if tail == head:
continue
g.add_edge(id_by_site[tail], id_by_site[head])
g = g.compressed()
#### Connect leads.
finalized_leads = []
lead_interfaces = []
lead_paddings = []
for lead_nr, lead in enumerate(builder.leads):
try:
with warnings.catch_warnings(record=True) as ws:
warnings.simplefilter("always")
# The following line is the whole "payload" of the entire
# try-block.
finalized_leads.append(lead.finalized())
for w in ws:
# Re-raise any warnings with an additional message and the
# proper stacklevel.
w = w.message
msg = 'When finalizing lead {0}:'.format(lead_nr)
warnings.warn(w.__class__(' '.join((msg,) + w.args)),
stacklevel=3)
except ValueError as e:
# Re-raise the exception with an additional message.
msg = 'Problem finalizing lead {0}:'.format(lead_nr)
e.args = (' '.join((msg,) + e.args),)
raise
try:
interface = [id_by_site[isite] for isite in lead.interface]
except KeyError as e:
msg = ("Lead {0} is attached to a site that does not "
"belong to the scattering region:\n {1}")
raise ValueError(msg.format(lead_nr, e.args[0]))
graph = _make_graph(builder.H, id_by_site)
lead_interfaces.append(np.array(interface))
padding = getattr(lead, 'padding', [])
# Some padding sites might have been removed after the lead was
# attached. Unlike in the case of the interface, this is not a
# problem.
finalized_padding = [
id_by_site[isite] for isite in padding if isite in id_by_site
]
lead_paddings.append(np.array(finalized_padding))
finalized_leads, lead_interfaces, lead_paddings =\
_finalize_leads(builder.leads, id_by_site)
# Because many onsites/hoppings share the same (value, parameter)
# pairs, we keep them in a cache so that we only store a given pair
......@@ -2063,7 +2282,7 @@ class FiniteSystem(_FinalizedBuilderMixin, system.FiniteSystem):
onsites = [cache(builder.H[site][1]) for site in sites]
cache = _value_params_pair_cache(2)
hoppings = [cache(builder._get_edge(sites[tail], sites[head]))
for tail, head in g]
for tail, head in graph]
# Compute the union of the parameters of onsites and hoppings. Here,
# 'onsites' and 'hoppings' are pairs whose second element is one of
......@@ -2083,7 +2302,7 @@ class FiniteSystem(_FinalizedBuilderMixin, system.FiniteSystem):
else:
parameters = frozenset(parameters)
self.graph = g
self.graph = graph
self.sites = sites
self.site_ranges = _site_ranges(sites)
self.id_by_site = id_by_site
......@@ -2096,8 +2315,502 @@ class FiniteSystem(_FinalizedBuilderMixin, system.FiniteSystem):
self.lead_paddings = lead_paddings
self._init_discrete_symmetries(builder)
def pos(self, i):
return self.sites[i].pos
def _make_site_arrays(sites):
tags_by_family = {}
for family, tag in sites: # Sites are tuples
tags_by_family.setdefault(family, []).append(tag)
site_arrays = []
for family, tags in sorted(tags_by_family.items()):
site_arrays.append(SiteArray(family, sorted(tags)))
return site_arrays
# Wrapper for accessing sites by their sequential number from a
# list of SiteArrays.
class _Sites:
def __init__(self, site_arrays):
self._site_arrays = site_arrays
lengths = [0] + [len(arr.tags) for arr in site_arrays]
self._start_idxs = np.cumsum(lengths)[:-1]
def __len__(self):
return sum(len(arr.tags) for arr in self._site_arrays)
def __getitem__(self, idx):
if not isinstance(idx, numbers.Integral):
raise TypeError("Only individual sites may be indexed")
if idx < 0 or idx >= len(self):
raise IndexError(idx)
which = bisect.bisect(self._start_idxs, idx) - 1
arr = self._site_arrays[which]
start = self._start_idxs[which]
fam = arr.family
tag = arr.tags[idx - start]
return Site(fam, tag)
def __iter__(self):
for arr in self._site_arrays:
for tag in arr.tags:
yield Site(arr.family, tag)
# Wrapper for accessing the sequential number of a site, given
# Site object, from a list of SiteArrays.
class _IdBySite:
def __init__(self, site_arrays):
self._site_arrays = site_arrays
lengths = [0] + [len(arr.tags) for arr in site_arrays]
self._start_idxs = np.cumsum(lengths)[:-1]
def __len__(self):
return sum(len(arr.tags) for arr in self._site_arrays)
def __getitem__(self, site):
# treat tags as 1D sequence by defining custom dtype
tag_dtype = np.asarray(site.tag).dtype
dtype = np.dtype([('f{}'.format(i), tag_dtype)
for i in range(len(site.tag))])
site_tag = np.asarray(site.tag).view(dtype)[0]
# This slightly convoluted logic is necessary because there
# may be >1 site_array with the same family for InfiniteSystems.
for start, arr in zip(self._start_idxs, self._site_arrays):
if arr.family != site.family:
continue
tags = arr.tags.view(dtype).reshape(-1)
idx_in_fam = np.recarray.searchsorted(tags, site_tag)
if idx_in_fam >= len(tags) or tags[idx_in_fam] != site_tag:
continue
return start + idx_in_fam
raise KeyError(site)
def _normalize_term_value(value, expected_n_values):
try:
value = np.asarray(value, dtype=complex)
except TypeError:
raise ValueError(
"Matrix elements declared with incompatible shapes."
) from None
# Upgrade to vector of matrices
if len(value.shape) == 1:
value = value[:, np.newaxis, np.newaxis]
if len(value.shape) != 3:
msg = (
"Vectorized value functions must return an array of"
"scalars or an array of matrices."
)
raise ValueError(msg)
if value.shape[0] != expected_n_values:
raise ValueError("Value functions must return a single value per "
"onsite/hopping.")
return value
def _sort_term(term, value):
term = np.asarray(term)
if not callable(value):
value = _normalize_term_value(value, len(term))
# Ensure that values still correspond to the correct
# sites in 'term' once the latter has been sorted.
value = value[term.argsort()]
term.sort()
return term, value
def _sort_hopping_term(term, value):
term = np.asarray(term)
# We want to sort the hopping terms in the same way
# as tuples (i, j), so we need to use a record array.
dtype = np.dtype([('f0', term.dtype), ('f1', term.dtype)])
term_prime = term.view(dtype=dtype).reshape(-1)
# _normalize_term will sort 'term' in-place via 'term_prime'
_, value = _sort_term(term_prime, value)
return term, value
def _make_onsite_terms(builder, sites, site_offsets, term_offset):
# Construct onsite terms.
#
# onsite_subgraphs
# Will contain tuples of the form described in
# kwant.system.VectorizedSystem.subgraphs, but during construction
# contains lists of the sites associated with each onsite term.
#
# onsite_term_values
# Contains a callable or array of constant values for each term.
#
# onsite_terms
# Constructed after the main loop. Contains a kwant.system.Term
# tuple for each onsite term.
#
# _onsite_term_by_site_id
# Maps the site ID to the number of the term that the site is
# a part of.
# lists the number of the
# Hamiltonian term associated with each site/hopping. For
# Hermitian conjugate hoppings "-term - 1" is stored instead.
onsite_subgraphs = []
onsite_term_values = []
onsite_term_parameters = []
# We just use the cache to handle non/callables and exceptions
cache = _value_params_pair_cache(1)
# maps onsite key -> onsite term number
onsite_to_term_nr = collections.OrderedDict()
_onsite_term_by_site_id = []
for site_id, site in enumerate(sites):
val = builder.H[builder.symmetry.to_fd(site)][1]
const_val = not callable(val)
which_site_array = bisect.bisect(site_offsets, site_id) - 1
# "key" uniquely identifies an onsite term.
if const_val:
key = (None, which_site_array)
else:
key = (id(val), which_site_array)
if key not in onsite_to_term_nr:
# Start a new term
onsite_to_term_nr[key] = len(onsite_subgraphs)
onsite_subgraphs.append([])
if const_val:
onsite_term_values.append([])
onsite_term_parameters.append(None)
else:
val, params = cache(val)
onsite_term_values.append(val)
onsite_term_parameters.append(params)
onsite_subgraphs[onsite_to_term_nr[key]].append(site_id)
_onsite_term_by_site_id.append(onsite_to_term_nr[key])
if const_val:
vals = onsite_term_values[onsite_to_term_nr[key]]
vals.append(val)
# Sort the sites in each term, and normalize any constant
# values to arrays of the correct dtype and shape.
onsite_subgraphs, onsite_term_values = zip(*(
_sort_term(term, value)
for term, value in
zip(onsite_subgraphs, onsite_term_values)))
# Store site array index and site offsets rather than sites themselves
tmp = []
for (_, which), s in zip(onsite_to_term_nr, onsite_subgraphs):
s = s - site_offsets[which]
tmp.append(((which, which), (s, s)))
onsite_subgraphs = tmp
# onsite_term_errors[i] contains an exception if the corresponding
# term has a value function with an illegal signature. We only raise
# the exception if we actually try to evaluate the offending term
# (to maintain backwards compatibility).
onsite_term_errors = [
err if isinstance(err, Exception) else None
for err in onsite_term_parameters
]
onsite_terms = [
system.Term(
subgraph=term_offset + i,
hermitian=False,
parameters=(
params if not isinstance(params, Exception) else None
),
)
for i, params in enumerate(onsite_term_parameters)
]
_onsite_term_by_site_id = np.array(_onsite_term_by_site_id)
return (onsite_subgraphs, onsite_terms, onsite_term_values,
onsite_term_parameters, onsite_term_errors, _onsite_term_by_site_id)
def _make_hopping_terms(builder, graph, sites, site_offsets, cell_size, term_offset):
# Construct the hopping terms
#
# The logic is the same as for the onsite terms, with the following
# differences.
#
# _hopping_term_by_edge_id
# Maps hopping edge IDs to the number of the term that the hopping
# is a part of. For Hermitian conjugate hoppings "-term_number -1"
# is stored instead.
hopping_subgraphs = []
hopping_term_values = []
hopping_term_parameters = []
# We just use the cache to handle non/callables and exceptions
cache = _value_params_pair_cache(2)
# maps hopping key -> hopping term number
hopping_to_term_nr = collections.OrderedDict()
_hopping_term_by_edge_id = []
for tail, head in graph:
tail_site, head_site = sites[tail], sites[head]
if tail >= cell_size:
# The tail belongs to the previous domain. Find the
# corresponding hopping with the tail in the fund. domain.
tail_site, head_site = builder.symmetry.to_fd(tail_site, head_site)
val = builder._get_edge(tail_site, head_site)
herm_conj = val is Other
if herm_conj:
val = builder._get_edge(*builder.symmetry.to_fd(head_site, tail_site))
const_val = not callable(val)
tail_site_array = bisect.bisect(site_offsets, tail) - 1
head_site_array = bisect.bisect(site_offsets, head) - 1
# "key" uniquely identifies a hopping term.
if const_val:
key = (None, tail_site_array, head_site_array)
else:
key = (id(val), tail_site_array, head_site_array)
if herm_conj:
key = (key[0], head_site_array, tail_site_array)
if key not in hopping_to_term_nr:
# start a new term
hopping_to_term_nr[key] = len(hopping_subgraphs)
hopping_subgraphs.append([])
if const_val:
hopping_term_values.append([])
hopping_term_parameters.append(None)
else:
val, params = cache(val)
hopping_term_values.append(val)
hopping_term_parameters.append(params)
if herm_conj:
# Hopping terms come after all onsite terms, so we need
# to offset the term ID by that many
term_id = -term_offset - hopping_to_term_nr[key] - 1
_hopping_term_by_edge_id.append(term_id)
else:
# Hopping terms come after all onsite terms, so we need
# to offset the term ID by that many
term_id = term_offset + hopping_to_term_nr[key]
_hopping_term_by_edge_id.append(term_id)
hopping_subgraphs[hopping_to_term_nr[key]].append((tail, head))
if const_val:
vals = hopping_term_values[hopping_to_term_nr[key]]
vals.append(val)
# Sort the hoppings in each term, and normalize any constant
# values to arrays of the correct dtype and shape.
if hopping_subgraphs:
hopping_subgraphs, hopping_term_values = zip(*(
_sort_hopping_term(term, value)
for term, value in
zip(hopping_subgraphs, hopping_term_values)))
# Store site array index and site offsets rather than sites themselves
tmp = []
for (_, tail_which, head_which), h in zip(hopping_to_term_nr,
hopping_subgraphs):
start = (site_offsets[tail_which], site_offsets[head_which])
# Transpose to get a pair of arrays rather than array of pairs
# We use the fact that the underlying array is stored in
# array-of-pairs order to search through it in 'hamiltonian'.
tmp.append(((tail_which, head_which), (h - start).transpose()))
hopping_subgraphs = tmp
# hopping_term_errors[i] contains an exception if the corresponding
# term has a value function with an illegal signature. We only raise
# the exception if this becomes a problem (to maintain backwards
# compatibility)
hopping_term_errors = [
err if isinstance(err, Exception) else None
for err in hopping_term_parameters
]
hopping_terms = [
system.Term(
subgraph=term_offset + i,
hermitian=True, # Builders are always Hermitian
parameters=(
params if not isinstance(params, Exception) else None
),
)
for i, params in enumerate(hopping_term_parameters)
]
_hopping_term_by_edge_id = np.array(_hopping_term_by_edge_id)
return (hopping_subgraphs, hopping_terms, hopping_term_values,
hopping_term_parameters, hopping_term_errors,
_hopping_term_by_edge_id)
class FiniteVectorizedSystem(_VectorizedFinalizedBuilderMixin, system.FiniteVectorizedSystem):
"""Finalized `Builder` with leads.
Usable as input for the solvers in `kwant.solvers`.
Attributes
----------
site_arrays : sequence of `~kwant.system.SiteArray`
The sites of the system.
sites : sequence
``sites[i]`` is the `~kwant.system.Site` instance that corresponds
to the integer-labeled site ``i`` of the low-level system. The sites
are ordered first by their family and then by their tag.
id_by_site : mapping
The inverse of ``sites``; maps high-level `~kwant.system.Site`
instances to their integer label.
Satisfies ``id_by_site[sites[i]] == i``.
"""
def __init__(self, builder):
assert builder.symmetry.num_directions == 0
assert builder.vectorize
sites = sorted(builder.H)
id_by_site = {site: site_id for site_id, site in enumerate(sites)}
if not all(s.family.norbs for s in sites):
raise ValueError("All sites must belong to families with norbs "
"specified for vectorized Builders. Specify "
"norbs when creating site families.")
graph = _make_graph(builder.H, id_by_site)
finalized_leads, lead_interfaces, lead_paddings =\
_finalize_leads(builder.leads, id_by_site)
del id_by_site # cleanup due to large size
site_arrays = _make_site_arrays(builder.H)
# We need this to efficiently find which array a given
# site belongs to
site_offsets = np.cumsum([0] + [len(arr) for arr in site_arrays])
(onsite_subgraphs, onsite_terms, onsite_term_values,
onsite_term_parameters, onsite_term_errors, _onsite_term_by_site_id) =\
_make_onsite_terms(builder, sites, site_offsets, term_offset=0)
(hopping_subgraphs, hopping_terms, hopping_term_values,
hopping_term_parameters, hopping_term_errors,
_hopping_term_by_edge_id) =\
_make_hopping_terms(builder, graph, sites, site_offsets,
len(sites), term_offset=len(onsite_terms))
# Construct the combined onsite/hopping term datastructures
subgraphs = tuple(onsite_subgraphs) + tuple(hopping_subgraphs)
terms = tuple(onsite_terms) + tuple(hopping_terms)
term_values = tuple(onsite_term_values) + tuple(hopping_term_values)
term_errors = tuple(onsite_term_errors) + tuple(hopping_term_errors)
# Construct system parameters
parameters = set()
for params in chain(onsite_term_parameters, hopping_term_parameters):
if params is not None:
parameters.update(params)
parameters = frozenset(parameters)
self.site_arrays = site_arrays
self.sites = _Sites(self.site_arrays)
self.id_by_site = _IdBySite(self.site_arrays)
self.graph = graph
self.subgraphs = subgraphs
self.terms = terms
self._term_values = term_values
self._term_errors = term_errors
self._hopping_term_by_edge_id = _hopping_term_by_edge_id
self._onsite_term_by_site_id = _onsite_term_by_site_id
self.parameters = parameters
self.symmetry = builder.symmetry
self.leads = finalized_leads
self.lead_interfaces = lead_interfaces
self.lead_paddings = lead_paddings
self._init_discrete_symmetries(builder)
def _make_lead_sites(builder, interface_order):
#### For each site of the fundamental domain, determine whether it has
#### neighbors in the previous domain or not.
sym = builder.symmetry
lsites_with = [] # Fund. domain sites with neighbors in prev. dom
lsites_without = [] # Remaining sites of the fundamental domain
for tail in builder.H: # Loop over all sites of the fund. domain.
for head in builder._out_neighbors(tail):
fd = sym.which(head)[0]
if fd == 1:
# Tail belongs to fund. domain, head to the next domain.
lsites_with.append(tail)
break
else:
# Tail is a fund. domain site not connected to prev. domain.
lsites_without.append(tail)
if not lsites_with:
warnings.warn('Infinite system with disconnected cells.',
RuntimeWarning, stacklevel=3)
### Create list of sites and a lookup table
minus_one = ta.array((-1,))
plus_one = ta.array((1,))
if interface_order is None:
# interface must be sorted
interface = [sym.act(minus_one, s) for s in lsites_with]
interface.sort()
else:
lsites_with_set = set(lsites_with)
lsites_with = []
interface = []
if interface_order:
shift = ta.array((-sym.which(interface_order[0])[0] - 1,))
for shifted_iface_site in interface_order:
# Shift the interface domain before the fundamental domain.
# That's the right place for the interface of a lead to be, but
# the sites of interface_order might live in a different
# domain.
iface_site = sym.act(shift, shifted_iface_site)
lsite = sym.act(plus_one, iface_site)
try:
lsites_with_set.remove(lsite)
except KeyError:
if (-sym.which(shifted_iface_site)[0] - 1,) != shift:
raise ValueError(
'The sites in interface_order do not all '
'belong to the same lead cell.')
else:
raise ValueError('A site in interface_order is not an '
'interface site:\n' + str(iface_site))
interface.append(iface_site)
lsites_with.append(lsite)
if lsites_with_set:
raise ValueError(
'interface_order did not contain all interface sites.')
# `interface_order` *must* be sorted, hence `interface` should also
if interface != sorted(interface):
raise ValueError('Interface sites must be sorted.')
del lsites_with_set
return sorted(lsites_with), sorted(lsites_without), interface
def _make_lead_graph(builder, sites, id_by_site, cell_size):
sym = builder.symmetry
g = graph.Graph()
g.num_nodes = len(sites) # Some sites could not appear in any edge.
for tail_id, tail in enumerate(sites[:cell_size]):
for head in builder._out_neighbors(tail):
head_id = id_by_site.get(head)
if head_id is None:
# Head belongs neither to the fundamental domain nor to the
# previous domain. Check that it belongs to the next
# domain and ignore it otherwise as an edge corresponding
# to this one has been added already or will be added.
fd = sym.which(head)[0]
if fd != 1:
msg = ('Further-than-nearest-neighbor cells '
'are connected by hopping\n{0}.')
raise ValueError(msg.format((tail, head)))
continue
if head_id >= cell_size:
# Head belongs to previous domain. The edge added here
# correspond to one left out just above.
g.add_edge(head_id, tail_id)
g.add_edge(tail_id, head_id)
return g.compressed()
class InfiniteSystem(_FinalizedBuilderMixin, system.InfiniteSystem):
......@@ -2106,10 +2819,10 @@ class InfiniteSystem(_FinalizedBuilderMixin, system.InfiniteSystem):
Attributes
----------
sites : sequence
``sites[i]`` is the `~kwant.builder.Site` instance that corresponds
``sites[i]`` is the `~kwant.system.Site` instance that corresponds
to the integer-labeled site ``i`` of the low-level system.
id_by_site : dict
The inverse of ``sites``; maps high-level `~kwant.builder.Site`
id_by_site : mapping
The inverse of ``sites``; maps high-level `~kwant.system.Site`
instances to their integer label.
Satisfies ``id_by_site[sites[i]] == i``.
......@@ -2133,69 +2846,12 @@ class InfiniteSystem(_FinalizedBuilderMixin, system.InfiniteSystem):
sym = builder.symmetry
assert sym.num_directions == 1
#### For each site of the fundamental domain, determine whether it has
#### neighbors in the previous domain or not.
lsites_with = [] # Fund. domain sites with neighbors in prev. dom
lsites_without = [] # Remaining sites of the fundamental domain
for tail in builder.H: # Loop over all sites of the fund. domain.
for head in builder._out_neighbors(tail):
fd = sym.which(head)[0]
if fd == 1:
# Tail belongs to fund. domain, head to the next domain.
lsites_with.append(tail)
break
else:
# Tail is a fund. domain site not connected to prev. domain.
lsites_without.append(tail)
lsites_with, lsites_without, interface =\
_make_lead_sites(builder, interface_order)
cell_size = len(lsites_with) + len(lsites_without)
if not lsites_with:
warnings.warn('Infinite system with disconnected cells.',
RuntimeWarning, stacklevel=3)
### Create list of sites and a lookup table
minus_one = ta.array((-1,))
plus_one = ta.array((1,))
if interface_order is None:
# interface must be sorted
interface = [sym.act(minus_one, s) for s in lsites_with]
interface.sort()
else:
lsites_with_set = set(lsites_with)
lsites_with = []
interface = []
if interface_order:
shift = ta.array((-sym.which(interface_order[0])[0] - 1,))
for shifted_iface_site in interface_order:
# Shift the interface domain before the fundamental domain.
# That's the right place for the interface of a lead to be, but
# the sites of interface_order might live in a different
# domain.
iface_site = sym.act(shift, shifted_iface_site)
lsite = sym.act(plus_one, iface_site)
try:
lsites_with_set.remove(lsite)
except KeyError:
if (-sym.which(shifted_iface_site)[0] - 1,) != shift:
raise ValueError(
'The sites in interface_order do not all '
'belong to the same lead cell.')
else:
raise ValueError('A site in interface_order is not an '
'interface site:\n' + str(iface_site))
interface.append(iface_site)
lsites_with.append(lsite)
if lsites_with_set:
raise ValueError(
'interface_order did not contain all interface sites.')
# `interface_order` *must* be sorted, hence `interface` should also
if interface != sorted(interface):
raise ValueError('Interface sites must be sorted.')
del lsites_with_set
# we previously sorted the interface, so don't sort it again
sites = sorted(lsites_with) + sorted(lsites_without) + interface
sites = lsites_with + lsites_without + interface
del lsites_with
del lsites_without
del interface
......@@ -2203,41 +2859,20 @@ class InfiniteSystem(_FinalizedBuilderMixin, system.InfiniteSystem):
for site_id, site in enumerate(sites):
id_by_site[site] = site_id
graph = _make_lead_graph(builder, sites, id_by_site, cell_size)
# In the following, because many onsites/hoppings share the same
# (value, parameter) pairs, we keep them in 'cache' so that we only
# store a given pair in memory *once*. This is like interning strings.
#### Make graph and extract onsite Hamiltonians.
#### Extract onsites
cache = _value_params_pair_cache(1)
g = graph.Graph()
g.num_nodes = len(sites) # Some sites could not appear in any edge.
onsites = []
for tail_id, tail in enumerate(sites[:cell_size]):
onsites.append(cache(builder.H[tail][1]))
for head in builder._out_neighbors(tail):
head_id = id_by_site.get(head)
if head_id is None:
# Head belongs neither to the fundamental domain nor to the
# previous domain. Check that it belongs to the next
# domain and ignore it otherwise as an edge corresponding
# to this one has been added already or will be added.
fd = sym.which(head)[0]
if fd != 1:
msg = ('Further-than-nearest-neighbor cells '
'are connected by hopping\n{0}.')
raise ValueError(msg.format((tail, head)))
continue
if head_id >= cell_size:
# Head belongs to previous domain. The edge added here
# correspond to one left out just above.
g.add_edge(head_id, tail_id)
g.add_edge(tail_id, head_id)
g = g.compressed()
onsites = [cache(builder.H[tail][1]) for tail in sites[:cell_size]]
#### Extract hoppings.
cache = _value_params_pair_cache(2)
hoppings = []
for tail_id, head_id in g:
for tail_id, head_id in graph:
tail = sites[tail_id]
head = sites[head_id]
if tail_id >= cell_size:
......@@ -2264,7 +2899,7 @@ class InfiniteSystem(_FinalizedBuilderMixin, system.InfiniteSystem):
else:
parameters = frozenset(parameters)
self.graph = g
self.graph = graph
self.sites = sites
self.site_ranges = _site_ranges(sites)
self.id_by_site = id_by_site
......@@ -2275,6 +2910,125 @@ class InfiniteSystem(_FinalizedBuilderMixin, system.InfiniteSystem):
self.cell_size = cell_size
self._init_discrete_symmetries(builder)
def hamiltonian(self, i, j, *args, params=None):
cs = self.cell_size
if i == j >= cs:
i -= cs
j -= cs
return super().hamiltonian(i, j, *args, params=params)
class InfiniteVectorizedSystem(_VectorizedFinalizedBuilderMixin, system.InfiniteVectorizedSystem):
"""Finalized infinite system, extracted from a `Builder`.
Attributes
----------
site_arrays : sequence of `~kwant.system.SiteArray`
The sites of the system.
sites : sequence
``sites[i]`` is the `~kwant.system.Site` instance that corresponds
to the integer-labeled site ``i`` of the low-level system.
id_by_site : mapping
The inverse of ``sites``; maps high-level `~kwant.system.Site`
instances to their integer label.
Satisfies ``id_by_site[sites[i]] == i``.
Notes
-----
In infinite systems ``sites`` and ``site_arrays`` consists of 3 parts:
sites in the fundamental domain (FD) with hoppings to neighboring cells,
sites in the FD with no hoppings to neighboring cells, and sites in FD+1
attached to the FD by hoppings. Each of these three subsequences is
individually sorted.
"""
def __init__(self, builder, interface_order=None):
"""
Finalize a builder instance which has to have exactly a single
symmetry direction.
If interface_order is not set, the order of the interface sites in the
finalized system will be arbitrary. If interface_order is set to a
sequence of interface sites, this order will be kept.
"""
sym = builder.symmetry
assert sym.num_directions == 1
assert builder.vectorize
lsites_with, lsites_without, interface =\
_make_lead_sites(builder, interface_order)
cell_size = len(lsites_with) + len(lsites_without)
sites = lsites_with + lsites_without + interface
id_by_site = {}
for site_id, site in enumerate(sites):
id_by_site[site] = site_id
# these are needed for constructing hoppings
lsites_with = frozenset(lsites_with)
lsites_without = frozenset(lsites_without)
interface = frozenset(interface)
if not all(s.family.norbs for s in sites):
raise ValueError("All sites must belong to families with norbs "
"specified for vectorized Builders. Specify "
"norbs when creating site families.")
graph = _make_lead_graph(builder, sites, id_by_site, cell_size)
del id_by_site # cleanup due to large size
# In order to conform to the kwant.system.InfiniteVectorizedSystem
# interface we need to put the sites that connect to the previous
# cell *first*, then the sites that do not couple to the previous
# cell, then the sites in the previous cell. Because sites in
# a SiteArray are sorted by tag this means that the sites in these
# 3 different sets need to be in different SiteArrays.
site_arrays = (
_make_site_arrays(lsites_with)
+ _make_site_arrays(lsites_without)
+ _make_site_arrays(interface)
)
site_offsets = np.cumsum([0] + [len(arr) for arr in site_arrays])
(onsite_subgraphs, onsite_terms, onsite_term_values,
onsite_term_parameters, onsite_term_errors, _onsite_term_by_site_id) =\
_make_onsite_terms(builder, sites, site_offsets, term_offset=0)
(hopping_subgraphs, hopping_terms, hopping_term_values,
hopping_term_parameters, hopping_term_errors,
_hopping_term_by_edge_id) =\
_make_hopping_terms(builder, graph, sites, site_offsets,
cell_size, term_offset=len(onsite_terms))
# Construct the combined onsite/hopping term datastructures
subgraphs = tuple(onsite_subgraphs) + tuple(hopping_subgraphs)
terms = tuple(onsite_terms) + tuple(hopping_terms)
term_values = tuple(onsite_term_values) + tuple(hopping_term_values)
term_errors = tuple(onsite_term_errors) + tuple(hopping_term_errors)
# Construct system parameters
parameters = set()
for params in chain(onsite_term_parameters, hopping_term_parameters):
if params is not None:
parameters.update(params)
parameters = frozenset(parameters)
self.site_arrays = site_arrays
self.sites = _Sites(self.site_arrays)
self.id_by_site = _IdBySite(self.site_arrays)
self.graph = graph
self.subgraphs = subgraphs
self.terms = terms
self._term_values = term_values
self._term_errors = term_errors
self._hopping_term_by_edge_id = _hopping_term_by_edge_id
self._onsite_term_by_site_id = _onsite_term_by_site_id
self.parameters = parameters
self.symmetry = builder.symmetry
self.cell_size = cell_size
self._init_discrete_symmetries(builder)
def hamiltonian(self, i, j, *args, params=None):
cs = self.cell_size
......@@ -2283,5 +3037,15 @@ class InfiniteSystem(_FinalizedBuilderMixin, system.InfiniteSystem):
j -= cs
return super().hamiltonian(i, j, *args, params=params)
def pos(self, i):
return self.sites[i].pos
def is_finite_system(syst):
return isinstance(syst, (FiniteSystem, FiniteVectorizedSystem))
def is_infinite_system(syst):
return isinstance(syst, (FiniteSystem, FiniteVectorizedSystem))
def is_system(syst):
return isinstance(syst, (FiniteSystem, FiniteVectorizedSystem,
InfiniteSystem, InfiniteVectorizedSystem))