diff --git a/Project.toml b/Project.toml index d39631a..7d2104d 100644 --- a/Project.toml +++ b/Project.toml @@ -14,17 +14,19 @@ LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e" OffsetArrays = "6fe1bfb0-de20-5000-8ca7-80f57d26f881" QuadGK = "1fd47b50-473d-5c70-9696-f719f8f3bcdc" Requires = "ae029012-a4dd-5104-9daa-d747884805df" +SplitApplyCombine = "03a91e81-4c3e-53e1-a0a4-9c0c8f19dd66" [compat] Compat = "3, 4" HDF5 = "0.15, 0.16" -ITensors = "0.3" +ITensors = "0.3.18" Infinities = "0.1" IterTools = "1" KrylovKit = "0.5" OffsetArrays = "1" QuadGK = "2" Requires = "1" +SplitApplyCombine = "1.2.2" julia = "1.6" [extras] diff --git a/examples/development/finite_mps_to_infinite_mps.jl b/examples/development/finite_mps_to_infinite_mps.jl new file mode 100644 index 0000000..01fe041 --- /dev/null +++ b/examples/development/finite_mps_to_infinite_mps.jl @@ -0,0 +1,63 @@ +using ITensors +using ITensorInfiniteMPS + +# H = -J Σⱼ XⱼXⱼ₊₁ - h Σⱼ Zⱼ +function ising_opsum_infinite(N; J, h) + os = OpSum() + for n in 1:N + os -= J, "X", n, "X", n + 1 + end + for n in 1:N + os -= h, "Z", n + end + return os +end + +function ising_opsum_finite(N; J, h) + return ising_opsum_infinite(N - 1; J, h) + (-h, "Z", N) +end + +function main(; N, J, h, nsites) + s = siteinds("S=1/2", N) + + H = MPO(ising_opsum_finite(N; J=J, h=h), s) + ψ0 = randomMPS(s) + + energy, ψ = dmrg(H, ψ0; nsweeps=10, cutoff=1e-10) + @show energy / N + + n1 = N ÷ 2 + n2 = n1 + 1 + + ψ = orthogonalize(ψ, n1) + + println("\nFinite MPS energy on bond (n1, n2) = $((n1, n2))") + + ϕ12 = ψ[n1] * ψ[n2] + ham12 = contract(MPO(ising_opsum_infinite(1; J=J, h=h), [s[n1], s[n2]])) + @show inner(ϕ12, apply(ham12, ϕ12)) + + # Get an infinite MPS that approximates + # the finite MPS in the bulk. + # `nsites` sets the number of sites in the unit + # cell of the infinite MPS. + ψ∞ = infinitemps_approx(ψ; nsites=nsites, nsweeps=10) + + println("\nInfinite MPS energy on a few different bonds") + println("Infinite MPS has unit cell of size $nsites") + + # Compute a few energies + for n1 in 1:3 + n2 = n1 + 1 + @show n1, n2 + ϕ12 = ψ∞.AL[n1] * ψ∞.AL[n2] * ψ∞.C[n2] + s1 = siteind(ψ∞, n1) + s2 = siteind(ψ∞, n2) + ham12 = contract(MPO(ising_opsum_infinite(1; J=J, h=h), [s1, s2])) + @show inner(ϕ12, apply(ham12, ϕ12)) / norm(ϕ12) + end + + return nothing +end + +main(; N=100, J=1.0, h=1.5, nsites=1) diff --git a/examples/development/orthogonalize_infinitemps.jl b/examples/development/orthogonalize_infinitemps.jl new file mode 100644 index 0000000..273b1cf --- /dev/null +++ b/examples/development/orthogonalize_infinitemps.jl @@ -0,0 +1,19 @@ +using ITensors +using ITensorInfiniteMPS +using Random + +Random.seed!(1234) + +# n-site unit cell +nsites = 4 +s = siteinds("S=1/2", nsites; conserve_szparity=true) +χ = 10 +@assert iseven(χ) +space = (("SzParity", 1, 2) => χ ÷ 2) ⊕ (("SzParity", 0, 2) => χ ÷ 2) +ψ = InfiniteMPS(ComplexF64, s; space=space) +for n in 1:nsites + ψ[n] = randomITensor(inds(ψ[n])) +end + +ψ = orthogonalize(ψ, :) +@show norm(contract(ψ.AL[1:nsites]) * ψ.C[nsites] - ψ.C[0] * contract(ψ.AR[1:nsites])) diff --git a/examples/development/vumps/transfer_matrix_spectrum_arpack.jl b/examples/development/vumps/transfer_matrix_spectrum_arpack.jl new file mode 100644 index 0000000..b7c65da --- /dev/null +++ b/examples/development/vumps/transfer_matrix_spectrum_arpack.jl @@ -0,0 +1,193 @@ +using ITensors +using ITensorInfiniteMPS + +############################################################################## +# VUMPS parameters +# + +maxdim = 10 # Maximum bond dimension +cutoff = 1e-6 # Singular value cutoff when increasing the bond dimension +max_vumps_iters = 20 # Maximum number of iterations of the VUMPS algorithm at a fixed bond dimension +tol = 1e-5 # Precision error tolerance for outer loop of VUMPS or TDVP +outer_iters = 5 # Number of times to increase the bond dimension +time_step = -Inf # -Inf corresponds to VUMPS +solver_tol = (x -> x / 100) # Tolerance for the local solver (eigsolve in VUMPS and exponentiate in TDVP) +multisite_update_alg = "parallel" # Choose between ["sequential", "parallel"] +conserve_qns = false +N = 2 # Number of sites in the unit cell (1 site unit cell is currently broken) + +# Parameters of the transverse field Ising model +model_params = (J=1.0, h=0.9) + +############################################################################## +# CODE BELOW HERE DOES NOT NEED TO BE MODIFIED +# + +model = Model("ising") + +initstate(n) = "↑" +s = infsiteinds("S=1/2", N; initstate, conserve_szparity=conserve_qns) +ψ = InfMPS(s, initstate) + +# Form the Hamiltonian +H = InfiniteSum{MPO}(model, s; model_params...) + +# Check translational invariance +@show norm(contract(ψ.AL[1:N]..., ψ.C[N]) - contract(ψ.C[0], ψ.AR[1:N]...)) + +vumps_kwargs = ( + tol=tol, + maxiter=max_vumps_iters, + solver_tol=solver_tol, + multisite_update_alg=multisite_update_alg, +) +subspace_expansion_kwargs = (cutoff=cutoff, maxdim=maxdim) +#ψ = tdvp(H, ψ; time_step=time_step, vumps_kwargs...) + +# Alternate steps of running VUMPS and increasing the bond dimension +@time for _ in 1:outer_iters + println("\nIncrease bond dimension") + global ψ = subspace_expansion(ψ, H; subspace_expansion_kwargs...) + println("Run VUMPS with new bond dimension") + global ψ = tdvp(H, ψ; time_step=time_step, vumps_kwargs...) +end + +# Check translational invariance +@show norm(contract(ψ.AL[1:N]..., ψ.C[N]) - contract(ψ.C[0], ψ.AR[1:N]...)) + +# +# Compare to DMRG +# + +Nfinite = 100 +sfinite = siteinds("S=1/2", Nfinite; conserve_szparity=true) +Hfinite = MPO(model, sfinite; model_params...) +ψfinite = randomMPS(sfinite, initstate) +@show flux(ψfinite) +sweeps = Sweeps(10) +setmaxdim!(sweeps, maxdim) +setcutoff!(sweeps, cutoff) +energy_finite_total, ψfinite = @time dmrg(Hfinite, ψfinite, sweeps) +@show energy_finite_total / Nfinite + +function energy_local(ψ1, ψ2, h) + ϕ = ψ1 * ψ2 + return (noprime(ϕ * h) * dag(ϕ))[] +end + +function ITensors.expect(ψ, o) + return (noprime(ψ * op(o, filterinds(ψ, "Site")...)) * dag(ψ))[] +end + +# Exact energy at criticality: 4/pi = 1.2732395447351628 + +nfinite = Nfinite ÷ 2 +orthogonalize!(ψfinite, nfinite) +hnfinite = ITensor(model, sfinite[nfinite], sfinite[nfinite + 1]; model_params...) +energy_finite = energy_local(ψfinite[nfinite], ψfinite[nfinite + 1], hnfinite) +energy_infinite = energy_local(ψ.AL[1], ψ.AL[2] * ψ.C[2], contract(H[(1, 2)])) +@show energy_finite, energy_infinite +@show abs(energy_finite - energy_infinite) + +energy_exact = reference(model, Observable("energy"); model_params...) +@show energy_exact + +Sz1_finite = expect(ψfinite[nfinite], "Sz") +orthogonalize!(ψfinite, nfinite + 1) +Sz2_finite = expect(ψfinite[nfinite + 1], "Sz") +Sz1_infinite = expect(ψ.AL[1] * ψ.C[1], "Sz") +Sz2_infinite = expect(ψ.AL[2] * ψ.C[2], "Sz") + +@show Sz1_finite, Sz2_finite +@show Sz1_infinite, Sz2_infinite + +############################################################################## +# Compute eigenspace of the transfer matrix +# + +using Arpack +using KrylovKit: eigsolve +using LinearAlgebra + +T = TransferMatrix(ψ.AL) +Tᵀ = transpose(T) +vⁱᴿ = randomITensor(dag(input_inds(T))) +vⁱᴸ = randomITensor(dag(input_inds(Tᵀ))) + +neigs = 10 +tol = 1e-10 +λ⃗ᴿ, v⃗ᴿ, right_info = eigsolve(T, vⁱᴿ, neigs, :LM; tol=tol) +λ⃗ᴸ, v⃗ᴸ, left_info = eigsolve(Tᵀ, vⁱᴸ, neigs, :LM; tol=tol) + +@show norm(T(v⃗ᴿ[1]) - λ⃗ᴿ[1] * v⃗ᴿ[1]) +@show norm(Tᵀ(v⃗ᴸ[1]) - λ⃗ᴸ[1] * v⃗ᴸ[1]) + +@show λ⃗ᴿ +@show λ⃗ᴸ +@show flux.(v⃗ᴿ) + +neigs = min(length(v⃗ᴸ), length(v⃗ᴿ)) +v⃗ᴸ = v⃗ᴸ[1:neigs] +v⃗ᴿ = v⃗ᴿ[1:neigs] + +# Normalize the vectors +N⃗ = [(translatecelltags(v⃗ᴸ[n], 1) * v⃗ᴿ[n])[] for n in 1:neigs] + +v⃗ᴿ ./= sqrt.(N⃗) +v⃗ᴸ ./= sqrt.(N⃗) + +# Form a second starting vector orthogonal to v⃗ᴿ[1] +# This doesn't work. TODO: project out v⃗ᴿ[1], v⃗ᴸ[1] from T +#λ⃗ᴿ², v⃗ᴿ², right_info_2 = eigsolve(T, vⁱᴿ², neigs, :LM; tol=tol) + +# Projector onto the n-th eigenstate +function proj(v⃗ᴸ, v⃗ᴿ, n) + Lⁿ = v⃗ᴸ[n] + Rⁿ = v⃗ᴿ[n] + return ITensorMap( + [translatecelltags(Lⁿ, 1), translatecelltags(Rⁿ, -1)]; + input_inds=inds(Rⁿ), + output_inds=inds(Lⁿ), + ) +end + +P⃗ = [proj(v⃗ᴸ, v⃗ᴿ, n) for n in 1:neigs] +T⁻P = T - sum(P⃗) + +#vⁱᴿ² = vⁱᴿ - (translatecelltags(v⃗ᴸ[1], 1) * vⁱᴿ)[] / norm(v⃗ᴿ[1]) * v⃗ᴿ[1] +#@show norm(dag(vⁱᴿ²) * v⃗ᴿ[1]) + +# XXX: Error with mismatched QN index arrows +# λ⃗ᴾᴿ, v⃗ᴾᴿ, right_info = eigsolve(T⁻P, vⁱᴿ, neigs, :LM; tol=tol) +# @show λ⃗ᴾᴿ + +vⁱᴿ⁻ᵈᵃᵗᵃ = vec(array(vⁱᴿ)) +λ⃗ᴿᴬ, v⃗ᴿ⁻ᵈᵃᵗᵃ = Arpack.eigs(T; v0=vⁱᴿ⁻ᵈᵃᵗᵃ, nev=neigs) + +## XXX: this is giving an error about trying to set the element of the wrong QN block for: +## maxdim = 5 +## cutoff = 1e-12 +## max_vumps_iters = 10 +## outer_iters = 10 +## model_params = (J=1.0, h=0.8) +## +## v⃗ᴿᴬ = [itensor(v⃗ᴿ⁻ᵈᵃᵗᵃ[:, n], input_inds(T); tol=1e-4) for n in 1:length(λ⃗ᴿᴬ)] +## @show flux.(v⃗ᴿᴬ) + +@show λ⃗ᴿᴬ + +# Full eigendecomposition + +Tfull = prod(T) +DV = eigen(Tfull, input_inds(T), output_inds(T)) + +@show norm(Tfull * DV.V - DV.Vt * DV.D) + +d = diag(array(DV.D)) + +p = sortperm(d; by=abs, rev=true) +@show p[1:neigs] +@show d[p[1:neigs]] + +println("Error if ED with Arpack") +@show d[p[1:neigs]] - λ⃗ᴿᴬ diff --git a/examples/vumps/vumps_from_dmrg.jl b/examples/development/vumps/vumps_from_dmrg.jl similarity index 100% rename from examples/vumps/vumps_from_dmrg.jl rename to examples/development/vumps/vumps_from_dmrg.jl diff --git a/examples/vumps/vumps_ising_noncontiguous.jl b/examples/development/vumps/vumps_ising_noncontiguous.jl similarity index 88% rename from examples/vumps/vumps_ising_noncontiguous.jl rename to examples/development/vumps/vumps_ising_noncontiguous.jl index 3fc08bf..032a9a2 100644 --- a/examples/vumps/vumps_ising_noncontiguous.jl +++ b/examples/development/vumps/vumps_ising_noncontiguous.jl @@ -18,13 +18,8 @@ N = 2# Number of sites in the unit cell J = 1.0 h = 1.0; -function space_shifted(q̃sz) - return [QN("SzParity", 1 - q̃sz, 2) => 1, QN("SzParity", 0 - q̃sz, 2) => 1] -end - -space_ = fill(space_shifted(1), N); -s = infsiteinds("S=1/2", N; space=space_) initstate(n) = "↑" +s = infsiteinds("S=1/2", N; initstate) ψ = InfMPS(s, initstate); model = Model("ising"); @@ -70,7 +65,7 @@ Sz_finite = expect(ψfinite, "Sz")[nfinite:(nfinite + N - 1)] @show ( energy_infinite, energy_finite_total / Nfinite, - reference(model, Observable("energy"); J=J, h=h, J₂=J₂), + reference(model, Observable("energy"); J=J, h=h), ) @show (Sz, Sz_finite) diff --git a/examples/vumps/vumps_localham.jl b/examples/development/vumps/vumps_localham.jl similarity index 97% rename from examples/vumps/vumps_localham.jl rename to examples/development/vumps/vumps_localham.jl index 55082c8..14ff15c 100644 --- a/examples/vumps/vumps_localham.jl +++ b/examples/development/vumps/vumps_localham.jl @@ -4,8 +4,6 @@ using Random Random.seed!(1234) -include("models.jl") - function _siteind(site_tag, n::Int; space) return addtags(Index(space, "Site,n=$n"), site_tag) end diff --git a/examples/vumps/vumps_localham_qns.jl b/examples/development/vumps/vumps_localham_qns.jl similarity index 100% rename from examples/vumps/vumps_localham_qns.jl rename to examples/development/vumps/vumps_localham_qns.jl diff --git a/examples/vumps/development/vumps_mpo.jl b/examples/development/vumps/vumps_mpo.jl similarity index 59% rename from examples/vumps/development/vumps_mpo.jl rename to examples/development/vumps/vumps_mpo.jl index 652a05c..acfa204 100644 --- a/examples/vumps/development/vumps_mpo.jl +++ b/examples/development/vumps/vumps_mpo.jl @@ -10,35 +10,30 @@ maxiter = 20 outer_iters = 3 N = 2 -model = Model"ising"() +model = Model("ising") -function space_shifted(::Model"ising", q̃sz) - return [QN("SzParity", 1 - q̃sz, 2) => 1, QN("SzParity", 0 - q̃sz, 2) => 1] -end - -space_ = fill(space_shifted(model, 1), N) -s = infsiteinds("S=1/2", N; space=space_) initstate(n) = "↑" +s = infsiteinds("S=1/2", N; initstate) ψ = InfMPS(s, initstate) model_params = (J=-1.0, h=0.9) vumps_kwargs = (multisite_update_alg="sequential", tol=tol, maxiter=maxiter) subspace_expansion_kwargs = (cutoff=cutoff, maxdim=maxdim) -H = InfiniteITensorSum(model, s; model_params...) +H = InfiniteSum{ITensor}(model, s; model_params...) # Alternate steps of running VUMPS and increasing the bond dimension ψ1 = vumps(H, ψ; vumps_kwargs...) for _ in 1:outer_iters - ψ1 = subspace_expansion(ψ1, H; subspace_expansion_kwargs...) - ψ1 = vumps(H, ψ1; vumps_kwargs...) + global ψ1 = subspace_expansion(ψ1, H; subspace_expansion_kwargs...) + global ψ1 = vumps(H, ψ1; vumps_kwargs...) end Hmpo = InfiniteMPOMatrix(model, s; model_params...) # Alternate steps of running VUMPS and increasing the bond dimension ψ2 = vumps(Hmpo, ψ; vumps_kwargs...) for _ in 1:outer_iters - ψ2 = subspace_expansion(ψ2, Hmpo; subspace_expansion_kwargs...) - ψ2 = vumps(Hmpo, ψ2; vumps_kwargs...) + global ψ2 = subspace_expansion(ψ2, Hmpo; subspace_expansion_kwargs...) + global ψ2 = vumps(Hmpo, ψ2; vumps_kwargs...) end SzSz = prod(op("Sz", s[1]) * op("Sz", s[2])) diff --git a/examples/finite_mps_to_infinite_mps.jl b/examples/finite_mps_to_infinite_mps.jl deleted file mode 100644 index e97e2c1..0000000 --- a/examples/finite_mps_to_infinite_mps.jl +++ /dev/null @@ -1,112 +0,0 @@ -using ITensors -using ITensorInfiniteMPS - -# H = -J Σⱼ XⱼXⱼ₊₁ - h Σⱼ Zⱼ -function ising_mpo(s; J, h) - N = length(s) - a = AutoMPO() - for n in 1:(N - 1) - a .+= -J, "X", n, "X", n + 1 - end - for n in 1:N - a .+= -h, "Z", n - end - return H = MPO(a, s) -end - -function ising_infinitempo_tensor(s⃗, l, r; J, h) - s = only(s⃗) - dₗ = 3 # The link dimension of the TFI - Hmat = fill(ITensor(s', dag(s)), dₗ, dₗ) - Hmat[1, 1] = op("Id", s) - Hmat[2, 2] = op("Id", s) - Hmat[3, 3] = op("Id", s) - Hmat[2, 1] = -J * op("X", s) - Hmat[3, 2] = -J * op("X", s) - Hmat[3, 1] = -h * op("Z", s) - H = emptyITensor() - for i in 1:dₗ, j in 1:dₗ - H += Hmat[i, j] * setelt(l => i) * setelt(r => j) - end - return H -end - -# H = -J Σⱼ XⱼXⱼ₊₁ - h Σⱼ Zⱼ -function ising_infinitempo(s; J, h) - N = length(s) - # Use a finite MPO to create the infinite MPO - H = InfiniteMPO(N) - # For the finite unit cell - s¹ = addtags.(s, "c=1") - l¹ = [Index(3, "Link,l=$n,c=1") for n in 1:N] - for n in 1:N - # TODO: use a CelledVector here to make the code logic simpler - l¹ₙ = n == 1 ? replacetags(dag(l¹[N]), "c=1" => "c=0") : l¹[n - 1] - r¹ₙ = l¹[n] - H[n] = ising_infinitempo_tensor(s¹[n], l¹ₙ, r¹ₙ; J=J, h=h) - end - return H -end - -# H = -J X₁X₂ - h Z₁ -function ising_localop(s1::Index, s2::Index; J, h) - return -J * op("X", s1) * op("X", s2) - h * op("Z", s1) * op("Id", s2) -end -# Version accepting IndexSet -ising_localop(s1, s2; kwargs...) = ising_localop(Index(s1), Index(s2); kwargs...) - -function ising_infinitesumlocalops(s; J, h) - N = length(s) - H = InfiniteSumLocalOps(N) - s∞ = CelledVector(s) - return InfiniteSumLocalOps([ising_localop(s∞[n], s∞[n + 1]; J=J, h=h) for n in 1:N]) -end - -let - N = 100 - s = siteinds("S=1/2", N) - - J = 1.0 - h = 1.5 - H = ising_mpo(s; J=J, h=h) - ψ0 = randomMPS(s) - - sweeps = Sweeps(10) - maxdim!(sweeps, 10) - cutoff!(sweeps, 1E-10) - energy, ψ = dmrg(H, ψ0, sweeps) - @show energy / N - - n1 = N ÷ 2 - n2 = n1 + 1 - - ψ = orthogonalize(ψ, n1) - - println("\nFinite MPS energy on bond (n1, n2) = $((n1, n2))") - - ϕ12 = ψ[n1] * ψ[n2] - ham12 = ising_localop(s[n1], s[n2]; J=J, h=h) - @show (ϕ12 * ham12 * prime(dag(ϕ12), commoninds(ϕ12, ham12)))[] - - # Get an infinite MPS that approximates - # the finite MPS in the bulk. - # `nsites` sets the number of sites in the unit - # cell of the infinite MPS. - nsites = 1 - ψ∞ = infinitemps_approx(ψ; nsites=nsites, nsweeps=10) - - println("\nInfinite MPS energy on a few different bonds") - println("Infinite MPS has unit cell of size $nsites") - - # Compute a few energies - for n1 in 1:3 - n2 = n1 + 1 - @show n1, n2 - ϕ12 = ψ∞.AL[n1] * ψ∞.AL[n2] * ψ∞.C[n2] - s1 = siteind(ψ∞, n1) - s2 = siteind(ψ∞, n2) - ham12 = ising_localop(s1, s2; J=J, h=h) - ϕ12ᴴ = prime(dag(ϕ12), commoninds(ϕ12, ham12)) - @show (ϕ12 * ham12 * ϕ12ᴴ)[] / norm(ϕ12) - end -end diff --git a/examples/infinitemps.jl b/examples/infinitemps.jl deleted file mode 100644 index 8c4afef..0000000 --- a/examples/infinitemps.jl +++ /dev/null @@ -1,20 +0,0 @@ -using ITensors -using ITensorInfiniteMPS -using Random - -Random.seed!(1234) - -# -# Example -# - -N = 10 -s = siteinds("S=1/2", N; conserve_szparity=true) -χ = 6 -@assert iseven(χ) -space = (("SzParity", 1, 2) => χ ÷ 2) ⊕ (("SzParity", 0, 2) => χ ÷ 2) -ψ = InfiniteMPS(ComplexF64, s; space=space) -randn!.(ψ) - -ψ = orthogonalize(ψ, :) -@show norm(prod(ψ.AL[1:N]) * ψ.C[N] - ψ.C[0] * prod(ψ.AR[1:N])) diff --git a/examples/vumps/src/vumps_subspace_expansion.jl b/examples/vumps/src/vumps_subspace_expansion.jl new file mode 100644 index 0000000..bca3068 --- /dev/null +++ b/examples/vumps/src/vumps_subspace_expansion.jl @@ -0,0 +1,28 @@ +using ITensors +using ITensorInfiniteMPS + +# Alternate steps of running TDVP and increasing the bond dimension. +# Attempt increasing the bond dimension `outer_iters` number of times. +function tdvp_subspace_expansion( + H, ψ; time_step, outer_iters, subspace_expansion_kwargs, vumps_kwargs +) + @time for _ in 1:outer_iters + println("\nIncrease bond dimension from $(maxlinkdim(ψ))") + println( + "cutoff = $(subspace_expansion_kwargs[:cutoff]), maxdim = $(subspace_expansion_kwargs[:maxdim])", + ) + ψ = subspace_expansion(ψ, H; subspace_expansion_kwargs...) + println("\nRun VUMPS with new bond dimension $(maxlinkdim(ψ))") + ψ = tdvp(H, ψ; time_step=time_step, vumps_kwargs...) + end + return ψ +end + +# Alternate steps of running VUMPS and increasing the bond dimension +function vumps_subspace_expansion( + H, ψ; outer_iters, subspace_expansion_kwargs, vumps_kwargs +) + return tdvp_subspace_expansion( + H, ψ; time_step=-Inf, outer_iters, subspace_expansion_kwargs, vumps_kwargs + ) +end diff --git a/examples/vumps/transfer_matrix_spectrum.jl b/examples/vumps/transfer_matrix_spectrum.jl index 9ff369a..ed6a4e4 100644 --- a/examples/vumps/transfer_matrix_spectrum.jl +++ b/examples/vumps/transfer_matrix_spectrum.jl @@ -1,6 +1,12 @@ using ITensors using ITensorInfiniteMPS +include( + joinpath( + pkgdir(ITensorInfiniteMPS), "examples", "vumps", "src", "vumps_subspace_expansion.jl" + ), +) + ############################################################################## # VUMPS parameters # @@ -23,19 +29,10 @@ model_params = (J=1.0, h=0.9) # CODE BELOW HERE DOES NOT NEED TO BE MODIFIED # -model = Model"ising"() +model = Model("ising") -function space_shifted(::Model"ising", q̃sz; conserve_qns=true) - if conserve_qns - return [QN("SzParity", 1 - q̃sz, 2) => 1, QN("SzParity", 0 - q̃sz, 2) => 1] - else - return [QN() => 2] - end -end - -space_ = fill(space_shifted(model, 0; conserve_qns=conserve_qns), N) -s = infsiteinds("S=1/2", N; space=space_) initstate(n) = "↑" +s = infsiteinds("S=1/2", N; initstate, conserve_szparity=conserve_qns) ψ = InfMPS(s, initstate) # Form the Hamiltonian @@ -51,15 +48,10 @@ vumps_kwargs = ( multisite_update_alg=multisite_update_alg, ) subspace_expansion_kwargs = (cutoff=cutoff, maxdim=maxdim) -#ψ = tdvp(H, ψ; time_step=time_step, vumps_kwargs...) - -# Alternate steps of running VUMPS and increasing the bond dimension -@time for _ in 1:outer_iters - println("\nIncrease bond dimension") - global ψ = subspace_expansion(ψ, H; subspace_expansion_kwargs...) - println("Run VUMPS with new bond dimension") - global ψ = tdvp(H, ψ; time_step=time_step, vumps_kwargs...) -end + +ψ = tdvp_subspace_expansion( + H, ψ; time_step, outer_iters, subspace_expansion_kwargs, vumps_kwargs +) # Check translational invariance @show norm(contract(ψ.AL[1:N]..., ψ.C[N]) - contract(ψ.C[0], ψ.AR[1:N]...)) @@ -69,7 +61,7 @@ end # Nfinite = 100 -sfinite = siteinds("S=1/2", Nfinite; conserve_szparity=true) +sfinite = siteinds("S=1/2", Nfinite; conserve_szparity=conserve_qns) Hfinite = MPO(model, sfinite; model_params...) ψfinite = randomMPS(sfinite, initstate) @show flux(ψfinite) @@ -81,11 +73,11 @@ energy_finite_total, ψfinite = @time dmrg(Hfinite, ψfinite, sweeps) function energy_local(ψ1, ψ2, h) ϕ = ψ1 * ψ2 - return (noprime(ϕ * h) * dag(ϕ))[] + return inner(ϕ, apply(h, ϕ)) end -function ITensors.expect(ψ, o) - return (noprime(ψ * op(o, filterinds(ψ, "Site")...)) * dag(ψ))[] +function ITensors.expect(ψ::ITensor, o::String) + return inner(ψ, apply(op(o, filterinds(ψ, "Site")...), ψ)) end # Exact energy at criticality: 4/pi = 1.2732395447351628 @@ -94,7 +86,7 @@ nfinite = Nfinite ÷ 2 orthogonalize!(ψfinite, nfinite) hnfinite = ITensor(model, sfinite[nfinite], sfinite[nfinite + 1]; model_params...) energy_finite = energy_local(ψfinite[nfinite], ψfinite[nfinite + 1], hnfinite) -energy_infinite = energy_local(ψ.AL[1], ψ.AL[2] * ψ.C[2], H[(1, 2)]) +energy_infinite = energy_local(ψ.AL[1], ψ.AL[2] * ψ.C[2], contract(H[(1, 2)])) @show energy_finite, energy_infinite @show abs(energy_finite - energy_infinite) @@ -114,8 +106,7 @@ Sz2_infinite = expect(ψ.AL[2] * ψ.C[2], "Sz") # Compute eigenspace of the transfer matrix # -using Arpack -using KrylovKit +using KrylovKit: eigsolve using LinearAlgebra T = TransferMatrix(ψ.AL) @@ -128,6 +119,8 @@ tol = 1e-10 λ⃗ᴿ, v⃗ᴿ, right_info = eigsolve(T, vⁱᴿ, neigs, :LM; tol=tol) λ⃗ᴸ, v⃗ᴸ, left_info = eigsolve(Tᵀ, vⁱᴸ, neigs, :LM; tol=tol) +println("\n##########################################") +println("Check transfer matrix left and right fixed points") @show norm(T(v⃗ᴿ[1]) - λ⃗ᴿ[1] * v⃗ᴿ[1]) @show norm(Tᵀ(v⃗ᴸ[1]) - λ⃗ᴸ[1] * v⃗ᴸ[1]) @@ -135,56 +128,25 @@ tol = 1e-10 @show λ⃗ᴸ @show flux.(v⃗ᴿ) -neigs = length(v⃗ᴿ) +neigs = min(length(v⃗ᴸ), length(v⃗ᴿ)) +v⃗ᴸ = v⃗ᴸ[1:neigs] +v⃗ᴿ = v⃗ᴿ[1:neigs] # Normalize the vectors -N⃗ = [(translatecell(v⃗ᴸ[n], 1) * v⃗ᴿ[n])[] for n in 1:neigs] +N⃗ = [(translatecelltags(v⃗ᴸ[n], 1) * v⃗ᴿ[n])[] for n in 1:neigs] v⃗ᴿ ./= sqrt.(N⃗) v⃗ᴸ ./= sqrt.(N⃗) -# Form a second starting vector orthogonal to v⃗ᴿ[1] -# This doesn't work. TODO: project out v⃗ᴿ[1], v⃗ᴸ[1] from T -#λ⃗ᴿ², v⃗ᴿ², right_info_2 = eigsolve(T, vⁱᴿ², neigs, :LM; tol=tol) - -# Projector onto the n-th eigenstate -function proj(v⃗ᴸ, v⃗ᴿ, n) - Lⁿ = v⃗ᴸ[n] - Rⁿ = v⃗ᴿ[n] - return ITensorMap( - [translatecell(Lⁿ, 1), translatecell(Rⁿ, -1)]; input_inds=inds(Rⁿ), output_inds=inds(Lⁿ) - ) -end - -P⃗ = [proj(v⃗ᴸ, v⃗ᴿ, n) for n in 1:neigs] -T⁻P = T - sum(P⃗) - -#vⁱᴿ² = vⁱᴿ - (translatecell(v⃗ᴸ[1], 1) * vⁱᴿ)[] / norm(v⃗ᴿ[1]) * v⃗ᴿ[1] -#@show norm(dag(vⁱᴿ²) * v⃗ᴿ[1]) - -λ⃗ᴾᴿ, v⃗ᴾᴿ, right_info = eigsolve(T⁻P, vⁱᴿ, neigs, :LM; tol=tol) -@show λ⃗ᴾᴿ - -vⁱᴿ⁻ᵈᵃᵗᵃ = vec(array(vⁱᴿ)) -λ⃗ᴿᴬ, v⃗ᴿ⁻ᵈᵃᵗᵃ = Arpack.eigs(T; v0=vⁱᴿ⁻ᵈᵃᵗᵃ, nev=neigs) - -## XXX: this is giving an error about trying to set the element of the wrong QN block for: -## maxdim = 5 -## cutoff = 1e-12 -## max_vumps_iters = 10 -## outer_iters = 10 -## model_params = (J=1.0, h=0.8) -## -## v⃗ᴿᴬ = [itensor(v⃗ᴿ⁻ᵈᵃᵗᵃ[:, n], input_inds(T); tol=1e-4) for n in 1:length(λ⃗ᴿᴬ)] -## @show flux.(v⃗ᴿᴬ) - -@show λ⃗ᴿᴬ - -# Full eigendecomposition - +# Compare to full eigendecomposition +# Note the this obtains eigenvectors from all QN +# sectors so will include vectors not found by +# `eigsolve` above, which only includes eigenvectors +# in the trivial QN sector. Tfull = prod(T) DV = eigen(Tfull, input_inds(T), output_inds(T)) +println("\nCheck full diagonalization on transfer matrix") @show norm(Tfull * DV.V - DV.Vt * DV.D) d = diag(array(DV.D)) @@ -192,6 +154,3 @@ d = diag(array(DV.D)) p = sortperm(d; by=abs, rev=true) @show p[1:neigs] @show d[p[1:neigs]] - -println("Error if ED with Arpack") -@show d[p[1:neigs]] - λ⃗ᴿᴬ diff --git a/examples/vumps/vumps_heisenberg.jl b/examples/vumps/vumps_heisenberg.jl index d042862..d2842a6 100644 --- a/examples/vumps/vumps_heisenberg.jl +++ b/examples/vumps/vumps_heisenberg.jl @@ -1,6 +1,12 @@ using ITensors using ITensorInfiniteMPS +include( + joinpath( + pkgdir(ITensorInfiniteMPS), "examples", "vumps", "src", "vumps_subspace_expansion.jl" + ), +) + ############################################################################## # VUMPS parameters # @@ -8,8 +14,9 @@ using ITensorInfiniteMPS maxdim = 100 # Maximum bond dimension cutoff = 1e-8 # Singular value cutoff when increasing the bond dimension max_vumps_iters = 100 # Maximum number of iterations of the VUMPS algorithm at each bond dimension -vumps_tol = 1e-6 -outer_iters = 8 # Number of times to increase the bond dimension +vumps_tol = 1e-5 +conserve_qns = true +outer_iters = 6 # Number of times to increase the bond dimension ############################################################################## # CODE BELOW HERE DOES NOT NEED TO BE MODIFIED @@ -17,11 +24,8 @@ outer_iters = 8 # Number of times to increase the bond dimension N = 2 # Number of sites in the unit cell -heisenberg_space_shift(q̃nf, q̃sz) = [QN("Sz", 1 - q̃sz) => 1, QN("Sz", -1 - q̃sz) => 1] - -heisenberg_space = fill(heisenberg_space_shift(1, 0), N) -s = infsiteinds("S=1/2", N; space=heisenberg_space) initstate(n) = isodd(n) ? "↑" : "↓" +s = infsiteinds("S=1/2", N; conserve_qns, initstate) ψ = InfMPS(s, initstate) model = Model("heisenberg") @@ -30,35 +34,26 @@ model = Model("heisenberg") H = InfiniteSum{MPO}(model, s) # Check translational invariance +println("\nCheck translation invariance of the initial VUMPS state") @show norm(contract(ψ.AL[1:N]..., ψ.C[N]) - contract(ψ.C[0], ψ.AR[1:N]...)) -vumps_kwargs = (tol=vumps_tol, maxiter=max_vumps_iters, eigsolve_tol=(x -> x / 1000)) +vumps_kwargs = (tol=vumps_tol, maxiter=max_vumps_iters, solver_tol=(x -> x / 1000)) subspace_expansion_kwargs = (cutoff=cutoff, maxdim=maxdim) -ψ = vumps(H, ψ; vumps_kwargs...) - -# Alternate steps of running VUMPS and increasing the bond dimension -for _ in 1:outer_iters - println("\nIncrease bond dimension") - global ψ = subspace_expansion(ψ, H; subspace_expansion_kwargs...) - println("Run VUMPS with new bond dimension") - global ψ = vumps(H, ψ; vumps_kwargs...) -end + +ψ = vumps_subspace_expansion(H, ψ; outer_iters, subspace_expansion_kwargs, vumps_kwargs) # Check translational invariance +println("\nCheck translation invariance of the final VUMPS state") @show norm(contract(ψ.AL[1:N]..., ψ.C[N]) - contract(ψ.C[0], ψ.AR[1:N]...)) -function ITensors.expect(ψ::InfiniteCanonicalMPS, o, n) - return (noprime(ψ.AL[n] * ψ.C[n] * op(o, s[n])) * dag(ψ.AL[n] * ψ.C[n]))[] -end - function expect_two_site(ψ::InfiniteCanonicalMPS, h::ITensor, n1n2) n1, n2 = n1n2 ϕ = ψ.AL[n1] * ψ.AL[n2] * ψ.C[n2] - return (noprime(ϕ * h) * dag(ϕ))[] + return inner(ϕ, apply(h, ϕ)) end function expect_two_site(ψ::InfiniteCanonicalMPS, h::MPO, n1n2) - return expect_two_site(ψ, prod(h), n1n2) + return expect_two_site(ψ, contract(h), n1n2) end Sz = [expect(ψ, "Sz", n) for n in 1:N] @@ -68,30 +63,24 @@ energy_infinite = map(b -> expect_two_site(ψ, H[b], b), bs) energy_exact = reference(model, Observable("energy")) -@show energy_infinite -@show energy_exact -@show Sz - # # Compare to DMRG # Nfinite = 100 -sfinite = siteinds("S=1/2", Nfinite; conserve_szparity=true) +sfinite = siteinds("S=1/2", Nfinite; conserve_qns) Hfinite = MPO(model, sfinite) ψfinite = randomMPS(sfinite, initstate; linkdims=10) @show flux(ψfinite) -sweeps = Sweeps(15) -setmaxdim!(sweeps, maxdim) -setcutoff!(sweeps, cutoff) -energy_finite_total, ψfinite = dmrg(Hfinite, ψfinite, sweeps) -@show energy_finite_total / Nfinite + +nsweeps = 10 +println("\nRunning finite DMRG for model $model on $Nfinite sites with $nsweeps") +energy_finite_total, ψfinite = dmrg(Hfinite, ψfinite; nsweeps, maxdim, cutoff) energy_exact_finite = reference(model, Observable("energy"); N=Nfinite) -@show energy_exact_finite -function ITensors.expect(ψ, o) - return (noprime(ψ * op(o, filterinds(ψ, "Site")...)) * dag(ψ))[] +function ITensors.expect(ψ::ITensor, o::String) + return inner(ψ, apply(op(o, filterinds(ψ, "Site")...), ψ)) end nfinite = Nfinite ÷ 2 @@ -100,4 +89,12 @@ Sz1_finite = expect(ψfinite[nfinite], "Sz") orthogonalize!(ψfinite, nfinite + 1) Sz2_finite = expect(ψfinite[nfinite + 1], "Sz") +println("\nCompare energy") +@show energy_finite_total / Nfinite +@show energy_infinite +@show energy_exact_finite +@show energy_exact + +println("\nCompare Sz") @show Sz1_finite, Sz2_finite +@show Sz diff --git a/examples/vumps/vumps_hubbard_extended.jl b/examples/vumps/vumps_hubbard_extended.jl index cffada3..cba2dd9 100644 --- a/examples/vumps/vumps_hubbard_extended.jl +++ b/examples/vumps/vumps_hubbard_extended.jl @@ -1,6 +1,12 @@ using ITensors using ITensorInfiniteMPS +include( + joinpath( + pkgdir(ITensorInfiniteMPS), "examples", "vumps", "src", "vumps_subspace_expansion.jl" + ), +) + ############################################################################## # VUMPS parameters # @@ -8,8 +14,10 @@ using ITensorInfiniteMPS maxdim = 50 # Maximum bond dimension cutoff = 1e-6 # Singular value cutoff when increasing the bond dimension max_vumps_iters = 200 # Maximum number of iterations of the VUMPS algorithm at each bond dimension +vumps_tol = 1e-5 outer_iters = 5 # Number of times to increase the bond dimension localham_type = MPO # or ITensor +conserve_qns = true eager = true model_params = (t=1.0, U=10.0, V=0.0) @@ -23,18 +31,8 @@ N = 2 # Unit cell size @show N @show localham_type -function electron_space_shift(q̃nf, q̃sz) - return [ - QN(("Nf", 0 - q̃nf, -1), ("Sz", 0 - q̃sz)) => 1, - QN(("Nf", 1 - q̃nf, -1), ("Sz", 1 - q̃sz)) => 1, - QN(("Nf", 1 - q̃nf, -1), ("Sz", -1 - q̃sz)) => 1, - QN(("Nf", 2 - q̃nf, -1), ("Sz", 0 - q̃sz)) => 1, - ] -end - -electron_space = fill(electron_space_shift(1, 0), N) -s = infsiteinds("Electron", N; space=electron_space) initstate(n) = isodd(n) ? "↑" : "↓" +s = infsiteinds("Electron", N; initstate, conserve_qns) ψ = InfMPS(s, initstate) model = Model"hubbard"() @@ -48,7 +46,7 @@ println("\nCheck translational invariance of initial infinite MPS") @show norm(contract(ψ.AL[1:N]..., ψ.C[N]) - contract(ψ.C[0], ψ.AR[1:N]...)) outputlevel = 1 -vumps_kwargs = (tol=1e-8, maxiter=max_vumps_iters, outputlevel=outputlevel, eager) +vumps_kwargs = (tol=vumps_tol, maxiter=max_vumps_iters, outputlevel, eager) subspace_expansion_kwargs = (cutoff=cutoff, maxdim=maxdim) # For now, to increase the bond dimension you must alternate @@ -57,23 +55,12 @@ subspace_expansion_kwargs = (cutoff=cutoff, maxdim=maxdim) # a larger bond dimension) println("\nRun VUMPS on initial product state, unit cell size $N") -ψ = @time vumps(H, ψ; vumps_kwargs...) - -@time for _ in 1:outer_iters - println("\nIncrease bond dimension") - global ψ = @time subspace_expansion(ψ, H; subspace_expansion_kwargs...) - println("Run VUMPS with new bond dimension") - global ψ = @time vumps(H, ψ; vumps_kwargs...) -end +ψ = vumps_subspace_expansion(H, ψ; outer_iters, subspace_expansion_kwargs, vumps_kwargs) # Check translational invariance println("\nCheck translational invariance of optimized infinite MPS") @show norm(contract(ψ.AL[1:N]..., ψ.C[N]) - contract(ψ.C[0], ψ.AR[1:N]...)) -function ITensors.expect(ψ::InfiniteCanonicalMPS, o, n) - return (noprime(ψ.AL[n] * ψ.C[n] * op(o, s[n])) * dag(ψ.AL[n] * ψ.C[n]))[] -end - function expect_two_site(ψ::InfiniteCanonicalMPS, h::ITensor, n1n2) n1, n2 = n1n2 ϕ = ψ.AL[n1] * ψ.AL[n2] * ψ.C[n2] @@ -103,19 +90,22 @@ energy_infinite = map(b -> expect_two_site(ψ, H[b], b), bs) # Nfinite = 100 -sfinite = siteinds("Electron", Nfinite; conserve_qns=true) +sfinite = siteinds("Electron", Nfinite; conserve_qns) Hfinite = MPO(model, sfinite; model_params...) ψfinite = randomMPS(sfinite, initstate; linkdims=10) println("\nQN sector of starting finite MPS") @show flux(ψfinite) -sweeps = Sweeps(10) + +nsweeps = 15 maxdims = min.(maxdim, [2, 2, 2, 2, 4, 4, 4, 4, 8, 8, 8, 8, 16, 16, 16, 16, 32, 32, 32, 32, 50]) @show maxdims -setmaxdim!(sweeps, maxdims...) -setcutoff!(sweeps, cutoff) + +## setmaxdim!(sweeps, maxdims...) +## setcutoff!(sweeps, cutoff) + println("\nRun DMRG on $Nfinite sites") -energy_finite_total, ψfinite = dmrg(Hfinite, ψfinite, sweeps) +energy_finite_total, ψfinite = dmrg(Hfinite, ψfinite; nsweeps, maxdims, cutoff) println("\nEnergy density") @show energy_finite_total / Nfinite diff --git a/examples/vumps/vumps_ising.jl b/examples/vumps/vumps_ising.jl index e54978c..6085045 100644 --- a/examples/vumps/vumps_ising.jl +++ b/examples/vumps/vumps_ising.jl @@ -1,5 +1,11 @@ -using ITensorInfiniteMPS using ITensors +using ITensorInfiniteMPS + +include( + joinpath( + pkgdir(ITensorInfiniteMPS), "examples", "vumps", "src", "vumps_subspace_expansion.jl" + ), +) ############################################################################## # VUMPS/TDVP parameters @@ -24,22 +30,10 @@ model_params = (J=1.0, h=0.9) # CODE BELOW HERE DOES NOT NEED TO BE MODIFIED # -model = Model"ising"() +model = Model("ising") -function space_shifted(::Model"ising", q̃sz; conserve_qns=true) - if conserve_qns - return [QN("SzParity", 1 - q̃sz, 2) => 1, QN("SzParity", 0 - q̃sz, 2) => 1] - else - return [QN() => 2] - end -end - -# Shift the QNs by 1 so that the flux of the unit cell is zero -# and therefore the flux density of the uniform state is zero -# to avoid diverging flux in the thermodynamic limit. -space_ = fill(space_shifted(model, 1; conserve_qns=conserve_qns), nsite) -s = infsiteinds("S=1/2", nsite; space=space_) initstate(n) = "↑" +s = infsiteinds("S=1/2", nsite; initstate, conserve_szparity=conserve_qns) ψ = InfMPS(s, initstate) # Form the Hamiltonian @@ -56,13 +50,7 @@ vumps_kwargs = ( ) subspace_expansion_kwargs = (cutoff=cutoff, maxdim=maxdim) -# Alternate steps of running VUMPS and increasing the bond dimension -@time for _ in 1:outer_iters - println("\nIncrease bond dimension") - global ψ = subspace_expansion(ψ, H; subspace_expansion_kwargs...) - println("Run VUMPS with new bond dimension") - global ψ = tdvp(H, ψ; time_step=time_step, vumps_kwargs...) -end +ψ = vumps_subspace_expansion(H, ψ; outer_iters, subspace_expansion_kwargs, vumps_kwargs) # Check translational invariance @show norm(contract(ψ.AL[1:nsite]..., ψ.C[nsite]) - contract(ψ.C[0], ψ.AR[1:nsite]...)) @@ -72,7 +60,7 @@ end # nsite_finite = 100 -s_finite = siteinds("S=1/2", nsite_finite; conserve_szparity=true) +s_finite = siteinds("S=1/2", nsite_finite; conserve_szparity=conserve_qns) H_finite = MPO(model, s_finite; model_params...) ψ_finite = randomMPS(s_finite, initstate) @show flux(ψ_finite) diff --git a/examples/vumps/vumps_ising_extended.jl b/examples/vumps/vumps_ising_extended.jl index 3061518..ed4af6f 100644 --- a/examples/vumps/vumps_ising_extended.jl +++ b/examples/vumps/vumps_ising_extended.jl @@ -1,32 +1,34 @@ using ITensors using ITensorInfiniteMPS +include( + joinpath( + pkgdir(ITensorInfiniteMPS), "examples", "vumps", "src", "vumps_subspace_expansion.jl" + ), +) + ############################################################################## # VUMPS parameters # maxdim = 64 # Maximum bond dimension cutoff = 1e-6 # Singular value cutoff when increasing the bond dimension -max_vumps_iters = 100 # Maximum number of iterations of the VUMPS algorithm at each bond dimension +max_vumps_iters = 10 # Maximum number of iterations of the VUMPS algorithm at each bond dimension vumps_tol = 1e-6 +conserve_qns = true outer_iters = 10 # Number of times to increase the bond dimension eager = true ############################################################################## # CODE BELOW HERE DOES NOT NEED TO BE MODIFIED # -N = 3# Number of sites in the unit cell +N = 2 # Number of sites in the unit cell J = -1.0 J₂ = -0.2 h = 1.0; -function space_shifted(q̃sz) - return [QN("SzParity", 1 - q̃sz, 2) => 1, QN("SzParity", 0 - q̃sz, 2) => 1] -end - -space_ = fill(space_shifted(1), N); -s = infsiteinds("S=1/2", N; space=space_) initstate(n) = "↑" +s = infsiteinds("S=1/2", N; initstate, conserve_szparity=conserve_qns) ψ = InfMPS(s, initstate); model = Model("ising_extended"); @@ -36,34 +38,28 @@ H = InfiniteSum{MPO}(model, s; J=J, J₂=J₂, h=h); vumps_kwargs = (tol=vumps_tol, maxiter=max_vumps_iters, eager) subspace_expansion_kwargs = (cutoff=cutoff, maxdim=maxdim) -ψ_0 = @time vumps(H, ψ; vumps_kwargs...) -@time for j in 1:outer_iters - println("\nIncrease bond dimension") - ψ_1 = @time subspace_expansion(ψ_0, H; subspace_expansion_kwargs...) - println("Run VUMPS with new bond dimension") - global ψ_0 = @time vumps(H, ψ_1; vumps_kwargs...) -end +ψ = vumps_subspace_expansion(H, ψ; outer_iters, subspace_expansion_kwargs, vumps_kwargs) -Sz = [expect(ψ_0, "Sz", n) for n in 1:N] -energy_infinite = expect(ψ_0, H) +Sz_infinite = [expect(ψ, "Sz", n) for n in 1:N] +energy_infinite = expect(ψ, H) Nfinite = 100 -sfinite = siteinds("S=1/2", Nfinite; conserve_szparity=true) +sfinite = siteinds("S=1/2", Nfinite; conserve_szparity=conserve_qns) Hfinite = MPO(model, sfinite; J=J, J₂=J₂, h=h) ψfinite = randomMPS(sfinite, initstate; linkdims=10) @show flux(ψfinite) -sweeps = Sweeps(10) -setmaxdim!(sweeps, maxdim) -setcutoff!(sweeps, cutoff) -energy_finite_total, ψfinite = dmrg(Hfinite, ψfinite, sweeps) +dmrg_kwargs = (nsweeps=10, maxdim, cutoff) +energy_finite_total, ψfinite = dmrg(Hfinite, ψfinite; dmrg_kwargs...) nfinite = Nfinite ÷ 2 Sz_finite = expect(ψfinite, "Sz")[nfinite:(nfinite + N - 1)] +println("\nEnergy") @show energy_infinite @show energy_finite_total / Nfinite @show reference(model, Observable("energy"); J=J, h=h, J₂=J₂) -@show Sz +println("\nSz") +@show Sz_infinite @show Sz_finite diff --git a/src/ITensorInfiniteMPS.jl b/src/ITensorInfiniteMPS.jl index a136574..996cfbb 100644 --- a/src/ITensorInfiniteMPS.jl +++ b/src/ITensorInfiniteMPS.jl @@ -16,13 +16,17 @@ using IterTools using HDF5 # For integration support when computing exact reference results using QuadGK +# For `groupreduce`, used when splitting up an OpSum +# with the unit cell terms into terms acting on each +# site. +using SplitApplyCombine using ITensors.NDTensors: eachdiagblock using KrylovKit: eigsolve, linsolve, exponentiate import Base: getindex, length, setindex!, +, -, * -import ITensors: AbstractMPS +import ITensors: AbstractMPS, ⊕ include("ITensors.jl") include("ITensorNetworks.jl") @@ -73,6 +77,7 @@ export Cell, reference, subspace_expansion, translatecell, + translatecelltags, translator, tdvp, vumps, diff --git a/src/abstractinfinitemps.jl b/src/abstractinfinitemps.jl index 2b32f16..a14c7c0 100644 --- a/src/abstractinfinitemps.jl +++ b/src/abstractinfinitemps.jl @@ -206,13 +206,23 @@ end # TODO: return a Dictionary or IndexSetNetwork? function ITensors.siteinds(ψ::AbstractInfiniteMPS, r::AbstractRange) - return [siteinds(ψ, n) for n in r] + return [siteind(ψ, n) for n in r] end -siterange(ψ::AbstractInfiniteMPS, c::Cell) = 1:(nsites(ψ) .+ (c.cell - 1)) +function siterange(ψ::AbstractInfiniteMPS, c::Cell) + #1:(nsites(ψ) .+ (c.cell - 1)) + start = nsites(ψ) * (c.cell - 1) + 1 + stop = start + nsites(ψ) - 1 + return start:stop +end ITensors.siteinds(ψ::AbstractInfiniteMPS, c::Cell) = siteinds(ψ, siterange(ψ, c)) +function ITensors.siteinds(ψ::AbstractInfiniteMPS) + return CelledVector(siteinds(ψ, Cell(1)), translator(ψ)) +end +infsiteinds(ψ::AbstractInfiniteMPS) = siteinds(ψ) + Base.getindex(ψ::AbstractInfiniteMPS, r::UnitRange{Int}) = MPS([ψ[n] for n in r]) struct UnitRangeToFunction{T<:Real,F<:Function} diff --git a/src/celledvectors.jl b/src/celledvectors.jl index ab2d6c2..16a8ad3 100644 --- a/src/celledvectors.jl +++ b/src/celledvectors.jl @@ -16,11 +16,9 @@ indextagprefix() = "n=" # translatecell # -# TODO: account for shifting by a tuple, for example: -# translatecell(ts"Site,c=1|2", (2, 3)) -> ts"Site,c=3|5" -# TODO: ts"c=10|12" has too many characters -# TODO: ts"c=1|2|3" has too many characters -# +# TODO: account for shifting by a tuple for multidimensional +# indexing, for example: +# translatecell(ts"Site,c=1×2", (2, 3)) -> ts"Site,c=3×5" # Determine the cell `n` from the tag `"c=n"` function getcell(ts::TagSet) @@ -67,8 +65,9 @@ function translatecell(translator::Function, is::Union{<:Tuple,<:Vector}, n::Int return translatecell.(translator, is, n) end -#Default behavior -#translatecell(T::ITensor, n::Integer) = ITensors.setinds(T, translatecell(inds(T), n)) +# Default behavior +translatecelltags(T::ITensor, n::Integer) = translatecell(translatecelltags, T, n) +translatecelltags(T::ITensors.Indices, n::Integer) = translatecell(translatecelltags, T, n) #translatecell(T::MPO, n::Integer) = translatecell.(T, n) #translatecell(T::Matrix{ITensor}, n::Integer) = translatecell.(T, n) @@ -99,11 +98,15 @@ function CelledVector{T}(v::Vector{T}, translator::Function) where {T} end """ - celllength(cv::CelledVector) + cell_length(cv::CelledVector) + + celllength(cv::CelledVector) # Deprecated The length of a unit cell of a CelledVector. """ -celllength(cv::CelledVector) = length(ITensors.data(cv)) +cell_length(cv::CelledVector) = length(ITensors.data(cv)) + +celllength(cv::CelledVector) = cell_length(cv) # For compatibility with Base Base.size(cv::CelledVector) = size(ITensors.data(cv)) diff --git a/src/infinitecanonicalmps.jl b/src/infinitecanonicalmps.jl index 3f1f04e..2831ba0 100644 --- a/src/infinitecanonicalmps.jl +++ b/src/infinitecanonicalmps.jl @@ -1,4 +1,22 @@ -# Get the promoted type of the Index objects in a collection +# TODO: Move to ITensors.jl +function setval(qnval::ITensors.QNVal, val::Int) + return ITensors.QNVal(ITensors.name(qnval), val, ITensors.modulus(qnval)) +end + +# TODO: Move to ITensors.jl +function Base.:/(qnval::ITensors.QNVal, n::Int) + div_val = ITensors.val(qnval) / n + if !isinteger(div_val) + error("Dividing $qnval by $n, the resulting QN value is not an integer") + end + return setval(qnval, Int(div_val)) +end + +# TODO: Move to ITensors.jl +function Base.:/(qn::QN, n::Int) + return QN(map(qnval -> qnval / n, qn.data)) +end + # of Index (Tuple, Vector, ITensor, etc.) indtype(i::Index) = typeof(i) indtype(T::Type{<:Index}) = T @@ -9,39 +27,42 @@ indtype(A::ITensor...) = indtype(inds.(A)) indtype(tn1, tn2) = promote_type(indtype(tn1), indtype(tn2)) indtype(tn) = mapreduce(indtype, promote_type, tn) -# More general siteind that allows specifying -# the space -function _siteind(site_tag, n::Int; space) - return addtags(Index(space, "Site,n=$n"), site_tag) +function infsiteinds(s::Vector{<:Index}, translator=translatecelltags) + return CelledVector(addtags(s, celltags(1)), translator) end -_siteinds(site_tag, N::Int; space) = __siteinds(site_tag, N, space) +shift_flux_to_zero(s::Vector{Index{Int}}, initestate::Function) = s +shift_flux_to_zero(s::Vector{Index{Int}}, flux_density::QN) = s -function __siteinds(site_tag, N::Int, space::Vector) - return [_siteind(site_tag, n; space=space[n]) for n in 1:N] +function shift_flux_to_zero(s::Vector{<:Index}, initstate::Function) + return shift_flux_to_zero(s, flux(MPS(s, initstate))) end -function __siteinds(site_tag, N::Int, space) - return [_siteind(site_tag, n; space=space) for n in 1:N] +function shift_flux(qnblock::Pair{QN,Int}, flux_density::QN) + return ((ITensors.qn(qnblock) - flux_density) => ITensors.blockdim(qnblock)) end - -infsiteinds(s::Vector{<:Index}) = CelledVector(addtags(s, celltags(1))) -function infsiteinds(s::Vector{<:Index}, translator) - return CelledVector(addtags(s, celltags(1)), translator) +function shift_flux(space::Vector{Pair{QN,Int}}, flux_density::QN) + return map(qnblock -> shift_flux(qnblock, flux_density), space) +end +function shift_flux(i::Index, flux_density::QN) + return ITensors.setspace(i, shift_flux(space(i), flux_density)) end -function infsiteinds(site_tag, N::Int; space=nothing, translator=nothing, kwargs...) - if !isnothing(space) - s = _siteinds(site_tag, N; space=space) - else - # TODO: add a shift option - s = siteinds(site_tag, N; kwargs...) - end - if !isnothing(translator) - return infsiteinds(s, translator) - else - return infsiteinds(s) +function shift_flux_to_zero(s::Vector{<:Index}, flux::QN) + if iszero(flux) + return s end + n = length(s) + flux_density = flux / n + return map(sₙ -> shift_flux(sₙ, flux_density), s) +end + +function infsiteinds( + site_tag, n::Int; translator=translatecelltags, initstate=nothing, kwargs... +) + s = siteinds(site_tag, n; kwargs...) + s = shift_flux_to_zero(s, initstate) + return infsiteinds(s, translator) end function ITensors.linkinds(ψ::InfiniteMPS) @@ -49,11 +70,7 @@ function ITensors.linkinds(ψ::InfiniteMPS) return CelledVector([linkinds(ψ, (n, n + 1)) for n in 1:N], translator(ψ)) end -function InfMPS(s::Vector, f::Function) - return InfMPS(infsiteinds(s), f) -end - -function InfMPS(s::Vector, f::Function, translator::Function) +function InfMPS(s::Vector, f::Function, translator::Function=translatecelltags) return InfMPS(infsiteinds(s, translator), f) end @@ -61,6 +78,8 @@ function indval(iv::Pair) return ind(iv) => val(iv) end +zero_qn(i::Index{Int}) = nothing + function zero_qn(i::Index) return zero(qn(first(space(i)))) end @@ -129,9 +148,11 @@ function InfMPS(s::CelledVector, f::Function) return ψ = InfiniteCanonicalMPS(ψL, ψC, ψR) end -function ITensors.expect(ψ::InfiniteCanonicalMPS, o, n) +function ITensors.expect(ψ::InfiniteCanonicalMPS, o::String, n::Int) s = siteinds(only, ψ.AL) - return (noprime(ψ.AL[n] * ψ.C[n] * op(o, s[n])) * dag(ψ.AL[n] * ψ.C[n]))[] + O = op(o, s[n]) + ϕ = ψ.AL[n] * ψ.C[n] + return inner(ϕ, apply(O, ϕ)) end function ITensors.expect(ψ::InfiniteCanonicalMPS, h::MPO) diff --git a/src/infinitempo.jl b/src/infinitempo.jl index 499967e..1c2c75f 100644 --- a/src/infinitempo.jl +++ b/src/infinitempo.jl @@ -10,13 +10,4 @@ mutable struct InfiniteMPO <: AbstractInfiniteMPS reverse::Bool end -## struct InfiniteSumLocalOps <: AbstractInfiniteMPS -## data::CelledVector{ITensor} -## end -## InfiniteSumLocalOps(N::Int) = InfiniteSumLocalOps(Vector{ITensor}(undef, N)) -## InfiniteSumLocalOps(data::Vector{ITensor}) = InfiniteSumLocalOps(CelledVector(data)) -## getindex(l::InfiniteSumLocalOps, n::Integer) = ITensors.data(l)[n] - -# TODO? Instead of having a big quasi empty ITensor, store only the non zero blocks - translator(mpo::InfiniteMPO) = mpo.data.translator diff --git a/src/infinitemps_approx.jl b/src/infinitemps_approx.jl index 44aa6d5..c2fa16d 100644 --- a/src/infinitemps_approx.jl +++ b/src/infinitemps_approx.jl @@ -1,4 +1,3 @@ - function siteind(ψ::InfiniteMPS, n::Integer) return uniqueind(ψ[n], ψ[n - 1], ψ[n + 1]) end @@ -37,7 +36,11 @@ function infinitemps_approx( A = InfiniteMPS( [settags(s[nrange[n]], addtags(site_tags, "n=$n")) for n in 1:length(nrange)]; space=χ ) - randn!.(A) + # XXX: Why doesn't this work? + # randn!.(A) + for n in 1:N + A[n] = randomITensor(inds(A[n])) + end ψ∞ = orthogonalize(A, :) # Site indices of the infinite MPS diff --git a/src/models/fqhe13.jl b/src/models/fqhe13.jl index 431a1ce..40d87e6 100644 --- a/src/models/fqhe13.jl +++ b/src/models/fqhe13.jl @@ -1,34 +1,14 @@ -# H = -J Σⱼ XⱼXⱼ₊₁ - h Σⱼ Zⱼ -function ITensors.OpSum( - ::Model"fqhe_2b_pot", n1, n2; Ly::Float64, Vs::Array{Float64,1}, prec::Float64 +function unit_cell_terms( + ::Model"fqhe_2b_pot"; Ly::Float64, Vs::Array{Float64,1}, prec::Float64 ) rough_N = round(Int64, 2 * Ly) coeff = build_two_body_coefficient_pseudopotential(; N_phi=rough_N, Ly=Ly, Vs=Vs) opt = optimize_coefficients(coeff; prec=prec) opt = filter_optimized_Hamiltonian_by_first_site(opt) #sorted_opt = sort_by_configuration(opt); - #println(opt) return generate_Hamiltonian(opt) end -function ITensors.MPO(::Model"fqhe_2b_pot", s::CelledVector, n::Int64; Ly, Vs, prec) - rough_N = round(Int64, 2 * Ly) - coeff = build_two_body_coefficient_pseudopotential(; N_phi=rough_N, Ly=Ly, Vs=Vs) - opt = optimize_coefficients(coeff; prec=prec) - opt = filter_optimized_Hamiltonian_by_first_site(opt) - range_model = check_max_range_optimized_Hamiltonian(opt) - while range_model >= rough_N - rough_N = 2 * (rough_N + 1) - coeff = build_two_body_coefficient_pseudopotential(; N_phi=rough_N, Ly=Ly, Vs=Vs) - opt = optimize_coefficients(coeff; prec=prec) - opt = filter_optimized_Hamiltonian_by_first_site(opt) - range_model = check_max_range_optimized_Hamiltonian(opt) - end - opsum = generate_Hamiltonian(opt) - - return splitblocks(linkinds, MPO(opsum, [s[x] for x in n:(n + range_model)])) -end - #Please contact Loic Herviou before using this part of the code for production # loic.herviou@epfl.ch ############################### diff --git a/src/models/heisenberg.jl b/src/models/heisenberg.jl index 87f38ab..2d9a338 100644 --- a/src/models/heisenberg.jl +++ b/src/models/heisenberg.jl @@ -1,23 +1,12 @@ -function ITensors.OpSum(::Model"heisenberg", n1, n2) +# H = Σⱼ (½ S⁺ⱼS⁻ⱼ₊₁ + ½ S⁻ⱼS⁺ⱼ₊₁ + SᶻⱼSᶻⱼ₊₁) +function unit_cell_terms(::Model"heisenberg") opsum = OpSum() - opsum += 0.5, "S+", n1, "S-", n2 - opsum += 0.5, "S-", n1, "S+", n2 - opsum += "Sz", n1, "Sz", n2 + opsum += 0.5, "S+", 1, "S-", 2 + opsum += 0.5, "S-", 1, "S+", 2 + opsum += "Sz", 1, "Sz", 2 return opsum end -# H = Σⱼ (½ S⁺ⱼS⁻ⱼ₊₁ + ½ S⁻ⱼS⁺ⱼ₊₁ + SᶻⱼSᶻⱼ₊₁) -function ITensors.MPO(::Model"heisenberg", s) - N = length(s) - os = OpSum() - for j in 1:(N - 1) - os .+= 0.5, "S+", j, "S-", j + 1 - os .+= 0.5, "S-", j, "S+", j + 1 - os .+= "Sz", j, "Sz", j + 1 - end - return splitblocks(linkinds, MPO(os, s)) -end - """ reference(::Model"heisenberg", ::Observable"energy"; N) diff --git a/src/models/hubbard.jl b/src/models/hubbard.jl index 8528371..77960e1 100644 --- a/src/models/hubbard.jl +++ b/src/models/hubbard.jl @@ -1,43 +1,14 @@ -function ITensors.OpSum(::Model"hubbard", n1, n2; t, U, V) +function unit_cell_terms(::Model"hubbard"; t, U, V=0.0) opsum = OpSum() - opsum += -t, "Cdagup", n1, "Cup", n2 - opsum += -t, "Cdagup", n2, "Cup", n1 - opsum += -t, "Cdagdn", n1, "Cdn", n2 - opsum += -t, "Cdagdn", n2, "Cdn", n1 - if U ≠ 0 - opsum += U, "Nupdn", n1 - end - if V ≠ 0 - opsum += V, "Ntot", n1, "Ntot", n2 - end + opsum += -t, "Cdagup", 1, "Cup", 2 + opsum += -t, "Cdagup", 2, "Cup", 1 + opsum += -t, "Cdagdn", 1, "Cdn", 2 + opsum += -t, "Cdagdn", 2, "Cdn", 1 + opsum += U, "Nupdn", 1 + opsum += V, "Ntot", 1, "Ntot", 2 return opsum end -function ITensors.MPO(::Model"hubbard", s; t, U, V) - N = length(s) - opsum = OpSum() - for n in 1:(N - 1) - n1, n2 = n, n + 1 - opsum .+= -t, "Cdagup", n1, "Cup", n2 - opsum .+= -t, "Cdagup", n2, "Cup", n1 - opsum .+= -t, "Cdagdn", n1, "Cdn", n2 - opsum .+= -t, "Cdagdn", n2, "Cdn", n1 - if V ≠ 0 - opsum .+= V, "Ntot", n1, "Ntot", n2 - end - end - if U ≠ 0 - for n in 1:N - opsum .+= U, "Nupdn", n - end - end - return splitblocks(linkinds, MPO(opsum, s)) -end - -function ITensors.ITensor(model::Model"hubbard", s; kwargs...) - return prod(MPO(model, s; kwargs...)) -end - """ @article{PhysRevB.6.930, title = {Magnetic Susceptibility at Zero Temperature for the One-Dimensional Hubbard Model}, diff --git a/src/models/ising.jl b/src/models/ising.jl index fdaada6..960638a 100644 --- a/src/models/ising.jl +++ b/src/models/ising.jl @@ -1,40 +1,11 @@ # H = -J Σⱼ XⱼXⱼ₊₁ - h Σⱼ Zⱼ -function ITensors.MPO(::Model"ising", s; J, h) - N = length(s) - a = OpSum() - for n in 1:(N - 1) - a .+= -J, "X", n, "X", n + 1 - end - for n in 1:N - a .+= -h, "Z", n - end - return splitblocks(linkinds, MPO(a, s)) -end - -# H = -J Σⱼ XⱼXⱼ₊₁ - h Σⱼ Zⱼ -function ITensors.OpSum(::Model"ising", n1, n2; J, h) +function unit_cell_terms(::Model"ising"; J=1.0, h=1.0) opsum = OpSum() - if J != 0 - opsum += -J, "X", n1, "X", n2 - end - if h != 0 - opsum += -h / 2, "Z", n1 - opsum += -h / 2, "Z", n2 - end + opsum += -J, "X", 1, "X", 2 + opsum += -h, "Z", 1 return opsum end -# H = -J X₁X₂ - h Z₁ -# XXX: use `op` instead of `ITensor` -function ITensors.ITensor(::Model"ising", s1::Index, s2::Index; J, h) - # -J * op("X", s1) * op("X", s2) - h * op("Z", s1) * op("Id", s2) - opsum = OpSum() - n = 1 - opsum += -J, "X", n, "X", n + 1 - opsum += -h, "Z", n - return prod(MPO(opsum, [s1, s2])) -end - # P. Pfeuty, The one-dimensional Ising model with a transverse field, Annals of Physics 57, p. 79 (1970) function reference(::Model"ising", ::Observable"energy"; h, J=1.0) (J == h) && return -4 / π diff --git a/src/models/ising_extended.jl b/src/models/ising_extended.jl index 929bd8e..e234efb 100644 --- a/src/models/ising_extended.jl +++ b/src/models/ising_extended.jl @@ -1,47 +1,13 @@ -nrange(::Model"ising_extended") = 3 +# The terms of the Hamiltonian in the first unit cell. # H = -J Σⱼ XⱼXⱼ₊₁ - h Σⱼ Zⱼ - J₂ Σⱼ XⱼZⱼ₊₁Xⱼ₊₂ -function ITensors.MPO(::Model"ising_extended", s; J=1.0, h=1.0, J₂=0.0) - N = length(s) - a = OpSum() - if J != 0 - for n in 1:(N - 1) - a += -J, "X", n, "X", n + 1 - end - end - if J₂ != 0 - for n in 1:(N - 2) - a += -J₂, "X", n, "Z", n + 1, "X", n + 2 - end - end - if h != 0 - for n in 1:N - a += -h, "Z", n - end - end - return splitblocks(linkinds, MPO(a, s)) -end - -# H = -J Σⱼ XⱼXⱼ₊₁ - h Σⱼ Zⱼ - J₂ Σⱼ XⱼZⱼ₊₁Xⱼ₊₂ -function ITensors.OpSum(::Model"ising_extended", n1, n2; J=1.0, h=1.0, J₂=0.0) +function unit_cell_terms(::Model"ising_extended"; J=1.0, h=1.0, J₂=0.0) opsum = OpSum() - if J != 0 - opsum += -J / 2, "X", n1, "X", n2 - opsum += -J / 2, "X", n2, "X", n2 + 1 - end - if J₂ != 0 - opsum += -J₂, "X", n1, "Z", n2, "X", n2 + 1 - end - opsum += -h / 3, "Z", n1 - opsum += -h / 3, "Z", n2 - opsum += -h / 3, "Z", n2 + 1 + opsum += -J, "X", 1, "X", 2 + opsum += -J₂, "X", 1, "Z", 2, "X", 3 + opsum += -h, "Z", 1 return opsum end -function ITensors.ITensor(::Model"ising_extended", s1, s2, s3; J=1.0, h=1.0, J₂=0.0) - opsum = OpSum(Model"ising_extended"(), 1, 2; J=J, h=h, J₂=J₂) - return prod(MPO(opsum, [s1, s2, s3])) -end - function reference(::Model"ising_extended", ::Observable"energy"; J=1.0, h=1.0, J₂=0.0) f(k) = sqrt((J * cos(k) + J₂ * cos(2k) - h)^2 + (J * sin(k) + J₂ * sin(2k))^2) return -1 / 2π * ITensorInfiniteMPS.∫(k -> f(k), -π, π) diff --git a/src/models/models.jl b/src/models/models.jl index 24baf98..75d5085 100644 --- a/src/models/models.jl +++ b/src/models/models.jl @@ -15,7 +15,6 @@ macro Observable_str(s) end ∫(f, a, b) = quadgk(f, a, b)[1] -ITensorInfiniteMPS.nrange(model::Model) = 2; #required to keep everything compatible with the current implementation for 2 band models # Create an infinite sum of Hamiltonian terms function InfiniteSum{T}(model::Model, s::Vector; kwargs...) where {T} @@ -23,29 +22,62 @@ function InfiniteSum{T}(model::Model, s::Vector; kwargs...) where {T} end function InfiniteSum{MPO}(model::Model, s::CelledVector; kwargs...) - N = length(s) - mpos = [splitblocks(linkinds, MPO(model, s, n; kwargs...)) for n in 1:N] #slightly improved version. Note: the current implementation does not really allow for staggered potentials for example - return InfiniteSum{MPO}(mpos, translator(s)) + return InfiniteSum{MPO}(opsum_infinite(model, cell_length(s); kwargs...), s) end function InfiniteSum{ITensor}(model::Model, s::CelledVector; kwargs...) N = length(s) - itensors = [ITensor(model, s, n; kwargs...) for n in 1:N] #slightly improved version. Note: the current implementation does not really allow for staggered potentials for example + itensors = [ITensor(model, s, n; kwargs...) for n in 1:N] return InfiniteSum{ITensor}(itensors, translator(s)) end -# MPO building version -function ITensors.MPO(model::Model, s::CelledVector, n::Int64; kwargs...) - n1, n2 = 1, 2 - opsum = OpSum(model, n1, n2; kwargs...) - return splitblocks(linkinds, MPO(opsum, [s[x] for x in n:(n + nrange(model) - 1)])) #modification to allow for more than two sites per term in the Hamiltonians +# Get the first site with nontrivial support of the OpSum +first_site(opsum::OpSum) = minimum(ITensors.sites(opsum)) +last_site(opsum::OpSum) = maximum(ITensors.sites(opsum)) + +function set_site(o::Op, s::Int) + return Op(ITensors.which_op(o), s; ITensors.params(o)...) end -# Version accepting IndexSet -function ITensors.ITensor(model::Model, s::CelledVector, n::Int64; kwargs...) - return prod(MPO(model, s, n; kwargs...)) #modification to allow for more than two sites per term in the Hamiltonians +function shift_site(o::Op, shift::Int) + return set_site(o, ITensors.site(o) + shift) +end + +function shift_sites(term::Scaled{C,Prod{Op}}, shift::Int) where {C} + shifted_term = ITensors.coefficient(term) + for o in ITensors.terms(term) + shifted_term *= shift_site(o, shift) + end + return shifted_term end +# Shift the sites of the terms of the OpSum by shift. +# By default, it shifts +function shift_sites(opsum::OpSum, shift::Int) + shifted_opsum = OpSum() + for o in ITensors.terms(opsum) + shifted_opsum += shift_sites(o, shift) + end + return shifted_opsum +end + +function InfiniteSum{MPO}(opsum::OpSum, s::CelledVector) + n = cell_length(s) + nrange = 0 # Maximum operator support + opsums = [OpSum() for _ in 1:n] + for o in ITensors.terms(opsum) + js = sort(ITensors.sites(o)) + j1 = first(js) + nrange = max(nrange, last(js) - j1 + 1) + opsums[j1] += o + end + shifted_opsums = [shift_sites(opsum, -first_site(opsum) + 1) for opsum in opsums] + mpos = [ + splitblocks(linkinds, MPO(shifted_opsums[j], [s[k] for k in j:(j + nrange - 1)])) for + j in 1:n + ] + return InfiniteSum{MPO}(mpos, translator(s)) +end # Helper function to make an MPO import ITensors: op op(::OpName"Zero", ::SiteType, s::Index) = ITensor(s', dag(s)) @@ -79,3 +111,118 @@ function InfiniteMPOMatrix(model::Model, s::CelledVector, translator::Function; #return mpos return InfiniteMPOMatrix(mpos, translator) end + +function ITensors.MPO(model::Model, s::Vector{<:Index}; kwargs...) + opsum = opsum_finite(model, length(s); kwargs...) + return splitblocks(linkinds, MPO(opsum, s)) +end + +translatecell(translator, opsum::OpSum, n::Integer) = translator(opsum, n) + +function infinite_terms(model::Model; kwargs...) + # An `OpSum` storing all of the terms in the + # first unit cell. + # TODO: Allow specifying the unit cell size + # explicitly. + return infinite_terms(unit_cell_terms(model; kwargs...)) +end + +function infinite_terms(opsum::OpSum; kwargs...) + # `Vector{OpSum}`, vector of length of number + # of sites in the unit cell where each element + # contains the terms with support starting on that + # site of the unit cell, i.e. `opsum_cell[i]` + # stores all terms starting on site `i`. + opsum_cell_dict = groupreduce(minimum ∘ ITensors.sites, +, opsum) + nsites = maximum(keys(opsum_cell_dict)) + # Assumes each site in the unit cell has a term + for j in 1:nsites + if !haskey(opsum_cell_dict, j) + error( + "The input unit cell terms for the $nsites-site unit cell doesn't have a term starting on site $j. Skipping sites in the unit cell is currently not supported. A workaround is to define a term in the unit cell with a coefficient of 0, for example `opsum += 0, \"I\", $j`.", + ) + end + end + opsum_cell = [opsum_cell_dict[j] for j in 1:nsites] + function _shift_cell(opsum::OpSum, cell::Int) + return shift_sites(opsum, nsites * cell) + end + return CelledVector(opsum_cell, _shift_cell) +end + +function opsum_infinite(model::Model, nsites::Int; kwargs...) + _infinite_terms = infinite_terms(model::Model; kwargs...) + # Desired unit cell size must be commensurate + # with the primitive unit cell of the model. + if !iszero(nsites % length(_infinite_terms)) + error( + "Desired unit cell size $nsites must be commensurate with the primitive unit cell size of the model $model, which is $(length(_infinite_terms))", + ) + end + opsum = OpSum() + for j in 1:nsites + opsum += _infinite_terms[j] + end + return opsum +end + +function filter_terms(f, opsum::OpSum; by=identity) + filtered_opsum = OpSum() + for t in ITensors.terms(opsum) + if f(by(t)) + filtered_opsum += t + end + end + return filtered_opsum +end + +function finite_terms(model::Model, n::Int; kwargs...) + _infite_terms = infinite_terms(model; kwargs...) + _finite_terms = OpSum[] + for j in 1:n + term_j = _infite_terms[j] + filtered_term_j = filter_terms(s -> all(≤(n), s), term_j; by=ITensors.sites) + push!(_finite_terms, filtered_term_j) + end + return _finite_terms +end + +# For finite Hamiltonian with open boundary conditions +# Obtain from infinite Hamiltonian, dropping terms +# that extend outside of the system. +function opsum_finite(model::Model, n::Int; kwargs...) + opsum = OpSum() + for term in finite_terms(model, n; kwargs...) + opsum += term + end + return opsum +end + +# The ITensor of a single term `n` of the model. +function ITensors.ITensor(model::Model, s::Vector{<:Index}, n::Int; kwargs...) + opsum = infinite_terms(model; kwargs...)[n] + opsum = shift_sites(opsum, -first_site(opsum) + 1) + return contract(MPO(opsum, [s...])) +end + +function ITensors.ITensor(model::Model, s::Vector{<:Index}; kwargs...) + return ITensor(model, s, 1; kwargs...) +end + +function ITensors.ITensor(model::Model, n::Int, s::Index...; kwargs...) + return ITensor(model, [s...], n; kwargs...) +end + +function ITensors.ITensor(model::Model, s::Index...; kwargs...) + return ITensor(model, 1, s...; kwargs...) +end + +function ITensors.ITensor(model::Model, s::CelledVector, n::Int64; kwargs...) + opsum = infinite_terms(model; kwargs...)[n] + opsum = shift_sites(opsum, -first_site(opsum) + 1) + site_range = n:(n + last_site(opsum) - 1) + return contract(MPO(opsum, [s[j] for j in site_range])) + + # Deprecated version + # return contract(MPO(model, s, n; kwargs...)) +end diff --git a/src/models/xx.jl b/src/models/xx.jl index 89a9407..061a49d 100644 --- a/src/models/xx.jl +++ b/src/models/xx.jl @@ -1,33 +1,9 @@ # H = Σⱼ XⱼXⱼ₊₁ + YⱼYⱼ₊₁ -function ITensors.MPO(::Model"xx", s) - N = length(s) +function unit_cell_terms(::Model"xx") os = OpSum() - for n in 1:(N - 1) - os .+= 1, "X", n, "X", n + 1 - os .+= 1, "Y", n, "Y", n + 1 - end - return splitblocks(linkinds, MPO(os, s)) -end - -# H = X₁X₂ + Y₁Y₂ -function ITensors.OpSum(::Model"xx", n1, n2) - opsum = OpSum() - for j in 1:(N - 1) - opsum .+= 0.5, "S+", j, "S-", j + 1 - opsum .+= 0.5, "S-", j, "S+", j + 1 - end - return opsum -end - -# H = X₁X₂ + Y₁Y₂ -# XXX: use `op` instead of `ITensor` -function ITensors.ITensor(::Model"xx", s1::Index, s2::Index) - # op("X", s1) * op("X", s2) + op("Y", s1) * op("Y", s2) - opsum = OpSum() - n = 1 - opsum += "X", n, "X", n + 1 - opsum += "Y", n, "Y", n + 1 - return prod(MPO(opsum, [s1, s2])) + os += "X", 1, "X", 2 + os += "Y", 1, "Y", 2 + return os end function reference(::Model"xx", ::Observable"energy"; N=∞) diff --git a/src/models/xx_extended.jl b/src/models/xx_extended.jl index cb94c7f..9e6939e 100644 --- a/src/models/xx_extended.jl +++ b/src/models/xx_extended.jl @@ -1,50 +1,18 @@ # H = JΣⱼ (½ S⁺ⱼS⁻ⱼ₊₁ + ½ S⁻ⱼS⁺ⱼ₊₁) + J₂Σⱼ (½ S⁺ⱼZⱼ₊₁S⁻ⱼ₊₂ + ½ S⁻ⱼZⱼ₊₁S⁺ⱼ₊₂) - h Σⱼ Sⱼᶻ -function ITensors.OpSum(::Model"xx_extended", n1, n2; J=1.0, J₂=1.0, h=0.0) +function unit_cell_terms(::Model"xx_extended"; J=1.0, J₂=1.0, h=0.0) opsum = OpSum() - if J != 0.0 - opsum += 0.25 * J, "S+", n1, "S-", n2 - opsum += 0.25 * J, "S-", n1, "S+", n2 - opsum += 0.25 * J, "S+", n2, "S-", n2 + 1 - opsum += 0.25 * J, "S-", n2, "S+", n2 + 1 - end - if J₂ != 0 - opsum += 0.5 * J₂, "S+", n1, "Sz", n2, "S-", n2 + 1 - opsum += 0.5 * J₂, "S-", n1, "Sz", n2, "S+", n2 + 1 - end - if h != 0 - opsum += -h / 3, "Sz", n1 - opsum += -h / 3, "Sz", n2 - opsum += -h / 3, "Sz", n2 + 1 - end + opsum += 0.25 * J, "S+", 1, "S-", 2 + opsum += 0.25 * J, "S-", 1, "S+", 2 + opsum += 0.25 * J, "S+", 2, "S-", 2 + 1 + opsum += 0.25 * J, "S-", 2, "S+", 2 + 1 + opsum += 0.5 * J₂, "S+", 1, "Sz", 2, "S-", 3 + opsum += 0.5 * J₂, "S-", 1, "Sz", 2, "S+", 3 + opsum += -h / 3, "Sz", 1 + opsum += -h / 3, "Sz", 2 + opsum += -h / 3, "Sz", 3 return opsum end -# H = JΣⱼ (½ S⁺ⱼS⁻ⱼ₊₁ + ½ S⁻ⱼS⁺ⱼ₊₁) + J₂Σⱼ (½ S⁺ⱼS⁻ⱼ₊₂ + ½ S⁻ⱼS⁺ⱼ₊₂) - h Σⱼ Sⱼᶻ -function ITensors.MPO(::Model"xx_extended", s; J=1.0, J₂=1.0, h=0.0) - N = length(s) - os = OpSum() - if h != 0 - for j in 1:N - os += -h, "Sz", j - end - end - if J != 0 - for j in 1:(N - 1) - os += 0.5 * J, "S+", j, "S-", j + 1 - os += 0.5 * J, "S-", j, "S+", j + 1 - end - end - if J₂ != 0 - for j in 1:(N - 2) - os += 0.5 * J₂, "S+", j, "Sz", j + 1, "S-", j + 2 - os += 0.5 * J₂, "S-", j, "Sz", j + 1, "S+", j + 2 - end - end - return splitblocks(linkinds, MPO(os, s)) -end - -nrange(::Model"xx_extended") = 3 - function ITensorInfiniteMPS.reference( ::Model"xx_extended", ::Observable"energy"; N=1000, J=1.0, J₂=1.0, h=0.0, filling=0.5 ) diff --git a/src/orthogonalize.jl b/src/orthogonalize.jl index 84660fd..5d1ba22 100644 --- a/src/orthogonalize.jl +++ b/src/orthogonalize.jl @@ -1,39 +1,18 @@ - # TODO: call as `orthogonalize(ψ, -∞)` # TODO: could use commontags(ψ) as a default for left_tags function right_orthogonalize( ψ::InfiniteMPS; left_tags=ts"Left", right_tags=ts"Right", tol::Real=1e-12 ) - # TODO: replace dag(ψ) with ψ'ᴴ? - ψᴴ = prime(linkinds, dag(ψ)) - - N = nsites(ψ) - - # The unit cell range - # Turn into function `eachcellindex(ψ::InfiniteMPS, cell::Integer = 1)` - cell₁ = 1:N - # A transfer matrix made from the 1st unit cell of - # the infinite MPS - # TODO: make a `Cell` Integer type and call as `ψ[Cell(1)]` - # TODO: make a TransferMatrix wrapper that automatically - # primes and daggers, so it can be called with: - # T = TransferMatrix(ψ[Cell(1)]) - ψ₁ = ψ[cell₁] - ψ₁ᴴ = ψᴴ[cell₁] - T₀₁ = ITensorMap( - ψ₁, - ψ₁ᴴ; - input_inds=unioninds(commoninds(ψ[N], ψ[N + 1]), commoninds(ψᴴ[N], ψᴴ[N + 1])), - output_inds=unioninds(commoninds(ψ[1], ψ[0]), commoninds(ψᴴ[1], ψᴴ[0])), - ) + # A transfer matrix made from the 1st unit cell of the infinite MPS + T = TransferMatrix(ψ) # TODO: make an optional initial state - v₁ᴿᴺ = randomITensor(dag(input_inds(T₀₁))) + v₁ᴿᴺ = randomITensor(dag(input_inds(T))) - # Start by getting the right eivenvector/eigenvalue of T₀₁ + # Start by getting the right eivenvector/eigenvalue of T # TODO: make a function `right_environments(::InfiniteMPS)` that computes # all of the right environments using `eigsolve` and shifting unit cells - λ⃗₁ᴿᴺ, v⃗₁ᴿᴺ, eigsolve_info = eigsolve(T₀₁, v₁ᴿᴺ, 1, :LM; tol=tol) + λ⃗₁ᴿᴺ, v⃗₁ᴿᴺ, eigsolve_info = eigsolve(T, v₁ᴿᴺ, 1, :LM; tol=tol) λ₁ᴿᴺ, v₁ᴿᴺ = λ⃗₁ᴿᴺ[1], v⃗₁ᴿᴺ[1] if imag(λ₁ᴿᴺ) / norm(λ₁ᴿᴺ) > 1e-15 @@ -51,7 +30,7 @@ function right_orthogonalize( @show norm(v₁ᴿᴺ - swapinds(dag(v₁ᴿᴺ), reverse(Pair(inds(v₁ᴿᴺ)...)))) error("v₁ᴿᴺ not hermitian") end - if norm(imag(v₁ᴿᴺ)) > 1e-15 + if norm(imag(v₁ᴿᴺ)) / norm(v₁ᴿᴺ) > 1e-13 println( "Norm of the imaginary part $(norm(imag(v₁ᴿᴺ))) is larger than the tolerance value 1e-15. Keeping as complex.", ) @@ -66,16 +45,6 @@ function right_orthogonalize( C₁ᴿᴺ = replacetags(C₁ᴿᴺ, left_tags => right_tags; plev=1) C₁ᴿᴺ = noprime(C₁ᴿᴺ, right_tags) - ## Flip the sign: - ## @show inds(ψ[N]) - ## @show inds(C₁ᴿᴺ) - ## @show C₁ᴿᴺ - ## @show rᴺ = dag(uniqueind(C₁ᴿᴺ, ψ[N])) - ## rᴺ⁻ = Index(-space(rᴺ); tags=tags(rᴺ), plev=plev(rᴺ), dir=dir(rᴺ)) - ## C₁ᴿᴺ = C₁ᴿᴺ * δ(rᴺ, rᴺ⁻) - ## @show C₁ᴿᴺ - ## @show inds(C₁ᴿᴺ) - # Normalize the center matrix normalize!(C₁ᴿᴺ) @@ -108,7 +77,7 @@ function right_orthogonalize_polar( λⁿ = norm(Cᴿ[n - 1]) Cᴿ[n - 1] /= λⁿ λ *= λⁿ - if !isapprox(ψ[n] * Cᴿ[n], λⁿ * Cᴿ[n - 1] * ψᴿ[n]; atol=1e-14) + if !isapprox(ψ[n] * Cᴿ[n], λⁿ * Cᴿ[n - 1] * ψᴿ[n]; rtol=1e-10) @show norm(ψ[n] * Cᴿ[n] - λⁿ * Cᴿ[n - 1] * ψᴿ[n]) error("ψ[n] * Cᴿ[n] ≠ λⁿ * Cᴿ[n-1] * ψᴿ[n]") end diff --git a/src/subspace_expansion.jl b/src/subspace_expansion.jl index a9786d1..eb33c09 100644 --- a/src/subspace_expansion.jl +++ b/src/subspace_expansion.jl @@ -213,10 +213,14 @@ function subspace_expansion( NR *= dag(V) ALⁿ¹, newl = ITensors.directsum( - ψ.AL[n1], dag(NL), uniqueinds(ψ.AL[n1], NL), uniqueinds(NL, ψ.AL[n1]); tags=("Left",) + ψ.AL[n1] => uniqueinds(ψ.AL[n1], NL), + dag(NL) => uniqueinds(NL, ψ.AL[n1]); + tags=("Left",), ) ARⁿ², newr = ITensors.directsum( - ψ.AR[n2], dag(NR), uniqueinds(ψ.AR[n2], NR), uniqueinds(NR, ψ.AR[n2]); tags=("Right",) + ψ.AR[n2] => uniqueinds(ψ.AR[n2], NR), + dag(NR) => uniqueinds(NR, ψ.AR[n2]); + tags=("Right",), ) C = ITensor(dag(newl)..., dag(newr)...) diff --git a/test/Project.toml b/test/Project.toml index 340ddf6..b9fb740 100644 --- a/test/Project.toml +++ b/test/Project.toml @@ -1,4 +1,7 @@ [deps] ITensors = "9136182c-28ba-11e9-034c-db9fb085ebd5" +KrylovKit = "0b1a1467-8014-51b9-945f-bf0ae24f4b77" +LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e" Random = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c" +Suppressor = "fd094767-a336-5f1f-9728-57cf17d0bbfb" Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40" diff --git a/test/test_examples.jl b/test/test_examples.jl new file mode 100644 index 0000000..a284669 --- /dev/null +++ b/test/test_examples.jl @@ -0,0 +1,12 @@ +using Test +using Suppressor + +examples_dir = joinpath(@__DIR__, "..", "examples", "vumps") +@testset "examples" begin + @testset "$file" for file in readdir(examples_dir; join=true) + if endswith(file, ".jl") + println("Running $file") + @suppress include(file) + end + end +end diff --git a/test/test_readwrite.jl b/test/test_readwrite.jl index 1119366..18115c9 100644 --- a/test/test_readwrite.jl +++ b/test/test_readwrite.jl @@ -4,15 +4,10 @@ using ITensorInfiniteMPS.HDF5 using Test @testset "HDF5 Read and Write" begin - function space_shifted(::Model"ising", q̃sz) - return [QN("SzParity", 1 - q̃sz, 2) => 1, QN("SzParity", 0 - q̃sz, 2) => 1] - end - N = 2 - model = Model"ising"() - space_ = fill(space_shifted(model, 0), N) - s = infsiteinds("S=1/2", N; space=space_) + model = Model("ising") initstate(n) = "↑" + s = infsiteinds("S=1/2", N; initstate) ψ = InfMPS(s, initstate) @testset "InfiniteCanonicalMPS" begin diff --git a/test/test_vumps.jl b/test/test_vumps.jl index 6718438..556540c 100644 --- a/test/test_vumps.jl +++ b/test/test_vumps.jl @@ -3,20 +3,23 @@ using ITensorInfiniteMPS using Test using Random +function expect_two_site(ψ1::ITensor, ψ2::ITensor, h::ITensor) + ϕ = ψ1 * ψ2 + return inner(ϕ, apply(h, ϕ)) +end + +expect_two_site(ψ1::ITensor, ψ2::ITensor, h::MPO) = expect_two_site(ψ1, ψ2, contract(h)) + +function expect_one_site(ψ::ITensor, o::String) + return inner(ψ, apply(op(o, filterinds(ψ, "Site")...), ψ)) +end + @testset "vumps" begin Random.seed!(1234) - model = Model"ising"() + model = Model("ising") model_kwargs = (J=1.0, h=1.1) - function space_shifted(::Model"ising", q̃sz; conserve_qns=true) - if conserve_qns - return [QN("SzParity", 1 - q̃sz, 2) => 1, QN("SzParity", 0 - q̃sz, 2) => 1] - else - return [QN() => 2] - end - end - # VUMPS arguments cutoff = 1e-8 maxdim = 20 @@ -31,10 +34,10 @@ using Random sfinite = siteinds("S=1/2", Nfinite; conserve_szparity=true) Hfinite = MPO(model, sfinite; model_kwargs...) ψfinite = randomMPS(sfinite, initstate) - sweeps = Sweeps(20) - setmaxdim!(sweeps, 10) - setcutoff!(sweeps, 1E-10) - energy_finite_total, ψfinite = dmrg(Hfinite, ψfinite, sweeps; outputlevel=0) + nsweeps = 20 + energy_finite_total, ψfinite = dmrg( + Hfinite, ψfinite; nsweeps, maxdims=10, cutoff=1e-10, outputlevel=0 + ) @testset "VUMPS/TDVP with: multisite_update_alg = $multisite_update_alg, conserve_qns = $conserve_qns, nsites = $nsites, time_step = $time_step, localham_type = $localham_type" for multisite_update_alg in [ @@ -45,20 +48,19 @@ using Random time_step in [-Inf, -0.5], localham_type in [ITensor, MPO] - # ITensor VUMPS currently broken for unit cells > 2 if (localham_type == ITensor) && (nsites > 2) + # ITensor VUMPS currently broken for unit cells > 2 continue end - @show localham_type + if nsites > 1 && isodd(nsites) && conserve_qns + # Odd site greater than 1 not commensurate with conserving parity + continue + end - #if conserve_qns - # Flux density is inconsistent for odd unit cells - # isodd(nsites) && continue - #end + @show localham_type - space_ = fill(space_shifted(model, 1; conserve_qns=conserve_qns), nsites) - s = infsiteinds("S=1/2", nsites; space=space_) + s = infsiteinds("S=1/2", nsites; initstate, conserve_szparity=conserve_qns) ψ = InfMPS(s, initstate) # Form the Hamiltonian @@ -86,38 +88,27 @@ using Random ## @test contract(ψ.AL[1:nsites]..., ψ.C[nsites]) ≈ contract(ψ.C[0], ψ.AR[1:nsites]...) rtol = ## 1e-6 - function energy(ψ1, ψ2, h::ITensor) - ϕ = ψ1 * ψ2 - return (noprime(ϕ * h) * dag(ϕ))[] - end - - energy(ψ1, ψ2, h::MPO) = energy(ψ1, ψ2, prod(h)) - - function expect(ψ, o) - return (noprime(ψ * op(o, filterinds(ψ, "Site")...)) * dag(ψ))[] - end - nfinite = Nfinite ÷ 2 hnfinite1 = ITensor(model, sfinite[nfinite], sfinite[nfinite + 1]; model_kwargs...) hnfinite2 = ITensor(model, sfinite[nfinite + 1], sfinite[nfinite + 2]; model_kwargs...) orthogonalize!(ψfinite, nfinite) - energy1_finite = energy(ψfinite[nfinite], ψfinite[nfinite + 1], hnfinite1) + energy1_finite = expect_two_site(ψfinite[nfinite], ψfinite[nfinite + 1], hnfinite1) orthogonalize!(ψfinite, nfinite + 1) - energy2_finite = energy(ψfinite[nfinite + 1], ψfinite[nfinite + 2], hnfinite2) + energy2_finite = expect_two_site(ψfinite[nfinite + 1], ψfinite[nfinite + 2], hnfinite2) - energy1_infinite = energy(ψ.AL[1], ψ.AL[2] * ψ.C[2], H[(1, 2)]) - energy2_infinite = energy(ψ.AL[2], ψ.AL[3] * ψ.C[3], H[(2, 3)]) + energy1_infinite = expect_two_site(ψ.AL[1], ψ.AL[2] * ψ.C[2], H[(1, 2)]) + energy2_infinite = expect_two_site(ψ.AL[2], ψ.AL[3] * ψ.C[3], H[(2, 3)]) orthogonalize!(ψfinite, nfinite) - Sz1_finite = expect(ψfinite[nfinite], "Sz") + Sz1_finite = expect_one_site(ψfinite[nfinite], "Sz") orthogonalize!(ψfinite, nfinite + 1) - Sz2_finite = expect(ψfinite[nfinite + 1], "Sz") + Sz2_finite = expect_one_site(ψfinite[nfinite + 1], "Sz") - Sz1_infinite = expect(ψ.AL[1] * ψ.C[1], "Sz") - Sz2_infinite = expect(ψ.AL[2] * ψ.C[2], "Sz") + Sz1_infinite = expect_one_site(ψ.AL[1] * ψ.C[1], "Sz") + Sz2_infinite = expect_one_site(ψ.AL[2] * ψ.C[2], "Sz") @test energy1_finite ≈ energy1_infinite rtol = 1e-4 @test energy2_finite ≈ energy2_infinite rtol = 1e-4 @@ -129,18 +120,10 @@ end @testset "vumps_ising_translator" begin Random.seed!(1234) - model = Model"ising"() + model = Model("ising") model_kwargs = (J=1.0, h=1.1) - function space_shifted(::Model"ising", q̃sz; conserve_qns=true) - if conserve_qns - return [QN("SzParity", 1 - q̃sz, 2) => 1, QN("SzParity", 0 - q̃sz, 2) => 1] - else - return [QN() => 2] - end - end - - #Not a correct and valid Hamiltonian #nned to think about a good test (or just import Laughlin 13) + #Not a correct and valid Hamiltonian need to think about a good test (or just import Laughlin 13) temp_translatecell(i::Index, n::Integer) = ITensorInfiniteMPS.translatecelltags(i, n) # VUMPS arguments @@ -157,19 +140,18 @@ end sfinite = siteinds("S=1/2", Nfinite; conserve_szparity=true) Hfinite = MPO(model, sfinite; model_kwargs...) ψfinite = randomMPS(sfinite, initstate) - sweeps = Sweeps(20) - setmaxdim!(sweeps, 10) - setcutoff!(sweeps, 1E-10) - energy_finite_total, ψfinite = dmrg(Hfinite, ψfinite, sweeps; outputlevel=0) + nsweeps = 20 + energy_finite_total, ψfinite = dmrg( + Hfinite, ψfinite; nsweeps, maxdims=10, cutoff=1e-10, outputlevel=0 + ) multisite_update_alg = "sequential" conserve_qns = true time_step = -Inf for nsite in 1:3 - space_ = fill(space_shifted(model, 1; conserve_qns=conserve_qns), nsite) - s_bis = infsiteinds("S=1/2", nsite; space=space_) - s = infsiteinds("S=1/2", nsite; space=space_, translator=temp_translatecell) + s_bis = infsiteinds("S=1/2", nsite; initstate) + s = infsiteinds("S=1/2", nsite; initstate, translator=temp_translatecell) ψ = InfMPS(s, initstate) # Form the Hamiltonian @@ -197,36 +179,27 @@ end ## @test contract(ψ.AL[1:nsites]..., ψ.C[nsites]) ≈ contract(ψ.C[0], ψ.AR[1:nsites]...) rtol = ## 1e-6 - function energy(ψ1, ψ2, h) - ϕ = ψ1 * ψ2 - return (noprime(ϕ * h) * dag(ϕ))[] - end - - function expect(ψ, o) - return (noprime(ψ * op(o, filterinds(ψ, "Site")...)) * dag(ψ))[] - end - nfinite = Nfinite ÷ 2 hnfinite1 = ITensor(model, sfinite[nfinite], sfinite[nfinite + 1]; model_kwargs...) hnfinite2 = ITensor(model, sfinite[nfinite + 1], sfinite[nfinite + 2]; model_kwargs...) orthogonalize!(ψfinite, nfinite) - energy1_finite = energy(ψfinite[nfinite], ψfinite[nfinite + 1], hnfinite1) + energy1_finite = expect_two_site(ψfinite[nfinite], ψfinite[nfinite + 1], hnfinite1) orthogonalize!(ψfinite, nfinite + 1) - energy2_finite = energy(ψfinite[nfinite + 1], ψfinite[nfinite + 2], hnfinite2) + energy2_finite = expect_two_site(ψfinite[nfinite + 1], ψfinite[nfinite + 2], hnfinite2) - energy1_infinite = energy(ψ.AL[1], ψ.AL[2] * ψ.C[2], prod(H[(1, 2)])) - energy2_infinite = energy(ψ.AL[2], ψ.AL[3] * ψ.C[3], prod(H[(2, 3)])) + energy1_infinite = expect_two_site(ψ.AL[1], ψ.AL[2] * ψ.C[2], prod(H[(1, 2)])) + energy2_infinite = expect_two_site(ψ.AL[2], ψ.AL[3] * ψ.C[3], prod(H[(2, 3)])) orthogonalize!(ψfinite, nfinite) - Sz1_finite = expect(ψfinite[nfinite], "Sz") + Sz1_finite = expect_one_site(ψfinite[nfinite], "Sz") orthogonalize!(ψfinite, nfinite + 1) - Sz2_finite = expect(ψfinite[nfinite + 1], "Sz") + Sz2_finite = expect_one_site(ψfinite[nfinite + 1], "Sz") - Sz1_infinite = expect(ψ.AL[1] * ψ.C[1], "Sz") - Sz2_infinite = expect(ψ.AL[2] * ψ.C[2], "Sz") + Sz1_infinite = expect_one_site(ψ.AL[1] * ψ.C[1], "Sz") + Sz2_infinite = expect_one_site(ψ.AL[2] * ψ.C[2], "Sz") @test energy1_finite ≈ energy1_infinite rtol = 1e-4 @test energy2_finite ≈ energy2_infinite rtol = 1e-4 diff --git a/test/test_vumps_extendedising.jl b/test/test_vumps_extendedising.jl index 9135bc1..7842cf2 100644 --- a/test/test_vumps_extendedising.jl +++ b/test/test_vumps_extendedising.jl @@ -3,43 +3,39 @@ using ITensorInfiniteMPS using Test using Random +function expect_three_site(ψ::MPS, h::ITensor, n::Int) + ϕ = ψ[n] * ψ[n + 1] * ψ[n + 2] + return inner(ϕ, apply(h, ϕ)) +end + @testset "vumps_extended_ising" begin Random.seed!(1234) - model = Model"ising_extended"() + model = Model("ising_extended") model_kwargs = (J=1.0, h=1.1, J₂=0.2) initstate(n) = "↑" - function space_shifted(::Model"ising_extended", q̃sz; conserve_qns=true) - if conserve_qns - return [QN("SzParity", 1 - q̃sz, 2) => 1, QN("SzParity", 0 - q̃sz, 2) => 1] - else - return [QN() => 2] - end - end - # Compare to DMRG Nfinite = 100 sfinite = siteinds("S=1/2", Nfinite; conserve_szparity=true) Hfinite = MPO(model, sfinite; model_kwargs...) ψfinite = randomMPS(sfinite, initstate) - sweeps = Sweeps(20) - setmaxdim!(sweeps, 30) - setcutoff!(sweeps, 1E-10) - energy_finite_total, ψfinite = dmrg(Hfinite, ψfinite, sweeps; outputlevel=0) + energy_finite_total, ψfinite = dmrg( + Hfinite, ψfinite; outputlevel=0, nsweeps=30, maxdim=30, cutoff=1e-10 + ) Szs_finite = expect(ψfinite, "Sz") - function energy(ψ, h, n) - ϕ = ψ[n] * ψ[n + 1] * ψ[n + 2] - return (noprime(ϕ * h) * dag(ϕ))[] - end - nfinite = Nfinite ÷ 2 hnfinite = ITensor( - model, sfinite[nfinite], sfinite[nfinite + 1], sfinite[nfinite + 2]; model_kwargs... + model, + nfinite, + sfinite[nfinite], + sfinite[nfinite + 1], + sfinite[nfinite + 2]; + model_kwargs..., ) orthogonalize!(ψfinite, nfinite) - energy_finite = energy(ψfinite, hnfinite, nfinite) + energy_finite = expect_three_site(ψfinite, hnfinite, nfinite) cutoff = 1e-8 maxdim = 100 @@ -60,8 +56,7 @@ using Random ) subspace_expansion_kwargs = (cutoff=cutoff, maxdim=maxdim) - space_ = fill(space_shifted(model, 1; conserve_qns=conserve_qns), nsites) - s = infsiteinds("S=1/2", nsites; space=space_) + s = infsiteinds("S=1/2", nsites; initstate) H = InfiniteSum{MPO}(model, s; model_kwargs...) ψ = InfMPS(s, initstate) @@ -86,18 +81,16 @@ using Random end end -# XXX: orthogonalize is broken right now -## @testset "ITensorInfiniteMPS.jl" begin -## @testset "Mixed canonical gauge" begin -## N = 10 -## s = siteinds("S=1/2", N; conserve_szparity=true) -## χ = 6 -## @test iseven(χ) -## space = (("SzParity", 1, 2) => χ ÷ 2) ⊕ (("SzParity", 0, 2) => χ ÷ 2) -## ψ = InfiniteMPS(ComplexF64, s; space=space) -## randn!.(ψ) -## -## ψ = orthogonalize(ψ, :) -## @test prod(ψ.AL[1:N]) * ψ.C[N] ≈ ψ.C[0] * prod(ψ.AR[1:N]) -## end -## end +@testset "Mixed canonical gauge" begin + N = 4 + s = siteinds("S=1/2", N; conserve_szparity=true) + χ = 6 + @test iseven(χ) + space = (("SzParity", 1, 2) => χ ÷ 2) ⊕ (("SzParity", 0, 2) => χ ÷ 2) + ψ = InfiniteMPS(ComplexF64, s; space=space) + for n in 1:N + ψ[n] = randomITensor(inds(ψ[n])) + end + ψ = orthogonalize(ψ, :) + @test contract(ψ.AL[1:N]) * ψ.C[N] ≈ ψ.C[0] * contract(ψ.AR[1:N]) +end diff --git a/test/test_vumpsmpo.jl b/test/test_vumpsmpo.jl index 51b3dec..98b1185 100644 --- a/test/test_vumpsmpo.jl +++ b/test/test_vumpsmpo.jl @@ -3,21 +3,18 @@ using ITensorInfiniteMPS using Test using Random +function expect_three_site(ψ::MPS, h::ITensor, n::Int) + ϕ = ψ[n] * ψ[n + 1] * ψ[n + 2] + return inner(ϕ, (apply(h, ϕ))) +end + #Ref time is 21.6s with negligible compilation time @time @testset "vumpsmpo_ising" begin Random.seed!(1234) - model = Model"ising"() + model = Model("ising") model_kwargs = (J=1.0, h=1.2) - function space_shifted(::Model"ising", q̃sz; conserve_qns=true) - if conserve_qns - return [QN("SzParity", 1 - q̃sz, 2) => 1, QN("SzParity", 0 - q̃sz, 2) => 1] - else - return [QN() => 2] - end - end - # VUMPS arguments cutoff = 1e-8 maxdim = 20 @@ -32,21 +29,17 @@ using Random sfinite = siteinds("S=1/2", Nfinite; conserve_szparity=true) Hfinite = MPO(model, sfinite; model_kwargs...) ψfinite = randomMPS(sfinite, initstate) + sweeps = Sweeps(20) setmaxdim!(sweeps, 10) setcutoff!(sweeps, 1E-10) energy_finite_total, ψfinite = dmrg(Hfinite, ψfinite, sweeps; outputlevel=0) Szs_finite = expect(ψfinite, "Sz") - function energy(ψ, h, n) - ϕ = ψ[n] * ψ[n + 1] * ψ[n + 2] - return (noprime(ϕ * h) * dag(ϕ))[] - end - nfinite = Nfinite ÷ 2 hnfinite = ITensor(model, sfinite[nfinite], sfinite[nfinite + 1]; model_kwargs...) orthogonalize!(ψfinite, nfinite) - energy_finite = energy(ψfinite, hnfinite, nfinite) + energy_finite = expect_three_site(ψfinite, hnfinite, nfinite) @testset "VUMPS/TDVP with: multisite_update_alg = $multisite_update_alg, conserve_qns = $conserve_qns, nsites = $nsites" for multisite_update_alg in [ @@ -65,8 +58,7 @@ using Random ) subspace_expansion_kwargs = (cutoff=cutoff, maxdim=maxdim) - space_ = fill(space_shifted(model, 1; conserve_qns=conserve_qns), nsites) - s = infsiteinds("S=1/2", nsites; space=space_) + s = infsiteinds("S=1/2", nsites; initstate, conserve_szparity=conserve_qns) ψ = InfMPS(s, initstate) Hmpo = InfiniteMPOMatrix(model, s; model_kwargs...) @@ -99,14 +91,6 @@ end model = Model"ising_extended"() model_kwargs = (J=1.0, h=1.1, J₂=0.2) - function space_shifted(::Model"ising_extended", q̃sz; conserve_qns=true) - if conserve_qns - return [QN("SzParity", 1 - q̃sz, 2) => 1, QN("SzParity", 0 - q̃sz, 2) => 1] - else - return [QN() => 2] - end - end - # VUMPS arguments cutoff = 1e-8 maxdim = 20 @@ -121,26 +105,21 @@ end sfinite = siteinds("S=1/2", Nfinite; conserve_szparity=true) Hfinite = MPO(model, sfinite; model_kwargs...) ψfinite = randomMPS(sfinite, initstate) - sweeps = Sweeps(20) - setmaxdim!(sweeps, 10) - setcutoff!(sweeps, 1E-10) - energy_finite_total, ψfinite = dmrg(Hfinite, ψfinite, sweeps; outputlevel=0) + nsweeps = 20 + energy_finite_total, ψfinite = dmrg( + Hfinite, ψfinite; nsweeps, maxdims=10, cutoff=1e-10, outputlevel=0 + ) Szs_finite = expect(ψfinite, "Sz") - function energy(ψ, h, n) - ϕ = ψ[n] * ψ[n + 1] * ψ[n + 2] - return (noprime(ϕ * h) * dag(ϕ))[] - end - nfinite = Nfinite ÷ 2 hnfinite = ITensor( model, sfinite[nfinite], sfinite[nfinite + 1], sfinite[nfinite + 2]; model_kwargs... ) orthogonalize!(ψfinite, nfinite) - energy_finite = energy(ψfinite, hnfinite, nfinite) + energy_finite = expect_three_site(ψfinite, hnfinite, nfinite) for multisite_update_alg in ["sequential"], - conserve_qns in [true], + conserve_qns in [true, false], nsites in [1, 2], time_step in [-Inf] @@ -153,8 +132,7 @@ end ) subspace_expansion_kwargs = (cutoff=cutoff, maxdim=maxdim) - space_ = fill(space_shifted(model, 1; conserve_qns=conserve_qns), nsites) - s = infsiteinds("S=1/2", nsites; space=space_) + s = infsiteinds("S=1/2", nsites; conserve_szparity=conserve_qns, initstate) ψ = InfMPS(s, initstate) Hmpo = InfiniteMPOMatrix(model, s; model_kwargs...) @@ -181,17 +159,9 @@ end @testset "vumpsmpo_extendedising_translator" begin Random.seed!(1234) - model = Model"ising_extended"() + model = Model("ising_extended") model_kwargs = (J=1.0, h=1.1, J₂=0.2) - function space_shifted(::Model"ising_extended", q̃sz; conserve_qns=true) - if conserve_qns - return [QN("SzParity", 1 - q̃sz, 2) => 1, QN("SzParity", 0 - q̃sz, 2) => 1] - else - return [QN() => 2] - end - end - # VUMPS arguments cutoff = 1e-8 maxdim = 20 @@ -212,25 +182,25 @@ end energy_finite_total, ψfinite = dmrg(Hfinite, ψfinite, sweeps; outputlevel=0) Szs_finite = expect(ψfinite, "Sz") - function energy(ψ, h, n) - ϕ = ψ[n] * ψ[n + 1] * ψ[n + 2] - return (noprime(ϕ * h) * dag(ϕ))[] - end - nfinite = Nfinite ÷ 2 hnfinite = ITensor( model, sfinite[nfinite], sfinite[nfinite + 1], sfinite[nfinite + 2]; model_kwargs... ) orthogonalize!(ψfinite, nfinite) - energy_finite = energy(ψfinite, hnfinite, nfinite) + energy_finite = expect_three_site(ψfinite, hnfinite, nfinite) temp_translatecell(i::Index, n::Integer) = ITensorInfiniteMPS.translatecelltags(i, n) for multisite_update_alg in ["sequential"], - conserve_qns in [true], + conserve_qns in [true, false], nsite in [1, 2, 3], time_step in [-Inf] + if nsite > 1 && isodd(nsite) && conserve_qns + # Parity conservation not commensurate with odd number of sites. + continue + end + vumps_kwargs = ( multisite_update_alg=multisite_update_alg, tol=tol, @@ -240,9 +210,14 @@ end ) subspace_expansion_kwargs = (cutoff=cutoff, maxdim=maxdim) - space_ = fill(space_shifted(model, 1; conserve_qns=conserve_qns), nsite) - s_bis = infsiteinds("S=1/2", nsite; space=space_) - s = infsiteinds("S=1/2", nsite; space=space_, translator=temp_translatecell) + s_bis = infsiteinds("S=1/2", nsite; conserve_szparity=conserve_qns, initstate) + s = infsiteinds( + "S=1/2", + nsite; + conserve_szparity=conserve_qns, + initstate, + translator=temp_translatecell, + ) ψ = InfMPS(s, initstate) Hmpo = InfiniteMPOMatrix(model, s; model_kwargs...) diff --git a/test/test_vumpsmpo_fqhe.jl b/test/test_vumpsmpo_fqhe.jl index 49b8b9a..76a9d22 100644 --- a/test/test_vumpsmpo_fqhe.jl +++ b/test/test_vumpsmpo_fqhe.jl @@ -3,22 +3,26 @@ using ITensorInfiniteMPS using Test using Random +function ITensors.space(::SiteType"FermionK", pos::Int; p=1, q=1, conserve_momentum=true) + if !conserve_momentum + return [QN("Nf", -p) => 1, QN("Nf", q - p) => 1] + else + return [ + QN(("Nf", -p), ("NfMom", -p * pos)) => 1, + QN(("Nf", q - p), ("NfMom", (q - p) * pos)) => 1, + ] + end +end + +# Forward all op definitions to Fermion +function ITensors.op!(Op::ITensor, opname::OpName, ::SiteType"FermionK", s::Index...) + return ITensors.op!(Op, opname, SiteType("Fermion"), s...) +end + #Currently, VUMPS cannot give the right result as the subspace expansion is too small #This is meant to test the generalized translation rules @time @testset "vumpsmpo_fqhe" begin Random.seed!(1234) - - function fermions_space_shift(pos; p=1, q=1, conserve_momentum=true, momentum_shift=1) - if !conserve_momentum - return [QN("Nf", -p) => 1, QN("Nf", q - p) => 1] - else - return [ - QN(("Nf", -p), ("NfMom", -p * pos + momentum_shift)) => 1, - QN(("Nf", q - p), ("NfMom", (q - p) * pos + momentum_shift)) => 1, - ] - end - end - function initstate(n) if mod(n, 3) == 1 return 2 @@ -35,7 +39,7 @@ using Random outer_iters = 1 #### Test 1/3 - model = Model"fqhe_2b_pot"() + model = Model("fqhe_2b_pot") model_kwargs = (Ly=3.0, Vs=[1.0], prec=1e-6) #parent Hamiltonian of Laughlin 13 p = 1 q = 3 @@ -52,14 +56,8 @@ using Random nsites in [6], time_step in [-Inf] - vumps_kwargs = ( - multisite_update_alg=multisite_update_alg, - tol=tol, - maxiter=maxiter, - outputlevel=0, - time_step=time_step, - ) - subspace_expansion_kwargs = (cutoff=cutoff, maxdim=maxdim) + vumps_kwargs = (; multisite_update_alg, tol, maxiter, outputlevel=0, time_step) + subspace_expansion_kwargs = (; cutoff, maxdim) function fermion_momentum_translator(i::Index, n::Integer; N=nsites) ts = tags(i) @@ -75,13 +73,14 @@ using Random return new_i end - fermionic_space = [ - fermions_space_shift( - x; p=p, q=q, conserve_momentum=conserve_momentum, momentum_shift=momentum_shift - ) for x in 1:nsites - ] s = infsiteinds( - "Fermion", nsites; space=fermionic_space, translator=fermion_momentum_translator + "FermionK", + nsites; + initstate, + translator=fermion_momentum_translator, + p, + q, + conserve_momentum, ) ψ = InfMPS(s, initstate)