Qsymm

Symmetry finder and symmetric Hamiltonian generator

qsymm is an open-source Python library that makes symmetry analysis simple.

It automatically generates model Hamiltonians from symmetry constraints and finds the full symmetry group of your Hamiltonian.

Check out the introductory tutorial to see examples of how to use qsymm.

Implemented algorithms

The two core concepts in qsymm are Hamiltonian families (Hamiltonians that may depend on free parameters) and symmetries. We provide powerful classes to handle these:

• Model is used to store symbolic Hamiltonians that may depend on momenta and other free parameters. We use sympy for symbolic manipulation, but our implementation utilizes numpy arrays for efficient calculations with matrix valued functions.

• PointGroupElement and ContinuousGroupGenerator are used to store symmetry operators. Besides the ability to combine symmetries, they can also be applied to a Model to transform it.

We implement algorithms that form a two-way connection between Hamiltonian families and symmetries.

• Symmetry finding is handled by symmetries, it takes a Model as input and finds all of its symmetries, including conserved quantities, time reversal, particle-hole, and spatial rotation symmetries. See the symmetry finder tutorial and the kekule tutorial for detailed examples.

• continuum_hamiltonian and bloch_family are used to generate k.p or lattice Hamiltonians from symmetry constraints. See the k-dot-p generator tutorial, the Bloch generator tutorial and the kekule tutorial for detailed examples.

Installation

qsymm works with Python 3.5 and is available on PyPI:

pip install qsymm

Some of the example notebooks also require Kwant.

Documentation

Qsymm’s documentation is hosted on Read the Docs

Citing

Check out CITING.md for instructions on how to cite Qsymm in your publications.

Development

qsymm is on Gitlab, visit there if you would like to to contribute, report issues, or get the latest development version.

Citing Qsymm

We provide Qsymm as free software under a BSD license. If you have used Qsymm for work that has lead to a scientific publication, please mention the fact that you used it explicitly in the text body.

In addition we ask you to cite the Qsymm paper:

Dániel Varjas, Tómas Ö Rosdahl, and Anton R Akhmerov
Qsymm: algorithmic symmetry finding and symmetric Hamiltonian generation
New J. Phys. 20 093026 (2018)

Changelog

All notable changes to this project will be documented in this file.

The format is based on Keep a Changelog, and this project adheres to Semantic Versioning.

[unreleased]

[1.2.7] - 2019-12-04

Added

• Added integration test: run the tests for Kwant’s qsymm module against the current version of qsymm.

• Added tests for fail cases in symmetry finder and Model with locals.

Fixed

• Bug when Model is initialized with numpy arrays in locals.

• Bug in symmetry finding when transformed Hamiltonian has very small coefficients, this could lead to not finding symmetries in continuum models.

• Bug in Model.copy with sparse arrays did not return a copy of the stored arrays.

• Move docstrings to class level from init functions and extend docstrings.

[1.2.6] - 2019-11-12

Fixed

• Downgraded scipy version requirement to 0.19, to maintain compatibility with Kwant

• Correctly set model format when reshaping models with sparse values

[1.2.5] - 2019-11-11

Added

• Added basic reference documentation

• Added citation instructions in CITING.md

• Added a tutorial with downloadable Python files and Jupyter notebooks

Fixed

• Allow ‘display_family’ to work when IPython is not installed

• Bug when multiplying a model by a sympy expression (undefined variable)

• Bug when multiplying BlochCoeff with sympy expression (sympy multiplication was being used, when BlochCoeff multiplication should have been used)

• Remove copying of Sympy symbols (the hash values for the original and copy is different in some cases (probably a bug in Sympy)), which broke comparison between Models.

• Remove ‘real’ assumption from Sympy symbols in ‘hamiltonian_from_family’; this makes all Sympy symbols in Qsymm consistently have no assumptions.

• Minor clarifications to various docstrings

Removed

• Removed the notebooks from the source repository

[1.2.4] - 2019-10-17

Fixed

• Set correct content-type on long description metadata

[1.2.3] - 2019-10-17

Fixed

• Add long description to package metadata

[1.2.2] - 2019-10-17

Fixed

• Adding zero to a model is now idempotent

[1.2.1] - 2019-09-02

Added

• Add CHANGELOG.md

Fixed

• Restore compatibility with Python 3.5

[1.2.0] - 2019-08-30

Added

• Factory functions for discrete symmetries: time_reversal, particle_hole, chiral, inversion, rotation, mirror.

• Better representation of PointGroupElements and ContinuousGroupGenerator, using _repr_pretty_ and _repr_latex_. Remove print_PG_elements and update example notebooks.

• Implemented new functionality in Model:

  • Implement __matmul__ (@).

  • Support of sparse matrices and LinearOperator as values.

  • Consistent support of scalar valued Models.

  • Add keep to only keep track of certain symbolic coefficients.

  • More options and more transparent initialization, allow string keys which are automatically symmpified by default.

  • Several new utility functions, such as trace, reshape, allclose etc.

Changed

• Slight change to the internal representation of PointGroupElements, allow mixing integer tinyarray with either sympy matrix or floating point tinyarray in the rotation part, but do not allow mixing the latter two representations. This removes the need to have two different representations for simple operators. Optimize multiplication by caching.

• Changes in the API of Model:

  • Add format attribute to keep track of the type of data entries (i.e scalar, dense or sparse array).

  • Change the behaviour of * for Model objects, matrix multiplication only works with @ from now on. This breaks backward compatibility, a few uses in the qsymm code base were fixed.

  • Stop supporting a sequence of ints as momenta, instead string keys (e.g. 'k_x') or sympy expressions can be used. Certain other, rarely used ways to initialize Model don’t work anymore, some tests were changed accordingly.

  • Stop rounding and removing small entries from Model objects. This means that testing that the Model is empty is not a good test anymore to see if it is approximately zero. Use Model.allclose instead.

• Optimizations in Model:

  • Remove unnecessary deep copying, which was slow.

  • Optimize the implementation of arithmetic operations to minimize number of loops and function calls.

  • Fast initialization by making restructuring optional when initialized with a dict.

  • Clean up the code of BlochModel to utilize improvements in Model.

  • Update symmetry finder to work with sparse models.

Deprecated

• Deprecate initializing empty Model without providing shape and format.

[1.1.0] - 2018-12-05

Many of these changes were made in anticipation of integrating qsymm with kwant, allowing conversion between the Hamiltonian formats used.

Added

• BlochModel as a subclass of Model in order to store Bloch Hamiltonians with floating point hopping vectors.

• bravais_point_group to find the point group of a Bravais lattice using the translation vectors only. This is intended for generating candidades for symmetries.

Changed

• Change the way equality of Models is tested, the current implementation treats tolerances more consistently, this fixes some bugs in symmetry finding.

• Allow using floating point rotation matrices in PointGroupElement.

Fixed

• Several bugs and internal code restructuring.

Qsymm Tutorial

Qsymm Basics

See also

The complete source code of this example can be found in basics.py. A Jupyter notebook can be found in basics.ipynb.

Getting started with Qsymm is as simple as importing it:

import qsymm

To make effective use of Qsymm we’ll also need a few other utilities: numpy handling numeric arrays, and sympy for symbolic mathematics:

import numpy as np
import sympy

In all the following tutorials we will use these standard imports, and they won’t be explicitly shown

Defining a Qsymm model

Let’s start by defining a 3D Rashba Hamiltonian symbolically as a Python string:

ham = ("hbar^2 / (2 * m) * (k_x**2 + k_y**2 + k_z**2) * eye(2) +" +
        "alpha * sigma_x * k_x + alpha * sigma_y * k_y + alpha * sigma_z * k_z")

We can then create a Qsymm Model directly from this symbolic Hamiltonian:

H = qsymm.Model(ham)

We can then directly inspect the contents by printing the Model:

print(H)

{alpha*k_z:
[[ 1.+0.j  0.+0.j]
 [ 0.+0.j -1.+0.j]],

hbar**2*k_x**2/m:
[[0.5+0.j 0. +0.j]
 [0. +0.j 0.5+0.j]],

hbar**2*k_y**2/m:
[[0.5+0.j 0. +0.j]
 [0. +0.j 0.5+0.j]],

hbar**2*k_z**2/m:
[[0.5+0.j 0. +0.j]
 [0. +0.j 0.5+0.j]],

alpha*k_x:
[[0.+0.j 1.+0.j]
 [1.+0.j 0.+0.j]],

alpha*k_y:
[[0.+0.j 0.-1.j]
 [0.+1.j 0.+0.j]],

}

We can also extract a more readable representation by using the tosympy method, which converts the Model to a sympy expression:

H.tosympy(nsimplify=True)

\[\begin{split}\displaystyle \left[\begin{matrix}\alpha k_{z} + \frac{\hbar^{2} k_{x}^{2}}{2 m} + \frac{\hbar^{2} k_{y}^{2}}{2 m} + \frac{\hbar^{2} k_{z}^{2}}{2 m} & \alpha k_{x} - i \alpha k_{y}\\\alpha k_{x} + i \alpha k_{y} & - \alpha k_{z} + \frac{\hbar^{2} k_{x}^{2}}{2 m} + \frac{\hbar^{2} k_{y}^{2}}{2 m} + \frac{\hbar^{2} k_{z}^{2}}{2 m}\end{matrix}\right]\end{split}\]

The argument nsimplify=True makes the output more readable by forcing sympy to elide factors of 1.0 that multiply each term. Note that Qsymm automatically interprets the symbols sigma_x, sigma_y and sigma_z as the Pauli matrices, and eye(2) as the 2x2 identity matrix.

Model as a momenta attribute that specifies which symbols are considered the momentum variables:

H.momenta

(k_x, k_y, k_z)

By default Qsymm assumes that your model is written in 3D (even if it does not include all 3 momenta). To define a lower-dimensional model you must explicitly specify the momentum variables, e.g:

ham2D = ("hbar^2 / (2 * m) * (k_x**2 + k_z**2) * eye(2) +" +
         "alpha * sigma_x * k_x + alpha * sigma_y * k_z")
H2D = qsymm.Model(ham2D, momenta=['k_x', 'k_z'])

H2D.tosympy(nsimplify=True)

\[\begin{split}\displaystyle \left[\begin{matrix}\frac{\hbar^{2} k_{x}^{2}}{2 m} + \frac{\hbar^{2} k_{z}^{2}}{2 m} & \alpha k_{x} - i \alpha k_{z}\\\alpha k_{x} + i \alpha k_{z} & \frac{\hbar^{2} k_{x}^{2}}{2 m} + \frac{\hbar^{2} k_{z}^{2}}{2 m}\end{matrix}\right]\end{split}\]

H2D.momenta

(k_x, k_z)

Defining group elements

Qsymm is all about finding and generating symmetries of models, so it is unsurprising that it contains utilities for defining group elements.

Below are a few examples of the sorts of things you can define with Qsymm:

# Identity in 3D
E = qsymm.identity(3)
# Inversion in 3D
I = qsymm.inversion(3)
# 4-fold rotation around the x-axis
C4 = qsymm.rotation(1/4, [1, 0, 0])
# 3-fold rotation around the [1, 1, 1] axis
C3 = qsymm.rotation(1/3, [1, 1, 1])
# Time reversal
TR = qsymm.time_reversal(3)
# Particle-hole
PH = qsymm.particle_hole(3)

The documentation page of the qsymm.groups module contains an exhaustive list of what can be generated.

As with other Qsymm objects we can get a readable representation of these group elements:

C4

\[R\left(\frac{\pi}{2}, [1 0 0]\right)\]

TR

\[ \mathcal{T}\]

Given a set of group generators we can also generate a group:

cubic_gens = {I, C4, C3, TR, PH}
cubic_group = qsymm.groups.generate_group(cubic_gens)
print(len(cubic_group))

192

Group elements can be multiplied and inverted, as we would expect:

C3 * C4

\[R\left(\pi, [1 1 0]\right)\]

C3**-1

\[R\left(-\frac{2\pi}{3}, [1 1 1]\right)\]

We can also apply group elements to the Model that we defined in the previous section:

H_with_TR = TR.apply(H)
H_with_TR.tosympy(nsimplify=True)

\[\begin{split}\displaystyle \left[\begin{matrix}- \alpha k_{z} + \frac{\hbar^{2} k_{x}^{2}}{2 m} + \frac{\hbar^{2} k_{y}^{2}}{2 m} + \frac{\hbar^{2} k_{z}^{2}}{2 m} & - \alpha k_{x} - i \alpha k_{y}\\- \alpha k_{x} + i \alpha k_{y} & \alpha k_{z} + \frac{\hbar^{2} k_{x}^{2}}{2 m} + \frac{\hbar^{2} k_{y}^{2}}{2 m} + \frac{\hbar^{2} k_{z}^{2}}{2 m}\end{matrix}\right]\end{split}\]

Defining continuous group generators

In addition to the group elements we can also define generators of continuous groups using qsymm.groups.ContinuousGroupGenerator:

sz = qsymm.ContinuousGroupGenerator(None, np.array([[1, 0], [0, -1]]))

The first argument to ContinuousGroupGenerator is the realspace rotation generator; by specifying None we indicate that we want the rotation part to be zero. The second argument is the unitary action of the generator on the Hilbert space as a Hermitian matrix.

Applying a ContinuousGroupGenerator to a Model calculates the commutator:

sz.apply(H).tosympy(nsimplify=True)

\[\begin{split}\displaystyle \left[\begin{matrix}0 & - 2 i \alpha k_{x} - 2 \alpha k_{y}\\2 i \alpha k_{x} - 2 \alpha k_{y} & 0\end{matrix}\right]\end{split}\]

For the 3D Rashba Hamiltonian we defined at the start of the tutorial spin-z is not conserved, hence the commutator is non-zero.

Finding symmetries

The function symmetries allows us to find the symmetries of a Model. Let us find out whether the 3D Rashba Hamiltonian defined earlier has cubic group symmetry:

discrete_symm, continuous_symm = qsymm.symmetries(H, cubic_group)
print(len(discrete_symm), len(continuous_symm))

48 0

It has 48 discrete symmetries (cubic group without inversion and time-reversal) and no continuous symmetries (conserved quantities).

For more detailed examples see Finding Symmetries and Finding symmetries of the Kekule-Y continuum model.

Generating Hamiltonians from symmetry constraints

The qsymm.hamiltonian_generator module contains algorithms for generating Hamiltonians from symmetry constraints.

For example let us generate all 2-band \(k \cdot p\) Hamiltonians with the same discrete symmetries as the Rashba Hamiltonian that we found in the previous section:

family = qsymm.continuum_hamiltonian(discrete_symm, dim=3, total_power=2, prettify=True)
qsymm.display_family(family)

\[\begin{split}\displaystyle \left[\begin{matrix}1 & 0\\0 & 1\end{matrix}\right]\end{split}\]

\[\begin{split}\displaystyle \left[\begin{matrix}- k_{z} & - k_{x} + i k_{y}\\- k_{x} - i k_{y} & k_{z}\end{matrix}\right]\end{split}\]

\[\begin{split}\displaystyle \left[\begin{matrix}k_{x}^{2} + k_{y}^{2} + k_{z}^{2} & 0\\0 & k_{x}^{2} + k_{y}^{2} + k_{z}^{2}\end{matrix}\right]\end{split}\]

It is exactly the Hamiltonian family we started with.

For more detailed examples see Generating k \cdot p models, Generating tight-binding models and Finding symmetries of the Kekule-Y continuum model.

Finding Symmetries

See also

The complete source code of this example can be found in symmetry_finder.py. A Jupyter notebook can be found in symmetry_finder.ipynb.

A basic example

Let’s start from the same 3D Rashba hamiltonian we used in the basic tutorial.

ham = ("hbar^2 / (2 * m) * (k_x**2 + k_y**2 + k_z**2) * eye(2) +" +
        "alpha * sigma_x * k_x + alpha * sigma_y * k_y + alpha * sigma_z * k_z")
# Convert to standard monomials form
H = qsymm.Model(ham)
H.tosympy()

\[\begin{split}\displaystyle \left[\begin{matrix}1.0 \alpha k_{z} + \frac{0.5 \hbar^{2} k_{x}^{2}}{m} + \frac{0.5 \hbar^{2} k_{y}^{2}}{m} + \frac{0.5 \hbar^{2} k_{z}^{2}}{m} & 1.0 \alpha k_{x} - 1.0 i \alpha k_{y}\\1.0 \alpha k_{x} + 1.0 i \alpha k_{y} & - 1.0 \alpha k_{z} + \frac{0.5 \hbar^{2} k_{x}^{2}}{m} + \frac{0.5 \hbar^{2} k_{y}^{2}}{m} + \frac{0.5 \hbar^{2} k_{z}^{2}}{m}\end{matrix}\right]\end{split}\]

We start from the cubic group as the set of candidates for point group symmetries:

cubic_group = qsymm.groups.cubic()

Then we can use symmetries to find point group symmetries, discrete onsite symmetries, and continuous symmetries:

sg, cg = qsymm.symmetries(H, cubic_group)

As before we see that we have 48 symmetries, but no continuous symmetries:

print(len(sg), len(cg))

48 0

If we take a look at a few of the symmetries, we see that it includes time reversal but not inversion:

from IPython.display import display

for i in range(1, 5):
  display(sg[-i])

\[R\left(-\frac{\pi}{2}, [1 0 0]\right) \mathcal{T}\]

\[R\left(\frac{\pi}{2}, [1 0 0]\right) \mathcal{T}\]

\[R\left(\pi, [1 0 0]\right) \mathcal{T}\]

\[R\left(\pi, [1 1 0]\right) \mathcal{T}\]

We can also get a more detailed representation of these group elements:

from IPython.display import Math

display(Math(qsymm.groups.pretty_print_pge(sg[28], latex=True, full=True)))

\[\begin{split}\displaystyle U H(\mathbf{k})^* U^{-1} = H(-R\mathbf{k}) \\R = R\left(\pi, [0 1 0]\right)\\U = \begin{bmatrix} 1.+0.j&-0.+0.j\\-0.+0.j& 1.+0.j\\\end{bmatrix}\end{split}\]

We can check that this is the same as the full cubic group without inversion, plus time-reversal

C4 = qsymm.rotation(1/4, [1, 0, 0])
C3 = qsymm.rotation(1/3, [1, 1, 1])
TR = qsymm.time_reversal(3)

set(sg) == qsymm.groups.generate_group({C3, C4, TR})

True

Adding more elements

Now we’re going to add a 2-fold degeneracy to the above hamiltonian by taking the kroneker product with the 2x2 identity:

ham_degenerate = "kron(eye(2), " + ham + ")"
H_degenerate = qsymm.Model(ham_degenerate)

If we look for the symmetries in the same way as previously we will see that the above modification adds a continuous SU(2) symmetry:

sg, cg = qsymm.symmetries(H_degenerate, cubic_group)

cg

\[\begin{split}\left[ H(\mathbf{{k}}), L \right] = 0 \\L = \begin{bmatrix}1.+0.j&0.+0.j&0.+0.j&0.+0.j\\0.+0.j&1.+0.j&0.+0.j&0.+0.j\\ 0.+0.j& 0.+0.j&-1.+0.j& 0.+0.j\\ 0.+0.j& 0.+0.j& 0.+0.j&-1.+0.j\\\end{bmatrix} \\\\\left[ H(\mathbf{{k}}), L \right] = 0 \\L = \begin{bmatrix}0.+0.j&0.+0.j&1.+0.j&0.+0.j\\0.+0.j&0.+0.j&0.+0.j&1.+0.j\\1.+0.j&0.+0.j&0.+0.j&0.+0.j\\0.+0.j&1.+0.j&0.+0.j&0.+0.j\\\end{bmatrix} \\\\\left[ H(\mathbf{{k}}), L \right] = 0 \\L = \begin{bmatrix}0.+0.j&0.+0.j&0.+1.j&0.+0.j\\0.+0.j&0.+0.j&0.+0.j&0.+1.j\\0.-1.j&0.+0.j&0.+0.j&0.+0.j\\0.+0.j&0.-1.j&0.+0.j&0.+0.j\\\end{bmatrix} \\\end{split}\]

If we instead add a particle-hole degree of freedom by using \(σ_z\) instead of the 2x2 identity then we instead get a continuous U(1) (charge conservation) symmetry:

ham_ph = "kron(sigma_z, " + ham + ")"
H_ph = qsymm.Model(ham_ph)
sg, cg = qsymm.symmetries(H_ph, cubic_group)
print(len(cg))

1

cg

\[\begin{split}\left[ H(\mathbf{{k}}), L \right] = 0 \\L = \begin{bmatrix}-0.5+0.j& 0. +0.j& 0. +0.j& 0. +0.j\\ 0. +0.j&-0.5+0.j& 0. +0.j& 0. +0.j\\0. +0.j&0. +0.j&0.5+0.j&0. +0.j\\0. +0.j&0. +0.j&0. +0.j&0.5+0.j\\\end{bmatrix} \\\end{split}\]

Rashba Hamiltonian with \(J = 3/2\)

J_x, J_y, J_z = qsymm.groups.spin_matrices(3/2)
ham_rashba = (
    "hbar^2 / (2 * m) * (k_x**2 + k_y**2 + k_z**2) * eye(4) +" +
    "alpha * J_x * k_x + alpha * J_y * k_y + alpha * J_z * k_z"
)
H_rashba = qsymm.Model(ham_rashba, locals=dict(J_x=J_x, J_y=J_y, J_z=J_z))
H_rashba.tosympy(nsimplify=True)

\[\begin{split}\displaystyle \left[\begin{matrix}\frac{3 \alpha k_{z}}{2} + \frac{\hbar^{2} k_{x}^{2}}{2 m} + \frac{\hbar^{2} k_{y}^{2}}{2 m} + \frac{\hbar^{2} k_{z}^{2}}{2 m} & \frac{\sqrt{3} \alpha k_{x}}{2} - \frac{\sqrt{3} i \alpha k_{y}}{2} & 0 & 0\\\frac{\sqrt{3} \alpha k_{x}}{2} + \frac{\sqrt{3} i \alpha k_{y}}{2} & \frac{\alpha k_{z}}{2} + \frac{\hbar^{2} k_{x}^{2}}{2 m} + \frac{\hbar^{2} k_{y}^{2}}{2 m} + \frac{\hbar^{2} k_{z}^{2}}{2 m} & \alpha k_{x} - i \alpha k_{y} & 0\\0 & \alpha k_{x} + i \alpha k_{y} & - \frac{\alpha k_{z}}{2} + \frac{\hbar^{2} k_{x}^{2}}{2 m} + \frac{\hbar^{2} k_{y}^{2}}{2 m} + \frac{\hbar^{2} k_{z}^{2}}{2 m} & \frac{\sqrt{3} \alpha k_{x}}{2} - \frac{\sqrt{3} i \alpha k_{y}}{2}\\0 & 0 & \frac{\sqrt{3} \alpha k_{x}}{2} + \frac{\sqrt{3} i \alpha k_{y}}{2} & - \frac{3 \alpha k_{z}}{2} + \frac{\hbar^{2} k_{x}^{2}}{2 m} + \frac{\hbar^{2} k_{y}^{2}}{2 m} + \frac{\hbar^{2} k_{z}^{2}}{2 m}\end{matrix}\right]\end{split}\]

