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