cleanup part II

Project.toml

@@ -18,8 +18,7 @@ TensorOperations = "3"
 julia = "1.6 - 1.9"
 
 [extras]
-Plots = "91a5bcdd-55d7-5caf-9e0b-520d859cae80"
 Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
 
 [targets]
-test = ["Test", "Plots"]
+test = ["Test"]

src/models.jl

@@ -71,45 +71,3 @@ function sixvertex(eltype, ::CU₁, lattice=InfiniteChain(1);
  blocks(mpo)[Irrep[CU₁](1, 2)] = reshape([a], (1, 1))
  return InfiniteMPO(permute(mpo, (1, 2), (4, 3)))
 end
-
-"""
- xxx(; J=1.0, spin=1)
-
-MPO for the hamiltonian of the xxx Heisenberg model, defined by
- ``H = J(∑_{<i,j>} X_i X_j + Y_i Y_j + Z_i Z_j)``
-"""
-function xxx(eltype=ComplexF64, symmetry=ℤ{1}, lattice=InfiniteChain(1);
- J=1.0, spin=1)
- SS = sigma_exchange(eltype, symmetry; spin=spin)
- return @mpoham sum(J * SS{i,j} for (i, j) in nearest_neighbours(lattice))
-end
-
-"""
- xxz(; J=1.0, Δ=1.0, spin=1)
-
-MPO for the hamiltonian of the xxz Heisenberg model, defined by
- ``H = J(∑_{<i,j>} X_i X_j + Y_i Y_j + Δ Z_i Z_j)``
-"""
-function xxz(eltype=ComplexF64, symmetry=ℤ{1}, lattice=InfiniteChain(1);
- J=1.0, Δ=1.0, spin=1)
- XX = sigma_xx(eltype, symmetry; spin=spin)
- YY = sigma_yy(eltype, symmetry; spin=spin)
- ZZ = sigma_zz(eltype, symmetry; spin=spin)
- return @mpoham sum(J * (XX{i,j} + YY{i,j} + Δ * ZZ{i,j})
- for (i, j) in nearest_neighbours(lattice))
-end
-
-"""
- xxz(; J=1.0, Δ=1.0, spin=1)
-
-MPO for the hamiltonian of the xyz Heisenberg model, defined by
- ``H = J(∑_{<i,j>} J_x X_{i}X_{j} + J_y Y_{i}Y_{j} + J_z Z_{i}Z_{j})``
-"""
-function xyz(eltype=ComplexF64, symmetry=ℤ{1}, lattice=InfiniteChain(1);
- Jx=1.0, Jy=1.0, Jz=1.0, spin=1)
- XX = sigma_xx(eltype, symmetry; spin=spin)
- YY = sigma_yy(eltype, symmetry; spin=spin)
- ZZ = sigma_zz(eltype, symmetry; spin=spin)
- return @mpoham sum(Jx * XX{i,j} + Jy * YY{i,j} + Jz * ZZ{i,j}
- for (i, j) in nearest_neighbours(lattice))
-end

src/models/hamiltonians.jl

@@ -1,3 +1,7 @@
+#===========================================================================================
+ Ising model
+===========================================================================================#
+
 """
  transverse_field_ising(; J=1.0, hx=1.0, hz=0.0, spin=1//2)
 
@@ -12,10 +16,56 @@ function transverse_field_ising(eltype=ComplexF64, symmetry=ℤ{1},
  if symmetry != ℤ{1}
  @assert hz == zero(hz) "parameters and symmetry incompatible"
  return @mpoham sum(-J * (ZZ{i,j} + hx * X{i})
- for (i, j) in nearest_neighbours(lattice))
+  for (i, j) in nearest_neighbours(lattice))
  else
  Z = sigma_z(eltype, symmetry; spin=spin)
  return @mpoham sum(-J * (ZZ{i,j} + hx * X{i} + hz * Z{i})
- for (i, j) in nearest_neighbours(lattice))
+  for (i, j) in nearest_neighbours(lattice))
  end
-end
+end
+
+#===========================================================================================
+ Heisenberg models
+===========================================================================================#
+
+"""
+ xxx(; J=1.0, spin=1)
+
+MPO for the hamiltonian of the xxx Heisenberg model, defined by
+ ``H = J(∑_{<i,j>} X_i X_j + Y_i Y_j + Z_i Z_j)``
+"""
+function xxx(eltype=ComplexF64, symmetry=ℤ{1}, lattice=InfiniteChain(1);
+ J=1.0, spin=1)
+ SS = sigma_exchange(eltype, symmetry; spin=spin)
+ return @mpoham sum(J * SS{i,j} for (i, j) in nearest_neighbours(lattice))
+end
+
+"""
+ xxz(; J=1.0, Δ=1.0, spin=1)
+
+MPO for the hamiltonian of the xxz Heisenberg model, defined by
+ ``H = J(∑_{<i,j>} X_i X_j + Y_i Y_j + Δ Z_i Z_j)``
+"""
+function xxz(eltype=ComplexF64, symmetry=ℤ{1}, lattice=InfiniteChain(1);
+ J=1.0, Δ=1.0, spin=1)
+ XX = sigma_xx(eltype, symmetry; spin=spin)
+ YY = sigma_yy(eltype, symmetry; spin=spin)
+ ZZ = sigma_zz(eltype, symmetry; spin=spin)
+ return @mpoham sum(J * (XX{i,j} + YY{i,j} + Δ * ZZ{i,j})
+ for (i, j) in nearest_neighbours(lattice))
+end
+
+"""
+ xyz(; Jx=1.0, Jy=1.0, Jz=1.0, spin=1)
+
+MPO for the hamiltonian of the xyz Heisenberg model, defined by
+ ``H = J(∑_{<i,j>} J_x X_{i}X_{j} + J_y Y_{i}Y_{j} + J_z Z_{i}Z_{j})``
+"""
+function xyz(eltype=ComplexF64, symmetry=ℤ{1}, lattice=InfiniteChain(1);
+ Jx=1.0, Jy=1.0, Jz=1.0, spin=1)
+ XX = sigma_xx(eltype, symmetry; spin=spin)
+ YY = sigma_yy(eltype, symmetry; spin=spin)
+ ZZ = sigma_zz(eltype, symmetry; spin=spin)
+ return @mpoham sum(Jx * XX{i,j} + Jy * YY{i,j} + Jz * ZZ{i,j}
+ for (i, j) in nearest_neighbours(lattice))
+end

src/mpoham.jl

@@ -78,8 +78,8 @@ function _is_free_channel(data, loc, channel)
  return all(ismissing.(data[mod1(loc, end), :, channel]))
 end
 
-function MPSKit.MPOHamiltonian(o::LocalOperator, unitcell=minimum(o.inds))
- return MPOHamiltonian(SumOfLocalOperators((o,)), unitcell)
+function MPSKit.MPOHamiltonian(o::LocalOperator)
+ return MPOHamiltonian(SumOfLocalOperators((o,)))
 end
 
 function MPSKit.MPOHamiltonian(opps::SumOfLocalOperators)
@@ -107,7 +107,7 @@ function MPSKit.MPOHamiltonian(opps::SumOfLocalOperators)
  hit, data = _find_free_channel(data, start)
 
  data[mod1(start, end), 1, hit] = mpo[1]
- for site in range(start + 1, stop - 1)
+ for site in (start + 1):(stop - 1)
  mpo_ind = findfirst(linds .== site)
  o = isnothing(mpo_ind) ? 1 : mpo[mpo_ind]
 

src/utility.jl

@@ -25,3 +25,7 @@ function nonsym_bosonictensors(cutoff::Int, elt=ComplexF64)
  a⁻ = TensorMap(collect(creadat'), ℂ^(cutoff + 1), ℂ^(cutoff + 1))
  return (a⁺, a⁻)
 end
+
+@static if !isdefined(Base, :allequal)
+ allequal(itr) = isempty(itr) ? true : all(isequal(first(itr)), itr)
+end

test/heisenberg.jl

@@ -0,0 +1,19 @@
+using MPSKitModels
+using MPSKit
+using TensorKit
+using Test
+
+alg = VUMPS()
+
+@testset "xxx" begin
+ E₀ = -1.401484014561
+ E₁ = 0.41047925
+ H = xxx()
+ Ψ = InfiniteMPS([ComplexSpace(3)], [ComplexSpace(48)])
+ @test sum(abs.(imag.(expectation_value(Ψ, H)))) ≈ 0 atol = 1e-10
+ Ψ, envs, δ = find_groundstate(Ψ, H, alg)
+ @test E₀ ≈ sum(expectation_value(Ψ, H, envs)) atol = 1e-4
+
+ ΔEs, qps = excitations(H, QuasiparticleAnsatz(), Float64(pi), Ψ, envs)
+ @test E₁ ≈ first(ΔEs) atol = 1e-4
+end

test/runtests.jl

@@ -8,14 +8,6 @@ end
  include("tfim.jl")
 end
 
-using Test, TensorKit, MPSKit, MPSKitModels
-
-@testset "xxz" begin
- @testset "nonsym_xxz_ham" begin
- th = nonsym_xxz_ham()
- ts = InfiniteMPS([ℂ^3], [ℂ^48])
- (ts, pars, _) = find_groundstate(ts, th, VUMPS(; maxiter=400, verbose=false))
- (energies, Bs) = excitations(th, QuasiparticleAnsatz(), Float64(pi), ts, pars)
- @test energies[1] ≈ 0.41047925 atol = 1e-4
- end
+@testset "heisenberg model" begin
+ include("heisenberg.jl")
 end

test/tfim.jl

@@ -1,32 +1,35 @@
 using MPSKitModels
 using MPSKit
+using TensorKit
 using Test
 
-const E₀ = -0.318309846883
+E₀ = -0.318309846883
+alg = VUMPS()
+
 
 @testset "single site" begin
  H = transverse_field_ising()
- Ψ₀ = InfiniteMPS(2, 16)
+ Ψ₀ = InfiniteMPS([ComplexSpace(2)], [ComplexSpace(16)])
  @test sum(abs.(imag.(expectation_value(Ψ₀, H)))) ≈ 0 atol = 1e-10
- Ψ, envs, δ = find_groundstate(Ψ₀, H; tol=1e-8, maxiter=400, verbose=false)
+ Ψ, envs, δ = find_groundstate(Ψ₀, H, alg)
  @test E₀ ≈ sum(expectation_value(Ψ, H, envs))
 end
 
 @testset "double sites" begin
  H2 = transverse_field_ising(ComplexF64, ℤ{1}, InfiniteChain(2))
- Ψ2 = InfiniteMPS([2, 2], [16, 16])
+ Ψ2 = InfiniteMPS(ComplexSpace.([2, 2]), ComplexSpace.([16, 16]))
 
  @test sum(abs.(imag.(expectation_value(Ψ2, H2)))) ≈ 0 atol = 1e-10
- Ψ2, envs, δ = find_groundstate(Ψ2, H2; tol=1e-8, maxiter=400, verbose=false)
+ Ψ2, envs, δ = find_groundstate(Ψ2, H2, alg)
  @test 2E₀ ≈ sum(expectation_value(Ψ2, H2, envs))
 end
 
 @testset "weird lattice" begin
  lattice = SnakePattern(InfiniteChain(2), i -> iseven(i) ? i - 1 : i + 1)
  H2 = transverse_field_ising(ComplexF64, ℤ{1}, lattice)
- Ψ2 = InfiniteMPS([2, 2], [16, 16])
+ Ψ2 = InfiniteMPS(ComplexSpace.([2, 2]), ComplexSpace.([16, 16]))
 
  @test sum(abs.(imag.(expectation_value(Ψ2, H2)))) ≈ 0 atol = 1e-10
- Ψ2, envs, δ = find_groundstate(Ψ2, H2; tol=1e-8, maxiter=400, verbose=false)
+ Ψ2, envs, δ = find_groundstate(Ψ2, H2, alg)
  @test 2E₀ ≈ sum(expectation_value(Ψ2, H2, envs))
 end