Note in the above that we had to tell the Model about the \(J_x\), \(J_y\) and \(J_z\) matrices by passing a dictionary locals, which maps symbols in ham_rashba to the corresponding matrices.

We again find the symmetries of this model, starting the search from the full cubic group:

sg, cg = qsymm.symmetries(H_rashba, cubic_group)

The symmetry group should be the same as the full cubic group without inversion, plus time-reversal:

set(sg) == qsymm.groups.generate_group({C4, C3, TR})

True

Bloch Hamiltonian

Let’s start by defining the hexagonal point group in 2D:

hex_group_2D = qsymm.groups.hexagonal()

Next we define a model for a single-band tight-binding model with nearest-neighbor hopping on a triangular lattice:

ham_tri = 'm * (cos(k_x) + cos(1/2*k_x + sqrt(3)/2*k_y) + cos(-1/2*k_x + sqrt(3)/2*k_y))'
display(qsymm.sympify(ham_tri))

\[\displaystyle m \left(\cos{\left(k_{x} \right)} + \cos{\left(- \frac{k_{x}}{2} + \frac{\sqrt{3} k_{y}}{2} \right)} + \cos{\left(\frac{k_{x}}{2} + \frac{\sqrt{3} k_{y}}{2} \right)}\right)\]

H_tri = qsymm.Model(ham_tri, momenta=['k_x', 'k_y'])

We find the symmetries, starting from the full hexagonal point group:

sg, cg = qsymm.symmetries(H_tri, hex_group_2D)

and verify that the point group of the model is the hexagonal group, without particle-hole symmetry:

set(sg) == qsymm.groups.hexagonal(ph=False)

True

Continuous rotation symmetry

The Rashba Hamiltonian we defined at the start of this tutorial actually has full rotation invariance:

display(qsymm.sympify(ham))
pg, cg = qsymm.symmetries(H, continuous_rotations=True, prettify=True)

\[\begin{split}\displaystyle \left[\begin{matrix}\alpha k_{z} + \frac{\hbar^{2} \left(k_{x}^{2} + k_{y}^{2} + k_{z}^{2}\right)}{2 m} & \alpha k_{x} - i \alpha k_{y}\\\alpha k_{x} + i \alpha k_{y} & - \alpha k_{z} + \frac{\hbar^{2} \left(k_{x}^{2} + k_{y}^{2} + k_{z}^{2}\right)}{2 m}\end{matrix}\right]\end{split}\]

cg

\[\begin{split}e^{-i\phi L}H(\mathbf{k})e^{i\phi L} = H(R_{\phi}\mathbf{k}) \\R_{\phi} = R\left(\phi, [0 0 1]\right)\\L = \begin{bmatrix}0.5+0.j&0. +0.j\\ 0. +0.j&-0.5+0.j\\\end{bmatrix} \\\\e^{-i\phi L}H(\mathbf{k})e^{i\phi L} = H(R_{\phi}\mathbf{k}) \\R_{\phi} = R\left(\phi, [ 0 -1 0]\right)\\L = \begin{bmatrix}0.+0.j &0.+0.5j\\0.-0.5j&0.+0.j \\\end{bmatrix} \\\\e^{-i\phi L}H(\mathbf{k})e^{i\phi L} = H(R_{\phi}\mathbf{k}) \\R_{\phi} = R\left(\phi, [1 0 0]\right)\\L = \begin{bmatrix}0. +0.j&0.5+0.j\\0.5+0.j&0. +0.j\\\end{bmatrix} \\\end{split}\]

\(k \cdot p\) Hamiltonian

Here we start from an 8x8 \(k \cdot p\) Hamiltonian that models a zinc-blende semiconductor:

ham_kp = 'Matrix([[hbar**2*k_x*gamma_0*k_x/(2*m_0)+hbar**2*k_y*gamma_0*k_y/(2*m_0)+hbar**2*k_z*gamma_0*k_z/(2*m_0)+E_0+E_v,0,-sqrt(2)*P*k_x/2-sqrt(2)*I*P*k_y/2,sqrt(6)*P*k_z/3,sqrt(6)*P*k_x/6-sqrt(6)*I*P*k_y/6,0,-sqrt(3)*P*k_z/3,-sqrt(3)*P*k_x/3+sqrt(3)*I*P*k_y/3],[0,hbar**2*k_x*gamma_0*k_x/(2*m_0)+hbar**2*k_y*gamma_0*k_y/(2*m_0)+hbar**2*k_z*gamma_0*k_z/(2*m_0)+E_0+E_v,0,-sqrt(6)*P*k_x/6-sqrt(6)*I*P*k_y/6,sqrt(6)*P*k_z/3,sqrt(2)*P*k_x/2-sqrt(2)*I*P*k_y/2,-sqrt(3)*P*k_x/3-sqrt(3)*I*P*k_y/3,sqrt(3)*P*k_z/3],[-sqrt(2)*k_x*P/2+sqrt(2)*I*k_y*P/2,0,-hbar**2*k_x*gamma_1*k_x/(2*m_0)-hbar**2*k_x*gamma_2*k_x/(2*m_0)-I*hbar**2*k_x*k_y/(2*m_0)-3*I*hbar**2*k_x*kappa*k_y/(2*m_0)-hbar**2*k_y*gamma_1*k_y/(2*m_0)-hbar**2*k_y*gamma_2*k_y/(2*m_0)+I*hbar**2*k_y*k_x/(2*m_0)+3*I*hbar**2*k_y*kappa*k_x/(2*m_0)-hbar**2*k_z*gamma_1*k_z/(2*m_0)+hbar**2*k_z*gamma_2*k_z/m_0+E_v,sqrt(3)*hbar**2*k_x*gamma_3*k_z/(2*m_0)+sqrt(3)*hbar**2*k_x*k_z/(6*m_0)+sqrt(3)*hbar**2*k_x*kappa*k_z/(2*m_0)-sqrt(3)*I*hbar**2*k_y*gamma_3*k_z/(2*m_0)-sqrt(3)*I*hbar**2*k_y*k_z/(6*m_0)-sqrt(3)*I*hbar**2*k_y*kappa*k_z/(2*m_0)+sqrt(3)*hbar**2*k_z*gamma_3*k_x/(2*m_0)-sqrt(3)*I*hbar**2*k_z*gamma_3*k_y/(2*m_0)-sqrt(3)*hbar**2*k_z*k_x/(6*m_0)+sqrt(3)*I*hbar**2*k_z*k_y/(6*m_0)-sqrt(3)*hbar**2*k_z*kappa*k_x/(2*m_0)+sqrt(3)*I*hbar**2*k_z*kappa*k_y/(2*m_0),sqrt(3)*hbar**2*k_x*gamma_2*k_x/(2*m_0)-sqrt(3)*I*hbar**2*k_x*gamma_3*k_y/(2*m_0)-sqrt(3)*hbar**2*k_y*gamma_2*k_y/(2*m_0)-sqrt(3)*I*hbar**2*k_y*gamma_3*k_x/(2*m_0),0,-sqrt(6)*hbar**2*k_x*gamma_3*k_z/(4*m_0)-sqrt(6)*hbar**2*k_x*k_z/(12*m_0)-sqrt(6)*hbar**2*k_x*kappa*k_z/(4*m_0)+sqrt(6)*I*hbar**2*k_y*gamma_3*k_z/(4*m_0)+sqrt(6)*I*hbar**2*k_y*k_z/(12*m_0)+sqrt(6)*I*hbar**2*k_y*kappa*k_z/(4*m_0)-sqrt(6)*hbar**2*k_z*gamma_3*k_x/(4*m_0)+sqrt(6)*I*hbar**2*k_z*gamma_3*k_y/(4*m_0)+sqrt(6)*hbar**2*k_z*k_x/(12*m_0)-sqrt(6)*I*hbar**2*k_z*k_y/(12*m_0)+sqrt(6)*hbar**2*k_z*kappa*k_x/(4*m_0)-sqrt(6)*I*hbar**2*k_z*kappa*k_y/(4*m_0),-sqrt(6)*hbar**2*k_x*gamma_2*k_x/(2*m_0)+sqrt(6)*I*hbar**2*k_x*gamma_3*k_y/(2*m_0)+sqrt(6)*hbar**2*k_y*gamma_2*k_y/(2*m_0)+sqrt(6)*I*hbar**2*k_y*gamma_3*k_x/(2*m_0)],[sqrt(6)*k_z*P/3,-sqrt(6)*k_x*P/6+sqrt(6)*I*k_y*P/6,sqrt(3)*hbar**2*k_x*gamma_3*k_z/(2*m_0)-sqrt(3)*hbar**2*k_x*k_z/(6*m_0)-sqrt(3)*hbar**2*k_x*kappa*k_z/(2*m_0)+sqrt(3)*I*hbar**2*k_y*gamma_3*k_z/(2*m_0)-sqrt(3)*I*hbar**2*k_y*k_z/(6*m_0)-sqrt(3)*I*hbar**2*k_y*kappa*k_z/(2*m_0)+sqrt(3)*hbar**2*k_z*gamma_3*k_x/(2*m_0)+sqrt(3)*I*hbar**2*k_z*gamma_3*k_y/(2*m_0)+sqrt(3)*hbar**2*k_z*k_x/(6*m_0)+sqrt(3)*I*hbar**2*k_z*k_y/(6*m_0)+sqrt(3)*hbar**2*k_z*kappa*k_x/(2*m_0)+sqrt(3)*I*hbar**2*k_z*kappa*k_y/(2*m_0),-hbar**2*k_x*gamma_1*k_x/(2*m_0)+hbar**2*k_x*gamma_2*k_x/(2*m_0)-I*hbar**2*k_x*k_y/(6*m_0)-I*hbar**2*k_x*kappa*k_y/(2*m_0)-hbar**2*k_y*gamma_1*k_y/(2*m_0)+hbar**2*k_y*gamma_2*k_y/(2*m_0)+I*hbar**2*k_y*k_x/(6*m_0)+I*hbar**2*k_y*kappa*k_x/(2*m_0)-hbar**2*k_z*gamma_1*k_z/(2*m_0)-hbar**2*k_z*gamma_2*k_z/m_0+E_v,hbar**2*k_x*k_z/(3*m_0)+hbar**2*k_x*kappa*k_z/m_0-I*hbar**2*k_y*k_z/(3*m_0)-I*hbar**2*k_y*kappa*k_z/m_0-hbar**2*k_z*k_x/(3*m_0)+I*hbar**2*k_z*k_y/(3*m_0)-hbar**2*k_z*kappa*k_x/m_0+I*hbar**2*k_z*kappa*k_y/m_0,sqrt(3)*hbar**2*k_x*gamma_2*k_x/(2*m_0)-sqrt(3)*I*hbar**2*k_x*gamma_3*k_y/(2*m_0)-sqrt(3)*hbar**2*k_y*gamma_2*k_y/(2*m_0)-sqrt(3)*I*hbar**2*k_y*gamma_3*k_x/(2*m_0),-sqrt(2)*hbar**2*k_x*gamma_2*k_x/(2*m_0)-sqrt(2)*I*hbar**2*k_x*k_y/(6*m_0)-sqrt(2)*I*hbar**2*k_x*kappa*k_y/(2*m_0)-sqrt(2)*hbar**2*k_y*gamma_2*k_y/(2*m_0)+sqrt(2)*I*hbar**2*k_y*k_x/(6*m_0)+sqrt(2)*I*hbar**2*k_y*kappa*k_x/(2*m_0)+sqrt(2)*hbar**2*k_z*gamma_2*k_z/m_0,3*sqrt(2)*hbar**2*k_x*gamma_3*k_z/(4*m_0)-sqrt(2)*hbar**2*k_x*k_z/(12*m_0)-sqrt(2)*hbar**2*k_x*kappa*k_z/(4*m_0)-3*sqrt(2)*I*hbar**2*k_y*gamma_3*k_z/(4*m_0)+sqrt(2)*I*hbar**2*k_y*k_z/(12*m_0)+sqrt(2)*I*hbar**2*k_y*kappa*k_z/(4*m_0)+3*sqrt(2)*hbar**2*k_z*gamma_3*k_x/(4*m_0)-3*sqrt(2)*I*hbar**2*k_z*gamma_3*k_y/(4*m_0)+sqrt(2)*hbar**2*k_z*k_x/(12*m_0)-sqrt(2)*I*hbar**2*k_z*k_y/(12*m_0)+sqrt(2)*hbar**2*k_z*kappa*k_x/(4*m_0)-sqrt(2)*I*hbar**2*k_z*kappa*k_y/(4*m_0)],[sqrt(6)*k_x*P/6+sqrt(6)*I*k_y*P/6,sqrt(6)*k_z*P/3,sqrt(3)*hbar**2*k_x*gamma_2*k_x/(2*m_0)+sqrt(3)*I*hbar**2*k_x*gamma_3*k_y/(2*m_0)-sqrt(3)*hbar**2*k_y*gamma_2*k_y/(2*m_0)+sqrt(3)*I*hbar**2*k_y*gamma_3*k_x/(2*m_0),-hbar**2*k_x*k_z/(3*m_0)-hbar**2*k_x*kappa*k_z/m_0-I*hbar**2*k_y*k_z/(3*m_0)-I*hbar**2*k_y*kappa*k_z/m_0+hbar**2*k_z*k_x/(3*m_0)+I*hbar**2*k_z*k_y/(3*m_0)+hbar**2*k_z*kappa*k_x/m_0+I*hbar**2*k_z*kappa*k_y/m_0,-hbar**2*k_x*gamma_1*k_x/(2*m_0)+hbar**2*k_x*gamma_2*k_x/(2*m_0)+I*hbar**2*k_x*k_y/(6*m_0)+I*hbar**2*k_x*kappa*k_y/(2*m_0)-hbar**2*k_y*gamma_1*k_y/(2*m_0)+hbar**2*k_y*gamma_2*k_y/(2*m_0)-I*hbar**2*k_y*k_x/(6*m_0)-I*hbar**2*k_y*kappa*k_x/(2*m_0)-hbar**2*k_z*gamma_1*k_z/(2*m_0)-hbar**2*k_z*gamma_2*k_z/m_0+E_v,-sqrt(3)*hbar**2*k_x*gamma_3*k_z/(2*m_0)+sqrt(3)*hbar**2*k_x*k_z/(6*m_0)+sqrt(3)*hbar**2*k_x*kappa*k_z/(2*m_0)+sqrt(3)*I*hbar**2*k_y*gamma_3*k_z/(2*m_0)-sqrt(3)*I*hbar**2*k_y*k_z/(6*m_0)-sqrt(3)*I*hbar**2*k_y*kappa*k_z/(2*m_0)-sqrt(3)*hbar**2*k_z*gamma_3*k_x/(2*m_0)+sqrt(3)*I*hbar**2*k_z*gamma_3*k_y/(2*m_0)-sqrt(3)*hbar**2*k_z*k_x/(6*m_0)+sqrt(3)*I*hbar**2*k_z*k_y/(6*m_0)-sqrt(3)*hbar**2*k_z*kappa*k_x/(2*m_0)+sqrt(3)*I*hbar**2*k_z*kappa*k_y/(2*m_0),3*sqrt(2)*hbar**2*k_x*gamma_3*k_z/(4*m_0)-sqrt(2)*hbar**2*k_x*k_z/(12*m_0)-sqrt(2)*hbar**2*k_x*kappa*k_z/(4*m_0)+3*sqrt(2)*I*hbar**2*k_y*gamma_3*k_z/(4*m_0)-sqrt(2)*I*hbar**2*k_y*k_z/(12*m_0)-sqrt(2)*I*hbar**2*k_y*kappa*k_z/(4*m_0)+3*sqrt(2)*hbar**2*k_z*gamma_3*k_x/(4*m_0)+3*sqrt(2)*I*hbar**2*k_z*gamma_3*k_y/(4*m_0)+sqrt(2)*hbar**2*k_z*k_x/(12*m_0)+sqrt(2)*I*hbar**2*k_z*k_y/(12*m_0)+sqrt(2)*hbar**2*k_z*kappa*k_x/(4*m_0)+sqrt(2)*I*hbar**2*k_z*kappa*k_y/(4*m_0),sqrt(2)*hbar**2*k_x*gamma_2*k_x/(2*m_0)-sqrt(2)*I*hbar**2*k_x*k_y/(6*m_0)-sqrt(2)*I*hbar**2*k_x*kappa*k_y/(2*m_0)+sqrt(2)*hbar**2*k_y*gamma_2*k_y/(2*m_0)+sqrt(2)*I*hbar**2*k_y*k_x/(6*m_0)+sqrt(2)*I*hbar**2*k_y*kappa*k_x/(2*m_0)-sqrt(2)*hbar**2*k_z*gamma_2*k_z/m_0],[0,sqrt(2)*k_x*P/2+sqrt(2)*I*k_y*P/2,0,sqrt(3)*hbar**2*k_x*gamma_2*k_x/(2*m_0)+sqrt(3)*I*hbar**2*k_x*gamma_3*k_y/(2*m_0)-sqrt(3)*hbar**2*k_y*gamma_2*k_y/(2*m_0)+sqrt(3)*I*hbar**2*k_y*gamma_3*k_x/(2*m_0),-sqrt(3)*hbar**2*k_x*gamma_3*k_z/(2*m_0)-sqrt(3)*hbar**2*k_x*k_z/(6*m_0)-sqrt(3)*hbar**2*k_x*kappa*k_z/(2*m_0)-sqrt(3)*I*hbar**2*k_y*gamma_3*k_z/(2*m_0)-sqrt(3)*I*hbar**2*k_y*k_z/(6*m_0)-sqrt(3)*I*hbar**2*k_y*kappa*k_z/(2*m_0)-sqrt(3)*hbar**2*k_z*gamma_3*k_x/(2*m_0)-sqrt(3)*I*hbar**2*k_z*gamma_3*k_y/(2*m_0)+sqrt(3)*hbar**2*k_z*k_x/(6*m_0)+sqrt(3)*I*hbar**2*k_z*k_y/(6*m_0)+sqrt(3)*hbar**2*k_z*kappa*k_x/(2*m_0)+sqrt(3)*I*hbar**2*k_z*kappa*k_y/(2*m_0),-hbar**2*k_x*gamma_1*k_x/(2*m_0)-hbar**2*k_x*gamma_2*k_x/(2*m_0)+I*hbar**2*k_x*k_y/(2*m_0)+3*I*hbar**2*k_x*kappa*k_y/(2*m_0)-hbar**2*k_y*gamma_1*k_y/(2*m_0)-hbar**2*k_y*gamma_2*k_y/(2*m_0)-I*hbar**2*k_y*k_x/(2*m_0)-3*I*hbar**2*k_y*kappa*k_x/(2*m_0)-hbar**2*k_z*gamma_1*k_z/(2*m_0)+hbar**2*k_z*gamma_2*k_z/m_0+E_v,sqrt(6)*hbar**2*k_x*gamma_2*k_x/(2*m_0)+sqrt(6)*I*hbar**2*k_x*gamma_3*k_y/(2*m_0)-sqrt(6)*hbar**2*k_y*gamma_2*k_y/(2*m_0)+sqrt(6)*I*hbar**2*k_y*gamma_3*k_x/(2*m_0),-sqrt(6)*hbar**2*k_x*gamma_3*k_z/(4*m_0)-sqrt(6)*hbar**2*k_x*k_z/(12*m_0)-sqrt(6)*hbar**2*k_x*kappa*k_z/(4*m_0)-sqrt(6)*I*hbar**2*k_y*gamma_3*k_z/(4*m_0)-sqrt(6)*I*hbar**2*k_y*k_z/(12*m_0)-sqrt(6)*I*hbar**2*k_y*kappa*k_z/(4*m_0)-sqrt(6)*hbar**2*k_z*gamma_3*k_x/(4*m_0)-sqrt(6)*I*hbar**2*k_z*gamma_3*k_y/(4*m_0)+sqrt(6)*hbar**2*k_z*k_x/(12*m_0)+sqrt(6)*I*hbar**2*k_z*k_y/(12*m_0)+sqrt(6)*hbar**2*k_z*kappa*k_x/(4*m_0)+sqrt(6)*I*hbar**2*k_z*kappa*k_y/(4*m_0)],[-sqrt(3)*k_z*P/3,-sqrt(3)*k_x*P/3+sqrt(3)*I*k_y*P/3,-sqrt(6)*hbar**2*k_x*gamma_3*k_z/(4*m_0)+sqrt(6)*hbar**2*k_x*k_z/(12*m_0)+sqrt(6)*hbar**2*k_x*kappa*k_z/(4*m_0)-sqrt(6)*I*hbar**2*k_y*gamma_3*k_z/(4*m_0)+sqrt(6)*I*hbar**2*k_y*k_z/(12*m_0)+sqrt(6)*I*hbar**2*k_y*kappa*k_z/(4*m_0)-sqrt(6)*hbar**2*k_z*gamma_3*k_x/(4*m_0)-sqrt(6)*I*hbar**2*k_z*gamma_3*k_y/(4*m_0)-sqrt(6)*hbar**2*k_z*k_x/(12*m_0)-sqrt(6)*I*hbar**2*k_z*k_y/(12*m_0)-sqrt(6)*hbar**2*k_z*kappa*k_x/(4*m_0)-sqrt(6)*I*hbar**2*k_z*kappa*k_y/(4*m_0),-sqrt(2)*hbar**2*k_x*gamma_2*k_x/(2*m_0)-sqrt(2)*I*hbar**2*k_x*k_y/(6*m_0)-sqrt(2)*I*hbar**2*k_x*kappa*k_y/(2*m_0)-sqrt(2)*hbar**2*k_y*gamma_2*k_y/(2*m_0)+sqrt(2)*I*hbar**2*k_y*k_x/(6*m_0)+sqrt(2)*I*hbar**2*k_y*kappa*k_x/(2*m_0)+sqrt(2)*hbar**2*k_z*gamma_2*k_z/m_0,3*sqrt(2)*hbar**2*k_x*gamma_3*k_z/(4*m_0)+sqrt(2)*hbar**2*k_x*k_z/(12*m_0)+sqrt(2)*hbar**2*k_x*kappa*k_z/(4*m_0)-3*sqrt(2)*I*hbar**2*k_y*gamma_3*k_z/(4*m_0)-sqrt(2)*I*hbar**2*k_y*k_z/(12*m_0)-sqrt(2)*I*hbar**2*k_y*kappa*k_z/(4*m_0)+3*sqrt(2)*hbar**2*k_z*gamma_3*k_x/(4*m_0)-3*sqrt(2)*I*hbar**2*k_z*gamma_3*k_y/(4*m_0)-sqrt(2)*hbar**2*k_z*k_x/(12*m_0)+sqrt(2)*I*hbar**2*k_z*k_y/(12*m_0)-sqrt(2)*hbar**2*k_z*kappa*k_x/(4*m_0)+sqrt(2)*I*hbar**2*k_z*kappa*k_y/(4*m_0),sqrt(6)*hbar**2*k_x*gamma_2*k_x/(2*m_0)-sqrt(6)*I*hbar**2*k_x*gamma_3*k_y/(2*m_0)-sqrt(6)*hbar**2*k_y*gamma_2*k_y/(2*m_0)-sqrt(6)*I*hbar**2*k_y*gamma_3*k_x/(2*m_0),-hbar**2*k_x*gamma_1*k_x/(2*m_0)-I*hbar**2*k_x*k_y/(3*m_0)-I*hbar**2*k_x*kappa*k_y/m_0-hbar**2*k_y*gamma_1*k_y/(2*m_0)+I*hbar**2*k_y*k_x/(3*m_0)+I*hbar**2*k_y*kappa*k_x/m_0-hbar**2*k_z*gamma_1*k_z/(2*m_0)-Delta+E_v,hbar**2*k_x*k_z/(3*m_0)+hbar**2*k_x*kappa*k_z/m_0-I*hbar**2*k_y*k_z/(3*m_0)-I*hbar**2*k_y*kappa*k_z/m_0-hbar**2*k_z*k_x/(3*m_0)+I*hbar**2*k_z*k_y/(3*m_0)-hbar**2*k_z*kappa*k_x/m_0+I*hbar**2*k_z*kappa*k_y/m_0],[-sqrt(3)*k_x*P/3-sqrt(3)*I*k_y*P/3,sqrt(3)*k_z*P/3,-sqrt(6)*hbar**2*k_x*gamma_2*k_x/(2*m_0)-sqrt(6)*I*hbar**2*k_x*gamma_3*k_y/(2*m_0)+sqrt(6)*hbar**2*k_y*gamma_2*k_y/(2*m_0)-sqrt(6)*I*hbar**2*k_y*gamma_3*k_x/(2*m_0),3*sqrt(2)*hbar**2*k_x*gamma_3*k_z/(4*m_0)+sqrt(2)*hbar**2*k_x*k_z/(12*m_0)+sqrt(2)*hbar**2*k_x*kappa*k_z/(4*m_0)+3*sqrt(2)*I*hbar**2*k_y*gamma_3*k_z/(4*m_0)+sqrt(2)*I*hbar**2*k_y*k_z/(12*m_0)+sqrt(2)*I*hbar**2*k_y*kappa*k_z/(4*m_0)+3*sqrt(2)*hbar**2*k_z*gamma_3*k_x/(4*m_0)+3*sqrt(2)*I*hbar**2*k_z*gamma_3*k_y/(4*m_0)-sqrt(2)*hbar**2*k_z*k_x/(12*m_0)-sqrt(2)*I*hbar**2*k_z*k_y/(12*m_0)-sqrt(2)*hbar**2*k_z*kappa*k_x/(4*m_0)-sqrt(2)*I*hbar**2*k_z*kappa*k_y/(4*m_0),sqrt(2)*hbar**2*k_x*gamma_2*k_x/(2*m_0)-sqrt(2)*I*hbar**2*k_x*k_y/(6*m_0)-sqrt(2)*I*hbar**2*k_x*kappa*k_y/(2*m_0)+sqrt(2)*hbar**2*k_y*gamma_2*k_y/(2*m_0)+sqrt(2)*I*hbar**2*k_y*k_x/(6*m_0)+sqrt(2)*I*hbar**2*k_y*kappa*k_x/(2*m_0)-sqrt(2)*hbar**2*k_z*gamma_2*k_z/m_0,-sqrt(6)*hbar**2*k_x*gamma_3*k_z/(4*m_0)+sqrt(6)*hbar**2*k_x*k_z/(12*m_0)+sqrt(6)*hbar**2*k_x*kappa*k_z/(4*m_0)+sqrt(6)*I*hbar**2*k_y*gamma_3*k_z/(4*m_0)-sqrt(6)*I*hbar**2*k_y*k_z/(12*m_0)-sqrt(6)*I*hbar**2*k_y*kappa*k_z/(4*m_0)-sqrt(6)*hbar**2*k_z*gamma_3*k_x/(4*m_0)+sqrt(6)*I*hbar**2*k_z*gamma_3*k_y/(4*m_0)-sqrt(6)*hbar**2*k_z*k_x/(12*m_0)+sqrt(6)*I*hbar**2*k_z*k_y/(12*m_0)-sqrt(6)*hbar**2*k_z*kappa*k_x/(4*m_0)+sqrt(6)*I*hbar**2*k_z*kappa*k_y/(4*m_0),-hbar**2*k_x*k_z/(3*m_0)-hbar**2*k_x*kappa*k_z/m_0-I*hbar**2*k_y*k_z/(3*m_0)-I*hbar**2*k_y*kappa*k_z/m_0+hbar**2*k_z*k_x/(3*m_0)+I*hbar**2*k_z*k_y/(3*m_0)+hbar**2*k_z*kappa*k_x/m_0+I*hbar**2*k_z*kappa*k_y/m_0,-hbar**2*k_x*gamma_1*k_x/(2*m_0)+I*hbar**2*k_x*k_y/(3*m_0)+I*hbar**2*k_x*kappa*k_y/m_0-hbar**2*k_y*gamma_1*k_y/(2*m_0)-I*hbar**2*k_y*k_x/(3*m_0)-I*hbar**2*k_y*kappa*k_x/m_0-hbar**2*k_z*gamma_1*k_z/(2*m_0)-Delta+E_v]])'
Hkp = qsymm.Model(ham_kp)

