1#####################################################################
  2#                   example 5 (DEPRECATED)                          #
  3#    In this script we demonstrate how to use QuSpin's to build     #
  4#    the Hamiltonian of the SSH model in real and momentum space.   #
  5#    Along the way, we showcase the block tools which allow the     #
  6#    user to create block-diagonal Hamiltonians. Last, we show      #
  7#    how to time-evolve free fermion states like the Fermi sea      #
  8#    and measure correlators.                                       #
  9#####################################################################
 10from quspin.operators import hamiltonian, exp_op  # Hamiltonians and operators
 11from quspin.basis import spinless_fermion_basis_1d  # Hilbert space fermion basis
 12from quspin.tools.block_tools import block_diag_hamiltonian  # block diagonalisation
 13import numpy as np  # generic math functions
 14import matplotlib.pyplot as plt  # plotting library
 15
 16try:  # import python 3 zip function in python 2 and pass if already using python 3
 17    import itertools.izip as zip
 18except ImportError:
 19    pass
 20##### define model parameters #####
 21L = 100  # system size
 22J = 1.0  # uniform hopping
 23deltaJ = 0.1  # bond dimerisation
 24Delta = 0.5  # staggered potential
 25beta = 100.0  # inverse temperature for Fermi-Dirac distribution
 26##### construct single-particle Hamiltonian #####
 27# define site-coupling lists
 28hop_pm = [[-J - deltaJ * (-1) ** i, i, (i + 1) % L] for i in range(L)]  # PBC
 29hop_mp = [[+J + deltaJ * (-1) ** i, i, (i + 1) % L] for i in range(L)]  # PBC
 30stagg_pot = [[Delta * (-1) ** i, i] for i in range(L)]
 31# define static and dynamic lists
 32static = [["+-", hop_pm], ["-+", hop_mp], ["n", stagg_pot]]
 33dynamic = []
 34# define basis
 35basis = spinless_fermion_basis_1d(L, Nf=1)
 36# build real-space Hamiltonian
 37H = hamiltonian(static, dynamic, basis=basis, dtype=np.float64)
 38# diagonalise real-space Hamiltonian
 39E, V = H.eigh()
 40##### compute Fourier transform and momentum-space Hamiltonian #####
 41# define momentm blocks and basis arguments
 42blocks = [
 43    dict(Nf=1, kblock=i, a=2) for i in range(L // 2)
 44]  # only L//2 distinct momenta
 45basis_args = (L,)
 46# construct block-diagonal Hamiltonian
 47FT, Hblock = block_diag_hamiltonian(
 48    blocks,
 49    static,
 50    dynamic,
 51    spinless_fermion_basis_1d,
 52    basis_args,
 53    np.complex128,
 54    get_proj_kwargs=dict(pcon=True),
 55)
 56# diagonalise momentum-space Hamiltonian
 57Eblock, Vblock = Hblock.eigh()
 58##### prepare the density observables and initial states #####
 59# grab single-particle states and treat them as initial states
 60psi0 = Vblock
 61# construct operator n_1 = $n_{j=0}$
 62n_1_static = [["n", [[1.0, 0]]]]
 63n_1 = hamiltonian(
 64    n_1_static, [], basis=basis, dtype=np.float64, check_herm=False, check_pcon=False
 65)
 66# construct operator n_2 = $n_{j=L/2}$
 67n_2_static = [["n", [[1.0, L // 2]]]]
 68n_2 = hamiltonian(
 69    n_2_static, [], basis=basis, dtype=np.float64, check_herm=False, check_pcon=False
 70)
 71# transform n_j operators to momentum space
 72n_1 = n_1.rotate_by(FT, generator=False)
 73n_2 = n_2.rotate_by(FT, generator=False)
 74##### evaluate nonequal time correlator <FS|n_2(t) n_1(0)|FS> #####
 75# define time vector
 76t = np.linspace(0.0, 90.0, 901)
 77# calcualte state acted on by n_1
 78n_psi0 = n_1.dot(psi0)
 79# construct time-evolution operator using exp_op class (sometimes faster)
 80U = exp_op(Hblock, a=-1j, start=t.min(), stop=t.max(), num=len(t), iterate=True)
 81# evolve states
 82psi_t = U.dot(psi0)
 83n_psi_t = U.dot(n_psi0)
 84# alternative method for time evolution using Hamiltonian class
 85# psi_t=Hblock.evolve(psi0,0.0,t,iterate=True)
 86# n_psi_t=Hblock.evolve(n_psi0,0.0,t,iterate=True)
 87# preallocate variable
 88correlators = np.zeros(t.shape + psi0.shape[1:])
 89# loop over the time-evolved states
 90for i, (psi, n_psi) in enumerate(zip(psi_t, n_psi_t)):
 91    correlators[i, :] = n_2.matrix_ele(psi, n_psi, diagonal=True).real
 92# evaluate correlator at finite temperature
 93n_FD = 1.0 / (np.exp(beta * E) + 1.0)
 94correlator = (n_FD * correlators).sum(axis=-1)
 95##### plot spectra
 96plt.plot(np.arange(H.Ns), E / L, marker="o", color="b", label="real space")
 97plt.plot(
 98    np.arange(Hblock.Ns),
 99    Eblock / L,
100    marker="x",
101    color="r",
102    markersize=2,
103    label="momentum space",
104)
105plt.xlabel("state number", fontsize=16)
106plt.ylabel("energy", fontsize=16)
107plt.xticks(fontsize=16)
108plt.yticks(fontsize=16)
109plt.legend(fontsize=16)
110plt.grid()
111plt.tight_layout()
112plt.savefig("example5a.pdf", bbox_inches="tight")
113# plt.show()
114plt.close()
115##### plot correlator
116plt.plot(t, correlator, linewidth=2)
117plt.xlabel("$t$", fontsize=16)
118plt.ylabel("$C_{0,L/2}(t,\\beta)$", fontsize=16)
119plt.xticks(fontsize=16)
120plt.yticks(fontsize=16)
121plt.grid()
122plt.tight_layout()
123plt.savefig("example5b.pdf", bbox_inches="tight")
124# plt.show()
125plt.close()