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