The symmetry is the full cubic group with time reversal; 96 elements:

sg, cg = qsymm.symmetries(Hkp, cubic_group)
print(len(sg))

96

If we restrict the parameters we can change the symmetry group:

kp_examples = [
    {'hbar': 1},
    {'hbar': 1, 'P': 0, 'Delta': 0},
    {'hbar': 1, 'P': 0, 'Delta': 0, 'gamma_3': 0},
    {'hbar': 1, 'gamma_2': 'gamma_1', 'gamma_3': 'gamma_1', 'k_z': 0},
    {'hbar': 1, 'P': 0, 'Delta': 0, 'gamma_3': 0},
    {'hbar': 1, 'Delta': 0, 'k_z': 0},
]

for subs in kp_examples:
    Hkp = qsymm.Model(ham_kp, locals=subs)
    sg, cg = qsymm.symmetries(Hkp, cubic_group)
    print(subs, len(sg), len(cg))

{'hbar': 1} 96 0

{'hbar': 1, 'P': 0, 'Delta': 0} 96 7

{'hbar': 1, 'P': 0, 'Delta': 0, 'gamma_3': 0} 96 15

{'hbar': 1, 'gamma_2': 'gamma_1', 'gamma_3': 'gamma_1', 'k_z': 0} 32 1

{'hbar': 1, 'P': 0, 'Delta': 0, 'gamma_3': 0} 96 15

{'hbar': 1, 'Delta': 0, 'k_z': 0} 32 7

Generating tight-binding models

See also

The complete source code of this example can be found in bloch_generator.py. A Jupyter notebook can be found in bloch_generator.ipynb.

In addition to finding the symmetry group of a given Hamiltonian, Qsymm can also generate a class of models that have a given symmetry. In this tutorial we will work through a few pertinent examples to show how we can use Qsymm to generate a tight-binding model from symmetry constraints.

Graphene

First we’re going to generate the spinless nearest-neighbor tight-binding Hamiltonian for graphene.

The generators of the symmetry group are time-reversal symmetry, sublattice (or chiral) symmetry, and threefold rotation symmetry.

# Time reversal
TR = qsymm.time_reversal(2, U=np.eye(2))

# Chiral symmetry
C = qsymm.chiral(2, U=np.array([[1, 0], [0, -1]]))

# Atom A rotates into A, B into B.
# We use sympy to get an exact representation of the multiplying factors
sphi = (2 / 3) * sympy.pi

C3 = qsymm.PointGroupElement(
    sympy.Matrix([
        [sympy.cos(sphi), -sympy.sin(sphi)],
        [sympy.sin(sphi), sympy.cos(sphi)]
    ]),
    U=np.eye(2),
)

symmetries = [C, TR, C3]

There are two carbon atoms per unit cell (A and B) with one orbital each. The lattice is triangular, and only includes hoppings between nearest neighbour atoms. This restricts hoppings to only those between atoms of different types, such that each atom couples to three neighbouring atoms.

Using the symmetrization strategy to generate the Hamiltonian, it is sufficient to specify hoppings to one such neighbour along with the symmetry generators, and we take the vector \((1,0)\) to connect this neighbouring pair of atoms.

norbs = [('A', 1), ('B', 1)]  # A and B atom per unit cell, one orbital each
hopping_vectors = [('A', 'B', [0, 1])]  # Hopping between neighbouring A and B atoms

Now we generate the hamiltonian using bloch_family:

family = qsymm.bloch_family(hopping_vectors, symmetries, norbs)
qsymm.display_family(family)

\[\begin{split}\displaystyle \left[\begin{matrix}0 & e^{i k_{y}} + e^{- \frac{i k_{y} \sin^{2}{\left(0.666666666666667 \pi \right)}}{\cos^{4}{\left(0.666666666666667 \pi \right)} + 2 \sin^{2}{\left(0.666666666666667 \pi \right)} \cos^{2}{\left(0.666666666666667 \pi \right)} + \sin^{4}{\left(0.666666666666667 \pi \right)}}} e^{\frac{i k_{y} \cos^{2}{\left(0.666666666666667 \pi \right)}}{\cos^{4}{\left(0.666666666666667 \pi \right)} + 2 \sin^{2}{\left(0.666666666666667 \pi \right)} \cos^{2}{\left(0.666666666666667 \pi \right)} + \sin^{4}{\left(0.666666666666667 \pi \right)}}} e^{- \frac{2 i k_{x} \sin{\left(0.666666666666667 \pi \right)} \cos{\left(0.666666666666667 \pi \right)}}{\cos^{4}{\left(0.666666666666667 \pi \right)} + 2 \sin^{2}{\left(0.666666666666667 \pi \right)} \cos^{2}{\left(0.666666666666667 \pi \right)} + \sin^{4}{\left(0.666666666666667 \pi \right)}}} + e^{- \frac{i k_{x} \sin{\left(0.666666666666667 \pi \right)}}{\cos^{2}{\left(0.666666666666667 \pi \right)} + \sin^{2}{\left(0.666666666666667 \pi \right)}}} e^{\frac{i k_{y} \cos{\left(0.666666666666667 \pi \right)}}{\cos^{2}{\left(0.666666666666667 \pi \right)} + \sin^{2}{\left(0.666666666666667 \pi \right)}}}\\e^{\frac{i k_{x} \sin{\left(0.666666666666667 \pi \right)}}{\cos^{2}{\left(0.666666666666667 \pi \right)} + \sin^{2}{\left(0.666666666666667 \pi \right)}}} e^{- \frac{i k_{y} \cos{\left(0.666666666666667 \pi \right)}}{\cos^{2}{\left(0.666666666666667 \pi \right)} + \sin^{2}{\left(0.666666666666667 \pi \right)}}} + e^{\frac{i k_{y} \sin^{2}{\left(0.666666666666667 \pi \right)}}{\cos^{4}{\left(0.666666666666667 \pi \right)} + 2 \sin^{2}{\left(0.666666666666667 \pi \right)} \cos^{2}{\left(0.666666666666667 \pi \right)} + \sin^{4}{\left(0.666666666666667 \pi \right)}}} e^{- \frac{i k_{y} \cos^{2}{\left(0.666666666666667 \pi \right)}}{\cos^{4}{\left(0.666666666666667 \pi \right)} + 2 \sin^{2}{\left(0.666666666666667 \pi \right)} \cos^{2}{\left(0.666666666666667 \pi \right)} + \sin^{4}{\left(0.666666666666667 \pi \right)}}} e^{\frac{2 i k_{x} \sin{\left(0.666666666666667 \pi \right)} \cos{\left(0.666666666666667 \pi \right)}}{\cos^{4}{\left(0.666666666666667 \pi \right)} + 2 \sin^{2}{\left(0.666666666666667 \pi \right)} \cos^{2}{\left(0.666666666666667 \pi \right)} + \sin^{4}{\left(0.666666666666667 \pi \right)}}} + e^{- i k_{y}} & 0\end{matrix}\right]\end{split}\]

Scale the bond length in terms of the graphene lattice constant, and have the function return a list of BlochModel objects. For this we can use a more convenient definition of the rotation

C3 = qsymm.rotation(1/3, U=np.eye(2))
symmetries = [C, TR, C3]

norbs = [('A', 1), ('B', 1)]
hopping_vectors = [('A', 'B', [0, 1/np.sqrt(3)])]

family = qsymm.bloch_family(hopping_vectors, symmetries, norbs, bloch_model=True)
qsymm.display_family(family)

\[\begin{split}\displaystyle \left[\begin{matrix}0 & e^{\frac{i k_{x}}{2}} e^{- \frac{\sqrt{3} i k_{y}}{6}} + e^{\frac{\sqrt{3} i k_{y}}{3}} + e^{- \frac{i k_{x}}{2}} e^{- \frac{\sqrt{3} i k_{y}}{6}}\\e^{\frac{i k_{x}}{2}} e^{\frac{\sqrt{3} i k_{y}}{6}} + e^{- \frac{\sqrt{3} i k_{y}}{3}} + e^{- \frac{i k_{x}}{2}} e^{\frac{\sqrt{3} i k_{y}}{6}} & 0\end{matrix}\right]\end{split}\]

Three-orbital tight-binding model for monolayer \(MX_2\)

We use the Hamiltonian generator to reproduce the tight binding model for monolayer \(MX_2\) published in Phys. Rev. B 88, 085433 (2013).

The generators of the symmetry group of the tight binding model are time reversal symmetry, mirror symmetry and threefold rotation symmetry.

# Time reversal
TR = qsymm.time_reversal(2, np.eye(3))

# Mirror symmetry
Mx = qsymm.mirror([1, 0], np.diag([1, -1, 1]))

# Threefold rotation on d_z^2, d_xy, d_x^2-y^2 states.
C3U = np.array([
    [1, 0, 0],
    [0, -0.5, -np.sqrt(3)/2],
    [0, np.sqrt(3)/2, -0.5]
])

# Could also use the predefined representation of rotations on d-orbitals
Ld = qsymm.groups.L_matrices(3, 2)
C3U2 = qsymm.groups.spin_rotation(2 * np.pi * np.array([0, 0, 1/3]), Ld)

# Restrict to d_z^2, d_xy, d_x^2-y^2 states
mask = np.array([1, 2 ,0])
C3U2 = C3U2[mask][:, mask]

assert np.allclose(C3U, C3U2)

C3 = qsymm.rotation(1/3, U=C3U)

symmetries = [TR, Mx, C3]

Next, we specify the hoppings to include. The tight binding model has a triangular lattice, three orbitals per M atom, and nearest neighbour hopping.

# One site per unit cell (M atom), with three orbitals
norbs = [('a', 3)]

Each atom has six nearest neighbour atoms at a distance of one primitive lattice vector. Since we use the symmetrization strategy to generate the Hamiltonian, it is sufficient to specify a hopping to one nearest neighbour atom along with the symmetry generators. We take the primitive vector connecting the pair of atoms to be \((1,0)\).

# Hopping to a neighbouring atom one primitive lattice vector away
hopping_vectors = [('a', 'a', [1, 0])]

We again use bloch_family to generate the tight-binding Hamiltonian:

family = qsymm.bloch_family(hopping_vectors, symmetries, norbs, bloch_model=True)

The Hamiltonian family should include 8 linearly independent components, including the onsite terms.

len(family)

8

qsymm.display_family(family)

\[\begin{split}\displaystyle \left[\begin{matrix}1 & 0 & 0\\0 & 0 & 0\\0 & 0 & 0\end{matrix}\right]\end{split}\]

\[\begin{split}\displaystyle \left[\begin{matrix}0 & 0 & 0\\0 & 1 & 0\\0 & 0 & 1\end{matrix}\right]\end{split}\]

\[\begin{split}\displaystyle \left[\begin{matrix}e^{\frac{i k_{x}}{2}} e^{\frac{\sqrt{3} i k_{y}}{2}} + e^{\frac{i k_{x}}{2}} e^{- \frac{\sqrt{3} i k_{y}}{2}} + e^{i k_{x}} + e^{- i k_{x}} + e^{- \frac{i k_{x}}{2}} e^{\frac{\sqrt{3} i k_{y}}{2}} + e^{- \frac{i k_{x}}{2}} e^{- \frac{\sqrt{3} i k_{y}}{2}} & 0 & 0\\0 & 0 & 0\\0 & 0 & 0\end{matrix}\right]\end{split}\]

\[\begin{split}\displaystyle \left[\begin{matrix}0 & \frac{e^{\frac{i k_{x}}{2}} e^{\frac{\sqrt{3} i k_{y}}{2}}}{2} + \frac{e^{\frac{i k_{x}}{2}} e^{- \frac{\sqrt{3} i k_{y}}{2}}}{2} + e^{i k_{x}} - e^{- i k_{x}} - \frac{e^{- \frac{i k_{x}}{2}} e^{\frac{\sqrt{3} i k_{y}}{2}}}{2} - \frac{e^{- \frac{i k_{x}}{2}} e^{- \frac{\sqrt{3} i k_{y}}{2}}}{2} & \frac{\sqrt{3} e^{\frac{i k_{x}}{2}} e^{\frac{\sqrt{3} i k_{y}}{2}}}{2} - \frac{\sqrt{3} e^{\frac{i k_{x}}{2}} e^{- \frac{\sqrt{3} i k_{y}}{2}}}{2} + \frac{\sqrt{3} e^{- \frac{i k_{x}}{2}} e^{\frac{\sqrt{3} i k_{y}}{2}}}{2} - \frac{\sqrt{3} e^{- \frac{i k_{x}}{2}} e^{- \frac{\sqrt{3} i k_{y}}{2}}}{2}\\- \frac{e^{\frac{i k_{x}}{2}} e^{\frac{\sqrt{3} i k_{y}}{2}}}{2} - \frac{e^{\frac{i k_{x}}{2}} e^{- \frac{\sqrt{3} i k_{y}}{2}}}{2} - e^{i k_{x}} + e^{- i k_{x}} + \frac{e^{- \frac{i k_{x}}{2}} e^{\frac{\sqrt{3} i k_{y}}{2}}}{2} + \frac{e^{- \frac{i k_{x}}{2}} e^{- \frac{\sqrt{3} i k_{y}}{2}}}{2} & 0 & 0\\- \frac{\sqrt{3} e^{\frac{i k_{x}}{2}} e^{\frac{\sqrt{3} i k_{y}}{2}}}{2} + \frac{\sqrt{3} e^{\frac{i k_{x}}{2}} e^{- \frac{\sqrt{3} i k_{y}}{2}}}{2} - \frac{\sqrt{3} e^{- \frac{i k_{x}}{2}} e^{\frac{\sqrt{3} i k_{y}}{2}}}{2} + \frac{\sqrt{3} e^{- \frac{i k_{x}}{2}} e^{- \frac{\sqrt{3} i k_{y}}{2}}}{2} & 0 & 0\end{matrix}\right]\end{split}\]

\[\begin{split}\displaystyle \left[\begin{matrix}0 & \frac{\sqrt{3} e^{\frac{i k_{x}}{2}} e^{\frac{\sqrt{3} i k_{y}}{2}}}{2} - \frac{\sqrt{3} e^{\frac{i k_{x}}{2}} e^{- \frac{\sqrt{3} i k_{y}}{2}}}{2} - \frac{\sqrt{3} e^{- \frac{i k_{x}}{2}} e^{\frac{\sqrt{3} i k_{y}}{2}}}{2} + \frac{\sqrt{3} e^{- \frac{i k_{x}}{2}} e^{- \frac{\sqrt{3} i k_{y}}{2}}}{2} & - \frac{e^{\frac{i k_{x}}{2}} e^{\frac{\sqrt{3} i k_{y}}{2}}}{2} - \frac{e^{\frac{i k_{x}}{2}} e^{- \frac{\sqrt{3} i k_{y}}{2}}}{2} + e^{i k_{x}} + e^{- i k_{x}} - \frac{e^{- \frac{i k_{x}}{2}} e^{\frac{\sqrt{3} i k_{y}}{2}}}{2} - \frac{e^{- \frac{i k_{x}}{2}} e^{- \frac{\sqrt{3} i k_{y}}{2}}}{2}\\\frac{\sqrt{3} e^{\frac{i k_{x}}{2}} e^{\frac{\sqrt{3} i k_{y}}{2}}}{2} - \frac{\sqrt{3} e^{\frac{i k_{x}}{2}} e^{- \frac{\sqrt{3} i k_{y}}{2}}}{2} - \frac{\sqrt{3} e^{- \frac{i k_{x}}{2}} e^{\frac{\sqrt{3} i k_{y}}{2}}}{2} + \frac{\sqrt{3} e^{- \frac{i k_{x}}{2}} e^{- \frac{\sqrt{3} i k_{y}}{2}}}{2} & 0 & 0\\- \frac{e^{\frac{i k_{x}}{2}} e^{\frac{\sqrt{3} i k_{y}}{2}}}{2} - \frac{e^{\frac{i k_{x}}{2}} e^{- \frac{\sqrt{3} i k_{y}}{2}}}{2} + e^{i k_{x}} + e^{- i k_{x}} - \frac{e^{- \frac{i k_{x}}{2}} e^{\frac{\sqrt{3} i k_{y}}{2}}}{2} - \frac{e^{- \frac{i k_{x}}{2}} e^{- \frac{\sqrt{3} i k_{y}}{2}}}{2} & 0 & 0\end{matrix}\right]\end{split}\]

