A tutorial on QuSpin’s user_basis¶

Spin-1/2, Fermions, Hardcore Bosons¶

Spin-\(1/2\), hardcore-boson, and fermion states all have in common that they allow at most one particle pre site. Thus, in the particle/z-spin basis, the basis states are given by strings of ones and zeros. The number of ones sets the particle number/magnetization of a given state.

Coincidentally, integers are stored in memory using the so-called binary representation: every integer is assigned a unique bit string. Therefore, a particularly efficient way of storing basis states is via this integer representation.

CONVENTION: in QuSpin, the integer corresponding to a given spin configuration is given by the integer representation of the Fock state bitstring. The lattice indices are defined in reversed order. For instance in a chain of four sites, we have:

• \(|0000\rangle \leftrightarrow 0\): empty lattice,

• \(|0001\rangle \leftrightarrow 1\): one particle at site 3,

• \(|0010\rangle \leftrightarrow 2\): one particle at site 2,

• \(|1000\rangle \leftrightarrow 8\): one particle at site 0,

In short, the bitstring is the same as the quantum state; the lattice site label is reversed relative to the bitstring (i.e. from left to right starting with 0).

Note: the user_basis allows the user to adopt their own convention. For consistency, this tutorial follows the above QuSpin convention.

>>> j = N - j - 1     # compute lattice site in reversed bit configuration (cf QuSpin convention for mapping from bits to sites)
>>> occ = (s>>j)&1    # occupation: either 0 or 1

>>> j = N - j - 1     # compute lattice site in reversed bit configuration (cf QuSpin convention for mapping from bits to sites)
>>> b = 1; b <<= j    # compute a "mask" integer b which is 1 on site j and zero elsewhere
>>> s ^= b            # flip occupation on site j

Local Hilbert space dimension larger than two: bosons, higher spins, etc.¶

When dealing with bosons or higher spins, the binary representation is no longer convenient, since the local on-site occupation can be larger than one.

