Gell-Mann Operators for Spin-1 systems
In this example, we show how to use the user_basis class to define the Gell-Mann operators to construct spin-1 Hamiltonians in QuSpin – the SU(3) equivalent of the Pauli operators for spin-1/2.
The eight generators of SU(3), are defined as
\[\begin{split}\lambda^1 = \begin{pmatrix} 0 & 1 & 0 \\ 1 & 0 & 0 \\ 0 & 0 & 0 \end{pmatrix},\quad \lambda^2 = \begin{pmatrix} 0 & -i & 0 \\ i & 0 & 0 \\ 0 & 0 & 0 \end{pmatrix}, \quad \lambda^3 = \begin{pmatrix} 1 & 0 & 0 \\ 0 & -1 & 0 \\ 0 & 0 & 0 \end{pmatrix},\end{split}\]
\[\begin{split}\lambda^4 = \begin{pmatrix} 0 & 0 & 1 \\ 0 & 0 & 0 \\ 1 & 0 & 0 \end{pmatrix},\quad \lambda^5 = \begin{pmatrix} 0 & 0 & -i \\ 0 & 0 & 0 \\ i & 0 & 0 \end{pmatrix},\end{split}\]
\[\begin{split}\lambda^6 = \begin{pmatrix} 0 & 0 & 0 \\ 0 & 0 & 1 \\ 0 & 1 & 0 \end{pmatrix},\quad \lambda^7 = \begin{pmatrix} 0 & 0 & 0 \\ 0 & 0 & -i \\ 0 & i & 0 \end{pmatrix}, \quad \lambda^8 = \frac{1}{\sqrt{3}} \begin{pmatrix} 1 & 0 & 0 \\ 0 & 1 & 0 \\ 0 & 0 & -2 \end{pmatrix}\end{split}\]
We now define them as operator strings for a custom spin-1 user_basis (please consult this post – user_basis tutorial – for more detailed explanations on using the user_basis class). In QuSpin, the basis constructor accepts operator strings of unit length; therefore, we use the operator strings 1, 2, 3, 4, 5, 6, 7, 8 to denote the operator \(\lambda^j\).
Script
1# -*- coding: utf-8 -*-
2from __future__ import unicode_literals
3
4
5#
6import sys, os
7
8os.environ["KMP_DUPLICATE_LIB_OK"] = (
9 "True" # uncomment this line if omp error occurs on OSX for python 3
10)
11os.environ["OMP_NUM_THREADS"] = "1" # set number of OpenMP threads to run in parallel
12os.environ["MKL_NUM_THREADS"] = "1" # set number of MKL threads to run in parallel
13#
14
15######################################################################
16# example 23 #
17# This example shows how to use the `user_basis` to define Gell-Mann #
18# operators. This allows to construct Hamiltonians using the 8 SU(3) #
19# generators. #
20######################################################################
21from quspin.operators import hamiltonian # Hamiltonians and operators
22from quspin.basis import boson_basis_1d # Hilbert space spin basis_1d
23from quspin.basis.user import user_basis # Hilbert space user basis
24from quspin.basis.user import (
25 next_state_sig_32,
26 op_sig_32,
27 map_sig_32,
28 count_particles_sig_32,
29) # user basis data types signatures
30from numba import carray, cfunc # numba helper functions
31from numba import uint32, int32, float64, complex128 # numba data types
32import numpy as np
33from scipy.special import comb
34
35#
36N = 2 # lattice sites
37sps = 3 # states per site: 3 for a spin-1 system
38Np = N // 2 # total number of spin-1/s
39
40
41#
42############ create boson user basis object #############
43#
44###### function to call when applying operators
45@cfunc(
46 op_sig_32,
47 locals=dict(
48 b=uint32, occ=int32, sps=uint32, me_offdiag=complex128, me_diag=float64
49 ),
50)
51def op(op_struct_ptr, op_str, site_ind, N, args):
52 # using struct pointer to pass op_struct_ptr back to C++ see numba Records
53 op_struct = carray(op_struct_ptr, 1)[0]
54 err = 0
55 sps = 3
56 me_offdiag = 1.0
57 me_diag = 1.0
58 #
59 site_ind = N - site_ind - 1 # convention for QuSpin for mapping from bits to sites.
60 occ = (op_struct.state // sps**site_ind) % sps # occupation
61 b = sps**site_ind
62 #
63 if op_str == 43: # "+" is integer value 43 = ord("+")
64 if (occ + 1) % sps != 1:
65 me_offdiag *= np.sqrt((occ + 1) % sps)
66 else:
67 me_offdiag *= np.sqrt(2.0)
68 op_struct.state += b if (occ + 1) < sps else 0
69
70 elif op_str == 45: # "-" is integer value 45 = ord("-")
71 if occ != 1:
72 me_offdiag *= np.sqrt(occ)
73 else:
74 me_offdiag *= np.sqrt(2.0)
75 op_struct.state -= b if occ > 0 else 0
76
77 elif op_str == 122: # "z" is integer value 122 = ord("n")
78 me_diag *= occ - 1
79
80 # define the 8 Gell-Mann matrices
81
82 elif op_str == 49: # "1" is integer value 49 = ord("1"): lambda_1 Gell-Mann matrix
83 if occ == 2:
84 op_struct.state -= b
85 elif occ == 1:
86 op_struct.state += b
87 else:
88 me_offdiag *= 0.0
89
90 elif op_str == 50: # "2" is integer value 50 = ord("2"): lambda_2 Gell-Mann matrix
91 if occ == 2:
92 op_struct.state -= b
93 me_offdiag *= 1.0j
94 elif occ == 1:
95 op_struct.state += b
96 me_offdiag *= -1.0j
97 else:
98 me_offdiag *= 0.0
99
100 elif op_str == 51: # "3" is integer value 51 = ord("3"): lambda_3 Gell-Mann matrix
101 if occ == 1:
102 me_diag *= -1.0
103 elif occ == 0:
104 me_diag *= 0.0
105
106 elif op_str == 52: # "4" is integer value 52 = ord("4"): lambda_4 Gell-Mann matrix
107 if occ == 2:
108 op_struct.state -= 2 * b
109 elif occ == 0:
110 op_struct.state += 2 * b
111 else:
112 me_offdiag *= 0.0
113
114 elif op_str == 53: # "5" is integer value 53 = ord("5"): lambda_5 Gell-Mann matrix
115 if occ == 2:
116 op_struct.state -= 2 * b
117 me_offdiag *= 1.0j
118 elif occ == 0:
119 op_struct.state += 2 * b
120 me_offdiag *= -1.0j
121 else:
122 me_offdiag *= 0.0
123
124 elif op_str == 54: # "6" is integer value 54 = ord("6"): lambda_6 Gell-Mann matrix
125 if occ == 1:
126 op_struct.state -= b
127 elif occ == 0:
128 op_struct.state += b
129 else:
130 me_offdiag *= 0.0
131
132 elif op_str == 55: # "7" is integer value 55 = ord("7"): lambda_7 Gell-Mann matrix
133 if occ == 1:
134 op_struct.state -= b
135 me_offdiag *= 1.0j
136 elif occ == 0:
137 op_struct.state += b
138 me_offdiag *= -1.0j
139 else:
140 me_offdiag *= 0.0
141
142 elif op_str == 56: # "8" is integer value 56 = ord("8"): lambda_8 Gell-Mann matrix
143 if occ == 2:
144 me_diag *= 1.0 / np.sqrt(3.0)
145 elif occ == 1:
146 me_diag *= 1.0 / np.sqrt(3.0)
147 else:
148 me_diag *= -2.0 / np.sqrt(3.0)
149
150 elif op_str == 73: # "I" is integer value 73 = ord("I")
151 pass
152 else:
153 me_diag = 0.0
154 err = -1
155 #
156 op_struct.matrix_ele *= me_diag * me_offdiag
157 #
158 return err
159
160
161#
162op_args = np.array([sps], dtype=np.uint32)
163
164
165###### function to implement magnetization/particle conservation
166#
167@cfunc(
168 next_state_sig_32,
169 locals=dict(
170 t=uint32,
171 i=int32,
172 j=int32,
173 n=int32,
174 sps=uint32,
175 b1=int32,
176 b2=int32,
177 l=int32,
178 n_left=int32,
179 ),
180)
181def next_state(s, counter, N, args):
182 """implements particle number conservation. Particle number set by initial state, cf `get_s0_pcon()` below."""
183 t = s
184 sps = args[1]
185 n = 0 # running total of number of particles
186 for i in range(N): # loop over lattices sites
187 b1 = (t // args[i]) % sps # get occupation at site i
188 if b1 > 0: # if there is a boson
189 n += b1
190 b2 = (t / args[i + 1]) % sps # get occupation st site ahead
191 if b2 < (sps - 1): # if I can move a boson to this site
192 n -= 1 # decrease one from the running total
193 t -= args[i] # remove one boson from site i
194 t += args[i + 1] # add one boson to site i+1
195 if n > 0: # if any bosons left
196 # so far: moved one boson forward;
197 # now: take rest of bosons and fill first l sites with maximum occupation
198 # to keep lexigraphic order
199 l = n // (
200 sps - 1
201 ) # how many sites can be fully occupied with n bosons
202 n_left = n % (
203 sps - 1
204 ) # leftover of particles on not maximally occupied sites
205 for j in range(i + 1):
206 t -= (t // args[j]) % sps * args[j]
207 if j < l: # fill in with maximal occupation
208 t += (sps - 1) * args[j]
209 elif j == l: # fill with leftover
210 t += n_left * args[j]
211 break # stop loop
212 return t
213
214
215next_state_args = np.array([sps**i for i in range(N)], dtype=np.uint32)
216
217
218# python function to calculate the starting state to generate the particle conserving basis
219def get_s0_pcon(N, Np):
220 sps = 3 # use as global variable
221 l = Np // (sps - 1)
222 s = sum((sps - 1) * sps**i for i in range(l))
223 s += (Np % (sps - 1)) * sps**l
224 return s
225
226
227# python function to calculate the size of the particle-conserved basis, i.e.
228# BEFORE applying pre_check_state and symmetry maps
229def get_Ns_pcon(N, Np):
230 Ns = 0
231 sps = 3
232 for r in range(Np // sps + 1):
233 r_2 = Np - r * sps
234 if r % 2 == 0:
235 Ns += comb(N, r, exact=True) * comb(N + r_2 - 1, r_2, exact=True)
236 else:
237 Ns += -comb(N, r, exact=True) * comb(N + r_2 - 1, r_2, exact=True)
238
239 return Ns
240
241
242#
243###### define symmetry maps
244#
245@cfunc(
246 map_sig_32,
247 locals=dict(
248 shift=uint32,
249 out=uint32,
250 sps=uint32,
251 i=int32,
252 j=int32,
253 ),
254)
255def translation(x, N, sign_ptr, args):
256 """works for all system sizes N."""
257 out = 0
258 shift = args[0]
259 sps = args[1]
260 for i in range(N):
261 j = (i + shift + N) % N
262 out += (x % sps) * sps**j
263 x //= sps
264 #
265 return out
266
267
268T_args = np.array([1, sps], dtype=np.uint32)
269
270
271#
272@cfunc(map_sig_32, locals=dict(out=uint32, sps=uint32, i=int32, j=int32))
273def parity(x, N, sign_ptr, args):
274 """works for all system sizes N."""
275 out = 0
276 sps = args[0]
277 for i in range(N):
278 j = (N - 1) - i
279 out += (x % sps) * (sps**j)
280 x //= sps
281 #
282 return out
283
284
285P_args = np.array([sps], dtype=np.uint32)
286
287
288#
289@cfunc(map_sig_32, locals=dict(out=uint32, sps=uint32, i=int32))
290def inversion(x, N, sign_ptr, args):
291 """works for all system sizes N."""
292 out = 0
293
294 sps = args[0]
295 for i in range(N):
296 out += (sps - x % sps - 1) * (sps**i)
297 x //= sps
298 #
299 return out
300
301
302Z_args = np.array([sps], dtype=np.uint32)
303
304
305#
306###### define function to count particles in bit representation
307#
308@cfunc(
309 count_particles_sig_32,
310 locals=dict(
311 s=uint32,
312 ),
313)
314def count_particles(x, p_number_ptr, args):
315 """Counts number of particles/spin-ups in a state stored in integer representation for up to N=32 sites"""
316 #
317 s = x # integer x cannot be changed
318 for i in range(args[0]):
319 p_number_ptr[0] += s % args[1]
320 s /= args[1]
321
322
323n_sectors = 1 # number of particle sectors
324count_particles_args = np.array([N, sps], dtype=np.int32)
325#
326###### construct user_basis
327# define maps dict
328maps = dict(
329 T_block=(translation, N, 0, T_args),
330 P_block=(parity, 2, 0, P_args),
331) # Z_block=(inversion,2,0,Z_args), )
332# define particle conservation and op dicts
333pcon_dict = dict(
334 Np=Np,
335 next_state=next_state,
336 next_state_args=next_state_args,
337 get_Ns_pcon=get_Ns_pcon,
338 get_s0_pcon=get_s0_pcon,
339 count_particles=count_particles,
340 count_particles_args=count_particles_args,
341 n_sectors=n_sectors,
342)
343op_dict = dict(op=op, op_args=op_args)
344# create user basiss
345basis = user_basis(
346 np.uint32,
347 N,
348 op_dict,
349 allowed_ops=set("+-zI12345678"),
350 sps=sps,
351 **maps,
352) # pcon_dict=pcon_dict)
353#
354print(basis)
355#
356############ create Hamiltonian #############
357#
358J = -1.0
359U = +np.sqrt(2.0)
360h = np.sqrt(7.0)
361#
362nn_int = [[+J, j, (j + 1) % N] for j in range(N)]
363onsite_int = [[U, j, j] for j in range(N)]
364onsite_field = [[h, j] for j in range(N)]
365#
366static = [
367 ["11", nn_int],
368 ["22", nn_int],
369 ["88", onsite_int],
370 ["5", onsite_field],
371]
372dynamic = []
373#
374no_checks = dict(check_symm=False, check_pcon=False, check_herm=False)
375H = hamiltonian(static, [], basis=basis, dtype=np.complex128, **no_checks)
376#
377np.set_printoptions(precision=2)
378print(H.toarray())