\[\begin{split}\displaystyle \left[\begin{matrix}0 & 0 & 0\\0 & \frac{e^{\frac{i k_{x}}{2}} e^{\frac{\sqrt{3} i k_{y}}{2}}}{4} + \frac{e^{\frac{i k_{x}}{2}} e^{- \frac{\sqrt{3} i k_{y}}{2}}}{4} + e^{i k_{x}} + e^{- i k_{x}} + \frac{e^{- \frac{i k_{x}}{2}} e^{\frac{\sqrt{3} i k_{y}}{2}}}{4} + \frac{e^{- \frac{i k_{x}}{2}} e^{- \frac{\sqrt{3} i k_{y}}{2}}}{4} & \frac{\sqrt{3} e^{\frac{i k_{x}}{2}} e^{\frac{\sqrt{3} i k_{y}}{2}}}{4} - \frac{\sqrt{3} e^{\frac{i k_{x}}{2}} e^{- \frac{\sqrt{3} i k_{y}}{2}}}{4} - \frac{\sqrt{3} e^{- \frac{i k_{x}}{2}} e^{\frac{\sqrt{3} i k_{y}}{2}}}{4} + \frac{\sqrt{3} e^{- \frac{i k_{x}}{2}} e^{- \frac{\sqrt{3} i k_{y}}{2}}}{4}\\0 & \frac{\sqrt{3} e^{\frac{i k_{x}}{2}} e^{\frac{\sqrt{3} i k_{y}}{2}}}{4} - \frac{\sqrt{3} e^{\frac{i k_{x}}{2}} e^{- \frac{\sqrt{3} i k_{y}}{2}}}{4} - \frac{\sqrt{3} e^{- \frac{i k_{x}}{2}} e^{\frac{\sqrt{3} i k_{y}}{2}}}{4} + \frac{\sqrt{3} e^{- \frac{i k_{x}}{2}} e^{- \frac{\sqrt{3} i k_{y}}{2}}}{4} & \frac{3 e^{\frac{i k_{x}}{2}} e^{\frac{\sqrt{3} i k_{y}}{2}}}{4} + \frac{3 e^{\frac{i k_{x}}{2}} e^{- \frac{\sqrt{3} i k_{y}}{2}}}{4} + \frac{3 e^{- \frac{i k_{x}}{2}} e^{\frac{\sqrt{3} i k_{y}}{2}}}{4} + \frac{3 e^{- \frac{i k_{x}}{2}} e^{- \frac{\sqrt{3} i k_{y}}{2}}}{4}\end{matrix}\right]\end{split}\]

\[\begin{split}\displaystyle \left[\begin{matrix}0 & 0 & 0\\0 & 0 & - e^{\frac{i k_{x}}{2}} e^{\frac{\sqrt{3} i k_{y}}{2}} - e^{\frac{i k_{x}}{2}} e^{- \frac{\sqrt{3} i k_{y}}{2}} + e^{i k_{x}} - e^{- i k_{x}} + e^{- \frac{i k_{x}}{2}} e^{\frac{\sqrt{3} i k_{y}}{2}} + e^{- \frac{i k_{x}}{2}} e^{- \frac{\sqrt{3} i k_{y}}{2}}\\0 & e^{\frac{i k_{x}}{2}} e^{\frac{\sqrt{3} i k_{y}}{2}} + e^{\frac{i k_{x}}{2}} e^{- \frac{\sqrt{3} i k_{y}}{2}} - e^{i k_{x}} + e^{- i k_{x}} - e^{- \frac{i k_{x}}{2}} e^{\frac{\sqrt{3} i k_{y}}{2}} - e^{- \frac{i k_{x}}{2}} e^{- \frac{\sqrt{3} i k_{y}}{2}} & 0\end{matrix}\right]\end{split}\]

\[\begin{split}\displaystyle \left[\begin{matrix}0 & 0 & 0\\0 & \frac{3 e^{\frac{i k_{x}}{2}} e^{\frac{\sqrt{3} i k_{y}}{2}}}{4} + \frac{3 e^{\frac{i k_{x}}{2}} e^{- \frac{\sqrt{3} i k_{y}}{2}}}{4} + \frac{3 e^{- \frac{i k_{x}}{2}} e^{\frac{\sqrt{3} i k_{y}}{2}}}{4} + \frac{3 e^{- \frac{i k_{x}}{2}} e^{- \frac{\sqrt{3} i k_{y}}{2}}}{4} & - \frac{\sqrt{3} e^{\frac{i k_{x}}{2}} e^{\frac{\sqrt{3} i k_{y}}{2}}}{4} + \frac{\sqrt{3} e^{\frac{i k_{x}}{2}} e^{- \frac{\sqrt{3} i k_{y}}{2}}}{4} + \frac{\sqrt{3} e^{- \frac{i k_{x}}{2}} e^{\frac{\sqrt{3} i k_{y}}{2}}}{4} - \frac{\sqrt{3} e^{- \frac{i k_{x}}{2}} e^{- \frac{\sqrt{3} i k_{y}}{2}}}{4}\\0 & - \frac{\sqrt{3} e^{\frac{i k_{x}}{2}} e^{\frac{\sqrt{3} i k_{y}}{2}}}{4} + \frac{\sqrt{3} e^{\frac{i k_{x}}{2}} e^{- \frac{\sqrt{3} i k_{y}}{2}}}{4} + \frac{\sqrt{3} e^{- \frac{i k_{x}}{2}} e^{\frac{\sqrt{3} i k_{y}}{2}}}{4} - \frac{\sqrt{3} e^{- \frac{i k_{x}}{2}} e^{- \frac{\sqrt{3} i k_{y}}{2}}}{4} & \frac{e^{\frac{i k_{x}}{2}} e^{\frac{\sqrt{3} i k_{y}}{2}}}{4} + \frac{e^{\frac{i k_{x}}{2}} e^{- \frac{\sqrt{3} i k_{y}}{2}}}{4} + e^{i k_{x}} + e^{- i k_{x}} + \frac{e^{- \frac{i k_{x}}{2}} e^{\frac{\sqrt{3} i k_{y}}{2}}}{4} + \frac{e^{- \frac{i k_{x}}{2}} e^{- \frac{\sqrt{3} i k_{y}}{2}}}{4}\end{matrix}\right]\end{split}\]

4-site model for monolayer \(WTe_2\)

We use the Hamiltonian generator to reproduce the tight binding model for monolayer WTe2 published in Phys. Rev. X 6, 041069 (2016).

The generators of the symmetry group of the tight binding model are time reversal symmetry, glide reflection and inversion symmetry.

# Define 4 sites with one orbital each
sites = ['Ad', 'Ap', 'Bd', 'Bp']
norbs = [(site, 1) for site in sites]

# Define symbolic coordinates for orbitals
rAp = qsymm.sympify('[x_Ap, y_Ap]')
rAd = qsymm.sympify('[x_Ad, y_Ad]')
rBp = qsymm.sympify('[x_Bp, y_Bp]')
rBd = qsymm.sympify('[x_Bd, y_Bd]')

# Define hoppings to include
hopping_vectors = [
    ('Bd', 'Bd', np.array([1, 0])),
    ('Ap', 'Ap', np.array([1, 0])),
    ('Bd', 'Ap', rAp - rBd),
    ('Ap', 'Bp', rBp - rAp),
    ('Ad', 'Bd', rBd - rAd),
]

# Inversion
perm_inv = {'Ad': 'Bd', 'Ap': 'Bp', 'Bd': 'Ad', 'Bp': 'Ap'}
onsite_inv = {site: (1 if site in ['Ad', 'Bd'] else -1) * np.eye(1) for site in sites}
inversion = qsymm.groups.symmetry_from_permutation(-np.eye(2), perm_inv, norbs, onsite_inv)

# Glide
perm_glide = {site: site for site in sites}
onsite_glide = {site: (1 if site in ['Ad', 'Bd'] else -1) * np.eye(1) for site in sites}
glide = qsymm.groups.symmetry_from_permutation(np.array([[-1, 0],[0, 1]]), perm_glide, norbs, onsite_glide)

# TR
time_reversal = qsymm.time_reversal(2, np.eye(4))

generators = {glide, inversion, time_reversal}
sg = qsymm.groups.generate_group(generators)

Again we generate the tight-binding Hamiltonian:

family = qsymm.bloch_family(hopping_vectors, generators, norbs=norbs)
qsymm.display_family(family)

\[\begin{split}\displaystyle \left[\begin{matrix}0 & 0 & 0 & e^{i k_{x} x_{Ap}} e^{- i k_{x} x_{Bd}} e^{- i k_{y} y_{Ap}} e^{i k_{y} y_{Bd}} - e^{- i k_{x} x_{Ap}} e^{i k_{x} x_{Bd}} e^{- i k_{y} y_{Ap}} e^{i k_{y} y_{Bd}}\\0 & 0 & - e^{i k_{x} x_{Ap}} e^{- i k_{x} x_{Bd}} e^{- i k_{y} y_{Ap}} e^{i k_{y} y_{Bd}} + e^{- i k_{x} x_{Ap}} e^{i k_{x} x_{Bd}} e^{- i k_{y} y_{Ap}} e^{i k_{y} y_{Bd}} & 0\\0 & e^{i k_{x} x_{Ap}} e^{- i k_{x} x_{Bd}} e^{i k_{y} y_{Ap}} e^{- i k_{y} y_{Bd}} - e^{- i k_{x} x_{Ap}} e^{i k_{x} x_{Bd}} e^{i k_{y} y_{Ap}} e^{- i k_{y} y_{Bd}} & 0 & 0\\- e^{i k_{x} x_{Ap}} e^{- i k_{x} x_{Bd}} e^{i k_{y} y_{Ap}} e^{- i k_{y} y_{Bd}} + e^{- i k_{x} x_{Ap}} e^{i k_{x} x_{Bd}} e^{i k_{y} y_{Ap}} e^{- i k_{y} y_{Bd}} & 0 & 0 & 0\end{matrix}\right]\end{split}\]

\[\begin{split}\displaystyle \left[\begin{matrix}1 & 0 & 0 & 0\\0 & 0 & 0 & 0\\0 & 0 & 1 & 0\\0 & 0 & 0 & 0\end{matrix}\right]\end{split}\]

\[\begin{split}\displaystyle \left[\begin{matrix}0 & 0 & 0 & 0\\0 & 1 & 0 & 0\\0 & 0 & 0 & 0\\0 & 0 & 0 & 1\end{matrix}\right]\end{split}\]

\[\begin{split}\displaystyle \left[\begin{matrix}0 & 0 & 0 & 0\\0 & 0 & 0 & e^{i k_{x} x_{Ap}} e^{- i k_{x} x_{Bp}} e^{- i k_{y} y_{Ap}} e^{i k_{y} y_{Bp}} + e^{- i k_{x} x_{Ap}} e^{i k_{x} x_{Bp}} e^{- i k_{y} y_{Ap}} e^{i k_{y} y_{Bp}}\\0 & 0 & 0 & 0\\0 & e^{i k_{x} x_{Ap}} e^{- i k_{x} x_{Bp}} e^{i k_{y} y_{Ap}} e^{- i k_{y} y_{Bp}} + e^{- i k_{x} x_{Ap}} e^{i k_{x} x_{Bp}} e^{i k_{y} y_{Ap}} e^{- i k_{y} y_{Bp}} & 0 & 0\end{matrix}\right]\end{split}\]

\[\begin{split}\displaystyle \left[\begin{matrix}e^{i k_{x}} + e^{- i k_{x}} & 0 & 0 & 0\\0 & 0 & 0 & 0\\0 & 0 & e^{i k_{x}} + e^{- i k_{x}} & 0\\0 & 0 & 0 & 0\end{matrix}\right]\end{split}\]

\[\begin{split}\displaystyle \left[\begin{matrix}0 & 0 & 0 & 0\\0 & e^{i k_{x}} + e^{- i k_{x}} & 0 & 0\\0 & 0 & 0 & 0\\0 & 0 & 0 & e^{i k_{x}} + e^{- i k_{x}}\end{matrix}\right]\end{split}\]

\[\begin{split}\displaystyle \left[\begin{matrix}0 & 0 & e^{i k_{x} x_{Ad}} e^{- i k_{x} x_{Bd}} e^{- i k_{y} y_{Ad}} e^{i k_{y} y_{Bd}} + e^{- i k_{x} x_{Ad}} e^{i k_{x} x_{Bd}} e^{- i k_{y} y_{Ad}} e^{i k_{y} y_{Bd}} & 0\\0 & 0 & 0 & 0\\e^{i k_{x} x_{Ad}} e^{- i k_{x} x_{Bd}} e^{i k_{y} y_{Ad}} e^{- i k_{y} y_{Bd}} + e^{- i k_{x} x_{Ad}} e^{i k_{x} x_{Bd}} e^{i k_{y} y_{Ad}} e^{- i k_{y} y_{Bd}} & 0 & 0 & 0\\0 & 0 & 0 & 0\end{matrix}\right]\end{split}\]

Square lattice with 4 sites in the unit cell

Now we’re going to make a model with square lattice that has 4 sites in the unit cell related by 4-fold rotation. Sites have spin-1/2 and we add time reversal and particle-hole symmetry.

hopping_vectors = [
    ('a', 'b', np.array([1, 0])),
    ('b', 'a', np.array([1, 0])),
    ('c', 'd', np.array([1, 0])),
    ('d', 'c', np.array([1, 0])),
    ('a', 'c', np.array([0, 1])),
    ('c', 'a', np.array([0, 1])),
    ('b', 'd', np.array([0, 1])),
    ('d', 'b', np.array([0, 1])),
]

# Define spin-1/2 operators
S = qsymm.groups.spin_matrices(1/2)
# Define real space rotation generators in 2D
L = qsymm.groups.L_matrices(d=2)

sites = ['a', 'b', 'c', 'd']
norbs = [(site, 2) for site in sites]

perm_C4 = {'a': 'b', 'b': 'd', 'd': 'c', 'c': 'a'}
onsite_C4 = {site: qsymm.groups.spin_rotation(2*np.pi * np.array([0, 0, 1/4]), S) for site in sites}
C4 = qsymm.groups.symmetry_from_permutation(
    qsymm.groups.spin_rotation(2*np.pi * np.array([1/4]), L, roundint=True),
    perm_C4,
    norbs,
    onsite_C4,
)

# Fermionic time-reversal
time_reversal = qsymm.time_reversal(
    2,
    np.kron(np.eye(4),
    qsymm.groups.spin_rotation(2*np.pi * np.array([0, 1/2, 0]), S)),
)

# define strange PH symmetry
particle_hole = qsymm.particle_hole(2, np.eye(8))

generators = {C4, time_reversal, particle_hole}
sg = qsymm.groups.generate_group(generators)

Once again we generate the tight-binding Hamiltonian:

family = qsymm.bloch_family(hopping_vectors, generators, norbs=norbs)
qsymm.display_family(family)

\[\begin{split}\displaystyle \left[\begin{matrix}0 & 0 & - i e^{- i k_{x}} & 0 & i e^{i k_{y}} & 0 & 0 & 0\\0 & 0 & 0 & i e^{- i k_{x}} & 0 & - i e^{i k_{y}} & 0 & 0\\i e^{i k_{x}} & 0 & 0 & 0 & 0 & 0 & - i e^{i k_{y}} & 0\\0 & - i e^{i k_{x}} & 0 & 0 & 0 & 0 & 0 & i e^{i k_{y}}\\- i e^{- i k_{y}} & 0 & 0 & 0 & 0 & 0 & i e^{- i k_{x}} & 0\\0 & i e^{- i k_{y}} & 0 & 0 & 0 & 0 & 0 & - i e^{- i k_{x}}\\0 & 0 & i e^{- i k_{y}} & 0 & - i e^{i k_{x}} & 0 & 0 & 0\\0 & 0 & 0 & - i e^{- i k_{y}} & 0 & i e^{i k_{x}} & 0 & 0\end{matrix}\right]\end{split}\]

\[\begin{split}\displaystyle \left[\begin{matrix}0 & 0 & i e^{i k_{x}} & 0 & - i e^{- i k_{y}} & 0 & 0 & 0\\0 & 0 & 0 & - i e^{i k_{x}} & 0 & i e^{- i k_{y}} & 0 & 0\\- i e^{- i k_{x}} & 0 & 0 & 0 & 0 & 0 & i e^{- i k_{y}} & 0\\0 & i e^{- i k_{x}} & 0 & 0 & 0 & 0 & 0 & - i e^{- i k_{y}}\\i e^{i k_{y}} & 0 & 0 & 0 & 0 & 0 & - i e^{i k_{x}} & 0\\0 & - i e^{i k_{y}} & 0 & 0 & 0 & 0 & 0 & i e^{i k_{x}}\\0 & 0 & - i e^{i k_{y}} & 0 & i e^{- i k_{x}} & 0 & 0 & 0\\0 & 0 & 0 & i e^{i k_{y}} & 0 & - i e^{- i k_{x}} & 0 & 0\end{matrix}\right]\end{split}\]

Generating \(k \cdot p\) models

See also

The complete source code of this example can be found in kdotp_generator.py. A Jupyter notebook can be found in kdotp_generator.ipynb.

In the Bloch generator tutorial we saw how to generate tight-binding Hamiltonians using Qsymm. In this tutorial we will explore how to generate \(k \cdot p\) Hamiltonians.

Hexagonal warping

We reproduce the effective Hamiltonian derived in Phys. Rev. Lett. 103, 266801 (2009) to explain hexagonal warping effects in the surface dispersion of a topological insulator.

The symmetry group is generated by threefold rotation symmetry, mirror symmetry and time-reversal symmetry.

# C3 rotational symmetry - invariant under 2pi/3
C3 = qsymm.rotation(1/3, spin=1/2)

# Time reversal
TR = qsymm.time_reversal(2, spin=1/2)

# Mirror symmetry in x
Mx = qsymm.mirror([1, 0], spin=1/2)

symmetries = [C3, TR, Mx]

The Hamiltonian includes the momenta \(k_x\) and \(k_y\), and contains terms up to third order in momenta.

dim = 2  # Momenta along x and y
total_power = 3  # Maximum power of momenta

We use continuum_hamiltonian to generate the Hamiltonian family:

family = qsymm.continuum_hamiltonian(symmetries, dim, total_power, prettify=True)
qsymm.display_family(family)

\[\begin{split}\displaystyle \left[\begin{matrix}1 & 0\\0 & 1\end{matrix}\right]\end{split}\]

\[\begin{split}\displaystyle \left[\begin{matrix}0 & i k_{x} + k_{y}\\- i k_{x} + k_{y} & 0\end{matrix}\right]\end{split}\]

\[\begin{split}\displaystyle \left[\begin{matrix}k_{x}^{2} + k_{y}^{2} & 0\\0 & k_{x}^{2} + k_{y}^{2}\end{matrix}\right]\end{split}\]

\[\begin{split}\displaystyle \left[\begin{matrix}k_{x}^{3} - 3 k_{x} k_{y}^{2} & 0\\0 & - k_{x}^{3} + 3 k_{x} k_{y}^{2}\end{matrix}\right]\end{split}\]

\[\begin{split}\displaystyle \left[\begin{matrix}0 & i k_{x}^{3} + k_{x}^{2} k_{y} + i k_{x} k_{y}^{2} + k_{y}^{3}\\- i k_{x}^{3} + k_{x}^{2} k_{y} - i k_{x} k_{y}^{2} + k_{y}^{3} & 0\end{matrix}\right]\end{split}\]

We can then verify that the Hamiltonian family satisfies the symmetries:

qsymm.hamiltonian_generator.check_symmetry(family, symmetries)

BHZ model

We reproduce the Hamiltonian for the quantum spin Hall effect derived in Science, 314, 1757 (2006).

The symmetry group is generated by spatial inversion symmetry, time-reversal symmetry and fourfold rotation symmetry.

# Spatial inversion
pU = np.array([
    [1.0, 0.0, 0.0, 0.0],
    [0.0, -1.0, 0.0, 0.0],
    [0.0, 0.0, 1.0, 0.0],
    [0.0, 0.0, 0.0, -1.0],
])
pS = qsymm.inversion(2, U=pU)

# Time reversal
trU = np.array([
    [0.0, 0.0, -1.0, 0.0],
    [0.0, 0.0, 0.0, -1.0],
    [1.0, 0.0, 0.0, 0.0],
    [0.0, 1.0, 0.0, 0.0],
])
trS = qsymm.time_reversal(2, U=trU)

# Rotation
phi = 2.0 * np.pi / 4.0  # Impose 4-fold rotational symmetry
rotU = np.array([
    [np.exp(-1j*phi/2), 0.0, 0.0, 0.0],
    [0.0, np.exp(-1j*3*phi/2), 0.0, 0.0],
    [0.0, 0.0, np.exp(1j*phi/2), 0.0],
    [0.0, 0.0, 0.0, np.exp(1j*3*phi/2)],
])
rotS = qsymm.rotation(1/4, U=rotU)

symmetries = [rotS, trS, pS]

The Hamiltonian includes the momenta \(k_x\) and \(k_y\), with terms up to second order.

dim = 2
total_power = 2

family = qsymm.continuum_hamiltonian(symmetries, dim, total_power, prettify=True)
qsymm.display_family(family)

\[\begin{split}\displaystyle \left[\begin{matrix}1 & 0 & 0 & 0\\0 & 0 & 0 & 0\\0 & 0 & 1 & 0\\0 & 0 & 0 & 0\end{matrix}\right]\end{split}\]

\[\begin{split}\displaystyle \left[\begin{matrix}0 & 0 & 0 & 0\\0 & 1 & 0 & 0\\0 & 0 & 0 & 0\\0 & 0 & 0 & 1\end{matrix}\right]\end{split}\]

\[\begin{split}\displaystyle \left[\begin{matrix}0 & - i k_{x} + k_{y} & 0 & 0\\i k_{x} + k_{y} & 0 & 0 & 0\\0 & 0 & 0 & - i k_{x} - k_{y}\\0 & 0 & i k_{x} - k_{y} & 0\end{matrix}\right]\end{split}\]

\[\begin{split}\displaystyle \left[\begin{matrix}0 & k_{x} + i k_{y} & 0 & 0\\k_{x} - i k_{y} & 0 & 0 & 0\\0 & 0 & 0 & - k_{x} + i k_{y}\\0 & 0 & - k_{x} - i k_{y} & 0\end{matrix}\right]\end{split}\]

\[\begin{split}\displaystyle \left[\begin{matrix}k_{x}^{2} + k_{y}^{2} & 0 & 0 & 0\\0 & 0 & 0 & 0\\0 & 0 & k_{x}^{2} + k_{y}^{2} & 0\\0 & 0 & 0 & 0\end{matrix}\right]\end{split}\]

\[\begin{split}\displaystyle \left[\begin{matrix}0 & 0 & 0 & 0\\0 & k_{x}^{2} + k_{y}^{2} & 0 & 0\\0 & 0 & 0 & 0\\0 & 0 & 0 & k_{x}^{2} + k_{y}^{2}\end{matrix}\right]\end{split}\]

Three-dimensional topological insulator

We reproduce the Hamiltonian for a three-dimensional topological insulator introduced in Nature Physics 5, 438–442 (2009).

The symmetry group is generated by threefold rotation symmetry around z, inversion symmetry and time-reversal symmetry.

# Spatial inversion
pU = np.diag([1, -1, 1, -1])
pS = qsymm.inversion(3, pU)