>>> j = N - j - 1            # compute lattice site in reversed bit configuration (cf QuSpin convention for mapping from bits to sites)
>>> occ = (s//(sps**j))%sps  # occupation: can be 0, 1, ..., sps-1

>>> j = N - j - 1            # compute lattice site in reversed bit configuration (cf QuSpin convention for mapping from bits to sites)
>>> b = sps**j               # obtain mask integer b
>>> occ = (s//b))%sps        # compute occupation on site j
>>> if (occ+1<sps): r += b   # increase occupation on site j by one

>>> j = N - j - 1            # compute lattice site in reversed bit configuration (cf QuSpin convention for mapping from bits to sites)
>>> b = sps**j               # obtain mask integer b
>>> occ = (s//b)%sps         # compute occupation on site j
>>> if (occ>0): r -= b       # decrease occupation on site j by one

Notes:¶

• even though in the case \(sps=2\), the above expressions reproduce the corresponding spin-1/2 expressions, they are not as efficient computationally.

• convenient quspin functions to transform between integer and quspin bit representations are basis.int_to_state() and basis.state_to_int().

• the attribute basis.states holds all states of the basis in their integer representation. The function basis.index() returns the position of a given state in the basis (useful for computing specific matrix elements or for defining states).

• printing a basis object print(basis) displays the states in their quantum mechanical notation, together with the order in which they appear in the basis, and the integer representation.

op(op_struct_ptr, op_str, site_ind, N, args)¶

This function method contains user-defined action of operators \(O\) on the integer states \(|s\rangle\) which produces the matrix elements \(\mathrm{me}\) via \(O|s\rangle = \mathrm{me}|s'\rangle\).

• op_struct_ptr: a C++ pointer to an object which, after being cast into an array using op_struct=carray(op_struct_ptr,1)[0], contains the attributes op_struct.state (which contains the quantum state in integer representation), and op_struct.matrix_ele (the value of the corresponding matrix element which defines the action of the operator \(O\).).

• op_str: holds the operator string (e.g. +, -, z, n, or any custom user-defined letter). Due to limitations in compiling python functions (see section on numba below), the charactors are passed in as integers using utf-8 unicode, e.g. + corresponds to the integer 43. Because of this one has to compare op_str to an integer representing the character of your choice in the body of op(). The integer, corresponding to any character str can be found in python using ord(str) or by looking up a utf-8 unicode table.

• N: the total number of lattice sites.

• args: optional arguments passed into the CFunc op; must be a np.ndarray of dtype basis_dtype.

The CFunc op returns an integer err which is used by QuSpin to throw different error messages. The following are reserved by QuSpin:

• err=0: the calculation was completed successfully.

• err=-1: no matching operator string was found.

• err=1: using real dtype for complex value.

If the error code returned is not one of these values QuSpin will raise a RunTimeError with the message: “user defined error code: <err>”, with <err> being the integer value returned by the user defined op CFunc. In this way the user can provide custom errors.

Notes

• this functionality will not support branching, i.e. no linear combination of multiple states in the basis, e.g. \(O|s\rangle = \mathrm{me}_1|s'_1\rangle + \mathrm{me}_2|s'_2\rangle + \dots\), is NOT allowed.

next_state(s, counter, N, args)¶

This function method provides a user-defined particle conservation rule, which constructs the basis in ascending order by numerical value. Given the initial state s0, next_state() generates all other states successively. Hence, if next_state() is set to conserve particle number then the particle number sector is defined by the initial state s0.

• s: quantum state in integer representation.

• counter: an integer which counts internally how many times the function has been called. The incrementation of counter will occur in the underlying C++ code, i.e. the user should not attempt to do this in the function body of next_state(). Can be used, e.g., to index an array passed in args, cf. Applying symmetries to reduce user-imported bases using QuSpin.

• args: a np.ndarray of the same data type as the user_basis. Can be used to pass optional arguments, e.g. to pass a precomputed basis into QuSpin in order to reduce it to a given symmetry sector: cf. Applying symmetries to reduce user-imported bases using QuSpin.

Two extra python functions required: they are not called inside next_state(), but are required by QuSpin to run next_state().

• get_s0_pcon(N,Np): given the total number of sites N and (the tuple of) particle sector Np this function computes the initial state, to be used by next_state() to construct the entire basis.

• get_Ns_pcon(N,Np): given the total number of sites N and (the tuple of) particle sector Np this function computes the Hilbert space dimension (i.e. the size of the basis) with particle number conservation only (in other words, get_Ns_pcon() should be equal to the number of iterations in next_state() required to exhaust the states search). get_Ns_pcon() returns an integer required to allocate memory for the particle-conserving basis. Note that get_Ns_pcon() ignores any possible reduction due to lattice symmetries (implemened via maps, see below), i.e. get_Ns_pcon() may not correspond to the final integer basis.Ns.

Notes

• there is no need to define next_state() if no particle number conservation use is intended, cf. Quantum scars Hamiltonian: spin-1/2 system with Hilbert space constraint.

• one can use this function, e.g., to implement sublattice particle number conservation, and similar features.

• next_state(), together with the entire set of related functions and variables is passed to the user_basis constructor via the pcon_dict dictionary.

• next_state() is a numba.CFunc object, but get_s0_pcon() and get_Ns_pcon() are regular python functions.

pre_check_state(s, N, args)¶

This optional function method provides user-defined extra filtering of basis states. The function body contains a boolean operation which, when applied to the basis states one at a time, determines whether to keep a state in the basis or not. This is independent of the lattice symmetries (implemented via maps, see below), and can be interpreted as a projection or a constraint on the Hilbert space.

A simple example of what pre_check_state() can be useful for is this: suppose you want a spinful_fermion_basis() without doubly occupied sites. This can be achieved by ajusting the body of pre_check_state() to eliminate such states. QuSpin will then first generate the basis with doble occupancies using next_state(), and subsequntly get rid of the doubly-occupied states using pre_check_state(). Another example is shown in Quantum scars Hamiltonian: spin-1/2 system with Hilbert space constraint.

• s: quantum state in integer representation.

• N: the total number of lattice sites.

• args: a np.ndarray of the same data type as the user_basis. Can be used to pass optional arguments.

count_particles(s, p_number_ptr, args)¶

This optional function method counts the total number of particles/magnetization in a given state.

• s: quantum state in integer representation.

• p_number_ptr: pointer of length n_sectors to fill in the number of particles. Each entry should correspond to the given particle sector in Np.

• args: a np.ndarray of data type np.integer. Can be used to pass optional arguments.

Notes

• this function does not return anything. Fill in the pointer p_number_ptr with the output instead.

• make sure that p_number_ptr[i] corresponds to the particle sector Np[i], etc.

symmetry_map(s, N, sign_ptr, args)¶

• s: quantum state in integer representation.

• N: the total number of lattice sites.

• sign_ptr: a pointer to a number which one can use to pass sign changes back to QuSpin; used for fermion systems.

• args: a np.ndarray of the same data type as the user_basis. Can be used to pass optional arguments, e.g. sps in the case of bosons.

Notes

• all four arguments must be present in the function, even if some are not used (this is required to keep the syntax general for all particle species).

System-size independent symmetries¶

System-size independent symmetries contain as a parameter the system size \(N\). As a result, they apply to all system sizes. Examples of such symmetries are

def parity(x,N,sign_ptr,args):
        """ works for all system sizes N, spin-1/2 only. """
        out = 0
        s = N-1
        #
        out ^= (x&1)
        x >>= 1
        while(x):
                out <<= 1
                out ^= (x&1)
                x >>= 1
                s -= 1
        #
        out <<= s
        return out

def translation(x,N,sign_ptr,args):
        """ works for all system sizes N, spin-1/2 only. """
        shift = 1 # translate state by shift sites
        period = N # periodicity/cyclicity of translation
        xmax = (1<<N)-1 # largest integer allowed to appear in the basis
        #
        l = (shift+period)%period
        x1 = (x >> (period - l))
        x2 = ((x << l) & xmax)
        #
        return (x2 | x1)

Symmetries for fixed system sizes using precomputed masks¶

The convenience to define symmetry maps which apply to all system sizes comes at a cost of speed. This can be circumvented by defining system-size specific maps, using integer masks to perform the bit operations. These masks also depend on the data type of the integer storing the state.

Luckily, there is a great tool to compute the symmetry maps, available at http://programming.sirrida.de/calcperm.php. All one needs to do is find the permutation of the lattice sites under the symmetry, and pass it to the tool to obain the symmetry map that acts on integers. Let us demonstrate how this works using two examples.

def parity(x,N,sign_ptr,args):
    """ works for N=10 sites and 32 bit-integers, spin-1/2 states only """
    return    (  ((x & 0x00004010) << 1)
               | ((x & 0x00002008) << 3)
               | ((x & 0x00001004) << 5)
               | ((x & 0x00000802) << 7)
               | ((x & 0x00000401) << 9)
               | ((x & 0x00080200) >> 9)
               | ((x & 0x00040100) >> 7)
               | ((x & 0x00020080) >> 5)
               | ((x & 0x00010040) >> 3)
               | ((x & 0x00008020) >> 1))

def translation(x,N,sign_ptr,args):
    """ works for N=10 sites and 32 bit-integers spin-1/2 states only. """
    return ((x & 0x0007fdff) << 1) | ((x & 0x00080200) >> 9)

Symmetry maps dictionary¶

In the user_basis, the functions encoding the symmetry action are referred to as maps. Every map has as its first argument the integer (state) to be tansformed, followed by the number of sites. For fermionic systems, the symmetry action can also modify the fermion sign of a given state. Therefore, the last argument is a sign_ptr.

Symmetries are passed to the user_basis constructor via a python dictionary, called maps. The keys are arbitrary strings which define a unique name for each map; the corresponding values are tuples of four entries: (map function, periodicity, quantum number, args). The symmetry periodicity (or cyclicity, or multiplicity) is the smallest integer \(m_Q\), such that \(Q^{m_Q} = 1\).

>>> maps = dict(T_block=(translation,10,0,T_args), P_block=(parity,2,0,P_args), )

Notes:

• all map functions need to be cast as decorated numba cfuncs (SEE below). The examples above are shown as python functions, so the user can test them in practice. Luckily, the same code can be used in numba.

• even though some arguments of the map functions are not used in the function bodies, the user is required to include them in the functin definition (and no more!). This allows to keep the code general. The names of these arguments are arbitrary, but their data typs are not.

Data types¶

Unlike python, C++ code requires the user to specify the data types of all variables (so-called strong typing). For this purpose, numba supports various data types, e.g. uint32, or int32. They are typically imported from numba in the beginning of the python script.

Function decorators¶

To indicate that the function we wrote in python should be compiled as a C++ code by numba, we use the @cfunc(signature,locals=dict()) decorator. The arguments of the decorator are the function variable signature (which contains the data types of all function variables), and locals which is a dictionary containing the data types of all other variables defined and used privately inside the function body.

In QuSpin, we provide the signatures next_state_sig_32, op_sig_32, map_sig_32, count_particles_32; next_state_sig_64, op_sig_64, map_sig_64, count_particles_64 which are compatible with the QuSpin base code. The name of the signature refers to the function type it is designed for, and the integer in the end specifies the data type the user_basis will be constructed with. These signatures can be imported from the user_basis.

As an example, consider the translation() python function defined above. To make this a numba.CFunc object, it suffices to place the decorator:

from quspin.basis.user import map_sig_32 # user basis data types
from numba import cfunc
from numba import uint32,int32 # numba data types
#
@cfunc(map_sig_32,
        locals=dict(shift=uint32,xmax=uint32,x1=uint32,x2=uint32,period=int32,l=int32,) )
def translation(x,N,sign_ptr,args):
        """ works for all system sizes N. """
        shift = 1 # translate state by shift sites
        period = N # periodicity/cyclicity of translation
        xmax = (1<<N-1)
        #
        l = (shift+period)%period
        x1 = (x >> (period - l))
        x2 = ((x << l) & xmax)
        #
        return (x2 | x1)

We use the signature map_sig_32 because it is designed to decorate symmetry map functions. Moreover, the local (private) variable data types are defined via locals=dict(shift=uint32,xmax=uint32,x1=uint32,x2=uint32,period=int32,l=int32,). These variables appear in the function body.

Notes

• because QuSpin provides predefined CFunc signatures, every CFunc (see function methods above) has a predefined, fixed number of arguments. Moreover, the data types of the arguments is also fixed. Even if some arguments are not used in the CFunc body, they have to appear in the function definition.

• if you mess up the data types, most likely you will receive a numba error. In such cases, we suggest that you remove the CFunc decorator and debug your function in python as you would normally do. Once you are confident that the function does its job, put back the decorator and pass it to the user_basis constructor.

• Unfortunately numba will not allow printing inside CFuncs because of complications dealing with Python’s Global Interpreter Lock (GIL). Because of this, debugging these functions can be a bit tedious. Always make sure that your code performs properly by running it purely with python before attemping to use it within the user_basis.

Scripts to construct spin, fermion, and boson bases¶

The following three examples demonstrate how the user_basis recovers the functionality of the basis_general classes:

• spin-1/2 Heisenberg model in 1d, download script

• spinless fermions with nearest-neighbor interactions in 1d, download script

• Bose-Hubbard model in 1d, download script

Scripts to demonstrate the additional functionality introduced by the user_basis¶

• Quantum scars Hamiltonian: spin-1/2 system with Hilbert space constraint, download script

• Spin-1/2 system with sublattice particle consevation, download script

• Applying symmetries to reduce user-imported bases using QuSpin, download script

• Gell-Mann Operators for Spin-1 systems, download script

• Majorana Fermions: Spinless Fermi-Hubbard Model, download script