Exact Diagonalisation of Spin Hamiltonians
This example shows how to code up the Heisenberg Hamiltonian:
\[H = \sum_{j=0}^{L-2}\frac{J_{xy}}{2}\left(S^+_{j+1}S^-_{j} + \mathrm{h.c.}\right) + J_{zz}S^z_{j+1}S^z_{j} + h_z\sum_{j=0}^{L-1}S^z_{j}.\]
Details about the code below can be found in SciPost Phys. 2, 003 (2017).
Script
1#####################################################################
2# example 0 #
3# In this script we demonstrate how to use QuSpin's exact #
4# diagonlization routines to solve for the eigenstates and #
5# energies of the XXZ chain. #
6#####################################################################
7from quspin.operators import hamiltonian # Hamiltonians and operators
8from quspin.basis import spin_basis_1d # Hilbert space spin basis
9import numpy as np # generic math functions
10
11#
12##### define model parameters #####
13L = 12 # system size
14Jxy = np.sqrt(2.0) # xy interaction
15Jzz_0 = 1.0 # zz interaction
16hz = 1.0 / np.sqrt(3.0) # z external field
17#
18##### set up Heisenberg Hamiltonian in an external z-field #####
19# compute spin-1/2 basis
20# basis = spin_basis_1d(L,pauli=False)
21# basis = spin_basis_1d(L,pauli=False,Nup=L//2) # zero magnetisation sector
22basis = spin_basis_1d(
23 L, pauli=False, Nup=L // 2, pblock=1
24) # and positive parity sector
25# define operators with OBC using site-coupling lists
26J_zz = [[Jzz_0, i, i + 1] for i in range(L - 1)] # OBC
27J_xy = [[Jxy / 2.0, i, i + 1] for i in range(L - 1)] # OBC
28h_z = [[hz, i] for i in range(L)]
29# static and dynamic lists
30static = [["+-", J_xy], ["-+", J_xy], ["zz", J_zz], ["z", h_z]]
31dynamic = []
32# compute the time-dependent Heisenberg Hamiltonian
33H_XXZ = hamiltonian(static, dynamic, basis=basis, dtype=np.float64)
34#
35##### various exact diagonalisation routines #####
36# calculate entire spectrum only
37E = H_XXZ.eigvalsh()
38# calculate full eigensystem
39E, V = H_XXZ.eigh()
40# calculate minimum and maximum energy only
41Emin, Emax = H_XXZ.eigsh(k=2, which="BE", maxiter=1e4, return_eigenvectors=False)
42# calculate the eigenstate closest to energy E_star
43E_star = 0.0
44E, psi_0 = H_XXZ.eigsh(k=1, sigma=E_star, maxiter=1e4)
45psi_0 = psi_0.reshape((-1,))