# Time reversal
trU = np.array([
    [0.0, 0.0, -1.0, 0.0],
    [0.0, 0.0, 0.0, -1.0],
    [1.0, 0.0, 0.0, 0.0],
    [0.0, 1.0, 0.0, 0.0],
])
trS = qsymm.time_reversal(3, trU)

# Rotation
phi = 2.0 * np.pi / 3.0  # Impose 3-fold rotational symmetry
rotU = np.array([
    [np.exp(1j*phi/2), 0.0, 0.0, 0.0],
    [0.0, np.exp(1j*phi/2), 0.0, 0.0],
    [0.0, 0.0, np.exp(-1j*phi/2), 0.0],
    [0.0, 0.0, 0.0, np.exp(-1j*phi/2)],
])
rotS = qsymm.rotation(1/3, axis=[0, 0, 1], U=rotU)

symmetries = [rotS, trS, pS]

The model includes the momenta \(k_x\), \(k_y\) and \(k_z\), up to second order.

dim = 3
total_power = 2

family = qsymm.continuum_hamiltonian(symmetries, dim, total_power, prettify=True)
qsymm.display_family(family)

\[\begin{split}\displaystyle \left[\begin{matrix}1 & 0 & 0 & 0\\0 & 0 & 0 & 0\\0 & 0 & 1 & 0\\0 & 0 & 0 & 0\end{matrix}\right]\end{split}\]

\[\begin{split}\displaystyle \left[\begin{matrix}0 & 0 & 0 & 0\\0 & 1 & 0 & 0\\0 & 0 & 0 & 0\\0 & 0 & 0 & 1\end{matrix}\right]\end{split}\]

\[\begin{split}\displaystyle \left[\begin{matrix}0 & 0 & 0 & - i k_{x} + k_{y}\\0 & 0 & - i k_{x} + k_{y} & 0\\0 & i k_{x} + k_{y} & 0 & 0\\i k_{x} + k_{y} & 0 & 0 & 0\end{matrix}\right]\end{split}\]

\[\begin{split}\displaystyle \left[\begin{matrix}0 & 0 & 0 & k_{x} + i k_{y}\\0 & 0 & k_{x} + i k_{y} & 0\\0 & k_{x} - i k_{y} & 0 & 0\\k_{x} - i k_{y} & 0 & 0 & 0\end{matrix}\right]\end{split}\]

\[\begin{split}\displaystyle \left[\begin{matrix}0 & k_{z} & 0 & 0\\k_{z} & 0 & 0 & 0\\0 & 0 & 0 & - k_{z}\\0 & 0 & - k_{z} & 0\end{matrix}\right]\end{split}\]

\[\begin{split}\displaystyle \left[\begin{matrix}0 & i k_{z} & 0 & 0\\- i k_{z} & 0 & 0 & 0\\0 & 0 & 0 & i k_{z}\\0 & 0 & - i k_{z} & 0\end{matrix}\right]\end{split}\]

\[\begin{split}\displaystyle \left[\begin{matrix}k_{x}^{2} + k_{y}^{2} & 0 & 0 & 0\\0 & 0 & 0 & 0\\0 & 0 & k_{x}^{2} + k_{y}^{2} & 0\\0 & 0 & 0 & 0\end{matrix}\right]\end{split}\]

\[\begin{split}\displaystyle \left[\begin{matrix}0 & 0 & 0 & 0\\0 & k_{x}^{2} + k_{y}^{2} & 0 & 0\\0 & 0 & 0 & 0\\0 & 0 & 0 & k_{x}^{2} + k_{y}^{2}\end{matrix}\right]\end{split}\]

\[\begin{split}\displaystyle \left[\begin{matrix}k_{z}^{2} & 0 & 0 & 0\\0 & 0 & 0 & 0\\0 & 0 & k_{z}^{2} & 0\\0 & 0 & 0 & 0\end{matrix}\right]\end{split}\]

\[\begin{split}\displaystyle \left[\begin{matrix}0 & 0 & 0 & 0\\0 & k_{z}^{2} & 0 & 0\\0 & 0 & 0 & 0\\0 & 0 & 0 & k_{z}^{2}\end{matrix}\right]\end{split}\]

Continuous rotations

Rotation in real and spin space with spin 1/2. Should have linear Rashba-term as there is no inversion symmetry.

# 3D real space rotation generators
L = qsymm.groups.L_matrices(3)
# Spin-1/2 matrices
S = qsymm.groups.spin_matrices(1/2)
# Spin-3/2 spin matrices
J = qsymm.groups.spin_matrices(3/2)
# Continuous rotation generators
symmetries = [qsymm.ContinuousGroupGenerator(l, s) for l, s in zip(L, S)]

dim = 3
total_power = 2

family = qsymm.continuum_hamiltonian(symmetries, dim, total_power, prettify=True)
qsymm.display_family(family)

\[\begin{split}\displaystyle \left[\begin{matrix}1 & 0\\0 & 1\end{matrix}\right]\end{split}\]

\[\begin{split}\displaystyle \left[\begin{matrix}- k_{z} & - k_{x} + i k_{y}\\- k_{x} - i k_{y} & k_{z}\end{matrix}\right]\end{split}\]

\[\begin{split}\displaystyle \left[\begin{matrix}k_{x}^{2} + k_{y}^{2} + k_{z}^{2} & 0\\0 & k_{x}^{2} + k_{y}^{2} + k_{z}^{2}\end{matrix}\right]\end{split}\]

Also impose inversion symmetry, which removes the linear Rashba term.

symmetries = [qsymm.ContinuousGroupGenerator(l, s) for l, s in zip(L, S)]
# Add inversion
symmetries.append(qsymm.PointGroupElement(-np.eye(3), False, False, np.eye(2)))
dim = 3
total_power = 2
family = qsymm.continuum_hamiltonian(symmetries, dim, total_power, prettify=True)
qsymm.display_family(family)

\[\begin{split}\displaystyle \left[\begin{matrix}1 & 0\\0 & 1\end{matrix}\right]\end{split}\]

\[\begin{split}\displaystyle \left[\begin{matrix}k_{x}^{2} + k_{y}^{2} + k_{z}^{2} & 0\\0 & k_{x}^{2} + k_{y}^{2} + k_{z}^{2}\end{matrix}\right]\end{split}\]

Rotations in real and spin space with spin 3/2.

symmetries = [qsymm.ContinuousGroupGenerator(l, s) for l, s in zip(L, J)]
dim = 3
total_power = 2
family = qsymm.continuum_hamiltonian(symmetries, dim, total_power, prettify=True)
qsymm.display_family(family)

\[\begin{split}\displaystyle \left[\begin{matrix}1 & 0 & 0 & 0\\0 & 1 & 0 & 0\\0 & 0 & 1 & 0\\0 & 0 & 0 & 1\end{matrix}\right]\end{split}\]

\[\begin{split}\displaystyle \left[\begin{matrix}- \sqrt{3} k_{z} & - k_{x} + i k_{y} & 0 & 0\\- k_{x} - i k_{y} & - \frac{\sqrt{3} k_{z}}{3} & - \frac{2 \sqrt{3} k_{x}}{3} + \frac{2 \sqrt{3} i k_{y}}{3} & 0\\0 & - \frac{2 \sqrt{3} k_{x}}{3} - \frac{2 \sqrt{3} i k_{y}}{3} & \frac{\sqrt{3} k_{z}}{3} & - k_{x} + i k_{y}\\0 & 0 & - k_{x} - i k_{y} & \sqrt{3} k_{z}\end{matrix}\right]\end{split}\]

\[\begin{split}\displaystyle \left[\begin{matrix}k_{x}^{2} + k_{y}^{2} + k_{z}^{2} & 0 & 0 & 0\\0 & k_{x}^{2} + k_{y}^{2} + k_{z}^{2} & 0 & 0\\0 & 0 & k_{x}^{2} + k_{y}^{2} + k_{z}^{2} & 0\\0 & 0 & 0 & k_{x}^{2} + k_{y}^{2} + k_{z}^{2}\end{matrix}\right]\end{split}\]

\[\begin{split}\displaystyle \left[\begin{matrix}\sqrt{3} k_{z}^{2} & 2 k_{x} k_{z} - 2 i k_{y} k_{z} & k_{x}^{2} - 2 i k_{x} k_{y} - k_{y}^{2} & 0\\2 k_{x} k_{z} + 2 i k_{y} k_{z} & \frac{2 \sqrt{3} k_{x}^{2}}{3} + \frac{2 \sqrt{3} k_{y}^{2}}{3} - \frac{\sqrt{3} k_{z}^{2}}{3} & 0 & k_{x}^{2} - 2 i k_{x} k_{y} - k_{y}^{2}\\k_{x}^{2} + 2 i k_{x} k_{y} - k_{y}^{2} & 0 & \frac{2 \sqrt{3} k_{x}^{2}}{3} + \frac{2 \sqrt{3} k_{y}^{2}}{3} - \frac{\sqrt{3} k_{z}^{2}}{3} & - 2 k_{x} k_{z} + 2 i k_{y} k_{z}\\0 & k_{x}^{2} + 2 i k_{x} k_{y} - k_{y}^{2} & - 2 k_{x} k_{z} - 2 i k_{y} k_{z} & \sqrt{3} k_{z}^{2}\end{matrix}\right]\end{split}\]

Distorted SnTe

Reproduce the \(k \cdot p\) model for SnTe used in Phys. Rev. Lett. 122, 186801 (2019)

The symmetry group is generated by three-fold rotation symmetry, mirror symmetry in \(x\), inversion symmetry, and time-reversal symmetry.

# Double spin-1/2 representation
spin = [np.kron(s, np.eye(2)) for s in qsymm.groups.spin_matrices(1/2)]
# C3 rotational symmetry
C3 = qsymm.rotation(1/3, axis=[0, 0, 1], spin=spin)

# Time reversal
TR = qsymm.time_reversal(3, spin=spin)

# Mirror x
Mx = qsymm.mirror([1, 0, 0], spin=spin)

# Inversion
IU = np.kron(np.eye(2), np.diag([1, -1]))
I = qsymm.inversion(3, IU)

dim = 3
total_power = 2
family = qsymm.continuum_hamiltonian([C3, TR, Mx, I], dim, total_power, prettify=True)
qsymm.display_family(family)

\[\begin{split}\displaystyle \left[\begin{matrix}1 & 0 & 0 & 0\\0 & 0 & 0 & 0\\0 & 0 & 1 & 0\\0 & 0 & 0 & 0\end{matrix}\right]\end{split}\]

\[\begin{split}\displaystyle \left[\begin{matrix}0 & 0 & 0 & 0\\0 & 1 & 0 & 0\\0 & 0 & 0 & 0\\0 & 0 & 0 & 1\end{matrix}\right]\end{split}\]

\[\begin{split}\displaystyle \left[\begin{matrix}0 & 0 & 0 & i k_{x} + k_{y}\\0 & 0 & i k_{x} + k_{y} & 0\\0 & - i k_{x} + k_{y} & 0 & 0\\- i k_{x} + k_{y} & 0 & 0 & 0\end{matrix}\right]\end{split}\]

\[\begin{split}\displaystyle \left[\begin{matrix}0 & i k_{z} & 0 & 0\\- i k_{z} & 0 & 0 & 0\\0 & 0 & 0 & i k_{z}\\0 & 0 & - i k_{z} & 0\end{matrix}\right]\end{split}\]

\[\begin{split}\displaystyle \left[\begin{matrix}k_{x}^{2} + k_{y}^{2} & 0 & 0 & 0\\0 & 0 & 0 & 0\\0 & 0 & k_{x}^{2} + k_{y}^{2} & 0\\0 & 0 & 0 & 0\end{matrix}\right]\end{split}\]

\[\begin{split}\displaystyle \left[\begin{matrix}0 & 0 & 0 & 0\\0 & k_{x}^{2} + k_{y}^{2} & 0 & 0\\0 & 0 & 0 & 0\\0 & 0 & 0 & k_{x}^{2} + k_{y}^{2}\end{matrix}\right]\end{split}\]

\[\begin{split}\displaystyle \left[\begin{matrix}k_{z}^{2} & 0 & 0 & 0\\0 & 0 & 0 & 0\\0 & 0 & k_{z}^{2} & 0\\0 & 0 & 0 & 0\end{matrix}\right]\end{split}\]

\[\begin{split}\displaystyle \left[\begin{matrix}0 & 0 & 0 & 0\\0 & k_{z}^{2} & 0 & 0\\0 & 0 & 0 & 0\\0 & 0 & 0 & k_{z}^{2}\end{matrix}\right]\end{split}\]

There are 3 combinations proportional to the identity. Generate all terms proportional to the identity and subtract them.

identity_terms = qsymm.continuum_hamiltonian([qsymm.groups.identity(dim, 1)], dim, total_power)
identity_terms = [
    qsymm.Model({
        key: np.kron(val, np.eye(4))
        for key, val in term.items()
    })
    for term in identity_terms
]

family = qsymm.hamiltonian_generator.subtract_family(family, identity_terms, prettify=True)
qsymm.display_family(family)

\[\begin{split}\displaystyle \left[\begin{matrix}1 & 0 & 0 & 0\\0 & -1 & 0 & 0\\0 & 0 & 1 & 0\\0 & 0 & 0 & -1\end{matrix}\right]\end{split}\]

\[\begin{split}\displaystyle \left[\begin{matrix}k_{x}^{2} + k_{y}^{2} & 0 & 0 & 0\\0 & - k_{x}^{2} - k_{y}^{2} & 0 & 0\\0 & 0 & k_{x}^{2} + k_{y}^{2} & 0\\0 & 0 & 0 & - k_{x}^{2} - k_{y}^{2}\end{matrix}\right]\end{split}\]

\[\begin{split}\displaystyle \left[\begin{matrix}k_{z}^{2} & 0 & 0 & 0\\0 & - k_{z}^{2} & 0 & 0\\0 & 0 & k_{z}^{2} & 0\\0 & 0 & 0 & - k_{z}^{2}\end{matrix}\right]\end{split}\]

\[\begin{split}\displaystyle \left[\begin{matrix}0 & 0 & 0 & i k_{x} + k_{y}\\0 & 0 & i k_{x} + k_{y} & 0\\0 & - i k_{x} + k_{y} & 0 & 0\\- i k_{x} + k_{y} & 0 & 0 & 0\end{matrix}\right]\end{split}\]

\[\begin{split}\displaystyle \left[\begin{matrix}0 & i k_{z} & 0 & 0\\- i k_{z} & 0 & 0 & 0\\0 & 0 & 0 & i k_{z}\\0 & 0 & - i k_{z} & 0\end{matrix}\right]\end{split}\]

Break C3 symmetry

Removing the C3 symmetry, we find 8 additional terms, or 11 after removing terms proportional to the identity.

dim = 3
total_power = 2
no_c3_family = qsymm.continuum_hamiltonian([TR, Mx, I], dim, total_power, prettify=True)
no_c3_family = qsymm.hamiltonian_generator.subtract_family(no_c3_family, identity_terms, prettify=True)

The new terms are given in the following.

new_terms = qsymm.hamiltonian_generator.subtract_family(no_c3_family, family, prettify=True)
qsymm.display_family(new_terms)

\[\begin{split}\displaystyle \left[\begin{matrix}k_{x}^{2} - k_{y}^{2} & 0 & 0 & 0\\0 & - k_{x}^{2} + k_{y}^{2} & 0 & 0\\0 & 0 & k_{x}^{2} - k_{y}^{2} & 0\\0 & 0 & 0 & - k_{x}^{2} + k_{y}^{2}\end{matrix}\right]\end{split}\]

\[\begin{split}\displaystyle \left[\begin{matrix}0 & 0 & 0 & - i k_{x} + k_{y}\\0 & 0 & - i k_{x} + k_{y} & 0\\0 & i k_{x} + k_{y} & 0 & 0\\i k_{x} + k_{y} & 0 & 0 & 0\end{matrix}\right]\end{split}\]

\[\begin{split}\displaystyle \left[\begin{matrix}0 & i k_{y} & 0 & 0\\- i k_{y} & 0 & 0 & 0\\0 & 0 & 0 & i k_{y}\\0 & 0 & - i k_{y} & 0\end{matrix}\right]\end{split}\]

\[\begin{split}\displaystyle \left[\begin{matrix}0 & k_{x} & 0 & 0\\k_{x} & 0 & 0 & 0\\0 & 0 & 0 & - k_{x}\\0 & 0 & - k_{x} & 0\end{matrix}\right]\end{split}\]

\[\begin{split}\displaystyle \left[\begin{matrix}k_{y} k_{z} & 0 & 0 & 0\\0 & - k_{y} k_{z} & 0 & 0\\0 & 0 & k_{y} k_{z} & 0\\0 & 0 & 0 & - k_{y} k_{z}\end{matrix}\right]\end{split}\]

\[\begin{split}\displaystyle \left[\begin{matrix}0 & 0 & 0 & k_{z}\\0 & 0 & k_{z} & 0\\0 & k_{z} & 0 & 0\\k_{z} & 0 & 0 & 0\end{matrix}\right]\end{split}\]

Break inversion symmetry

Breaking inversion symmetry as well yields 22 additional terms.

dim = 3
total_power = 2
broken_family = qsymm.continuum_hamiltonian([TR, Mx], dim, total_power, prettify=True)
broken_family = qsymm.hamiltonian_generator.subtract_family(broken_family, identity_terms, prettify=True)
new_terms = qsymm.hamiltonian_generator.subtract_family(broken_family, no_c3_family, prettify=True)
qsymm.display_family(new_terms)

\[\begin{split}\displaystyle \left[\begin{matrix}0 & 1 & 0 & 0\\1 & 0 & 0 & 0\\0 & 0 & 0 & 1\\0 & 0 & 1 & 0\end{matrix}\right]\end{split}\]

\[\begin{split}\displaystyle \left[\begin{matrix}0 & 0 & 0 & i\\0 & 0 & - i & 0\\0 & i & 0 & 0\\- i & 0 & 0 & 0\end{matrix}\right]\end{split}\]

\[\begin{split}\displaystyle \left[\begin{matrix}0 & k_{x}^{2} & 0 & 0\\k_{x}^{2} & 0 & 0 & 0\\0 & 0 & 0 & k_{x}^{2}\\0 & 0 & k_{x}^{2} & 0\end{matrix}\right]\end{split}\]

\[\begin{split}\displaystyle \left[\begin{matrix}0 & 0 & 0 & i k_{x}^{2}\\0 & 0 & - i k_{x}^{2} & 0\\0 & i k_{x}^{2} & 0 & 0\\- i k_{x}^{2} & 0 & 0 & 0\end{matrix}\right]\end{split}\]

\[\begin{split}\displaystyle \left[\begin{matrix}0 & k_{z}^{2} & 0 & 0\\k_{z}^{2} & 0 & 0 & 0\\0 & 0 & 0 & k_{z}^{2}\\0 & 0 & k_{z}^{2} & 0\end{matrix}\right]\end{split}\]

\[\begin{split}\displaystyle \left[\begin{matrix}0 & 0 & 0 & i k_{z}^{2}\\0 & 0 & - i k_{z}^{2} & 0\\0 & i k_{z}^{2} & 0 & 0\\- i k_{z}^{2} & 0 & 0 & 0\end{matrix}\right]\end{split}\]

\[\begin{split}\displaystyle \left[\begin{matrix}0 & k_{y}^{2} & 0 & 0\\k_{y}^{2} & 0 & 0 & 0\\0 & 0 & 0 & k_{y}^{2}\\0 & 0 & k_{y}^{2} & 0\end{matrix}\right]\end{split}\]

\[\begin{split}\displaystyle \left[\begin{matrix}0 & 0 & 0 & i k_{y}^{2}\\0 & 0 & - i k_{y}^{2} & 0\\0 & i k_{y}^{2} & 0 & 0\\- i k_{y}^{2} & 0 & 0 & 0\end{matrix}\right]\end{split}\]

\[\begin{split}\displaystyle \left[\begin{matrix}0 & 0 & 0 & k_{x} k_{y}\\0 & 0 & - k_{x} k_{y} & 0\\0 & - k_{x} k_{y} & 0 & 0\\k_{x} k_{y} & 0 & 0 & 0\end{matrix}\right]\end{split}\]

\[\begin{split}\displaystyle \left[\begin{matrix}0 & i k_{x} k_{y} & 0 & 0\\- i k_{x} k_{y} & 0 & 0 & 0\\0 & 0 & 0 & - i k_{x} k_{y}\\0 & 0 & i k_{x} k_{y} & 0\end{matrix}\right]\end{split}\]

\[\begin{split}\displaystyle \left[\begin{matrix}0 & 0 & k_{y} & 0\\0 & 0 & 0 & 0\\k_{y} & 0 & 0 & 0\\0 & 0 & 0 & 0\end{matrix}\right]\end{split}\]

\[\begin{split}\displaystyle \left[\begin{matrix}0 & 0 & 0 & 0\\0 & 0 & 0 & k_{y}\\0 & 0 & 0 & 0\\0 & k_{y} & 0 & 0\end{matrix}\right]\end{split}\]

\[\begin{split}\displaystyle \left[\begin{matrix}0 & 0 & 0 & k_{x} k_{z}\\0 & 0 & - k_{x} k_{z} & 0\\0 & - k_{x} k_{z} & 0 & 0\\k_{x} k_{z} & 0 & 0 & 0\end{matrix}\right]\end{split}\]

\[\begin{split}\displaystyle \left[\begin{matrix}0 & i k_{x} k_{z} & 0 & 0\\- i k_{x} k_{z} & 0 & 0 & 0\\0 & 0 & 0 & - i k_{x} k_{z}\\0 & 0 & i k_{x} k_{z} & 0\end{matrix}\right]\end{split}\]

\[\begin{split}\displaystyle \left[\begin{matrix}k_{x} & 0 & 0 & 0\\0 & 0 & 0 & 0\\0 & 0 & - k_{x} & 0\\0 & 0 & 0 & 0\end{matrix}\right]\end{split}\]

\[\begin{split}\displaystyle \left[\begin{matrix}0 & 0 & 0 & 0\\0 & k_{x} & 0 & 0\\0 & 0 & 0 & 0\\0 & 0 & 0 & - k_{x}\end{matrix}\right]\end{split}\]

\[\begin{split}\displaystyle \left[\begin{matrix}0 & 0 & i k_{x} & 0\\0 & 0 & 0 & 0\\- i k_{x} & 0 & 0 & 0\\0 & 0 & 0 & 0\end{matrix}\right]\end{split}\]

\[\begin{split}\displaystyle \left[\begin{matrix}0 & 0 & 0 & 0\\0 & 0 & 0 & i k_{x}\\0 & 0 & 0 & 0\\0 & - i k_{x} & 0 & 0\end{matrix}\right]\end{split}\]

