diff --git a/src/algorithms/ED.jl b/src/algorithms/ED.jl index c4ccd2e6a..36c112630 100644 --- a/src/algorithms/ED.jl +++ b/src/algorithms/ED.jl @@ -31,47 +31,51 @@ equivalent to dense eigenvectors. """ function exact_diagonalization(H::FiniteMPOHamiltonian; - sector=first(sectors(oneunit(physicalspace(H, 1)))), - len::Int=length(H), num::Int=1, which::Symbol=:SR, + sector=one(sectortype(H)), + num::Int=1, which::Symbol=:SR, alg=Defaults.alg_eigsolve(; dynamic_tols=false)) - left = ℂ[typeof(sector)](sector => 1) - right = oneunit(left) + L = length(H) + @assert L > 1 "FiniteMPOHamiltonian must have length > 1" + middle_site = (L >> 1) + 1 - middle_site = Int(round(len / 2)) - - Ot = eltype(H[1]) - - mpst_type = tensormaptype(spacetype(Ot), 2, 1, storagetype(Ot)) - mpsb_type = tensormaptype(spacetype(Ot), 1, 1, storagetype(Ot)) - Cs = Vector{Union{Missing,mpsb_type}}(missing, len + 1) - ALs = Vector{Union{Missing,mpst_type}}(missing, len) - ARs = Vector{Union{Missing,mpst_type}}(missing, len) - ACs = Vector{Union{Missing,mpst_type}}(missing, len) + T = storagetype(eltype(H)) + TA = tensormaptype(spacetype(H), 2, 1, T) + # fuse from left to right + ALs = Vector{Union{Missing,TA}}(missing, L) + left = oneunit(spacetype(H)) for i in 1:(middle_site - 1) - ALs[i] = isomorphism(storagetype(Ot), left * physicalspace(H, i), - fuse(left * physicalspace(H, i))) - left = _lastspace(ALs[i])' + P = physicalspace(H, i) + ALs[i] = isomorphism(T, left ⊗ P ← fuse(left ⊗ P)) + left = right_virtualspace(ALs[i]) end - for i in len:-1:(middle_site + 1) - ARs[i] = _transpose_front(isomorphism(storagetype(Ot), - fuse(right * physicalspace(H, i)'), - right * physicalspace(H, i)')) - right = _firstspace(ARs[i]) + + # fuse from right to left + ARs = Vector{Union{Missing,TA}}(missing, L) + right = spacetype(H)(sector => 1) + for i in reverse((middle_site + 1):L) + P = physicalspace(H, i) + ARs[i] = _transpose_front(isomorphism(T, fuse(right ⊗ P') ← right ⊗ P')) + right = left_virtualspace(ARs[i]) end - ACs[middle_site] = randomize!(similar(H[1][1, 1, 1, 1], - left * physicalspace(H, middle_site) ← right)) - norm(ACs[middle_site]) == 0 && throw(ArgumentError("invalid sector")) - normalize!(ACs[middle_site]) - #construct the largest possible finite mps of that length + # center + ACs = Vector{Union{Missing,TA}}(missing, L) + P = physicalspace(H, middle_site) + ACs[middle_site] = rand!(similar(ALs[1], left ⊗ P ← right)) + + TB = tensormaptype(spacetype(H), 1, 1, T) + Cs = Vector{Union{Missing,TB}}(missing, L + 1) state = FiniteMPS(ALs, ARs, ACs, Cs) envs = environments(state, H) - #optimize the middle site. Because there is no truncation, this single site captures the entire possible hilbert space - H_ac = ∂∂AC(middle_site, state, H, envs) # this linear operator is now the actual full hamiltonian! - (vals, vecs, convhist) = eigsolve(H_ac, state.AC[middle_site], num, which, alg) + # optimize the middle site + # Because the MPS is full rank - this is equivalent to the full Hamiltonian + AC₀ = state.AC[middle_site] + H_ac = ∂∂AC(middle_site, state, H, envs) + vals, vecs, convhist = eigsolve(H_ac, AC₀, num, which, alg) + # repack eigenstates state_vecs = map(vecs) do v cs = copy(state) cs.AC[middle_site] = v diff --git a/src/operators/abstractmpo.jl b/src/operators/abstractmpo.jl index 21adff0be..fa3a1f9ef 100644 --- a/src/operators/abstractmpo.jl +++ b/src/operators/abstractmpo.jl @@ -261,7 +261,7 @@ end function add_physical_charge(O::MPOTensor, charge::Sector) sectortype(O) === typeof(charge) || throw(SectorMismatch()) - auxspace = Vect[typeof(charge)](charge => 1) + auxspace = Vect[typeof(charge)](charge => 1)' F = fuser(scalartype(O), physicalspace(O), auxspace) @plansor O_charged[-1 -2; -3 -4] := F[-2; 1 2] * O[-1 1; 4 3] * @@ -270,16 +270,14 @@ function add_physical_charge(O::MPOTensor, charge::Sector) end function add_physical_charge(O::BraidingTensor, charge::Sector) sectortype(O) === typeof(charge) || throw(SectorMismatch()) - auxspace = Vect[typeof(charge)](charge => 1) + auxspace = Vect[typeof(charge)](charge => 1)' V = left_virtualspace(O) ⊗ fuse(physicalspace(O), auxspace) ← fuse(physicalspace(O), auxspace) ⊗ right_virtualspace(O) return BraidingTensor{scalartype(O)}(V) end function add_physical_charge(O::AbstractBlockTensorMap{<:Any,<:Any,2,2}, charge::Sector) sectortype(O) == typeof(charge) || throw(SectorMismatch()) - - auxspace = Vect[typeof(charge)](charge => 1) - + auxspace = Vect[typeof(charge)](charge => 1)' Odst = similar(O, left_virtualspace(O) ⊗ fuse(physicalspace(O), auxspace) ← fuse(physicalspace(O), auxspace) ⊗ right_virtualspace(O)) diff --git a/src/states/finitemps.jl b/src/states/finitemps.jl index b1944516e..d6c309638 100644 --- a/src/states/finitemps.jl +++ b/src/states/finitemps.jl @@ -234,13 +234,13 @@ function FiniteMPS(f, elt, Pspaces::Vector{<:Union{S,CompositeSpace{S}}}, Vspaces[1] = left for k in 2:N - Vspaces[k] = infimum(fuse(Vspaces[k - 1], fusedPspaces[k]), maxVspaces[k - 1]) + Vspaces[k] = infimum(fuse(Vspaces[k - 1], fusedPspaces[k - 1]), maxVspaces[k - 1]) dim(Vspaces[k]) > 0 || @warn "no fusion channels available at site $k" end Vspaces[end] = right - for k in N:-1:2 - Vspaces[k] = infimum(maxVspaces[k - 1], fuse(Vspaces[k + 1], dual(fusedPspaces[k]))) + for k in reverse(2:N) + Vspaces[k] = infimum(Vspaces[k], fuse(Vspaces[k + 1], dual(fusedPspaces[k]))) dim(Vspaces[k]) > 0 || @warn "no fusion channels available at site $k" end @@ -383,18 +383,36 @@ function TensorKit.space(ψ::FiniteMPS{<:GenericMPSTensor}, n::Integer) return ProductSpace{S}(space.(Ref(t), Base.front(Base.tail(TensorKit.allind(t))))) end +""" + max_virtualspaces(ψ::FiniteMPS) + max_virtualspaces(Ps::Vector{<:Union{S,CompositeSpace{S}}}; left=oneunit(S), right=oneunit(S)) + +Compute the maximal virtual spaces of a given finite MPS or its physical spaces. +""" +function max_virtualspaces(Ps::Vector{<:Union{S,CompositeSpace{S}}}; left=oneunit(S), + right=oneunit(S)) where {S<:ElementarySpace} + Vs = similar(Ps, length(Ps) + 1) + Vs[1] = left + Vs[end] = right + for k in 2:length(Ps) + Vs[k] = fuse(Vs[k - 1], fuse(Ps[k - 1])) + end + for k in reverse(2:length(Ps)) + Vs[k] = infimum(Vs[k], fuse(Vs[k + 1], dual(fuse(Ps[k])))) + end + return Vs +end +function max_virtualspaces(ψ::FiniteMPS) + return max_virtualspaces(physicalspace(ψ); left=left_virtualspace(ψ, 1), + right=right_virtualspace(ψ, length(ψ))) +end + """ max_Ds(ψ::FiniteMPS) -> Vector{Float64} Compute the dimension of the maximal virtual space at a given site. """ -function max_Ds(ψ::FiniteMPS) - N = length(ψ) - physicaldims = dim.(space.(Ref(ψ)), 1:N) - D_left = cumprod(vcat(dim(left_virtualspace(ψ, 1)), physicaldims)) - D_right = cumprod(vcat(dim(right_virtualspace(ψ, N)), reverse(physicaldims))) - return min.(D_left, D_right) -end +max_Ds(ψ::FiniteMPS) = dim.(max_virtualspaces(ψ)) function Base.summary(io::IO, ψ::FiniteMPS) return print(io, "$(length(ψ))-site FiniteMPS ($(scalartype(ψ)), $(spacetype(ψ)))") diff --git a/src/states/quasiparticle_state.jl b/src/states/quasiparticle_state.jl index 58fabbb85..ed1b81d9e 100644 --- a/src/states/quasiparticle_state.jl +++ b/src/states/quasiparticle_state.jl @@ -41,10 +41,11 @@ function LeftGaugedQP(datfun, left_gs, right_gs=left_gs; Xs = map(enumerate(VLs)) do (loc, vl) x = similar(vl, right_virtualspace(vl) ← - excitation_space' ⊗ right_virtualspace(right_gs, loc)) + excitation_space ⊗ right_virtualspace(right_gs, loc)) fill_data!(x, datfun) return x end + sum(dim, Xs) == 0 && @warn "LeftGaugedQP: No possible fusion channels" left_gs isa InfiniteMPS || momentum == zero(momentum) || @warn "momentum is ignored for finite quasiparticles" @@ -63,13 +64,15 @@ end function RightGaugedQP(datfun, left_gs, right_gs=left_gs; sector=one(sectortype(left_gs)), momentum=0.0) - #find the left null spaces for the TNS + # find the left null spaces for the TNS excitation_space = Vect[typeof(sector)](sector => 1) - VRs = [adjoint(leftnull(adjoint(v))) for v in _transpose_tail.(right_gs.AR)] - Xs = [TensorMap{scalartype(left_gs)}(undef, left_virtualspace(right_gs, loc)', - excitation_space' * _firstspace(VRs[loc])) - for loc in 1:length(left_gs)] - fill_data!.(Xs, datfun) + VRs = convert(Vector, map(rightnull! ∘ _transpose_tail, right_gs.AR)) + Xs = map(enumerate(VRs)) do (i, vr) + x = similar(vr, + left_virtualspace(left_gs, i)' ← + excitation_space ⊗ _firstspace(vr)) + return fill_data!(x, datfun) + end left_gs isa InfiniteMPS || momentum == zero(momentum) || @warn "momentum is ignored for finite quasiparticles" @@ -84,10 +87,10 @@ function Base.similar(v::RightGaugedQP, ::Type{T}=scalartype(v)) where {T<:Numbe return RightGaugedQP(v.left_gs, v.right_gs, similar.(v.Xs, T), v.VRs, v.momentum) end -Base.getindex(v::LeftGaugedQP, i::Int) = v.VLs[mod1(i, end)] * v.Xs[mod1(i, end)]; +Base.getindex(v::LeftGaugedQP, i::Int) = v.VLs[mod1(i, end)] * v.Xs[mod1(i, end)] function Base.getindex(v::RightGaugedQP, i::Int) @plansor t[-1 -2; -3 -4] := v.Xs[mod1(i, end)][-1; -3 1] * v.VRs[mod1(i, end)][1; -4 -2] -end; +end function Base.setindex!(v::LeftGaugedQP, B, i::Int) v.Xs[mod1(i, end)] = v.VLs[mod1(i, end)]' * B diff --git a/test/algorithms.jl b/test/algorithms.jl index f8a8245c8..00a5d3b8a 100644 --- a/test/algorithms.jl +++ b/test/algorithms.jl @@ -11,6 +11,7 @@ using MPSKit using MPSKit: fuse_mul_mpo using TensorKit using TensorKit: ℙ +using LinearAlgebra: eigvals verbosity_full = 5 verbosity_conv = 1 @@ -917,4 +918,50 @@ end @test Z_tdvp ≈ Z_dense_2 atol = 1e-2 end +@testset "Sector conventions" begin + L = 4 + H = TestSetup.XY_model(U1Irrep; L) + + H_dense = convert(TensorMap, H) + vals_dense = eigvals(H_dense) + + tol = 1e-18 # tolerance required to separate degenerate eigenvalues + alg = MPSKit.Defaults.alg_eigsolve(; dynamic_tols=false, tol) + + maxVspaces = MPSKit.max_virtualspaces(physicalspace(H)) + gs, = find_groundstate(FiniteMPS(physicalspace(H), maxVspaces[2:(end - 1)]), H; + verbosity=0) + E₀ = expectation_value(gs, H) + @test E₀ ≈ first(vals_dense[one(U1Irrep)]) + + for (sector, vals) in vals_dense + # ED tests + num = length(vals) + E₀s, ψ₀s, info = exact_diagonalization(H; num, sector, alg) + @test E₀s[1:num] ≈ vals[1:num] + # this is a trick to make the mps full-rank again, which is not guaranteed by ED + ψ₀ = changebonds(first(ψ₀s), SvdCut(; trscheme=notrunc())) + Vspaces = left_virtualspace.(Ref(ψ₀), 1:L) + push!(Vspaces, right_virtualspace(ψ₀, L)) + @test all(splat(==), zip(Vspaces, MPSKit.max_virtualspaces(ψ₀))) + + # Quasiparticle tests + Es, Bs = excitations(H, QuasiparticleAnsatz(; tol), gs; sector, num=1) + Es = Es .+ E₀ + # first excited state is second eigenvalue if sector is trivial + @test Es[1] ≈ vals[isone(sector) ? 2 : 1] atol = 1e-8 + end + + # shifted charges tests + # targeting states with Sz = 1 => vals_shift_dense[0] == vals_dense[1] + # so effectively shifting the charges by -1 + H_shift = MPSKit.add_physical_charge(H, U1Irrep.([1, 0, 0, 0])) + H_shift_dense = convert(TensorMap, H_shift) + vals_shift_dense = eigvals(H_shift_dense) + for (sector, vals) in vals_dense + sector′ = only(sector ⊗ U1Irrep(-1)) + @test vals ≈ vals_shift_dense[sector′] + end +end + end diff --git a/test/excis.jl b/test/excis.jl new file mode 100644 index 000000000..3c0c87ab3 --- /dev/null +++ b/test/excis.jl @@ -0,0 +1,100 @@ +using MPSKit, TensorKit +import LinearAlgebra.eigvals +using MPSKit: max_virtualspaces + +V = U1Space(i => 1 for i in 0:1) +O = randn(V^2 ← V^2) +O += O' # Hermitian + +N = 3 +H = FiniteMPOHamiltonian(fill(V, N), (i, i + 1) => O for i in 1:(N - 1)); + +h = convert(TensorMap, H); + +sectors = collect(blocksectors(h)) +sec = U1Irrep(1) +@assert sec in sectors +num = 10 +vals1 = eigvals(block(h, sec)) +vals2 = eigvals(h)[sec] +@assert vals1 ≈ vals2 +vals3, vecs = exact_diagonalization(H; sector=(sec), num); +vals1 +vals3 + +right = U1Space(sec => 1) +max_virtualspaces(physicalspace(H); right) + +psi_full = rand(oneunit(V) * V^N ← right); +dim(space(psi_full)) +psi = MPSKit.decompose_localmps(psi_full); +left_virtualspace.(psi) + +left_virtualspace.(vecs[1].AL) +right_virtualspace.(vecs[1].AL) +max_virtualspaces(physicalspace(H); right) +psi1 = FiniteMPS(physicalspace(H), max_virtualspaces(physicalspace(H))[2:(end - 1)]) +psi1, = find_groundstate(psi1, H); +psi2 = FiniteMPS(physicalspace(H), max_virtualspaces(physicalspace(H); right)[2:(end - 1)]; + right) +left_virtualspace.(psi2.AL) +right_virtualspace.(psi2.AL) + +psi2, = find_groundstate(psi2, H); +expectation_value(psi2, H) + +Es, Bs = excitations(H, QuasiparticleAnsatz(), FiniteMPS(psi); sector=sec); +@inferred excitations(H, QuasiparticleAnsatz(), FiniteMPS(psi); sector=sec); + +Es .+ expectation_value(psi1, H) + +vals1 +vals3 + +psi2 +using TestEnv +using Test +TestEnv.activate() +include("setup.jl") +using .TestSetup +using TensorKit, MPSKit +using MPSKit: Multiline +using KrylovKit + +H = repeat(TestSetup.sixvertex(), 2) +ψ = InfiniteMPS([ℂ^2, ℂ^2], [ℂ^10, ℂ^10]) +ψ, envs, _ = leading_boundary(ψ, H, + VUMPS(; maxiter=400, tol=1e-10)) +energies, ϕs = @inferred excitations(H, QuasiparticleAnsatz(), + [0.0, Float64(pi / 2)], ψ, + envs; verbosity=0) +@test abs(energies[1]) > abs(energies[2]) # has a minimum at pi/2 +alg = QuasiparticleAnsatz() +ps = [0.0, Float64(pi / 2)] +excitations(H, alg, ps, ψ, envs; verbosity=0); +using Cthulhu +Hm = convert(MultilineMPO, H); +psim = convert(MultilineMPS, ψ); +envs = environments(psim, Hm); +excitations(Hm, alg, ps[1], psim, envs, psim, envs; verbosity=0); + +@descend excitations(Hm, alg, ps[1], psim, envs, psim, envs); + +qp = Multiline([LeftGaugedQP(rand, psim[1], psim[1]; sector=one(sectortype(psim[1])), + momentum=ps[1])]); +excitations(Hm, alg, qp, envs, envs; num=1); + +@descend excitations(Hm, alg, qp, envs, envs; num=1); +@code_warntype excitations(Hm, alg, qp, envs, envs; num=1); + +Heff = MPSKit.EffectiveExcitationHamiltonian(H, envs[1], envs[1], fill(1.0, length(H))); +Heffs = Multiline([Heff]); +@code_warntype KrylovKit.apply(Heff, qp[1]) +@descend KrylovKit.apply(Heff, qp[1]) + +KrylovKit.apply(Multiline([Heff]), qp); +@code_warntype KrylovKit.apply(Heffs, qp); + +@descend eigsolve(Heffs, qp, 1, :LM, Lanczos()); +@code_warntype eigsolve(Heffs, qp, 1, :LM); +@descend eigsolve(Heffs, qp, 1, :LM); diff --git a/test/setup.jl b/test/setup.jl index 804ca6e3a..dd2e928b4 100644 --- a/test/setup.jl +++ b/test/setup.jl @@ -14,7 +14,7 @@ using Combinatorics: permutations export S_xx, S_yy, S_zz, S_x, S_y, S_z export force_planar export symm_mul_mpo -export transverse_field_ising, heisenberg_XXX, bilinear_biquadratic_model +export transverse_field_ising, heisenberg_XXX, bilinear_biquadratic_model, XY_model export classical_ising, finite_classical_ising, sixvertex # using TensorOperations @@ -173,6 +173,27 @@ function heisenberg_XXX(; spin=1, L=Inf) end end +function XY_model(::Type{U1Irrep}; g=1 / 2, L=Inf) + h = zeros(ComplexF64, + U1Space(-1 // 2 => 1, 1 // 2 => 1)^2 ← U1Space(-1 // 2 => 1, 1 // 2 => 1)^2) + h[U1Irrep.((-1 // 2, 1 // 2, -1 // 2, 1 // 2))] .= 1 + h[U1Irrep.((1 // 2, -1 // 2, 1 // 2, -1 // 2))] .= 1 + Sz = zeros(ComplexF64, space(h, 1) ← space(h, 1)) + Sz[U1Irrep.((-1 // 2, 1 // 2))] .= -1 // 2 + Sz[U1Irrep.((1 // 2, -1 // 2))] .= 1 // 2 + if L == Inf + lattice = PeriodicArray([space(h, 1)]) + H = InfiniteMPOHamiltonian(lattice, (i, i + 1) => -h for i in 1:1) + iszero(g) && return H + return H + g * InfiniteMPOHamiltonian(lattice, (i,) => -Sz for i in 1:1) + else + lattice = fill(space(h, 1), L) + H = FiniteMPOHamiltonian(lattice, (i, i + 1) => -h for i in 1:(L - 1)) + iszero(g) && return H + return H + g * FiniteMPOHamiltonian(lattice, (i,) => -Sz for i in 1:L) + end +end + function bilinear_biquadratic_model(::Type{SU2Irrep}; θ=atan(1 / 3), L=Inf) h1 = ones(ComplexF64, SU2Space(1 => 1)^2 ← SU2Space(1 => 1)^2) for (c, b) in blocks(h1)