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) h_.tosympy()
$$ \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, continuous_rotations=True, prettify=True, sparse_linalg=True ) 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( [herm,P_1,M_y,rot], dim=2, hermitian = False, total_power=1, ) len(family)
0
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], dim=2, hermitian = False, total_power=1, ) h_=qsymm.hamiltonian_generator.hamiltonian_from_family(family, tosympy=False) sg, cg = qsymm.symmetries( h_, candidates=candidates, continuous_rotations=True, prettify=True, sparse_linalg=True ) 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(): rounded_model[a]=np.round(b,2) model = rounded_model qsymm.Model(model).tosympy()
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.
Encoding
import json from astropy.utils.misc import JsonCustomEncoder def pickle_this(model, name): '''model is a qsymm.Model ''' #name 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)
Decoding
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' ] symbolic.append( 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))] else: # no power, proceed normally for term in split_key: if symbolic == []: symbolic = [Symbol(term)] else: 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]) array.append(lineski) 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( site, omodel, model_parameters **kwargs ): """ 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( site1, site2, hmodel, hmodel_paramters, **kwargs, ): """ 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) else: 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'] h_os=qsymm.hamiltonian_generator.hamiltonian_from_family( family, coeffs = onsite_coeffs, tosympy=False ) h_hop=qsymm.hamiltonian_generator.hamiltonian_from_family( family, coeffs = hopping_coeffs, tosympy=False ) h_=h_os+h_hop
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.