\[\begin{split}\displaystyle \left[\begin{matrix}0 & k_{y} k_{z} & 0 & 0\\k_{y} k_{z} & 0 & 0 & 0\\0 & 0 & 0 & k_{y} k_{z}\\0 & 0 & k_{y} k_{z} & 0\end{matrix}\right]\end{split}\]

\[\begin{split}\displaystyle \left[\begin{matrix}0 & 0 & 0 & i k_{y} k_{z}\\0 & 0 & - i k_{y} k_{z} & 0\\0 & i k_{y} k_{z} & 0 & 0\\- i k_{y} k_{z} & 0 & 0 & 0\end{matrix}\right]\end{split}\]

\[\begin{split}\displaystyle \left[\begin{matrix}0 & 0 & k_{z} & 0\\0 & 0 & 0 & 0\\k_{z} & 0 & 0 & 0\\0 & 0 & 0 & 0\end{matrix}\right]\end{split}\]

\[\begin{split}\displaystyle \left[\begin{matrix}0 & 0 & 0 & 0\\0 & 0 & 0 & k_{z}\\0 & 0 & 0 & 0\\0 & k_{z} & 0 & 0\end{matrix}\right]\end{split}\]

Finding symmetries of the Kekule-Y continuum model

See also

The complete source code of this example can be found in kekule.py. A Jupyter notebook can be found in kekule.ipynb.

Kekule-Y

Find symmetries of Kekule-Y effective model, starting from the full hexagonal symmetry group with time reversal and particle-hole symmetry.

ham_kek_y = """
    vs * kron((k_x * sigma_x + k_y * sigma_y), eye(2))
    + vt * kron(eye(2), (k_x * sigma_x + k_y * sigma_y))
"""
display(qsymm.sympify(ham_kek_y))

\[\begin{split}\displaystyle \left[\begin{matrix}0 & vt \left(k_{x} - i k_{y}\right) & vs \left(k_{x} - i k_{y}\right) & 0\\vt \left(k_{x} + i k_{y}\right) & 0 & 0 & vs \left(k_{x} - i k_{y}\right)\\vs \left(k_{x} + i k_{y}\right) & 0 & 0 & vt \left(k_{x} - i k_{y}\right)\\0 & vs \left(k_{x} + i k_{y}\right) & vt \left(k_{x} + i k_{y}\right) & 0\end{matrix}\right]\end{split}\]

H_kek_y= qsymm.Model(ham_kek_y, momenta=['k_x', 'k_y'])
candidates = qsymm.groups.hexagonal()
sgy, cg = qsymm.symmetries(H_kek_y, candidates)

Enumerate subgroups that protect double Dirac cone:

all_subgroups = qsymm.groups.generate_subgroups(sgy)

# Terms proportional to the identity matrix don't open a gap.
identity_terms = [qsymm.Model({qsymm.sympify(1): np.eye(4)})]

dd_sg = []
for subgroup, gens in all_subgroups.items():
    family = qsymm.continuum_hamiltonian(gens, 2, 0)
    family = qsymm.hamiltonian_generator.subtract_family(family, identity_terms)
    if family == []:
        dd_sg.append(subgroup)

dd_sg.sort(key=len)

dd_sg[0]

frozenset({1, R(π) C})

Twofold rotation combined with sublattice symmetry is enough. However 2-fold rotation is not a physical symmetry of the lattice model. Let’s get rid of 2-fold rotations:

rdd_sg = []
for subgroup, gens in all_subgroups.items():
    invs = {
        qsymm.inversion(2) * g
        for g in [
            qsymm.identity(2),
            qsymm.time_reversal(2),
            qsymm.particle_hole(2),
            qsymm.chiral(2),
        ]
    }
    # Skip subgroups that have inversion
    if subgroup & invs:
        continue
    family = qsymm.continuum_hamiltonian(gens, 2, 0)
    family = qsymm.hamiltonian_generator.subtract_family(family, identity_terms)
    if family == []:
        rdd_sg.append(subgroup)

rdd_sg.sort(key=len)

rdd_sg[0]

frozenset({1, R(2π/3), R(-2π/3), C, R(2π/3) C, R(-2π/3) C})

Sublattice symmetry with 3-fold rotations is the minimal symmetry group required. Let’s test if any of the solutions has no sublattice symmetry:

any(g.antisymmetry for sg in rdd_sg for g in sg)

True

Sublattice symmetry is required for the protection of the double Dirac cone.

Continuous rotation symmetries of continuum models

Kekule-Y

_, cg = qsymm.symmetries(H_kek_y, candidates=[], continuous_rotations=True, prettify=True)
cg

\[\begin{split}e^{-i\phi L}H(\mathbf{k})e^{i\phi L} = H(R_{\phi}\mathbf{k}) \\R_{{\phi}} = R\left(\phi\right)\\L = \begin{bmatrix}1.+0.j&0.+0.j&0.+0.j&0.+0.j\\0.+0.j&0.+0.j&0.+0.j&0.+0.j\\0.+0.j&0.+0.j&0.+0.j&0.+0.j\\ 0.+0.j& 0.+0.j& 0.+0.j&-1.+0.j\\\end{bmatrix} \\\end{split}\]

Kekule-O

ham_kek_o= """
    vs * (
        k_x * kron(eye(2), sigma_x)
        + k_y * kron(sigma_z, sigma_y)
    )
    + Delta * kron(sigma_y, sigma_y)
"""
display(qsymm.sympify(ham_kek_o))

\[\begin{split}\displaystyle \left[\begin{matrix}0 & vs \left(k_{x} - i k_{y}\right) & 0 & - \Delta\\vs \left(k_{x} + i k_{y}\right) & 0 & \Delta & 0\\0 & \Delta & 0 & vs \left(k_{x} + i k_{y}\right)\\- \Delta & 0 & vs \left(k_{x} - i k_{y}\right) & 0\end{matrix}\right]\end{split}\]

H_kek_o= qsymm.Model(ham_kek_o, momenta=['k_x', 'k_y'])
_, cg = qsymm.symmetries(H_kek_o, candidates=[], continuous_rotations=True, prettify=True)
cg

\[\begin{split}\left[ H(\mathbf{{k}}), L \right] = 0 \\L = \begin{bmatrix} 0. +0.j&-0. +0.j&-0. +0.j& 0.5+0.j\\-0. -0.j& 0. +0.j& 0.5-0.j&-0. -0.j\\-0. -0.j& 0.5+0.j& 0. +0.j&-0. +0.j\\ 0.5-0.j&-0. +0.j&-0. -0.j& 0. +0.j\\\end{bmatrix} \\\\e^{-i\phi L}H(\mathbf{k})e^{i\phi L} = H(R_{\phi}\mathbf{k}) \\R_{{\phi}} = R\left(\phi\right)\\L = \begin{bmatrix} 0.5+0.j&-0. -0.j&-0. -0.j&-0. +0.j\\-0. +0.j&-0.5-0.j& 0. -0.j&-0. +0.j\\-0. +0.j& 0. +0.j&-0.5+0.j&-0. +0.j\\-0. -0.j&-0. -0.j&-0. -0.j& 0.5+0.j\\\end{bmatrix} \\\end{split}\]

The rotation generators are different in the two cases

Lattice models

Generate Kekule-O model

6 sites (3 A and 3 B) per unit cell, spinless orbitals. Time-reversal, sublattice and 6-fold rotation and mirror symmetry.

sites = ['A1', 'B1', 'A2', 'B2', 'A3', 'B3']
norbs = [(site, 1) for site in sites]

# Time reversal
TR = qsymm.PointGroupElement(sympy.eye(2), True, False, np.eye(6))

# Chiral symmetry
C = qsymm.PointGroupElement(sympy.eye(2), False, True, np.kron(np.eye(3), ([[1, 0], [0, -1]])))

# Atom A rotates into B, B into A.
sphi = 2*sympy.pi/6
RC6 = sympy.Matrix([[sympy.cos(sphi), -sympy.sin(sphi)],
                  [sympy.sin(sphi), sympy.cos(sphi)]])
permC6 = {'A1': 'B1', 'B1': 'A2', 'A2': 'B2', 'B2': 'A3', 'A3': 'B3', 'B3': 'A1'}
C6 = qsymm.groups.symmetry_from_permutation(RC6, permC6, norbs)

RMx = sympy.Matrix([[-1, 0], [0, 1]])
permMx = {'A1': 'A3', 'B1': 'B2', 'A2': 'A2', 'B2': 'B1', 'A3': 'A1', 'B3': 'B3'}
Mx = qsymm.groups.symmetry_from_permutation(RMx, permMx, norbs)

symmetries = [C, TR, C6, Mx]

# Only need the unique hoppings, one weak and one strong
hopping_vectors = [
    ('A1', 'B1', np.array([0, -1])), # around ring
    ('A2', 'B3', np.array([0, -1])), # next UC
]

# construct Hamiltonian terms to symmetry-allowed terms
family = qsymm.bloch_family(hopping_vectors, symmetries, norbs=norbs)
qsymm.display_family(family)

\[\begin{split}\displaystyle \left[\begin{matrix}0 & e^{- i k_{y}} & 0 & 0 & 0 & e^{\frac{i k_{y}}{2}} e^{- \frac{\sqrt{3} i k_{x}}{2}}\\e^{i k_{y}} & 0 & e^{- \frac{i k_{y}}{2}} e^{- \frac{\sqrt{3} i k_{x}}{2}} & 0 & 0 & 0\\0 & e^{\frac{i k_{y}}{2}} e^{\frac{\sqrt{3} i k_{x}}{2}} & 0 & e^{\frac{i k_{y}}{2}} e^{- \frac{\sqrt{3} i k_{x}}{2}} & 0 & 0\\0 & 0 & e^{- \frac{i k_{y}}{2}} e^{\frac{\sqrt{3} i k_{x}}{2}} & 0 & e^{i k_{y}} & 0\\0 & 0 & 0 & e^{- i k_{y}} & 0 & e^{\frac{i k_{y}}{2}} e^{\frac{\sqrt{3} i k_{x}}{2}}\\e^{- \frac{i k_{y}}{2}} e^{\frac{\sqrt{3} i k_{x}}{2}} & 0 & 0 & 0 & e^{- \frac{i k_{y}}{2}} e^{- \frac{\sqrt{3} i k_{x}}{2}} & 0\end{matrix}\right]\end{split}\]

\[\begin{split}\displaystyle \left[\begin{matrix}0 & 0 & 0 & e^{\frac{i k_{y}}{2}} e^{\frac{\sqrt{3} i k_{x}}{2}} & 0 & 0\\0 & 0 & 0 & 0 & e^{- \frac{i k_{y}}{2}} e^{\frac{\sqrt{3} i k_{x}}{2}} & 0\\0 & 0 & 0 & 0 & 0 & e^{- i k_{y}}\\e^{- \frac{i k_{y}}{2}} e^{- \frac{\sqrt{3} i k_{x}}{2}} & 0 & 0 & 0 & 0 & 0\\0 & e^{\frac{i k_{y}}{2}} e^{- \frac{\sqrt{3} i k_{x}}{2}} & 0 & 0 & 0 & 0\\0 & 0 & e^{i k_{y}} & 0 & 0 & 0\end{matrix}\right]\end{split}\]

Make terms at k=0 and check degeneracy for some linear combination.

Gamma_terms = [term.subs(dict(k_x=0, k_y=0)) for term in family]
evals, evecs = np.linalg.eigh(
    np.sum([
        list(term.values())[0] * coeff
        for term, coeff in zip(Gamma_terms, [1, 2])
    ],
    axis=0)
)

evals

array([-4., -1., -1.,  1.,  1.,  4.])

Find symmetry representation in the low energy subspace

# projector to low energy subspace
proj = evecs[:, 1:5]
UC3 = (C6 * C6).U

# projected symmetry operator
pC3 = proj.T.conj() @ UC3 @ proj

assert np.allclose(pC3 @ pC3.T.conj(), np.eye(4))

evals2, evecs2 = np.linalg.eig(pC3)

# Check angles of the eigenvalues
np.angle(evals2) / (2 * np.pi)

array([ 0.33333333, -0.33333333,  0.33333333, -0.33333333])

Generate Kekule-Y model

6 sites (3 A and 3 B) per unit cell, spinless orbitals. Time-reversal, sublattice and 3-fold rotation and mirror symmetry.

sites = ['A1', 'B1', 'A2', 'B2', 'A3', 'B3']
norbs = [(site, 1) for site in sites]

# Time reversal
TR = qsymm.time_reversal(2, U=np.eye(6))

# Chiral symmetry
C = qsymm.chiral(2, U=np.kron(np.eye(3), ([[1, 0], [0, -1]])))

# Atom A rotates into B, B into A.
sphi = 2*sympy.pi/3
RC3 = sympy.Matrix([
    [sympy.cos(sphi), -sympy.sin(sphi)],
    [sympy.sin(sphi), sympy.cos(sphi)],
])
permC3 = {'A1': 'A1', 'B1': 'B3', 'A2': 'A2', 'B2': 'B1', 'A3': 'A3', 'B3': 'B2'}
C3 = qsymm.groups.symmetry_from_permutation(RC3, permC3, norbs)

RMx = sympy.Matrix([[-1, 0], [0, 1]])
permMx = {'A1': 'A1', 'B1': 'B1', 'A2': 'A3', 'B2': 'B3', 'A3': 'A2', 'B3': 'B2'}
Mx = qsymm.groups.symmetry_from_permutation(RMx, permMx, norbs)

symmetries = [C, TR, C3, Mx]

# Only need the unique hoppings, one weak and one strong
hopping_vectors = [
    ('A1', 'B1', np.array([0, -1])), # Y
    ('A2', 'B3', np.array([0, -1])), # other
]

# construct Hamiltonian terms to symmetry-allowed terms
family = qsymm.bloch_family(hopping_vectors, symmetries, norbs=norbs)
qsymm.display_family(family)

\[\begin{split}\displaystyle \left[\begin{matrix}0 & e^{- i k_{y}} & 0 & e^{\frac{i k_{y}}{2}} e^{\frac{\sqrt{3} i k_{x}}{2}} & 0 & e^{\frac{i k_{y}}{2}} e^{- \frac{\sqrt{3} i k_{x}}{2}}\\e^{i k_{y}} & 0 & 0 & 0 & 0 & 0\\0 & 0 & 0 & 0 & 0 & 0\\e^{- \frac{i k_{y}}{2}} e^{- \frac{\sqrt{3} i k_{x}}{2}} & 0 & 0 & 0 & 0 & 0\\0 & 0 & 0 & 0 & 0 & 0\\e^{- \frac{i k_{y}}{2}} e^{\frac{\sqrt{3} i k_{x}}{2}} & 0 & 0 & 0 & 0 & 0\end{matrix}\right]\end{split}\]

\[\begin{split}\displaystyle \left[\begin{matrix}0 & 0 & 0 & 0 & 0 & 0\\0 & 0 & e^{- \frac{i k_{y}}{2}} e^{- \frac{\sqrt{3} i k_{x}}{2}} & 0 & e^{- \frac{i k_{y}}{2}} e^{\frac{\sqrt{3} i k_{x}}{2}} & 0\\0 & e^{\frac{i k_{y}}{2}} e^{\frac{\sqrt{3} i k_{x}}{2}} & 0 & e^{\frac{i k_{y}}{2}} e^{- \frac{\sqrt{3} i k_{x}}{2}} & 0 & e^{- i k_{y}}\\0 & 0 & e^{- \frac{i k_{y}}{2}} e^{\frac{\sqrt{3} i k_{x}}{2}} & 0 & e^{i k_{y}} & 0\\0 & e^{\frac{i k_{y}}{2}} e^{- \frac{\sqrt{3} i k_{x}}{2}} & 0 & e^{- i k_{y}} & 0 & e^{\frac{i k_{y}}{2}} e^{\frac{\sqrt{3} i k_{x}}{2}}\\0 & 0 & e^{i k_{y}} & 0 & e^{- \frac{i k_{y}}{2}} e^{- \frac{\sqrt{3} i k_{x}}{2}} & 0\end{matrix}\right]\end{split}\]

Make terms at k=0 and check degeneracy for some linear combination.

Gamma_terms = [term.subs(dict(k_x=0, k_y=0)) for term in family]
evals, evecs = np.linalg.eigh(
    np.sum([
        list(term.values())[0] * coeff
        for term, coeff in zip(Gamma_terms, [1, 2])
    ],
    axis=0)
)

evals

array([-5.19615242e+00, -8.69969719e-18,  0.00000000e+00,  0.00000000e+00,
        0.00000000e+00,  5.19615242e+00])

Find symmetry representation of 3-fold rotations in the low energy subspace

# projector to low energy subspace
proj = evecs[:, 1:5]
UC3 = C3.U

# projected symmetry operator
pC3 = proj.T.conj() @ UC3 @ proj

assert np.allclose(pC3 @ pC3.T.conj(), np.eye(4))

evals2, evecs2 = np.linalg.eig(pC3)

# Check angles of the eigenvalues
np.angle(evals2) / (2 * np.pi)

array([ 0.        ,  0.        ,  0.33333333, -0.33333333])

The symmetry representation in the low energy subspace near the Gamma point is different in the two cases.

Qsymm Reference Documentation

Submodules

qsymm.model module

class qsymm.model.BlochCoeff

Bases: tuple

Container for Bloch coefficient in BlochModel, in the form of (hop, coeff), equivalent to coeff * exp(I * hop.dot(k)).

copy()

tosympy(momenta, nsimplify=False)

class qsymm.model.BlochModel(hamiltonian=None, locals=None, momenta=('k_x', 'k_y', 'k_z'), keep=None, shape=None, format=None)

Bases: qsymm.model.Model

A Model where coefficients are periodic functions of momenta.

Internally it is a sum(BlochCoeff * value), where BlochCoeff is a symbolic representation of coefficients and a periodic function of k. value can be scalar, array (both dense and sparse) or LinearOperator. This is accessible as a dict {BlochCoeff: value}.

Parameters

• hamiltonian (Model, str, SymPy expression, dict or None (default)) – Symbolic representation of a Hamiltonian. If a string, it is converted to a SymPy expression using kwant_continuum.sympify. If a dict is provided, it should have the form {symbol: array} with all arrays the same size (dense or sparse). If symbol is not a BlochCoeff, it is passed through sympy.sympify, and should consist purely of a product of symbolic coefficients, no constant factors other than 1. symbol is then converted to BlochCoeff. None initializes a zero BlochModel.

• locals (dict or None (default)) – Additional namespace entries for sympify. May be used to simplify input of matrices or modify input before proceeding further. For example: locals={'k': 'k_x + I * k_y'} or locals={'sigma_plus': [[0, 2], [0, 0]]}.

• momenta (iterable of strings or Sympy symbols) – Names of momentum variables, default ('k_x', 'k_y', 'k_z') or corresponding sympy symbols. Momenta are treated the same as other keys for the purpose of keep. Ignored when initialized with Model.

• keep (iterable of BlochCoeff (optional)) – Set of symbolic coefficients that are kept, anything that does not appear here is discarded. Useful for perturbative calculations where only terms to a given order are needed. By default all keys are kept. Ignored when initialized with Model.

• shape (tuple or None (default)) – Shape of the Model, must match the shape of all the values. If not provided, it is automatically found based on the shape of the input. Must be provided is hamiltonian is None or {}. Empty tuple corresponds to scalar values. Ignored when initialized with Model.

• format (class or None (default)) – Type of the values in the model. Supported types are np.complex128, np.ndarray, scipy.sparse.spmatrix and scipy.sparse.linalg.LinearOperator. If hamiltonian is provided as a dict, all values must be of this type, except for scalar values, which are recast to np.complex128. If format is not provided, it is inferred from the type of the values. If hamiltonian is not a dictionary, format is ignored and set to np.ndarray or hamiltonian.format if it is a Model.

conj()

Complex conjugation.

rotate_momenta(R)

Rotate momenta with rotation matrix R.

subs(*args, **kwargs)

Substitute symbolic expressions. See Model.subs.

tomodel(nsimplify=False)

Convert to Model.

tosympy(nsimplify=False)

Return sympy representation of the Model. If nsimplify=True, attempt to rewrite numerical coefficients as exact formulas.

transform_symbolic(func)

Transform keys by applying func to all of them. Useful for symbolic substitutions, differentiation, etc.

class qsymm.model.Model(hamiltonian=None, locals=None, momenta=('k_x', 'k_y', 'k_z'), keep=None, symbol_normalizer=None, normalize=False, shape=None, format=None)

Bases: collections.UserDict

Symbolic matrix-valued function that depends on momenta and other parameters.

Implements the algebra of matrix valued functions. Implements many sympy and numpy methods and overrides arithmetic operators. Internally it represents sum(symbol * value), where symbol is a symbolic expression, and value can be scalar, array (both dense and sparse) or LinearOperator. This is accessible as a dict {symbol: value}.

Parameters

• hamiltonian (str, SymPy expression, dict or None (default)) – Symbolic representation of a Hamiltonian. If a string, it is first converted to a SymPy expression using kwant_continuum.sympify. If a dict is provided, it should have the form {symbol: array} with all arrays the same size (dense or sparse). symbol by default is passed through sympy.sympify, and should consist purely of a product of symbolic coefficients, no constant factors other than 1, except if normalize=True. None initializes a zero Model.

• locals (dict or None (default)) – Additional namespace entries for sympify. May be used to simplify input of matrices or modify input before proceeding further. For example: locals={'k': 'k_x + I * k_y'} or locals={'sigma_plus': [[0, 2], [0, 0]]}.

• keep (iterable of expressions (optional)) – Set of symbolic coefficients that are kept, anything that does not appear here is discarded. Useful for perturbative calculations where only terms to a given order are needed. By default all keys are kept.

• momenta (iterable of strings or Sympy symbols) – Names of momentum variables, default ('k_x', 'k_y', 'k_z') or corresponding sympy symbols. Momenta are treated the same as other keys for the purpose of keep.

• symbol_normalizer (callable (optional)) – Function applied to symbols when initializing the internal dict. By default the keys are passed through sympy.sympify and sympy.expand_power_exp. Keys when accessing a term and keys in keep are also passed through symbol_normalizer.

• normalize (bool, default False) – Whether to clean input dict by splitting summands in symbols, moving numerical factors in the symbols to values, removing entries with values allclose to zero. Ignored if hamiltonian is not a dict.

• shape (tuple or None (default)) – Shape of the Model, must match the shape of all the values. If not provided, it is automatically found based on the shape of the input. Must be provided if hamiltonian is None or {}. Empty tuple corresponds to scalar values.

• format (class or None (default)) – Type of the values in the model. Supported types are np.complex128, scipy.sparse.linalg.LinearOperator, np.ndarray, and subclasses of scipy.sparse.spmatrix . If hamiltonian is provided as a dict, all values must be of this type, except for scalar values, which are recast to np.complex128. If format is not provided, it is inferred from the type of the values. Must be provided if hamiltonian is None or {}. If hamiltonian is not a dictionary, format is ignored an set to np.ndarray.

