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 spin_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  # numba data types
 23import numpy as np
 24from scipy.special import comb
 25
 26#
 27N = 6  # lattice sites
 28Np = N // 2  # total number of spin ups
 29
 30
 31#
 32############   create spin-1/2 user basis object   #############
 33#
 34######  function to call when applying operators
 35@cfunc(
 36    op_sig_32,
 37    locals=dict(s=int32, n=int32, b=uint32),
 38)
 39def op(op_struct_ptr, op_str, site_ind, N, args):
 40    # using struct pointer to pass op_struct_ptr back to C++ see numba Records
 41    op_struct = carray(op_struct_ptr, 1)[0]
 42    err = 0
 43    #
 44    site_ind = N - site_ind - 1  # convention for QuSpin for mapping from bits to sites.
 45    n = (op_struct.state >> site_ind) & 1  # either 0 or 1
 46    s = (((op_struct.state >> site_ind) & 1) << 1) - 1  # either -1 or 1
 47    b = 1 << site_ind
 48    #
 49    if op_str == 120:  # "x" is integer value 120 = ord("x")
 50        op_struct.state ^= b
 51
 52    elif op_str == 121:  # "y" is integer value 120 = ord("y")
 53        op_struct.state ^= b
 54        op_struct.matrix_ele *= 1.0j * s
 55
 56    elif op_str == 43:  # "+" is integer value 43 = ord("+")
 57        if n:
 58            op_struct.matrix_ele = 0
 59        else:
 60            op_struct.state ^= b  # create spin
 61
 62    elif op_str == 45:  # "-" is integer value 45 = ord("-")
 63        if n:
 64            op_struct.state ^= b  # destroy spin
 65        else:
 66            op_struct.matrix_ele = 0
 67
 68    elif op_str == 122:  # "z" is integer value 120 = ord("z")
 69        op_struct.matrix_ele *= s
 70
 71    elif op_str == 110:  # "n" is integer value 110 = ord("n")
 72        op_struct.matrix_ele *= n
 73
 74    elif op_str == 73:  # "I" is integer value 73 = ord("I")
 75        pass
 76
 77    else:
 78        op_struct.matrix_ele = 0
 79        err = -1
 80    #
 81    return err
 82
 83
 84op_args = np.array([], dtype=np.uint32)
 85
 86
 87#
 88######  function to implement magnetization/particle conservation
 89#
 90@cfunc(
 91    next_state_sig_32,
 92    locals=dict(t=uint32),
 93)
 94def next_state(s, counter, N, args):
 95    """implements magnetization conservation."""
 96    if s == 0:
 97        return s
 98    #
 99    t = (s | (s - 1)) + 1
100    return t | ((((t & (0 - t)) // (s & (0 - s))) >> 1) - 1)
101
102
103next_state_args = np.array([], dtype=np.uint32)  # compulsory, even if empty
104
105
106# python function to calculate the starting state to generate the particle conserving basis
107def get_s0_pcon(N, Np):
108    return sum(1 << i for i in range(Np))
109
110
111# python function to calculate the size of the particle-conserved basis,
112# i.e. BEFORE applying pre_check_state and symmetry maps
113def get_Ns_pcon(N, Np):
114    return comb(N, Np, exact=True)
115
116
117#
118######  define symmetry maps
119#
120@cfunc(
121    map_sig_32,
122    locals=dict(
123        shift=uint32,
124        xmax=uint32,
125        x1=uint32,
126        x2=uint32,
127        period=int32,
128        l=int32,
129    ),
130)
131def translation(x, N, sign_ptr, args):
132    """works for all system sizes N."""
133    shift = args[0]  # translate state by shift sites
134    period = N  # periodicity/cyclicity of translation
135    xmax = args[1]
136    #
137    l = (shift + period) % period
138    x1 = x >> (period - l)
139    x2 = (x << l) & xmax
140    #
141    return x2 | x1
142
143
144T_args = np.array([1, (1 << N) - 1], dtype=np.uint32)
145
146
147#
148@cfunc(
149    map_sig_32,
150    locals=dict(
151        out=uint32,
152        s=int32,
153    ),
154)
155def parity(x, N, sign_ptr, args):
156    """works for all system sizes N."""
157    out = 0
158    s = args[0]  # N-1
159    #
160    out ^= x & 1
161    x >>= 1
162    while x:
163        out <<= 1
164        out ^= x & 1
165        x >>= 1
166        s -= 1
167    #
168    out <<= s
169    return out
170
171
172P_args = np.array([N - 1], dtype=np.uint32)
173
174
175#
176@cfunc(
177    map_sig_32,
178    locals=dict(
179        xmax=uint32,
180    ),
181)
182def spin_inversion(x, N, sign_ptr, args):
183    """works for all system sizes N."""
184    xmax = args[0]  # maximum integer
185    return x ^ xmax
186
187
188Z_args = np.array([(1 << N) - 1], dtype=np.uint32)
189
190
191#
192######  define function to count particles in bit representation
193#
194@cfunc(count_particles_sig_32, locals=dict(s_count=uint32))
195def count_particles(x, p_number_ptr, args):
196    """Counts number of particles/spin-ups in a state stored in integer representation for up to N=32 sites"""
197    #
198    s_count = x & ((0x7FFFFFFF) >> (31 - args[0]))
199    s_count = s_count - ((s_count >> 1) & 0x55555555)
200    s_count = (s_count & 0x33333333) + ((s_count >> 2) & 0x33333333)
201    s_count = (((s_count + (s_count >> 4)) & 0x0F0F0F0F) * 0x01010101) >> 24
202    #
203    p_number_ptr[0] = s_count
204
205
206n_sectors = 1  # number of particle sectors
207count_particles_args = np.array([N], dtype=np.int32)
208#
209######  construct user_basis
210# define maps dict
211maps = dict(
212    T_block=(translation, N, 0, T_args),
213    P_block=(parity, 2, 0, P_args),
214    Z_block=(spin_inversion, 2, 0, Z_args),
215)
216# define particle conservation and op dicts
217pcon_dict = dict(
218    Np=Np,
219    next_state=next_state,
220    next_state_args=next_state_args,
221    get_Ns_pcon=get_Ns_pcon,
222    get_s0_pcon=get_s0_pcon,
223    count_particles=count_particles,
224    count_particles_args=count_particles_args,
225    n_sectors=n_sectors,
226)
227op_dict = dict(op=op, op_args=op_args)
228# create user basis
229basis = user_basis(
230    np.uint32,
231    N,
232    op_dict,
233    allowed_ops=set("+-xyznI"),
234    sps=2,
235    pcon_dict=pcon_dict,
236    **maps,
237)
238#
239#
240#
241############   create same spin-1/2 basis_1d object   #############
242basis_1d = spin_basis_1d(N, Nup=Np, pauli=True, kblock=0, pblock=1, zblock=1)  #
243#
244#
245print(basis)
246print(basis_1d)
247#
248############   create and compare Hamiltonians   #############
249#
250J = 1.0
251spin_spin = [[J, j, (j + 1) % N] for j in range(N)]
252static = [["xx", spin_spin], ["yy", spin_spin], ["zz", spin_spin]]
253dynamic = []
254#
255no_checks = dict(check_symm=False, check_pcon=False, check_herm=False)
256H = hamiltonian(static, [], basis=basis, dtype=np.float64, **no_checks)
257H_1d = hamiltonian(static, [], basis=basis_1d, dtype=np.float64)
258print(H.toarray())
259print(H_1d.toarray())
260print(np.linalg.norm((H - H_1d).toarray()))