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