Notes

Sympy symbols are immutable and references to the same symbols is stored in different Models. Be warned that setting any assumptions for symbols (such as real) will result in an identically named, but different symbol, and these are not handled properly. Model assumes that all sympy symbols are real without any assumptions explicitly set.

T()

Transpose

allclose(other, rtol=1e-05, atol=1e-08, equal_nan=False)

Test whether two Models are approximately equal

around(decimals=3)

Return Model with matrices rounded to given number of decimals.

conj()

Complex conjugation

copy()

Return a copy.

eliminate_zeros(rtol=1e-05, atol=1e-08)

Return a model with small terms removed. Tolerances are in Frobenius matrix norm, relative tolerance compares to the value with largest norm.

evalf(subs=None)

Evaluate using parameter values in subs.

lambdify(nsimplify=False, *, onsite=False, hopping=False)

Return a callable object for the model, with sympy symbols as parameters.

Parameters

• nsimplify (bool, default False) – Whether or not to attempt to rewrite numerical coefficients as exact symbols in sympification.

• onsite (bool, default False) – If True, adds ‘site’ as the first argument to the callable object. Helpful for passing Model objects to kwant Builder objects as onsite functions.

• hopping (bool, default False) – If True, adds ‘site1’ and ‘site2’ as the first two arguments to the callable object. Helpful for passing Model objects to kwant Builder objects as hopping functions.

• Notes –

• and hopping are mutually exclusive. If both are set to True, (onsite) –

• error is thrown. (an) –

reshape(*args, **kwargs)

Reshape, see numpy.reshape.

rotate_momenta(R)

Rotate momenta with rotation matrix R.

subs(*args, **kwargs)

Substitute symbolic expressions. See documentation of sympy.Expr.subs() for details.

Allows for the replacement of momenta in the Model object. Replacing a momentum k with a sympy.Symbol object p replaces the momentum k with p in the Model. Replacing a momentum k with a number removes the momentum k from the Model momenta. Replacing a momentum k with a sympy expression that does not contain any of the Model.momenta also removes the momentum k from the momenta.

toarray()

Convert to dense numpy ndarray format.

tocsr()

Convert to sparse csr format.

tosympy(nsimplify=False)

Return sympy representation of the Model. If nsimplify=True, attempt to rewrite numerical coefficients as exact formulas.

trace()

Take trace of the matrix and return a scalar valued Model.

transform_symbolic(func)

Transform keys by applying func to all of them. Useful for symbolic substitutions, differentiation, etc.

value_list(key_list)

Return a list of the matrix coefficients corresponding to the keys in key_list.

zeros_like()

Return an empty model object that inherits the other properties

qsymm.model.copy(a)

qsymm.model.substitute_exponents(expr)

Substitute trignometric functions with exp.

sin(X) -> (e^(I * X) - e^(-I * X)) / (2 * I) cos(X) -> (e^(I * X) + e^(-I * X)) / 2 exp(X) -> e^X

qsymm.symmetry_finder module

qsymm.symmetry_finder.bravais_group_2d(neighbors, num_eq, sets_eq, verbose=False)

qsymm.symmetry_finder.bravais_group_3d(neighbors, num_eq, sets_eq, verbose=False)

qsymm.symmetry_finder.bravais_point_group(periods, tr=False, ph=False, generators=False, verbose=False)

Find the point group of the Bravais-lattice defined by periods.

Parameters

• periods (array) – Translation vectors as row vectors, arranged into a 2D array.

• ph (tr,) – Whether to return time reversal and particle-hole operators. If False, only pure point-group operators are returned.

• generators (bool (default False)) – If True only a (not necessarily minimal) generator set of the symmetry group is returned.

• verbose (bool (default False)) – If True, the name of the Bravais lattice is printed.

Returns

Return type

set of PointGroupElements

qsymm.symmetry_finder.casimir(gens)

qsymm.symmetry_finder.check_bravais_symmetry(neighbors, group)

qsymm.symmetry_finder.conserved_quantities(Ps, prettify=False, num_digits=3)

Construct a full set of conserved quantities from the projectors.

Parameters

• Ps (list fo 3 index ndarrays) – projectors ‘Ps’ returned ‘_reduce_hamiltonian()’

• prettify (bool) – If true it finds a nice sparse basis of the conserved quantities, they are generally returned in a random basis, but any linear combination is also conserved.

Returns

conserved quantities that all commute with the family of Hamiltonians. The identity matrix is excluded.

Return type

list of ContinuousGroupGenerators

qsymm.symmetry_finder.continuous_symmetries(model, Ps=None, prettify=True, num_digits=8, sparse_linalg=False)

Find continuous rotation symmetries of Hamiltonian family represented by model. Hamiltonian is reduced, so on-site continuous symmetries are factored out.

Parameters

• model (Model) – symbolic representation of the Hamiltonian family

• Ps (ndarray, optional) – Projectors as returned by _reduce_hamiltonian

• prettify (bool) – Whether to carry out sparsification of the results, in general an arbitrary linear combination of the symmetry generators is returned.

• num_digits (float) – Absolute precision when prettifying the result.

• sparse_linalg (bool) – Whether to use sparse linear algebra in the calculation. Can give large performance gain in large systems.

Returns

symmetries – List of linearly independent symmetry generators.

Return type

list of ContinuousGroupGenerator

qsymm.symmetry_finder.discrete_symmetries(model, candidates, Ps=None, generators=False, verbose=False, sparse_linalg=False)

Find point group symmetries of Hamiltonians family. Optimized version to reduce number of tests, uses sympy exact rotation matrices

Parameters

• model (Model) – Model which represents family of Hamiltonians

• candidates (set of PointGroupElements) – Set of candidate PointGroupElements. Must have .U attribute set to None.

• Ps (ndarray, optional) – Projectors as returned by _reduce_hamiltonian.

• generators (bool) – If true, only a set of generators are returned, otherwise the full symmetry group is returned.

• sparse_linalg (bool) – Whether to use sparse linear algebra in the calculation. Can give large performance gain in large systems.

• verbose (bool) –

Returns

• genset or symset (set of PointGroupElement) – Symmetries of model.

• ### TODO (remove Ps from return)

• Ps (ndarray) – Projectors as returned by _reduce_hamiltonian.

qsymm.symmetry_finder.equals(vectors)

qsymm.symmetry_finder.killing(gens)

qsymm.symmetry_finder.pick_perp(vectors, n, other_vectors=None)

qsymm.symmetry_finder.separate_lie_algebra(gens)

Separate Lie-algebra defined by the list of linearly independent generators ‘gens’ into its center and semisimple parts.

Parameters

gens (ndarray) – 3 index array of shape (m, N, N), Lie-algebra generators

Returns

• gensc (ndarray) – generators of the center, shape (mc, N, N)

• genss (ndarray) – generators of semisimple part, shape (ms, N, N)

qsymm.symmetry_finder.struct_const(gens)

qsymm.symmetry_finder.symmetries(model, candidates=None, continuous_rotations=False, generators=False, prettify=False, num_digits=8, verbose=False, sparse_linalg=False)

Find symmetries of the Hamiltonian family described by model.

Parameters

• model (Model) – Model or BlochModel which represents family of Hamiltonians. Every symbolic prefactor is treated as a free parameter, and model.momenta as independent momentum variables.

• candidates (iterable of PointGroupElements or None) – Set of candidate PointGroupElements used for finding discrete symmetries. Must have .U attribute set to None. It is advised that candidates forms a group, as combinations of elements in candidates are not tested. If model is a Model that describes a Bloch Hamiltonian, the rotation matrices must be either integer or Sympy matrices, as exact arythmetic is assumed. Floating point rotation matrices are acceptable for k.p Models and BlochModels. If None (default), only discrete onsite symmetries (time reversal, particle-hole, chiral) are found.

• continuous_rotations (bool (default False)) – Whether to search for continuous rotation symmetries.

• generators (bool (default False)) – If True, only a set of generators are returned, otherwise the full discrete symmetry group is returned.

• prettify (bool (default False)) – Whether to carry out sparsification of the continuous symmetry generators, in general an arbitrary linear combination of the symmetry generators is returned.

• num_digits (float) – Absolute precision when deciding whether symmetry leaves Hamiltonian invariant and prettifying the result with prettify=True.

• verbose (bool) – Whether to print additional information.

• sparse_linalg (bool) – Whether to use sparse linear algebra in the calculation. Can give large performance gain in large systems.

Returns

• disc_sym (list of PointGroupElement) – Discrete symmetries of model, always a group generated by a subset of candidates. If generators=False, it is a closed group, if generators=True it is a (not necessarily minimal) generator set, may be empty.

• cont_sym (list of ContinuousGroupGenerator) – List of linearly independent continuous symmetry generators, may be empty. Onsite conserved quantities are listed first, then continuous rotation generators. The trivial conserved quantity proportional to identity is not included.

qsymm.symmetry_finder.symmetry_adapted_sun(gens, check=False)

Find symmetry adapted basis of the simple ‘su(d)’ Lie-algebra representation defined by generators ‘gens’. It is assumed that the representation is the direct sum of ‘n’ identical ‘su(d)’ representations, i.e. in the symmetry adapted basis the generators have the form ‘L otimes 1_{nxn}’ where ‘L’ runs over a generator set of all ‘d*d’ Hermitian matrices.

Parameters

gens (ndarray) – 3 index array of shape ‘(d**2-1, n*d, n*d)’ with ‘n’ and ‘d’ integers, list of ‘su(d)’ generators.

Returns

• Ps (ndarray) – 3 index array of shape ‘(d, n*d, n)’, list of projectors onto the symmetry adapted subspaces. ‘np.einsum(‘aij,ab,bkj->ik’, Ps, L, Ps.conjugate()))’ is a symmetry spanned by ‘gens’ for any ‘d*d’ Hermitian matrix ‘L’. Stacking these produces the unitary transformation ‘U = np.hstack(Ps)’ to the symmetry adapted basis.

• check (bool) – Whether to check the final result.

qsymm.symmetry_finder.threefold_axis(vectors, neighbors)

qsymm.symmetry_finder.twofold_axis(vectors, neighbors)

qsymm.hamiltonian_generator module

qsymm.hamiltonian_generator.bloch_family(hopping_vectors, symmetries, norbs, onsites=True, momenta=[k_x, k_y, k_z], symmetrize=True, prettify=True, num_digits=10, bloch_model=False)

Generate a family of symmetric Bloch-Hamiltonians.

Parameters

• hopping_vectors (list of tuples (a, b, vec)) – a and b are identifiers for the different sites (e.g. strings) of the unit cell, vec is the real space hopping vector. Vec may contain contain integers, sympy symbols, or floating point numbers.

• symmetries (list of PointGroupElement or ContinuousGroupGenerator) – Generators of the symmetry group. ContinuousGroupGenerators can only have onsite action as a lattice system cannot have continuous rotation invariance. It is assumed that the block structure of the unitary action is consistent with norbs, as returned by symmetry_from_permutation.

• norbs (OrderedDict : {site : norbs_site} or tuple of tuples ((site, norbs_site), )) – sites are ordered in the order specified, with blocks of size norbs_site corresponding to each site.

• onsites (bool, default True) – Whether to include on-site terms consistent with the symmetry.

• momenta (iterable of strings or Sympy symbols) – Names of momentum variables, default ['k_x', 'k_y', 'k_z'] or corresponding sympy symbols.

• symmetrize (bool, default True) – Whether to use the symmetrization strategy. This does not require a full set of hoppings to start, all symmetry related hoppings are generated. Otherwise the constraining strategy is used, this does not generate any new hoppings beyond the ones specified and constrains coefficients to enforce symmetry.

• prettify (bool) – Whether to prettify the result. This step may be numerically unstable.

• num_digits (int, default 10) – Number of significant digits kept when prettifying.

• bloch_model (bool, default False) – Determines the return format of this function. If set to False, returns a list of Model objects. If True, returns a list of BlochModel objects. BlochModel objects are more suitable than Model objects if the hopping vectors include floating point numbers.

Returns

• family (list of Model or BlochModel objects) – A list of Model or BlochModel objects representing the family that satisfies the symmetries specified. Each object satisfies all the symmetries by construction.

• Notes

• ——

• There is no need to specify the translation vectors, all necessary information

• is encoded in the symmetries and hopping vectors.

• In the generic case the Bloch-Hamiltonian produced is not Brillouin-zone periodic,

• instead it acquires a unitary transformation W_G = delta_{ab} exp(i G (r_a - r_b))

• where G is a reciprocal lattice vector, a and b are sites and r_a is the

• real space position of site a. If the lattice is primitive (there is only one

• site per unit cell), the hopping vectors are also translation vectors and the

• resulting Hamiltonian is BZ periodic.

• If symmetrize=True, all onsite unitary symmetries need to be explicitely

• specified as ContinuousGroupGenerators. Onsite PointGroupSymmetries (ones

• with R=identity) are ignored.

• If floating point numbers are used in the argument hopping_vectors, it is

• recommended to have this function return BlochModel objects instead of Model

• objects, by setting the bloch_model flag to True.

qsymm.hamiltonian_generator.check_symmetry(family, symmetries, num_digits=None)

Check that a family satisfies symmetries. A symmetry is satisfied if all members of the family satisfy it.

If the input family was rounded before hand, it is necessary to use specify the number of significant digits using num_digits, otherwise this check might fail.

Parameters

• family (iterable of Model or BlochModel objects representing) – a family.

• symmetries (iterable representing the symmetries to check.) – If the family is a Hamiltonian family, symmetries is an iterable of PointGroupElement objects representing the symmetries to check.

• num_digits (integer) – In the case that the input family has been rounded, num_digits should be the number of significant digits to which the family was rounded.

qsymm.hamiltonian_generator.constrain_family(symmetries, family, sparse_linalg=False)

Apply symmetry constraints to a family.

Parameters

• symmetries (iterable of PointGroupElement objects, representing the symmetries) – that are used to constrain the Hamiltonian family.

• family (iterable of Model or BlochModel objects, representing the Hamiltonian) – family to which the symmetry constraints are applied.

• sparse_linalg (bool) – Whether to use sparse linear algebra. Using sparse solver can result in performance increase for large, highly constrained families, but not as well tested as the default dense version.

Returns

family – family with the symmetry constraints applied.

Return type

iterable of Model or BlochModel objects, that represents the

qsymm.hamiltonian_generator.continuum_hamiltonian(symmetries, dim, total_power, all_powers=None, momenta=[k_x, k_y, k_z], sparse_linalg=False, prettify=False, num_digits=10)

Generate a family of continuum Hamiltonians that satisfy symmetry constraints.

Parameters

• symmetries (iterable of PointGroupElement objects.) – An iterable of PointGroupElement objects, each describing a symmetry that the family should possess.

• dim (integer) – The number of spatial dimensions along which the Hamiltonian family is translationally invariant. Only the first dim entries in all_powers and momenta are used.

• total_power (integer or list of integers) – Allowed total powers of the momentum variables in the continuum Hamiltonian. If an integer is given, all powers below it are included as well.

• all_powers (list of integer or list of list of integers) – Allowed powers of the momentum variables in the continuum Hamiltonian listed for each spatial direction. If an integer is given, all powers below it are included as well. If a list of integers is given, only these powers are used.

• momenta (iterable of strings or Sympy symbols) – Names of momentum variables, default ['k_x', 'k_y', 'k_z'] or corresponding sympy symbols.

• sparse_linalg (bool) – Whether to use sparse linear algebra. Using sparse solver can result in performance increase for large, highly constrained families, but not as well tested as the default dense version.

• prettify (boolean, default False) – Whether to make the basis pretty by rounding and basis change. For details see docstring of make_basis_pretty. May be numerically unstable.

• num_digits (integer, default 10) – Number of significant digits to which the basis is rounded using prettify. This argument is ignored if prettify = False.

Returns

family – A list of Model objects representing the family that satisfies the symmetries specified. Each Model object satisfies all the symmetries by construction.

Return type

list

qsymm.hamiltonian_generator.continuum_pairing(symmetries, dim, total_power, all_powers=None, momenta=[k_x, k_y, k_z], phases=None, ph_square=1, sparse_linalg=False, prettify=False, num_digits=10)

Generate a family of continuum superconducting pairing functions that satisfy symmetry constraints.

The specified symmetry operators should act on the normal state Hamiltonian, not in particle-hole space.

Parameters

• symmetries (iterable of PointGroupElement objects.) – An iterable of PointGroupElement objects, each describing a symmetry that the family should possess.

• dim (integer) – The number of spatial dimensions along which the Hamiltonian family is translationally invariant. Only the first dim entries in all_powers and momenta are used.

• total_power (integer or list of integers) – Allowed total powers of the momentum variables in the continuum Hamiltonian. If an integer is given, all powers below it are included as well.

• all_powers (list of integer or list of list of integers) – Allowed powers of the momentum variables in the continuum Hamiltonian listed for each spatial direction. If an integer is given, all powers below it are included as well. If a list of integers is given, only these powers are used.

• momenta (list of int or list of Sympy objects) – Indices of momenta from [‘k_x’, ‘k_y’, ‘k_z’] or a list of names for the momentum variables. Default is [‘k_x’, ‘k_y’, ‘k_z’].

• phases (iterable of numbers) – Phase factors to multiply the hole block of the symmetry operators in particle-hole space. By default, all phase factors are 1.

• ph_square (integer, either 1 or -1.) – Specifies whether the particle-hole operator squares to +1 or -1.

• sparse_linalg (bool) – Whether to use sparse linear algebra. Using sparse solver can result in performance increase for large, highly constrained families, but not as well tested as the default dense version.

• prettify (boolean, default False) – Whether to make the basis pretty by rounding and basis change. For details see docstring of make_basis_pretty. May be numerically unstable.

• num_digits (integer, default 10) – Number of significant digits to which the basis is rounded using prettify. This argument is ignored if prettify = False.

Returns

family – A list of Model objects representing the family that satisfies the symmetries specified. Each Model object satisfies all the symmetries by construction.

Return type

list

qsymm.hamiltonian_generator.continuum_variables(dim, total_power, all_powers=None, momenta=[k_x, k_y, k_z])

Make a list of all linearly independent combinations of momentum variables with given total power.

Parameters

• dim (integer) – The number of spatial dimensions along which the Hamiltonian family is translationally invariant. Only the first dim entries in all_powers and momenta are used.

• total_power (integer) – Allowed total power of the momentum variables in the continuum Hamiltonian.

• all_powers (list of integer or list of list of integers) – Allowed powers of the momentum variables in the continuum Hamiltonian listed for each spatial direction. If an integer is given, all powers below it are included as well. If a list of integers is given, only these powers are used.

• momenta (list of int or list of Sympy objects) – Indices of momenta from [‘k_x’, ‘k_y’, ‘k_z’] or a list of names for the momentum variables. Default is [‘k_x’, ‘k_y’, ‘k_z’].

Returns

Return type

A list of Sympy objects, representing the momentum variables that enter the Hamiltonian.

qsymm.hamiltonian_generator.display_family(family, summed=False, coeffs=None, nsimplify=True)

Display a Hamiltonian family in a Jupyter notebook

If this function is used from a Jupyter notebook then it uses the notebook’s rich LaTeX display features. If used from a console or script, then this function just uses print().

Parameters

• family (iterable of Model or BlochModel objects) – List of terms in a Hamiltonian family.

• summed (boolean) – Whether to display the Hamiltonian family by individual member (False), or as a sum over all members with expansion coefficients (True).

• coeffs (list of sympy objects, optional) – Coefficients used when combining terms in the family if summed is True.

• nsimplify (boolean) – Whether to use sympy.nsimplify on the output or not, which attempts to replace floating point numbers with simpler expressions, e.g. fractions.

qsymm.hamiltonian_generator.hamiltonian_from_family(family, coeffs=None, nsimplify=True, tosympy=True)

Form a Hamiltonian from a Hamiltonian family by taking a linear combination of its elements.

Parameters

• family (iterable of Model or BlochModel objects) – List of terms in the Hamiltonian family.

• coeffs (list of sympy objects, optional) – Coefficients used to form the linear combination of terms in the family. Element n of coeffs multiplies member n of family. The default choice of the coefficients is c_n.

• nsimplify (bool) – Whether to use sympy.nsimplify on the output or not, which attempts to replace floating point numbers with simpler expressions, e.g. fractions.

• tosympy (bool) – Whether to convert the Hamiltonian to a sympy expression. If False, a Model or BlochModel object is returned instead, depending on the type of the Hamiltonian family.

Returns

ham – The Hamiltonian, i.e. the linear combination of entries in family.

Return type

sympy.Matrix or Model/BlochModel object.

qsymm.hamiltonian_generator.make_basis_pretty(family, num_digits=2)

Attempt to make a family more legible by reducing the number of nonzero entries in the matrix coefficients.

Parameters

• family (iterable of Model or BlochModel objects representing) – a family.

• num_digits (positive integer) – Number of significant digits that matrix coefficients are rounded to.

• attempts to make the family more legible by prettifying a matrix, (This) –

• is done by bringing it to reduced row-echelon form. This (which) –

• may be numerically unstable, so this function should be used (procedure) –

• caution. (with) –

qsymm.hamiltonian_generator.remove_duplicates(family, tol=1e-08)

Remove linearly dependent terms in Hamiltonian family using SVD.

Parameters

• family (iterable of Model or BlochModel objects representing) – a family.

• tol (float) – tolerance used in SVD when finding the span.

Returns

rfamily – the family with only linearly independent terms.

Return type

list of Model or BlochModel objects representing

qsymm.hamiltonian_generator.round_family(family, num_digits=3)

Round the matrix coefficients of a family to specified significant digits.

Parameters

• family (iterable of Model objects that represents) – a family.

• num_digits (integer) – Number if significant digits to which the matrix coefficients are rounded.

Returns

• A list of Model objects representing the family, with

• the matrix coefficients rounded to num_digits significant digits.

qsymm.hamiltonian_generator.subtract_family(family1, family2, tol=1e-08, prettify=False)

Remove the linear span of family2 from the span of family1 using SVD. family2 must be a span of terms that are either inside the span of family1 or orthogonal to it. This guarantees that projecting out family2 from family1 results in a subfamily of family1.

Parameters

• family2 (family1,) – Hamiltonian families.

• tol (float) – tolerance used in SVD when finding the span.

Returns

rfamily – family1 with the span of family2 removed.

Return type

list of Model or BlochModel objects representing

qsymm.hamiltonian_generator.symmetrize_monomial(monomial, symmetries)

Symmetrize monomial by averaging over all symmetry images under symmetries.

Parameters

• monomial (Model or BlochModel object) – Hamiltonian term to be symmetrized

• symmetries (iterable of PointGroupElement objects) – Symmetries to use for symmetrization. symmetries must form a closed group.

Returns

