Applying symmetries to reduce user-imported bases using QuSpin

This example makes use of the user_basis class to take a user-imported basis consisting of ordered integers, and apply QuSpin to reduce it to a desired symmetry sector.

In the example, we manually construct the basis for a two-legged ladder (and handle it as a user-imported array of basis states). Then we use QuSpin to apply translation and parity (reflection) symmetries to reduce the Hilbert space dimension.

Please consult this post – A tutorial on QuSpin’s user_basis – for more detailed explanations on using the user_basis class.

Script

download script

  1
  2
  3
  4
  5
  6
  7
  8
  9
 10
 11
 12
 13
 14
 15
 16
 17
 18
 19
 20
 21
 22
 23
 24
 25
 26
 27
 28
 29
 30
 31
 32
 33
 34
 35
 36
 37
 38
 39
 40
 41
 42
 43
 44
 45
 46
 47
 48
 49
 50
 51
 52
 53
 54
 55
 56
 57
 58
 59
 60
 61
 62
 63
 64
 65
 66
 67
 68
 69
 70
 71
 72
 73
 74
 75
 76
 77
 78
 79
 80
 81
 82
 83
 84
 85
 86
 87
 88
 89
 90
 91
 92
 93
 94
 95
 96
 97
 98
 99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
from __future__ import print_function, division
#
import sys,os
os.environ['KMP_DUPLICATE_LIB_OK']='True' # uncomment this line if omp error occurs on OSX for python 3
os.environ['OMP_NUM_THREADS']='1' # set number of OpenMP threads to run in parallel
os.environ['MKL_NUM_THREADS']='1' # set number of MKL threads to run in parallel
#
quspin_path = os.path.join(os.getcwd(),"../../")
sys.path.insert(0,quspin_path)
###########################################################################
#                            example 16                                   #
#  In this script we demonstrate how to apply the user_basis to reduce    #
#  user-imported arrays of bases states (in integer representation)       #
#  to user-defined symmetry-reduced subspaces.                            #
###########################################################################
from quspin.basis import spin_basis_1d,spin_basis_general # Hilbert space spin basis_1d
from quspin.basis.user import user_basis # Hilbert space user basis
from quspin.basis.user import next_state_sig_32,op_sig_32,map_sig_32 # user basis data types
from numba import carray,cfunc # numba helper functions
from numba import uint32,int32 # numba data types
import numpy as np
from scipy.special import comb
#
#####
N_half = 10 # number of sites for each leg of the ladder
N = 2*N_half # total number of lattice sites
#
def make_basis(N_half):
    """ Generates a list of integers to represent external, user-imported basis """
    old_basis = spin_basis_general(N_half,m=0)
    #
    states = old_basis.states
    shift_states = np.left_shift(states,N_half)
    #
    shape=states.shape+states.shape
    #
    states_b = np.broadcast_to(states,shape)
    shift_states_b = np.broadcast_to(shift_states,shape)
    # this does the kronecker sum in a more memory efficient way. 
    return (states_b+shift_states_b.T).ravel()
#
external_basis = make_basis(N_half)
#
Np = () # dummy argument, could be any value (particle conservation should've been 
# taken into account when constructing the basis object)
#
######  function to call when applying operators
@cfunc(op_sig_32, locals=dict(s=int32,b=uint32))
def op(op_struct_ptr,op_str,ind,N,args):
    # using struct pointer to pass op_struct_ptr back to C++ see numba Records
    op_struct = carray(op_struct_ptr,1)[0]
    err = 0
    ind = N - ind - 1 # convention for QuSpin for mapping from bits to sites.
    s = (((op_struct.state>>ind)&1)<<1)-1
    b = (1<<ind)
    #
    if op_str==120: # "x" is integer value 120 (check with ord("x"))
        op_struct.state ^= b
    elif op_str==121: # "y" is integer value 120 (check with ord("y"))
        op_struct.state ^= b
        op_struct.matrix_ele *= 1.0j*s
    elif op_str==122: # "z" is integer value 120 (check with ord("z"))
        op_struct.matrix_ele *= s
    else:
        op_struct.matrix_ele = 0
        err = -1
    #
    return err
#
op_args=np.array([],dtype=np.uint32)
#
######  function to read user-imported basis into QuSpin 
#
# function to call when generating next_state
@cfunc(next_state_sig_32)
def next_state(s,counter,N,args):
    # return pre-calculated basis state.
    # add one to counter because the first state is already checked.
    return args[counter+1] # = basis
#
next_state_args = external_basis # this has to be an array of same dtype as the user_basis
#
class function_wrapper(object):
    """
    This class provides a wrapper for the user-imported basis,
    as well as the functions required for the `user_basis` functionality.
    #
    This is needed to easily pass parameters (defined as class attributes) to the
    functions `get_so_pcon()` and `get_Ns_pcon`.
    """
    def __init__(self,basis):
        self.basis = basis
    #
    # python function to calculate the starting state to generate the particle conserving basis
    def get_s0_pcon(self,N,Np):
        """ calculates the starting state to generate the particle conserving basis. """
        # ignore input arguments as basis is already calculated.
        return self.basis[0]
    # 
    # python function to calculate the size of the particle-conserved basis, 
    # i.e. BEFORE applying pre_check_state and symmetry maps
    def get_Ns_pcon(self,N,Np):
        """ calculates the size of the particle conservation basis (ignoring symmetries at this stage). """
        # ignore input arguments as basis is already calculated.
        return self.basis.size
#
######  define symmetry maps
#
if N_half!=10:
    print ("symmetry masks are hard-coded and work only for N=10; \n\
To do a different system size, it is required to update the masks accordingly.\n\
exiting...")
    exit()
#
@cfunc(map_sig_32)
def translation(x,N,sign_ptr,args):
    # bit permutation target bits: 1 2 3 4 5 6 7 8 9 0 11 12 13 14 15 16 17 18 19 10
    # code generated here: http://programming.sirrida.de/calcperm.php
    # only works for N_half=10 
    return ((x & 0x0007fdff) << 1) | ((x & 0x00080200) >> 9)
T_args=np.array([],dtype=np.uint32)
#
@cfunc(map_sig_32)
def parity(x,N,sign_ptr,args):
    # bit permutation target bits: 9 8 7 6 5 4 3 2 1 0 19 18 17 16 15 14 13 12 11 10
    # code generated here: http://programming.sirrida.de/calcperm.php
    # only works for N_half=10
    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))
P_args=np.array([],dtype=np.uint32)
#
######  construct user_basis 
# define maps dict
maps = dict(T_block=(translation,N_half,0,T_args), P_block=(parity,2,0,P_args), )
# define particle conservation and op dicts
FW = function_wrapper(external_basis)
pcon_dict = dict(Np=Np,next_state=next_state,next_state_args=next_state_args,
                 get_Ns_pcon=FW.get_Ns_pcon,get_s0_pcon=FW.get_s0_pcon)
op_dict = dict(op=op,op_args=op_args)
# create user basis
basis = user_basis(np.uint32,N,op_dict,allowed_ops=set("xyz"),sps=2,pcon_dict=pcon_dict,**maps)
# print basis
print(basis)