Spin-1/2 system with sublattice particle consevation
This example makes use of the user_basis class to define the Hamiltonian
\[H = \sum_{j=0}^{N/2-1} t (\tau^+_{j+1}\tau^-_j + \sigma^+_{j+1}\sigma^-_j + \mathrm{h.c.}) + U \sigma^z_j\tau^z_j\]
where \(\sigma\) and \(\tau\) describe hardcore bosons on the two legs of a ladder geometry. Note that particles cannot be exchanged between the legs of the ladder, which allows to further reduce the Hilbert space dimension using a Hilbert space constraint.
Please consult this post – user_basis tutorial – for more detailed explanations on 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###########################################################################
12# example 15 #
13# In this script we demonstrate how to apply the user_basis to #
14# construct a spin-1/2 model with sublattice particle conservation. #
15###########################################################################
16from quspin.operators import hamiltonian
17from quspin.basis import spin_basis_1d # Hilbert space spin basis_1d
18from quspin.basis.user import user_basis # Hilbert space user basis
19from quspin.basis.user import (
20 next_state_sig_32,
21 op_sig_32,
22 map_sig_32,
23) # user basis data types
24from numba import carray, cfunc # numba helper functions
25from numba import uint32, int32 # numba data types
26import numpy as np
27from scipy.special import comb
28
29#
30N_half = 4 # sublattice total number of sites
31N = 2 * N_half # total number of sites
32Np = (N_half // 2, N_half // 2) # sublattice magnetizations
33
34
35#
36###### function to call when applying operators
37@cfunc(op_sig_32, locals=dict(n=int32, b=uint32))
38def op(op_struct_ptr, op_str, ind, N, args):
39 # using struct pointer to pass op_struct_ptr back to C++ see numba Records
40 op_struct = carray(op_struct_ptr, 1)[0]
41 err = 0
42 ind = N - ind - 1 # convention for QuSpin for mapping from bits to sites.
43 n = (op_struct.state >> ind) & 1 # either 0 or 1
44 b = 1 << ind
45 #
46 if op_str == 110: # "n" is integer value 110 (check with ord("n"))
47 op_struct.matrix_ele *= n
48 elif op_str == 43: # "+" is integer value 43 (check with ord("+"))
49 if n:
50 op_struct.matrix_ele = 0
51 else:
52 op_struct.state ^= b # create hcb
53 elif op_str == 45: # "-" is integer value 45 (check with ord("-"))
54 if n:
55 op_struct.state ^= b # destroy hcb
56 else:
57 op_struct.matrix_ele = 0
58 else:
59 op_struct.matrix_ele = 0
60 err = -1
61 #
62 return err
63
64
65#
66op_args = np.array([], dtype=np.uint32)
67
68
69#
70###### function to read user-imported basis into QuSpin
71#
72# function to call when generating next_state
73@cfunc(
74 next_state_sig_32,
75 locals=dict(N_half=int32, t=uint32, s_right=uint32, s_left=uint32),
76)
77def next_state(s, counter, N, args):
78 # unpack args
79 mask = args[0]
80 s_right_min = args[1]
81 s_right_max = args[2]
82 N_half = args[3] # = (N>>1), sublattice system size
83 #
84 # split sublattice
85 s_left = s >> N_half
86 s_right = s & mask
87 # increment s_right unless it has reached the last state,
88 if s_right < s_right_max:
89 if s_right > 0:
90 t = (s_right | (s_right - 1)) + 1
91 s_right = t | ((((t & (0 - t)) // (s_right & (0 - s_right))) >> 1) - 1)
92 # otherwise op_structet s_right to first state and increment s_left.
93 else:
94 s_right = s_right_min
95 if s_left > 0:
96 t = (s_left | (s_left - 1)) + 1
97 s_left = t | ((((t & (0 - t)) // (s_left & (0 - s_left))) >> 1) - 1)
98 # combine and return next state.
99 return (s_left << N_half) + s_right
100
101
102#
103### optional arguments to pass into next_state
104s_right_min = sum(1 << i for i in range(Np[1])) # fill first bits
105s_right_max = sum(1 << (N_half - i - 1) for i in range(Np[1])) # fill last bits
106mask = 2**N_half - 1 # fill all bits
107next_state_args = np.array([mask, s_right_min, s_right_max, N >> 1], dtype=np.uint32)
108
109
110#
111# python function to calculate the starting state to generate the particle conserving basis
112def get_s0_pcon(N, Np):
113 """calculates the starting state to generate the particle conserving basis."""
114 N_half = N >> 1
115 Np_left, Np_right = Np
116
117 s_left = sum(1 << i for i in range(Np_left))
118 s_right = sum(1 << i for i in range(Np_right))
119 return (s_left << N_half) + s_right
120
121
122#
123# python function to calculate the size of the particle-conserved basis,
124# i.e. BEFORE applying pre_check_state and symmetry maps
125def get_Ns_pcon(N, Np):
126 """calculates the size of the particle conservation basis (ignoring symmetries at this stage)."""
127 N_half = N >> 1
128 return comb(N_half, Np[0], exact=True) * comb(N_half, Np[1], exact=True)
129
130
131#
132###### construct user_basis
133# define maps dict
134maps = dict() # no symmetries
135# define particle conservation and op dicts
136pcon_dict = dict(
137 Np=Np,
138 next_state=next_state,
139 next_state_args=next_state_args,
140 get_Ns_pcon=get_Ns_pcon,
141 get_s0_pcon=get_s0_pcon,
142)
143op_dict = dict(op=op, op_args=op_args)
144# create user basis
145basis = user_basis(
146 np.uint32, N, op_dict, allowed_ops=set("n+-"), sps=2, pcon_dict=pcon_dict, **maps
147)
148# print basis
149print(basis)
150#
151###### construct Hamiltonian
152# site-coupling lists
153t_list = [
154 [1.0, i, (i + 1) % N_half] for i in range(N_half)
155] # first sublattice/leg of the ladder
156t_list += [
157 [t, N_half + i, N_half + j] for t, i, j in t_list
158] # second sublattice/leg of the ladder
159U_list = [[1.0, i, i + N_half] for i in range(N_half)]
160# operator string lists
161static = [["+-", t_list], ["-+", t_list], ["nn", U_list]]
162# compute Hamiltonian, no checks have been implemented
163no_checks = dict(check_symm=False, check_pcon=False, check_herm=False)
164H = hamiltonian(static, [], basis=basis, dtype=np.float64, **no_checks)