Symmetrized term.

Return type

Model or BlochModel object

qsymm.groups module

class qsymm.groups.ContinuousGroupGenerator(R=None, U=None)

Bases: object

A generator of a continuous group.

Generates a family of symmetry operators that act on the Hamiltonian as:

\[H(k) → \exp{-iλU} H(\exp{iλR} k) \exp{iλU}\]

with λ a real parameter.

Parameters

• R (ndarray (optional)) – Real space rotation generator: Hermitian antisymmetric. If not provided, the zero matrix is used.

• U (ndarray (optional)) – Hilbert space unitary rotation generator: Hermitian. If not provided, the zero matrix is used.

R

U

apply(model)

Return copy of model H(k) with applied infinitesimal generator. 1j * (H(k) U - U H(k)) + 1j * dH(k)/dk_i R_{ij} k_j

qsymm.groups.L_matrices(d=3, l=1)

Construct real space rotation generator matrices in d=2 or 3 dimensions. Can also be used to get angular momentum operators for real atomic orbitals in 3 dimensions, for p-orbitals use l=1, for d-orbitals l=2. The basis of p-orbitals is p_x, p_y, p_z, for d-orbitals d_{x^2 - y^2}, d_{3 z^2 - r^2}, d_{xy}, d_{yz}, d_{zx}. The matrices are all purely imaginary and antisymmetric. To generate finite rotations, use ‘spin_rotation(n, L)’.

class qsymm.groups.PointGroupElement(R, conjugate=False, antisymmetry=False, U=None, _strict_eq=False)

Bases: object

An element of a point group.

Parameters

• R (sympy.ImmutableMatrix or array) – Real space rotation action of the operator. Square matrix with size of the number of spatial dimensions.

• conjugate (boolean (default False)) – Whether the operation includes conplex conjugation (antiunitary operator)

• antisymmetry (boolean (default False)) – Whether the operator flips the sign of the Hamiltonian (antisymmetry)

• U (ndarray (optional)) – The unitary action on the Hilbert space. May be None, to be able to treat symmetry candidates

• _strict_eq (boolean (default False)) – Whether to test the equality of the unitary parts when comparing with other PointGroupElements. By default the unitary parts are ignored. If True, PointGroupElements are considered equal, if the unitary parts are proportional, an overall phase difference is still allowed.

Notes

As U is floating point and has a phase ambiguity at least, hence it is ignored when comparing objects by default.

R is the real space rotation acion. Do not include minus sign for the k-space action of antiunitary operators, such as time reversal. This minus sign will be included automatically if ‘conjugate=True’.

For most uses R can be provided as a floating point array. It is necessary to use exact sympy matrix representation if the PGE has to act on Models with complicated momentum dependence (not polynomial), as the function parts of models are compared exactly. If the momentum dependence is periodic (sine, cosine and exponential), use BlochModel, this works with floating point rotations.

R

U

antisymmetry

apply(model)

Return copy of model with applied symmetry operation.

if unitary: (+/-) U H(inv(R) k) U^dagger if antiunitary: (+/-) U H(- inv(R) k).conj() U^dagger

(+/-) stands for (symmetry / antisymmetry)

If self.U is None, U is taken as the identity.

conjugate

identity()

Return identity element with the same structure as self.

inv()

Invert PointGroupElement

class qsymm.groups.PrettyList

Bases: list

Subclass of list that displays its elements with latex.

qsymm.groups.chiral(realspace_dim, U=None)

Return a chiral symmetry operator

Parameters

• realspace_dim (int) – Realspace dimension

• U (ndarray (optional)) – The unitary action on the Hilbert space. May be None, to be able to treat symmetry candidates

Returns

P

Return type

PointGroupElement

qsymm.groups.cubic(tr=True, ph=True, generators=False, spin=None)

Generate cubic point group in standard basis.

Parameters

• ph (tr,) – Whether to include time-reversal and particle-hole symmetry.

• generators (bool (default false)) – Only return the group generators if True.

• spin (float or sequence of arrays (optional)) – Spin representation to use for the unitary action of the operator. If not provided, the PointGroupElements have the unitary action set to None. If float is provided, it should be integer or half-integer specifying the spin representation in the standard basis, see spin_matrices. Otherwise a sequence of 3 arrays of identical square size must be provided representing 3 components of the angular momentum operator. The unitary action of rotation operator is U = exp(-i n⋅s). If ph is True, spin may not be provided, as it is not possible to deduce the unitary representation of particle-hole symmetry from spin alone. In this case construct the particle-hole operator manually.

Returns

• set of PointGroupElement objects with integer rotations

• Notes

• —— – We assume inversion acts trivially in spin space.

qsymm.groups.generate_group(gens)

Generate group from gens

Parameters

gens (iterable of PointGroupElement) – generator set of the group

Returns

group – group generated by gens, closed under multiplication

Return type

set of PointGroupElement

qsymm.groups.generate_subgroups(group)

Generate all subgroups of group, including the trivial group and itself.

Parameters

group (set of PointGroupElement) – A closed group as set of its elements.

Returns

subgroups – frozesets are the subgroups, sets are a generator set of the subgroup.

Return type

dict of forzenset: set

qsymm.groups.hexagonal(dim=2, tr=True, ph=True, generators=False, sympy_R=True, spin=None)

Generate hexagonal point group in standard basis in 2 or 3 dimensions. Mirror symmetries with the main coordinate axes as normals are included. In 3D the hexagonal axis is the z axis.

Parameters

• dim (int (default 2)) – Real sapce dimensionality, 2 or 3.

• ph (tr,) – Whether to include time-reversal and particle-hole symmetry.

• generators (bool (default True)) – Only return the group generators if True.

• sympy_R (bool (default True)) – Whether the rotation matrices should be exact sympy representations.

• spin (float or sequence of arrays (optional)) – Spin representation to use for the unitary action of the operator. If not provided, the PointGroupElements have the unitary action set to None. If float is provided, it should be integer or half-integer specifying the spin representation in the standard basis, see spin_matrices. Otherwise a sequence of 3 arrays of identical square size must be provided representing 3 components of the angular momentum operator. The unitary action of rotation operator is U = exp(-i n⋅s). In 2D the z axis is assumed to be the rotation axis. If ph is True, spin may not be provided, as it is not possible to deduce the unitary representation of particle-hole symmetry from spin alone. In this case construct the particle-hole operator manually.

Returns

• set of PointGroupElements

• Notes

• —— – Warning: in 2D the real space action of a mirror and and a 2-fold rotation around an axis in the plane is identical, however the action on angular momentum is different. Here we consider the action of the mirror, which is the same as the action of a 2-fold rotation around the mirror axis, assuming inversion acts trivially.

qsymm.groups.identity(dim, shape=None)

Return identity operator with appropriate shape.

Parameters

• dim (int) – Dimension of real space.

• shape (int (optional)) – Size of the unitary part of the operator. If not provided, U is set to None.

Returns

id

Return type

PointGroupElement

qsymm.groups.inversion(realspace_dim, U=None)

Return an inversion operator

Parameters

• realspace_dim (int) – Realspace dimension

• U (ndarray (optional)) – The unitary action on the Hilbert space. May be None, to be able to treat symmetry candidates

Returns

P

Return type

PointGroupElement

qsymm.groups.is_sympy_matrix(R)

qsymm.groups.mirror(axis, U=None, spin=None)

Return a mirror operator

Parameters

• axis (ndarray) – Normal of the mirror. The dimensionality of the operator is the same as the length of axis.

• U (ndarray (optional)) – The unitary action on the Hilbert space. May be None, to be able to treat symmetry candidates

• spin (float or sequence of arrays (optional)) – Spin representation to use for the unitary action of the operator. If float is provided, it should be integer or half-integer specifying the spin representation in the standard basis, see spin_matrices. Otherwise a sequence of 3 arrays of identical square size must be provided representing 3 components of the angular momentum operator. The unitary action of mirror operator is U = exp(-i π n⋅s) where n is normalized to 1. In 2D the axis is treated as x and y coordinates. Only one of U and spin may be provided.

Returns

• P (PointGroupElement)

• Notes

• —— – Warning: in 2D the real space action of a mirror and and a 2-fold rotation around an axis in the plane is identical, however the action on angular momentum is different. Here we consider the action of the mirror, which is the same as the action of a 2-fold rotation around the mirror axis.

qsymm.groups.particle_hole(realspace_dim, U=None)

Return a particle-hole symmetry operator

Parameters

• realspace_dim (int) – Realspace dimension

• U (ndarray (optional)) – The unitary action on the Hilbert space. May be None, to be able to treat symmetry candidates

Returns

P

Return type

PointGroupElement

qsymm.groups.pretty_print_cgg(g, latex=False)

Return a human readable string representation of ContinuousGroupGenerator

Parameters

• g (ContinuousGroupGenerator) – Point group element to be represented.

• latex (bool (default False)) – Whether to output LateX formatted string.

Returns

name – The first line is the action on the Hamiltonian, the following lines display the real space rotation as R(ϕ) for 2D rotation and R(ϕ, axis) for 3D rotation (axis is not normalized), and the conserved quantity as a matrix. If either of these is trivial, it is omitted. Note that a ContinuousGroupGenerator defines a continuous group of symmetries, so there is always a free parameter ϕ.

Return type

string

qsymm.groups.pretty_print_pge(g, full=False, latex=False)

Return a human readable string representation of PointGroupElement

Parameters

• g (PointGroupElement) – Point group element to be represented.

• full (bool (default False)) – Whether to return a full representation. The default short representation only contains the real space action and the symbol of the Altland-Zirnbauer part of the symmetry (see below). The full representation presents the symmetry action on the Hamiltonian and the unitary Hilbert-space action if set.

• latex (bool (default False)) – Whether to output LateX formatted string.

Returns

name – In the short representation it is a sting rot_name + az_name. In the long representation the first line is the action on the Hamiltonian, the second line is rot_name and the third line is the unitary action as a matrix, if set.

rot_name can be: - 1 for identity - I for inversion (in 1D mirror is the same as inversion) - R(angle) for 2D rotation - R(angle, axis) for 3D rotation (axis is not normalized) - M(normal) for mirror - S(angle, axis) for 3D rotoinversion (axis is not normalized)

az_name can be: - T for time-reversal (antiunitary symmetry) - P for particle-hole (antiunitary antisymmetry) - C for chiral (unitary antisymmetry) - missing if the symmetry is unitary

Return type

string

qsymm.groups.rotation(angle, axis=None, inversion=False, U=None, spin=None)

Return a rotation operator

Parameters

• angle (float) – Rotation angle in units of 2 pi.

• axis (ndarray or None (default)) – Rotation axis, optional. If not provided, a 2D rotation is generated around the axis normal to the plane. If a 3D vector is provided, a 3D rotation is generated around this axis. Does not need to be normalized to 1.

• inversion (bool (default False)) – Whether to generate a rotoinversion. By default a proper rotation is returned. Only valid in 3D.

• U (ndarray (optional)) – The unitary action on the Hilbert space. May be None, to be able to treat symmetry candidates

• spin (float or sequence of arrays (optional)) – Spin representation to use for the unitary action of the operator. If float is provided, it should be integer or half-integer specifying the spin representation in the standard basis, see spin_matrices. Otherwise a sequence of 3 arrays of identical square size must be provided representing 3 components of the angular momentum operator. The unitary action of rotation operator is U = exp(-i n⋅s). In 2D the z axis is assumed to be the rotation axis. Only one of U and spin may be provided.

Returns

P

Return type

PointGroupElement

qsymm.groups.rotation_to_angle(R)

qsymm.groups.set_multiply(G, H)

qsymm.groups.spin_matrices(s, include_0=False)

Construct spin-s matrices for any half-integer spin.

Parameters

• s (float or int) – Spin representation to use, must be integer or half-integer.

• include_0 (bool (default False)) – If include_0 is True, S[0] is the identity, indices 1, 2, 3 correspond to x, y, z. Otherwise indices 0, 1, 2 are x, y, z.

Returns

Sequence of spin-s operators in the standard spin-z basis. Array of shape (3, 2*s + 1, 2*s + 1), or if include_0 is True (4, 2*s + 1, 2*s + 1).

Return type

ndarray

qsymm.groups.spin_rotation(n, s, roundint=False)

Construct the unitary spin rotation matrix for rotation specified by the vector n (in radian units) with angular momentum s, given by U = exp(-i n⋅s).

Parameters

• n (iterable) – Rotation vector. Its norm is the rotation angle in radians.

• s (float or sequence of arrays) – Spin representation to use for the unitary action of the operator. If float is provided, it should be integer or half-integer specifying the spin representation in the standard basis, see spin_matrices. Otherwise a sequence of 3 arrays of identical square size must be provided representing 3 components of the angular momentum operator.

• roundint (bool (default False)) – If roundint is True, result is converted to integer tinyarray if possible.

Returns

U – Unitary spin rotation matrix of the same shape as the spin matrices or (2*s + 1, 2*s + 1).

Return type

ndarray

qsymm.groups.square(tr=True, ph=True, generators=False, spin=None)

Generate square point group in standard basis.

Parameters

• ph (tr,) – Whether to include time-reversal and particle-hole symmetry.

• generators (bool (default false)) – Only return the group generators if True.

• spin (float or sequence of arrays (optional)) – Spin representation to use for the unitary action of the operator. If not provided, the PointGroupElements have the unitary action set to None. If float is provided, it should be integer or half-integer specifying the spin representation in the standard basis, see spin_matrices. Otherwise a sequence of 3 arrays of identical square size must be provided representing 3 components of the angular momentum operator. The unitary action of rotation operator is U = exp(-i n⋅s). In 2D the z axis is assumed to be the rotation axis. If ph is True, spin may not be provided, as it is not possible to deduce the unitary representation of particle-hole symmetry from spin alone. In this case construct the particle-hole operator manually.

Returns

• set of PointGroupElement objects with integer rotations

• Notes

• —— – Warning: in 2D the real space action of a mirror and and a 2-fold rotation around an axis in the plane is identical, however the action on angular momentum is different. Here we consider the action of the mirror, which is the same as the action of a 2-fold rotation around the mirror axis, assuming inversion acts trivially.

qsymm.groups.symmetry_from_permutation(R, perm, norbs, onsite_action=None, antiunitary=False, antisymmetry=False)

Construct symmetry operator for lattice systems with multiple sites.

Parameters

• R (real space rotation) –

• perm (dict : {site: image_site}) – permutation of the sites under the symmetry action

• norbs (OrderedDict : {site : norbs_site} or tuple of tuples ((site, norbs_site), )) – sites are ordered in the order specified, with blocks of size norbs_site corresponding to each site.

• onsite_action (dict : {site: ndarray} or ndarray or None) – onsite symmetry action, such as spin rotation for each site. If only one array is specified, it is used for every site. If None (default), identity is used on every site. Size of the arrays must be consistent with norbs.

• antisymmetry (antiunitary,) –

Returns

• g (PointGroupElement) – PointGroupElement corresponding to the operation.

• Notes

• ——

• Sites can be indexed by any hashable identifiers, such as integers or stings.

qsymm.groups.time_reversal(realspace_dim, U=None, spin=None)

Return a time-reversal symmetry operator

Parameters

• realspace_dim (int) – Realspace dimension

• U (ndarray (optional)) – The unitary action on the Hilbert space. May be None, to be able to treat symmetry candidates.

• spin (float or sequence of arrays (optional)) – Spin representation to use for the unitary action of the time reversal operator. If float is provided, it should be integer or half-integer specifying the spin representation in the standard basis, see spin_matrices. Otherwise a sequence of 3 arrays of identical square size must be provided representing 3 components of the angular momentum operator. The unitary action of time-reversal operator is U = exp(-i π s_y). Only one of U and spin may be provided.

Returns

T

Return type

PointGroupElement

qsymm.linalg module

qsymm.linalg.allclose(a, b, rtol=1e-05, atol=1e-08, equal_nan=False)

Alternative to numpy.allclose to test that two ndarrays are elementwise close. Unlike the numpy version, the relative tolerance is not elementwise, but it is relative to the largest entry in the arrays.

qsymm.linalg.commutator(A, B)

qsymm.linalg.family_to_vectors(family, all_keys=None)

Convert a list of Model to standard vector representation.

Parameters

• family (iterable of Model objects) – The Model to convert to vector representation.

• all_keys (iterable or None) – Iterable that must contain all the keys in the family. An optional parameter that allows to fix the ordering of the returned vectors to the same as in all_keys. Otherwise, the keys are ordered arbitrarily (all_keys = None).

Returns

• A matrix representing the family, in which the family members form the

• rows.

qsymm.linalg.matrix_basis(dim, traceless=False, antihermitian=False, real=False, sparse=False)

“Construct a basis for the vector space of dim x dim matrices, that may in addition be traceless, and/or composed of Hermitian or skew-Hermitian basis elements.

Parameters

• dim (positive integer) – The dimension of the matrices (dim x dim).

• traceless (boolean) – Whether the vector space only includes traceless matrices (True) or not (False).

• antihermitian (boolean) – Whether to use a basis composed of Hermitian matrices (False), or skew-Hermitian matrices (True).

• real – Whether the basis only spans real matrices. Depending on ‘antihermitian’ it spans symmetric (False) or antisymmetric (True) matrices only.

• sparse – If True, scipy.sparse.csr_matrix is returned, otherwise np.ndarray.

Returns

basis – A generator that returns matrices that span the vector space.

Return type

generator

qsymm.linalg.mtm(a, B, c)

qsymm.linalg.nullspace(A, atol=1e-06, return_complement=False, sparse=None, k_max=-10)

Compute an approximate basis for the nullspace of matrix A.

Parameters

• A (ndarray or spmatrix) – A should be 2-D.

• atol (real) – Absolute tolerance when deciding whether an eigen/singular value is zero, so the corresponding vector is included in the null space.

• return_complement (bool) – Whether to return the basis of the orthocomplement ot the null space as well. Treated as False if sparse is True.

• sparse (bool or None (default)) – Whether to use sparse linear algebra to find the null space. The default value is inferred from the type of the input array.

• k_max (int (default -10)) – If k_max<0, the full null-space is found by increasing the number of eigenvectors found in each step by abs(k_max). If k_max>0, k_max is the maximum number of null space vectors requested. If the null space dimension is higher than k_max, basis of a random subspace is returned. Can be used to increase performance if the maximum null space dimensionality is known. Ignored unless sparse is True.

Returns

• ns (ndarray) – Basis of the null space of A as row vectors.

• nsc (ndarray) – Basis of the complement of the null space as row vectors, the othonormal basis of the row span of A. Only returned if return_complement=True and sparse=False.

qsymm.linalg.prop_to_id(A)

qsymm.linalg.rref(m, rtol=0.001, return_S=False)

Bring a matrix to reduced row echelon form.

Parameters

• m (2D numpy.ndarray) – The matrix on which to perform row reduction.

• rtol (float) – The tolerance relative to the largest matrix element for what is considered a zero entry.

• return_S (bool) – Whether to return the matrix of linear operations that brings the input matrix to reduced row echelon form.

Returns

• red (2D numpy.ndarray) – The input matrix in reduced row echelon form.

• S (2D numpy.ndarray) – The matrix of operations that brings the input matrix to reduced row echelon form. Only returned if return_S = True.

• Note (If r is the reduced row echelon form of the input matrix m, then) – m = S.dot(r).

qsymm.linalg.simult_diag(mats, tol=1e-06, checks=0)

Simultaneously diagonalize commuting normal square matrices. Recursive algorithm, first diagonalize first matrix, express the rest of the matrices in its eigenbasis, and repeat the procedure for the rest of the matrices in each of the degenerate eigensubspaces. Recursion depth is limited by the number of nontrivial nested eigensubspaces, which is at most the number of matrices or the size of the matrices.

Parameters

• mats (ndarray) – List of commuting matrices, 3 index array.

• tol (float) – Tolerance when deciding degeneracy of eigenvalues.

• checks (int) – Amount of checks to do. 0: no checks 1: check the final result 2: check all the intermediate results

Returns

• Ps (list of ndarrays) – List of rectangular matrices of common eigenvectors spanning each of the common eigensubspaces. Stacking them gives unitary matrix ‘U = np.hstack(Ps)’ such that ‘U.T.conj() @ mats[i] @ U’ is diagonal for every ‘i’.

• Notes

• It is recommended to order matrices such that those with eigenvalues that

• are far spaced or exactly degenerate (such as symmetries) are first, and

• those which may have accidental near degeneracies (such as Hamiltonians)

• are last.

• Tolerance is divided by the size of the block when deciding which eigenvalues

• are degenerate. The default value seems to work well, but not extensively

• tested, numerical instabilities are possible.

qsymm.linalg.solve_mat_eqn(HL, HR=None, hermitian=False, traceless=False, conjugate=False, sparse=None, k_max=-10)

Solve for X the simultaneous matrix equatioins X HL[i] = HR[i] X for every i. It is mapped to a system of linear equations, the null space of which gives a basis for all sulutions.

Parameters

• HL (ndarray or list of arrays or sparse matrices) – Coefficient matrices of identical square shape, list of arrays of shape (n, n) or one array of shape (m, n, n).

• HR (ndarray or list of ndarrays or None) – Same as HL. If HR = None, HR = HL is used.

• traceless (hermitian,) – If X has to be hermitian, use hermitian = True. If X has to be traceless, use traceless = True.

• conjugate (boolean or list of booleans) – If True, solve X HL[i] = HR[i] X^* instead. If a list with the same length as HL and HR is provided, conjugation is applied to the equations with index corresponding to the True entries.

• sparse (bool or None (default)) – Whether to use sparse linear algebra to find the null space. The default value is inferred from the type of the input array.

• k_max (int) – If k_max<0, all solutions are found by increasing the number of soulutions sought in each step by abs(k_max). If k_max>0, k_max is the maximum number of solutions requested. If the solution space dimension is higher than k_max, basis of a random subspace is returned. Can be used to increase performance if the maximum number of solutions is known. Ignored unless sparse is True.

Returns

List of linearly independent square matrices that span all solutions of the eaquations.

Return type

ndarray of shape (l, n, n), or list of csr_matrix

qsymm.linalg.sparse_basis(bas, num_digits=3, reals=False)

Reduce the number of nonzero entries in a (basis) matrix by rounding the matrix and bringing it to reduced row echelon form.

Parameters

• bas (numpy.ndarray) – The (basis) matrix to prettify and round.

• num_digits (positive integer) – Number of significant digits to which to round.

• reals (boolean) – If True, the real and imaginary parts of the matrix are treated separately.

Returns

• A rounded, sparsified matrix with the same row span as the input matrix.

• This is an attempt at matrix sparsification, i.e. we attempt to construct

• the sparsest matrix that has the same row rank. Note that the reduction to

• row echelon form is numerically unstable, so this function should be used

• with caution.

qsymm.linalg.split_list(vals, tol=1e-06)