1from quspin.operators import hamiltonian, exp_op, quantum_operator  # operators
  2from quspin.basis import spinful_fermion_basis_1d  # Hilbert space basis
  3from quspin.tools.measurements import obs_vs_time  # calculating dynamics
  4import numpy as np  # general math functions
  5from numpy.random import uniform, choice  # tools for doing random sampling
  6from time import time  # tool for calculating computation time
  7import matplotlib.pyplot as plt  # plotting library
  8
  9#####################################################################
 10#                            example 6                              #
 11#   In this script we demonstrate how to use QuSpin's to create	    #
 12#   a disordered Fermi-Hubbard model with a parameter-dependent     #
 13#   Hamiltonian, and measure the imbalance on different lattice     #
 14#   sites (see arXiv:1501.05661). We also show how to prepare       #
 15#   fermion Fock states, and do disorder averaging.                 #
 16#####################################################################
 17##### setting parameters for simulation
 18# simulation parameters
 19n_real = 100  # number of realizations
 20n_boot = 100  # number of bootstrap samples to calculate error
 21# physical parameters
 22L = 8  # system size
 23N = L // 2  # number of particles
 24N_up = N // 2 + N % 2  # number of fermions with spin up
 25N_down = N // 2  # number of fermions with spin down
 26w_list = [1.0, 4.0, 10.0]  # disorder strength
 27J = 1.0  # hopping strength
 28U = 5.0  # interaction strength
 29# range in time to evolve system
 30start, stop, num = 0.0, 35.0, 101
 31t = np.linspace(start, stop, num=num, endpoint=True)
 32#
 33###### create the basis
 34# build spinful fermions basis
 35basis = spinful_fermion_basis_1d(L, Nf=(N_up, N_down))
 36#
 37##### create model
 38# define site-coupling lists
 39hop_right = [[J, i, i + 1] for i in range(L - 1)]  # hopping to the right OBC
 40hop_left = [[+J, i, i + 1] for i in range(L - 1)]  # hopping to the left OBC
 41int_list = [[U, i, i] for i in range(L)]  # onsite interaction
 42# site-coupling list to create the sublattice imbalance observable
 43sublat_list = [[(-1.0) ** i / N, i] for i in range(0, L)]
 44# create static lists
 45operator_list_0 = [
 46    ["+-|", hop_left],  # up hop left
 47    ["-+|", hop_right],  # up hop right
 48    ["|+-", hop_left],  # down hop left
 49    ["|-+", hop_right],  # down hop right
 50    ["n|n", int_list],  # onsite interaction
 51]
 52imbalance_list = [["n|", sublat_list], ["|n", sublat_list]]
 53# create operator dictionary for quantum_operator class
 54# add key for Hubbard hamiltonian
 55operator_dict = dict(H0=operator_list_0)
 56# add keys for local potential in each site
 57for i in range(L):
 58    # add to dictioanry keys h0,h1,h2,...,hL with local potential operator
 59    operator_dict["n" + str(i)] = [["n|", [[1.0, i]]], ["|n", [[1.0, i]]]]
 60#
 61###### setting up operators
 62# set up hamiltonian dictionary and observable (imbalance I)
 63no_checks = dict(check_pcon=False, check_symm=False, check_herm=False)
 64H_dict = quantum_operator(operator_dict, basis=basis, **no_checks)
 65I = hamiltonian(imbalance_list, [], basis=basis, **no_checks)
 66# strings which represent the initial state
 67s_up = "".join("1000" for i in range(N_up))
 68s_down = "".join("0010" for i in range(N_down))
 69# basis.index accepts strings and returns the index
 70# which corresponds to that state in the basis list
 71i_0 = basis.index(s_up, s_down)  # find index of product state
 72psi_0 = np.zeros(basis.Ns)  # allocate space for state
 73psi_0[i_0] = 1.0  # set MB state to be the given product state
 74print(
 75    "H-space size: {:d}, initial state: |{:s}>(x)|{:s}>".format(basis.Ns, s_up, s_down)
 76)
 77
 78
 79#
 80# define function to do dynamics for different disorder realizations.
 81def real(H_dict, I, psi_0, w, t, i):
 82    # body of function goes below
 83    ti = time()  # start timing function for duration of reach realisation
 84    # create a parameter list which specifies the onsite potential with disorder
 85    params_dict = dict(H0=1.0)
 86    for j in range(L):
 87        params_dict["n" + str(j)] = uniform(-w, w)
 88    # using the parameters dictionary construct a hamiltonian object with those
 89    # parameters defined in the list
 90    H = H_dict.tohamiltonian(params_dict)
 91    # use exp_op to get the evolution operator
 92    U = exp_op(H, a=-1j, start=t.min(), stop=t.max(), num=len(t), iterate=True)
 93    psi_t = U.dot(psi_0)  # get generator psi_t for time evolved state
 94    # use obs_vs_time to evaluate the dynamics
 95    t = U.grid  # extract time grid stored in U, and defined in exp_op
 96    obs_t = obs_vs_time(psi_t, t, dict(I=I))
 97    # print reporting the computation time for realization
 98    print("realization {}/{} completed in {:.2f} s".format(i + 1, n_real, time() - ti))
 99    # return observable values
100    return obs_t["I"].real
101
102
103#
104###### looping over different disorder strengths
105for w in w_list:
106    I_data = np.vstack([real(H_dict, I, psi_0, w, t, i) for i in range(n_real)])
107    ##### averaging and error estimation
108    I_avg = I_data.mean(axis=0)  # get mean value of I for all time points
109    # generate bootstrap samples
110    bootstrap_gen = (
111        I_data[choice(n_real, size=n_real)].mean(axis=0) for i in range(n_boot)
112    )
113    # generate the fluctuations about the mean of I
114    sq_fluc_gen = ((bootstrap - I_avg) ** 2 for bootstrap in bootstrap_gen)
115    I_error = np.sqrt(sum(sq_fluc_gen) / n_boot)
116    ##### plotting results
117    plt.errorbar(t, I_avg, I_error, marker=".", label="w={:.2f}".format(w))
118# configuring plots
119plt.xlabel("$Jt$", fontsize=18)
120plt.ylabel("$\\mathcal{I}$", fontsize=18)
121plt.grid(True)
122plt.tick_params(labelsize=16)
123plt.legend(loc=0)
124plt.tight_layout()
125plt.savefig("fermion_MBL.pdf", bbox_inches="tight")
126# plt.show()
127plt.close()