Qsymm workflow

Qsymm is a great tool, developed here in the group, but not everyone is very familiar with it. This is a Qsymm tutorial designed to change that. Qsymm already has extensive documentation and a few tutorials on read the docs, but this tutorial shows a start-to-finish workflow. It covers using the nonhermitian branch of qsymm to create real-space Hamiltonians, saving and loading models, how to test our models for various symmetries, and how to use our models in kwant systems.

import qsymm
import numpy as np

In order to use the nonhermitian branch of qsymm, clone that branch into your environment, navigate into the directory and do pip install .

Making a real-space model for kwant to use

We want to make a hamiltonian matrix that we can supply to hopping and onsite functions when constructing our system in kwant. We want this matrix to obey certain symmetries that we're interested in. This hamiltonian will be defined in real space, and will vary based on the bond length between sites d rather than the wave vector k.

Step 1: defining symmetries and generating a model

For nonhermitian models, we need to modify the actions of our symmetries accordingly.

  • particle-hole and time-reversal symmetries send $H(k)$ to $-H^\ast(-k)$ or $H^\ast(-k)$ respectively. In real space, we write them such that we end up with $-H^\ast(k)$ and $H^\ast(k)$ respectively.
  • we have to supply a modified hermiticity condition. Here reversing the bond means hopping in the opposite direction, so it is equal to the hermitian conjugate of the hopping. This translates to $H(-d)=H^+(d)$.
