1#####################################################################
  2#                   example 7  (DEPRECATED)                         #
  3#   In this script we demonstrate how to use QuSpin's to create	    #
  4#   a ladder geometry and study quanch dynamics in the Bose-Hubbard #
  5#   model. We also show how to compute the entanglement entropy of  #
  6#   bosonic systems. Last, we demonstrate how to use the block_ops  #
  7#   tool to decompose a state in the symemtry sectors of a          #
  8#   Hamiltonian, evolve the separate parts, and put back the state  #
  9#   in the end.                                                     #
 10#####################################################################
 11from quspin.operators import hamiltonian  # Hamiltonians and operators
 12from quspin.basis import boson_basis_1d  # bosonic Hilbert space
 13from quspin.tools.block_tools import block_ops  # dynamics in symmetry blocks
 14import numpy as np  # general math functions
 15import matplotlib.pyplot as plt  # plotting library
 16import matplotlib.animation as animation  # animating movie of dynamics
 17
 18#
 19##### define model parameters
 20# initial seed for random number generator
 21np.random.seed(0)  # seed is 0 to produce plots from QuSpin2 paper
 22# setting up parameters of simulation
 23L = 6  # length of chain
 24N = 2 * L  # number of sites
 25nb = 0.5  # density of bosons
 26sps = 3  # number of states per site
 27J_par_1 = 1.0  # top side of ladder hopping
 28J_par_2 = 1.0  # bottom side of ladder hopping
 29J_perp = 0.5  # rung hopping
 30U = 20.0  # Hubbard interaction
 31#
 32##### set up Hamiltonian and observables
 33# define site-coupling lists
 34int_list_1 = [[-0.5 * U, i] for i in range(N)]  # interaction $-U/2 \sum_i n_i$
 35int_list_2 = [[0.5 * U, i, i] for i in range(N)]  # interaction: $U/2 \num_i n_i^2$
 36# setting up hopping lists
 37hop_list = [[-J_par_1, i, (i + 2) % N] for i in range(0, N, 2)]  # PBC bottom leg
 38hop_list.extend([[-J_par_2, i, (i + 2) % N] for i in range(1, N, 2)])  # PBC top leg
 39hop_list.extend([[-J_perp, i, i + 1] for i in range(0, N, 2)])  # perp/rung hopping
 40hop_list_hc = [[J.conjugate(), i, j] for J, i, j in hop_list]  # add h.c. terms
 41# set up static and dynamic lists
 42static = [
 43    ["+-", hop_list],  # hopping
 44    ["-+", hop_list_hc],  # hopping h.c.
 45    ["nn", int_list_2],  # U n_i^2
 46    ["n", int_list_1],  # -U n_i
 47]
 48dynamic = []  # no dynamic operators
 49# create block_ops object
 50blocks = [dict(kblock=kblock) for kblock in range(L)]  # blocks to project on to
 51baisis_args = (N,)  # boson_basis_1d manditory arguments
 52basis_kwargs = dict(nb=nb, sps=sps, a=2)  # boson_basis_1d optional args
 53get_proj_kwargs = dict(pcon=True)  # set projection to full particle basis
 54H_block = block_ops(
 55    blocks,
 56    static,
 57    dynamic,
 58    boson_basis_1d,
 59    baisis_args,
 60    np.complex128,
 61    basis_kwargs=basis_kwargs,
 62    get_proj_kwargs=get_proj_kwargs,
 63)
 64# setting up local Fock basis
 65basis = boson_basis_1d(N, nb=nb, sps=sps)
 66# setting up observables
 67no_checks = dict(check_herm=False, check_symm=False, check_pcon=False)
 68n_list = [
 69    hamiltonian([["n", [[1.0, i]]]], [], basis=basis, dtype=np.float64, **no_checks)
 70    for i in range(N)
 71]
 72##### do time evolution
 73# set up initial state
 74i0 = np.random.randint(basis.Ns)  # pick random state from basis set
 75psi = np.zeros(basis.Ns, dtype=np.float64)
 76psi[i0] = 1.0
 77# print info about setup
 78state_str = "".join(
 79    str(int((basis[i0] // basis.sps ** (L - i - 1))) % basis.sps) for i in range(N)
 80)
 81print("total H-space size: {}, initial state: |{}>".format(basis.Ns, state_str))
 82# setting up parameters for evolution
 83start, stop, num = 0, 30, 301  # 0.1 equally spaced points
 84times = np.linspace(start, stop, num)
 85# calculating the evolved states
 86n_jobs = 1  # paralelisation: increase to see if calculation runs faster!
 87psi_t = H_block.expm(
 88    psi, a=-1j, start=start, stop=stop, num=num, block_diag=False, n_jobs=n_jobs
 89)
 90# calculating the local densities as a function of time
 91expt_n_t = np.vstack([n.expt_value(psi_t).real for n in n_list]).T
 92# reshape data for plotting
 93n_t = np.zeros((num, 2, L))
 94n_t[:, 0, :] = expt_n_t[:, 0::2]
 95n_t[:, 1, :] = expt_n_t[:, 1::2]
 96# calculating entanglement entropy
 97sub_sys_A = range(0, N, 2)  # bottom side of ladder
 98gen = (basis.ent_entropy(psi, sub_sys_A=sub_sys_A)["Sent_A"] for psi in psi_t.T[:])
 99ent_t = np.fromiter(gen, dtype=np.float64, count=num)
100# plotting static figures
101# """
102fig, ax = plt.subplots(nrows=5, ncols=1)
103im = []
104im_ind = []
105for i, t in enumerate(np.logspace(-1, np.log10(stop - 1), 5, base=10)):
106    j = times.searchsorted(t)
107    im_ind.append(j)
108    im.append(ax[i].imshow(n_t[j], cmap="hot", vmax=n_t.max(), vmin=0))
109    ax[i].tick_params(labelbottom=False, labelleft=False)
110cax = fig.add_axes([0.85, 0.1, 0.03, 0.8])
111fig.colorbar(im[2], cax)
112plt.savefig("boson_density.pdf")
113plt.figure()
114plt.plot(times, ent_t, lw=2)
115plt.plot(times[im_ind], ent_t[im_ind], marker="o", linestyle="", color="red")
116plt.xlabel("$Jt$", fontsize=20)
117plt.ylabel("$s_\\mathrm{ent}(t)$", fontsize=20)
118plt.grid()
119plt.savefig("boson_entropy.pdf")
120# plt.show()
121plt.close()
122# """
123# setting up two plots to animate side by side
124fig, (ax1, ax2) = plt.subplots(1, 2)
125fig.set_size_inches(10, 5)
126ax1.set_xlabel(r"$Jt$", fontsize=18)
127ax1.set_ylabel(r"$s_\mathrm{ent}$", fontsize=18)
128ax1.grid()
129(line1,) = ax1.plot(times, ent_t, lw=2)
130line1.set_data([], [])
131im = ax2.matshow(n_t[0], cmap="hot")
132fig.colorbar(im)
133
134
135def run(i):  # function to update frame
136    # set new data for plots
137    if i == num - 1:
138        exit()  # comment this line to retain last plot
139    line1.set_data(times[:i], ent_t[:i])
140    im.set_data(n_t[i])
141    return im, line1
142
143
144# define and display animation
145ani = animation.FuncAnimation(fig, run, range(num), interval=50, repeat=False)
146plt.show()
147plt.close()
148#
149""" 
150###### ladder lattice 
151# hopping coupling parameters:
152# - : J_par_1
153# = : J_par_2
154# | : J_perp
155#
156# lattice graph
157#
158 = 1 = 3 = 5 = 7 = 9 =
159   |   |   |   |   |
160 - 0 - 2 - 4 - 6 - 8 -
161#
162# translations along leg-direction (i -> i+2):
163#
164 = 9 = 1 = 3 = 5 = 7 =
165 |   |   |   |   | 
166 - 8 - 0 - 2 - 4 - 6 -
167#
168# if J_par_1=J_par_2, one can use regular chain parity (i -> N - i) as combination 
169# of the two ladder parity operators:
170#
171 - 8 - 6 - 4 - 2 - 0 - 
172 |   |   |   |   | 
173 - 9 - 7 - 5 - 3 - 1 -
174"""