Example 1: S=1/2 Heisenberg Chain

using SparseArrays
using LinearAlgebra
using Arpack
# using Plots

using LatticeTools
using QuantumHamiltonian

println("# S=1/2 Heisenberg Chain")
n_sites = 8;
(hs, σ) = QuantumHamiltonian.Toolkit.spin_half_system(n_sites)
println("Quantum number sectors (2Sz) : ", quantum_number_sectors(hs))

S = Dict(μ => sum(σ(i, μ) for i in 1:n_sites) for μ in [:x, :y, :z])
spin_squared = simplify(S[:x]^2 + S[:y]^2 + S[:z]^2)
j1 = simplify(sum(σ(i, j) * σ(mod(i, n_sites)+1 , j) for i in 1:n_sites for j in [:x, :y, :z]))
hs_rep = represent(hs);     # Use FrozenSortedArrayIndex{UInt} for basis lookup

println("## All sectors at once")
# plt = plot(size=(400, 400))
begin
  j1_rep = represent(hs_rep, j1)
  eigenvalues, eigenvectors = eigs(j1_rep; nev=32, which=:SR, )
  eigenvalues = real.(eigenvalues)
  println("E : ", eigenvalues[1:5])
  # scatter!(plt,
  #          zeros(size(eigenvalues)), eigenvalues,
  #          markershape=:hline,
  #          markersize=10,
  #          markerstrokecolor=:red,
  #          legend=:none)
end

println("## Sz=0, each momentum sectors")
hs_sector = HilbertSpaceSector(hs, 0)
hs_rep = represent_dict(hs_sector) # Use Dict{UInt, Int} for basis lookup

unitcell = makeunitcell(1.0; SiteType=String)
addsite!(unitcell, "Spin", FractCoord([0], [0.0]))
lattice = make_lattice(unitcell, 8)
tsymbed = translation_symmetry_embedding(lattice)

#translation_group = TranslationGroup([Permutation([ mod(i, n_sites)+1 for i in 1:n_sites])])
ks = symmetry(tsymbed).fractional_momenta

for tsic in get_irrep_components(tsymbed)
  k = ks[tsic.irrep_index]
  hs_redrep = symmetry_reduce(hs_rep, tsic)
  j1_redrep = represent(hs_redrep, j1)
  j1_redrep_dense = Matrix(j1_redrep)   # Make a dense matrix
  eigenvalues = eigvals(Hermitian(j1_redrep_dense))
  println("E(k=", join(string.(k), ","), ") : ", eigenvalues[1:5])
  # scatter!(plt,
  #          ones(size(eigenvalues)).*tsic.irrep_index, eigenvalues,
  #          markershape=:hline,
  #          markersize=10,
  #          markerstrokecolor=:blue,
  #          legend=:none)
end
# xticks!(plt, collect(0:n_sites), ["All", string.(1:n_sites)...])
# xlims!(plt, -1, n_sites+1)
# xlabel!(plt, "Momentum index")
# ylabel!(plt, "Energy")
# savefig(plt, "spinchain.svg"); nothing

# S=1/2 Heisenberg Chain
Quantum number sectors (2Sz) : [(-8,), (-6,), (-4,), (-2,), (0,), (2,), (4,), (6,), (8,)]
## All sectors at once
E : [-14.604373635748718, -12.513676255378321, -12.513676255378291, -12.51367625537829, -10.798512593101258]
## Sz=0, each momentum sectors
┌ Warning: Use of HilbertSpaceSector deprecated
└ @ QuantumHamiltonian ~/work/QuantumHamiltonian.jl/QuantumHamiltonian.jl/src/HilbertSpace/hilbert_space_sector.jl:35
┌ Warning: represent_dict deprecated. Use reprent(hs, BR, BT) instead.
└ @ QuantumHamiltonian ~/work/QuantumHamiltonian.jl/QuantumHamiltonian.jl/src/Representation/hilbert_space_representation.jl:178
┌ Warning: Assignment to `eigenvalues` in soft scope is ambiguous because a global variable by the same name exists: `eigenvalues` will be treated as a new local. Disambiguate by using `local eigenvalues` to suppress this warning or `global eigenvalues` to assign to the existing global variable.
└ @ spinhalf.md:54
┌ Warning: Use of IrrepComponent with symmetry_reduce deprecated. Use simpler symmetry_reduce function that takes symops and amplitudes as input
└ @ QuantumHamiltonian ~/work/QuantumHamiltonian.jl/QuantumHamiltonian.jl/src/Symmetry/deprecated_symmetry_reduce.jl:13
E(k=0//1) : [-14.604373635748695, -7.207750943219352, -3.999999999999999, -2.9044377801431303, -1.7801674716505156]
E(k=1//8) : [-9.83495403557933, -4.576491222541476, -2.5135137075322476, -2.0507328122554997, -1.1715728752538104]
E(k=1//4) : [-7.418550718873845, -6.4721359549995805, -5.068140393445468, -2.3878731329492617, -1.034608090016609]
E(k=3//8) : [-8.580593495682868, -6.828427124746193, -3.435694198841033, -1.748064097795283, 1.061810846301981]
E(k=1//2) : [-12.513676255378314, -10.798512593101274, -4.806558702893607, -3.0435068869737827, -8.388394747161498e-16]
E(k=5//8) : [-8.580593495682864, -6.828427124746194, -3.4356941988410292, -1.748064097795285, 1.0618108463019813]
E(k=3//4) : [-7.418550718873845, -6.4721359549995805, -5.068140393445468, -2.3878731329492617, -1.034608090016609]
E(k=7//8) : [-9.834954035579335, -4.576491222541476, -2.5135137075322476, -2.0507328122554997, -1.17157287525381]