herm = qsymm.PointGroupElement(
    R = -np.eye(2), # 2 spatial dimensions
    conjugate = True, 
    antisymmetry = False, 
    U = np.eye(2), # unitary part of the symmetry
    transpose = True, #for hemitian conjugate, have to conjugate and transpose

family = qsymm.continuum_hamiltonian(
    [herm], #the generating symmetries
    dim=2, #spatial dimension
    hermitian = False, #raises an error if you're not on the nonhermitian branch of qsymm
    total_power=1, #maximum power of k we want to appear in our model
h_=qsymm.hamiltonian_generator.hamiltonian_from_family(family,  tosympy=False)

$$ \left[\begin{matrix} 1.0 c_{2} + 1.0 i c_{4} k_{x} + 1.0 i c_{8} k_{y} & 1.0 i c_{0} - 1.0 c_{1} + 1.0 i c_{10} k_{y} - 1.0 c_{11} k_{y} + 1.0 i c_{6} k_{x} - 1.0 c_{7} k_{x} \\ - 1.0 i c_{0} - 1.0 c_{1} + 1.0 i c_{10} k_{y} + 1.0 c_{11} k_{y} + 1.0 i c_{6} k_{x} + 1.0 c_{7} k_{x} & - 1.0 c_{3} + 1.0 i c_{5} k_{x} + 1.0 i c_{9} k_{y} \end{matrix}\right] $$

We can check symmetries of our Hamiltonian model

M1 = qsymm.groups.mirror([0, 1])
M2 = qsymm.groups.mirror([1, 0])
TR = qsymm.time_reversal(2)
PH = qsymm.particle_hole(2)
C = qsymm.groups.chiral(2)
gens = {
    M1, M2, C, TR, PH}

candidates = qsymm.groups.generate_group(gens)
sg, cg = qsymm.symmetries(
    h_, candidates=candidates,
sg, cg
([1], [])

No symmetries, since we didn't specify any. Let's now specify some.

When supplying symmetries to qsymm to build a model, we have to supply the unitary part. Part of it is easy to deduce, for example if the symmetry squares to -1, it will contain sigma_y. Other than that, you have to ensure that the symmetry operation does what you expect it to in the basis you choose.

S = 2*qsymm.groups.spin_matrices(1/2)

P_1 = qsymm.PointGroupElement(np.eye(2), conjugate=True, antisymmetry=True, U = S[0]) #particle-hole symmetry
M_y = qsymm.PointGroupElement(np.diag([1,-1]), False, False, S[0]) #mirror symmetry in y
rot = qsymm.ContinuousGroupGenerator(R = S[1], U = 0.5*S[2]) #continuous rotation symmetry

family = qsymm.continuum_hamiltonian(
    hermitian = False,

No terms allowed it seems. This means either we have to go to 4x4, or that the unitary parts of our symmetries aren't consistent within the basis.

family = qsymm.continuum_hamiltonian(
    [herm, M_y, rot], 
    hermitian = False,
h_=qsymm.hamiltonian_generator.hamiltonian_from_family(family,  tosympy=False)
sg, cg = qsymm.symmetries(
    h_, candidates=candidates,
sg, cg
([1, R(π), M([-1  0]), M([ 0 -1]), T, R(π) T, M([-1  0]) T, M([ 0 -1]) T],
 [exp(-i ϕ L)⋅H(k)⋅exp(i ϕ L) = H(R_ϕ⋅k)
R_ϕ = R(ϕ)
L = [[ 0.5+0.j  0. +0.j]
     [ 0. +0.j -0.5+0.j]]

You can also get additional symmetries you didn't intend in your model, like the time-reversal symmetry we see above, even though we didn't ask for it. You have to try out different unitary parts to get the result you want.

Once you have a hamiltonian, you make it into a model with:

model = qsymm.Model(h_, momenta=["k_x", "k_y"])
# sometimes the coefficients get freaky so we can round them.
rounded_model = dict()
for a,b in model.items():
model = rounded_model

Step 2: saving and loading a model

Here we show how we can move between a qsymm Model, a json serialization, then to a set of numpy arrays.

import json
from astropy.utils.misc import JsonCustomEncoder

def pickle_this(model, name):
    '''model is a qsymm.Model
    name = name + '.json'
    model_dict = { str(key) : model[key].tolist() for key in model.keys() }
    with open(name, 'w') as fp:
        json.dump(model_dict, fp, cls = JsonCustomEncoder)

To load it, you have to use a custom decoder, that's not part of a package.

class JsonCustomDecoder(json.JSONDecoder):
    A decoder function tailored specifically to our purposes: take a dictionary of strings,
    output a dictionary of symbolic keys (Mul types etc.) and arrays.

    def __init__(self, *args, **kwargs):
        json.JSONDecoder.__init__(self, object_hook=self.object_hook, *args, **kwargs)

    def object_hook(self, obj):
        new_dict = dict()
        for key in obj.keys():
            # have to back-convert the key from string to symbols
            symbolic = []
            split_key = key.split("*")

            if str() in set(split_key):
                # there's a power: proceed with caution. This procedure is also valid for
                # only one single power being present.
                power_location = np.where(np.array(split_key) == str())[0][0]
                split_key_power = [
                    split_key[power_location - 1],
                    split_key[power_location + 1],
                key_1 = split_key[: power_location - 1]
                key_2 = split_key[power_location + 2 :]
                split_key = key_1 + key_2  # key terms without the power
                # split key power is of the form ['k', 'x' ]
                    sympy.Pow(Symbol(split_key_power[0]), int(split_key_power[1]))
                # now add the rest
                for term in split_key:
                    symbolic = [sympy.Mul(symbolic[-1], Symbol(term))]
                # no power, proceed normally
                for term in split_key:
                    if symbolic == []:
                        symbolic = [Symbol(term)]
                        symbolic = [sympy.Mul(symbolic[-1], Symbol(term))]
            symbolic_key = symbolic[0]
            array = []
            for line in obj[key]:
                lineski = []
                for entry in line:
                    lineski.append(entry[0] + 1j * entry[1])
            new_dict[symbolic_key] = np.array(array)
        return new_dict

To unpickle, do something like

with open('filename.json', 'r') as fp:
    unpickled_dictionary = json.load(fp, cls = JsonCustomDecoder)

Step 3: using our model as a value function in a kwant system

While building a kwant system, you do something like

syst[lat.shape(shape_edge, (0, 0))] = onsite
syst[lat.neighbors(n=1)] = hopping

for example.

onsite() and hopping() and value functions that take each site and each pair of sites connected by a bond in the system respectively, and return either an onsite or a hopping value for this site or bond.

We can use load our models into the function that builds our kwant system before calling onsite and hopping, and define the latter as:

def onsite(
    Value function used in the creation of the kwant FiniteSystem instances.
    Determines the onsite values.

    :param Site site: A kwant Site instance
    :param Model omodel: A qsymm lambdified Model
    :param dict omodel_parameters: the parameters to be used by omodel

    :rtype: valueFunction
    k_and_r = dict(k_x=0, k_y=0, r=0, os=1)

    return omodel(site, **omodel_parameters, **k_and_r)

def hopping(
    Value function used in the creation of the kwant FiniteSystem instances.
    Determines the hopping values based on the distance between the two sites
    to be connected. Periodicity determines how distances are treated.

    :param Site site1,2: kwant Site instances
    :param Model hmodel: A qsymm lambdified Model
    :param dict hmodel_parameters: the parameters to be used by hmodel

    :rtype: valueFunction
    #extracting the bond length
    delta = site2.pos - site1.pos
    if la.norm(delta) != 0:
        d = delta / la.norm(delta)
        d = delta

    k_and_r = dict(k_x=d[0], k_y=d[1], r=1, os=0)  # set onsite to 0

    return hmodel(site1, site2, **hmodel_parameters, **k_and_r)

You see that we have to be able to turn off the onsite part of the model when we return the hopping value function, and vice versa for the onsite. We can do this by either adding additional coefficients to our models like so:

onsite_coeffs = ['mu1','lam1','lam2','mu2',0,0,0,0]
hopping_coeffs = [0,0,0,0,'t1','o1','o2','t2']

    family, coeffs = onsite_coeffs, tosympy=False
    family, coeffs = hopping_coeffs, tosympy=False

Or by saving the two different models h_os and h_hop separately.

You also notice that we provide "lambdified" models to kwant. What this means is that we are preparing the model to behave like a value function. It's important to note that this lambdification takes a bit of time, which adds up over all the sites and bonds of a lattice, so best to do it once at the beginning.

You'll also notice that onsite() gets omodel and hopping gets hmodel. These aren't two distinct qsymm models, they're just the same qsymm model lambdified two different ways:

with open("filename.json", "r") as fp:
    model = json.load(fp, cls=JsonCustomDecoder)
model = qsymm.Model(model, momenta=["k_x", "k_y"])
omodel = model.lambdify(onsite=True)
hmodel = model.lambdify(hopping=True)

And there you have it.