From fbf3e22ff62bc588667e8930dc79c8afc7c23ba3 Mon Sep 17 00:00:00 2001 From: Matthew Fishman Date: Wed, 15 Jun 2022 13:32:41 -0400 Subject: [PATCH 01/32] Start simplifying the necessary OpSum definitions for Models --- examples/vumps/vumps_ising_extended.jl | 32 +++++---- src/celledvectors.jl | 8 ++- src/models/ising_extended.jl | 65 +++++++++--------- src/models/models.jl | 95 +++++++++++++++++++++----- 4 files changed, 133 insertions(+), 67 deletions(-) diff --git a/examples/vumps/vumps_ising_extended.jl b/examples/vumps/vumps_ising_extended.jl index 6f03e48..1cf766a 100644 --- a/examples/vumps/vumps_ising_extended.jl +++ b/examples/vumps/vumps_ising_extended.jl @@ -7,14 +7,14 @@ using ITensorInfiniteMPS 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 -outer_iters = 4 # Number of times to increase the bond dimension +outer_iters = 10 # Number of times to increase the bond dimension ############################################################################## # CODE BELOW HERE DOES NOT NEED TO BE MODIFIED # -N = 3# Number of sites in the unit cell +N = 3 # Number of sites in the unit cell J = -1.0 J₂ = -0.2 h = 1.0; @@ -39,9 +39,9 @@ subspace_expansion_kwargs = (cutoff=cutoff, maxdim=maxdim) for j in 1:outer_iters println("\nIncrease bond dimension") - ψ_1 = subspace_expansion(ψ_0, H; subspace_expansion_kwargs...) + ψ_1 = @time subspace_expansion(ψ_0, H; subspace_expansion_kwargs...) println("Run VUMPS with new bond dimension") - global ψ_0 = vumps(H, ψ_1; vumps_kwargs...) + global ψ_0 = @time vumps(H, ψ_1; vumps_kwargs...) end Sz = [expect(ψ_0, "Sz", n) for n in 1:N] @@ -52,18 +52,20 @@ sfinite = siteinds("S=1/2", Nfinite; conserve_szparity=true) Hfinite = MPO(model, sfinite; J=J, J₂=J₂, h=h) ψ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) +# sweeps = Sweeps(10) +# setmaxdim!(sweeps, maxdim) +# setcutoff!(sweeps, cutoff) +dmrg_kwargs = (nsweeps=10, maxdim, mindim=64, cutoff) +energy_finite_total, ψfinite = dmrg(Hfinite, ψfinite; dmrg_kwargs...) nfinite = Nfinite ÷ 2 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₂), -) +println("\nEnergy") +@show energy_infinite +@show energy_finite_total / Nfinite +@show reference(model, Observable("energy"); J=J, h=h, J₂=J₂) -@show (Sz, Sz_finite) +println("\nSz") +@show Sz +@show Sz_finite diff --git a/src/celledvectors.jl b/src/celledvectors.jl index ab2d6c2..c026510 100644 --- a/src/celledvectors.jl +++ b/src/celledvectors.jl @@ -99,11 +99,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/models/ising_extended.jl b/src/models/ising_extended.jl index 6a16475..d8b2be0 100644 --- a/src/models/ising_extended.jl +++ b/src/models/ising_extended.jl @@ -1,48 +1,45 @@ -nrange(::Model"ising_extended") = 3 +# For infinite Hamiltonian # 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 +function opsum_infinite(::Model"ising_extended", n; J=1.0, h=1.0, J₂=0.0) + opsum = OpSum() + for j in 1:n + opsum += -J, "X", j, "X", j + 1 + opsum += -J₂, "X", j, "Z", j + 1, "X", j + 2 + opsum -= h, "Z", j end - return MPO(a, s) + return opsum end +# For finite Hamiltonian with open boundary conditions # 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 opsum_finite(::Model"ising_extended", n; J=1.0, h=1.0, J₂=0.0) + # TODO: define using opsum_infinite, filtering sites + # that extend over the boundary. 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 + for j in 1:n + if j ≤ (n - 1) + opsum += -J, "X", j, "X", j + 1 + end + if j ≤ (n - 2) + opsum += -J₂, "X", j, "Z", j + 1, "X", j + 2 + end + opsum -= h, "Z", j end - opsum += -h / 3, "Z", n1 - opsum += -h / 3, "Z", n2 - opsum += -h / 3, "Z", n2 + 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), -π, π) end + +# +# Deprecated +# + +# nrange(::Model"ising_extended") = 3 + +# function ITensors.ITensor(::Model"ising_extended", s1::Index, s2::Index, s3::Index; 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 diff --git a/src/models/models.jl b/src/models/models.jl index 87a2bcf..a1020db 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,52 @@ function InfiniteSum{T}(model::Model, s::Vector; kwargs...) where {T} end function InfiniteSum{MPO}(model::Model, s::CelledVector; kwargs...) - N = length(s) - mpos = [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 - return InfiniteSum{ITensor}(itensors, translator(s)) +# Get the first site with nontrivial support of the OpSum +first_site(opsum::OpSum) = minimum(ITensors.sites(opsum)) + +function set_site(o::Op, s::Int) + return Op(ITensors.which_op(o), s; ITensors.params(o)...) 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 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 +function shift_site(o::Op, shift::Int) + return set_site(o, ITensors.site(o) + shift) 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_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 = [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 +101,44 @@ 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 MPO(opsum, s) +end + +# +# Deprecated +# + +#ITensorInfiniteMPS.nrange(model::Model) = 2; #required to keep everything compatible with the current implementation for 2 band models + +# 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 +# end + +# function InfiniteSum{ITensor}(model::Model, s::CelledVector; kwargs...) +# N = length(s) +# 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; kwargs...) +# opsum = OpSum(model, nsites(s); kwargs...) +# +# # @show nrange(ops) +# # @show n:(n + nrange(opsum) - 1) +# # @show opsum +# +# # return MPO(opsum, [s[x] for x in n:(n + nrange(opsum) - 1)]) #modification to allow for more than two sites per term in the Hamiltonians +# return MPO(opsum, [s[x] for x in ITensors.sites(opsum)]) #modification to allow for more than two sites per term in the Hamiltonians +# end + +# function InfiniteSum{MPO}(model::Model, s::CelledVector; kwargs...) +# N = length(s) +# # mpos = [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 +# mpo = MPO(model, s; kwargs...) +# return InfiniteSum{MPO}([mpo], translator(s)) +# end From 712c9fbb221cffc96316e279a893a6b07da24249 Mon Sep 17 00:00:00 2001 From: Matthew Fishman Date: Tue, 21 Jun 2022 09:18:32 -0400 Subject: [PATCH 02/32] Format --- examples/vumps/vumps_ising_extended.jl | 2 +- src/models/models.jl | 5 ++++- 2 files changed, 5 insertions(+), 2 deletions(-) diff --git a/examples/vumps/vumps_ising_extended.jl b/examples/vumps/vumps_ising_extended.jl index 8dd26df..2f9cf36 100644 --- a/examples/vumps/vumps_ising_extended.jl +++ b/examples/vumps/vumps_ising_extended.jl @@ -53,7 +53,7 @@ sfinite = siteinds("S=1/2", Nfinite; conserve_szparity=true) Hfinite = MPO(model, sfinite; J=J, J₂=J₂, h=h) ψfinite = randomMPS(sfinite, initstate; linkdims=10) @show flux(ψfinite) -dmrg_kwargs = (nsweeps=10, maxdim cutoff) +dmrg_kwargs = (nsweeps=10, maxdim, cutoff) energy_finite_total, ψfinite = dmrg(Hfinite, ψfinite; dmrg_kwargs...) energy_finite_total, ψfinite = dmrg(Hfinite, ψfinite, sweeps) diff --git a/src/models/models.jl b/src/models/models.jl index e9fa70e..57128d6 100644 --- a/src/models/models.jl +++ b/src/models/models.jl @@ -65,7 +65,10 @@ function InfiniteSum{MPO}(opsum::OpSum, s::CelledVector) 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] + 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 From 469b8e1ee3cdc9c5535ec9bb240fe5a346f46c1c Mon Sep 17 00:00:00 2001 From: Matthew Fishman Date: Tue, 21 Jun 2022 19:46:03 -0400 Subject: [PATCH 03/32] Fix finite to infinite MPS conversion --- .../finite_mps_to_infinite_mps.jl | 63 ++++++++++ .../experimental/orthogonalize_infinitemps.jl | 19 +++ examples/finite_mps_to_infinite_mps.jl | 112 ------------------ examples/infinitemps.jl | 20 ---- src/infinitemps_approx.jl | 7 +- src/models/models.jl | 6 + src/orthogonalize.jl | 45 ++----- test/test_vumps_extendedising.jl | 33 +++--- 8 files changed, 114 insertions(+), 191 deletions(-) create mode 100644 examples/experimental/finite_mps_to_infinite_mps.jl create mode 100644 examples/experimental/orthogonalize_infinitemps.jl delete mode 100644 examples/finite_mps_to_infinite_mps.jl delete mode 100644 examples/infinitemps.jl diff --git a/examples/experimental/finite_mps_to_infinite_mps.jl b/examples/experimental/finite_mps_to_infinite_mps.jl new file mode 100644 index 0000000..01fe041 --- /dev/null +++ b/examples/experimental/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/experimental/orthogonalize_infinitemps.jl b/examples/experimental/orthogonalize_infinitemps.jl new file mode 100644 index 0000000..273b1cf --- /dev/null +++ b/examples/experimental/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/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/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/models.jl b/src/models/models.jl index 57128d6..77a63ed 100644 --- a/src/models/models.jl +++ b/src/models/models.jl @@ -110,6 +110,12 @@ function ITensors.MPO(model::Model, s::Vector{<:Index}; kwargs...) return splitblocks(linkinds, MPO(opsum, s)) end +# The ITensor of a single term of the model +function ITensors.ITensor(model::Model, s::Index...; kwargs...) + opsum = opsum_infinite(model, 1; kwargs...) + return contract(MPO(opsum, [s...])) +end + # # Deprecated # 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/test/test_vumps_extendedising.jl b/test/test_vumps_extendedising.jl index 9135bc1..607de13 100644 --- a/test/test_vumps_extendedising.jl +++ b/test/test_vumps_extendedising.jl @@ -23,10 +23,7 @@ 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, 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) @@ -86,18 +83,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 From b2b43ba9662b0820cd8fe096c315c64198ec6655 Mon Sep 17 00:00:00 2001 From: Matthew Fishman Date: Wed, 22 Jun 2022 17:09:22 -0400 Subject: [PATCH 04/32] New interface for models in terms of unit_cell_terms, defining the terms on each site of the primite unit cell in a Vector{OpSum} --- src/models/ising_extended.jl | 51 +++++++++++++------------ src/models/models.jl | 64 ++++++++++++++++++++++++++++++-- test/test_vumps_extendedising.jl | 2 +- 3 files changed, 87 insertions(+), 30 deletions(-) diff --git a/src/models/ising_extended.jl b/src/models/ising_extended.jl index d8b2be0..91fc707 100644 --- a/src/models/ising_extended.jl +++ b/src/models/ising_extended.jl @@ -1,31 +1,12 @@ -# For infinite Hamiltonian +# A vector of the terms associated with each site of the +# unit cell. # H = -J Σⱼ XⱼXⱼ₊₁ - h Σⱼ Zⱼ - J₂ Σⱼ XⱼZⱼ₊₁Xⱼ₊₂ -function opsum_infinite(::Model"ising_extended", n; 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() - for j in 1:n - opsum += -J, "X", j, "X", j + 1 - opsum += -J₂, "X", j, "Z", j + 1, "X", j + 2 - opsum -= h, "Z", j - end - return opsum -end - -# For finite Hamiltonian with open boundary conditions -# H = -J Σⱼ XⱼXⱼ₊₁ - h Σⱼ Zⱼ - J₂ Σⱼ XⱼZⱼ₊₁Xⱼ₊₂ -function opsum_finite(::Model"ising_extended", n; J=1.0, h=1.0, J₂=0.0) - # TODO: define using opsum_infinite, filtering sites - # that extend over the boundary. - opsum = OpSum() - for j in 1:n - if j ≤ (n - 1) - opsum += -J, "X", j, "X", j + 1 - end - if j ≤ (n - 2) - opsum += -J₂, "X", j, "Z", j + 1, "X", j + 2 - end - opsum -= h, "Z", j - end - return opsum + opsum -= J, "X", 1, "X", 2 + opsum -= J₂, "X", 1, "Z", 2, "X", 3 + opsum -= h, "Z", 1 + return [opsum] end function reference(::Model"ising_extended", ::Observable"energy"; J=1.0, h=1.0, J₂=0.0) @@ -43,3 +24,21 @@ end # opsum = OpSum(Model"ising_extended"(), 1, 2; J=J, h=h, J₂=J₂) # return prod(MPO(opsum, [s1, s2, s3])) # end + +# # For finite Hamiltonian with open boundary conditions +# # H = -J Σⱼ XⱼXⱼ₊₁ - h Σⱼ Zⱼ - J₂ Σⱼ XⱼZⱼ₊₁Xⱼ₊₂ +# function opsum_finite(::Model"ising_extended", n; J=1.0, h=1.0, J₂=0.0) +# # TODO: define using opsum_infinite, filtering sites +# # that extend over the boundary. +# opsum = OpSum() +# for j in 1:n +# if j ≤ (n - 1) +# opsum += -J, "X", j, "X", j + 1 +# end +# if j ≤ (n - 2) +# opsum += -J₂, "X", j, "Z", j + 1, "X", j + 2 +# end +# opsum -= h, "Z", j +# end +# return opsum +# end diff --git a/src/models/models.jl b/src/models/models.jl index 77a63ed..802faa3 100644 --- a/src/models/models.jl +++ b/src/models/models.jl @@ -110,9 +110,67 @@ function ITensors.MPO(model::Model, s::Vector{<:Index}; kwargs...) return splitblocks(linkinds, MPO(opsum, s)) end -# The ITensor of a single term of the model -function ITensors.ITensor(model::Model, s::Index...; kwargs...) - opsum = opsum_infinite(model, 1; kwargs...) +translatecell(translator, opsum::OpSum, n::Integer) = translator(opsum, n) + +function infinite_terms(model::Model; kwargs...) + opsum_cell = unit_cell_terms(model; kwargs...) + nsites = length(opsum_cell) + 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, n::Int, s::Index...; kwargs...) + opsum = infinite_terms(model; kwargs...)[n] + opsum = shift_sites(opsum, -first_site(opsum) + 1) return contract(MPO(opsum, [s...])) end diff --git a/test/test_vumps_extendedising.jl b/test/test_vumps_extendedising.jl index 607de13..9c5ece7 100644 --- a/test/test_vumps_extendedising.jl +++ b/test/test_vumps_extendedising.jl @@ -33,7 +33,7 @@ using Random 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) From b89c072b36ff8564b2d3121623c92248bf413fa6 Mon Sep 17 00:00:00 2001 From: Matthew Fishman Date: Wed, 22 Jun 2022 18:53:19 -0400 Subject: [PATCH 05/32] Rewrite models in terms of new interface --- .../finite_mps_to_infinite_mps.jl | 0 .../orthogonalize_infinitemps.jl | 0 .../{ => development}/vumps_from_dmrg.jl | 0 src/models/fqhe13.jl | 55 ++++++----- src/models/heisenberg.jl | 42 ++++---- src/models/hubbard.jl | 83 +++++++++------- src/models/ising.jl | 73 +++++++------- src/models/models.jl | 37 +++++-- src/models/xx.jl | 63 ++++++------ src/models/xx_extended.jl | 98 +++++++++++-------- test/test_vumps.jl | 4 +- test/test_vumps_extendedising.jl | 11 ++- 12 files changed, 276 insertions(+), 190 deletions(-) rename examples/{experimental => development}/finite_mps_to_infinite_mps.jl (100%) rename examples/{experimental => development}/orthogonalize_infinitemps.jl (100%) rename examples/vumps/{ => development}/vumps_from_dmrg.jl (100%) diff --git a/examples/experimental/finite_mps_to_infinite_mps.jl b/examples/development/finite_mps_to_infinite_mps.jl similarity index 100% rename from examples/experimental/finite_mps_to_infinite_mps.jl rename to examples/development/finite_mps_to_infinite_mps.jl diff --git a/examples/experimental/orthogonalize_infinitemps.jl b/examples/development/orthogonalize_infinitemps.jl similarity index 100% rename from examples/experimental/orthogonalize_infinitemps.jl rename to examples/development/orthogonalize_infinitemps.jl diff --git a/examples/vumps/vumps_from_dmrg.jl b/examples/vumps/development/vumps_from_dmrg.jl similarity index 100% rename from examples/vumps/vumps_from_dmrg.jl rename to examples/vumps/development/vumps_from_dmrg.jl diff --git a/src/models/fqhe13.jl b/src/models/fqhe13.jl index 431a1ce..10b792c 100644 --- a/src/models/fqhe13.jl +++ b/src/models/fqhe13.jl @@ -1,6 +1,5 @@ -# 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) @@ -8,25 +7,7 @@ function ITensors.OpSum( 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)])) + return [generate_Hamiltonian(opt)] end #Please contact Loic Herviou before using this part of the code for production @@ -322,3 +303,33 @@ function generate_Hamiltonian(coeff::Dict; global_factor=1, prec=1e-12) mpo = OpSum() return generate_Hamiltonian(mpo, coeff; global_factor=global_factor, prec=prec) end + +## function ITensors.OpSum( +## ::Model"fqhe_2b_pot", n1, n2; 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 diff --git a/src/models/heisenberg.jl b/src/models/heisenberg.jl index 87f38ab..2388b99 100644 --- a/src/models/heisenberg.jl +++ b/src/models/heisenberg.jl @@ -1,21 +1,10 @@ -function ITensors.OpSum(::Model"heisenberg", n1, n2) - opsum = OpSum() - opsum += 0.5, "S+", n1, "S-", n2 - opsum += 0.5, "S-", n1, "S+", n2 - opsum += "Sz", n1, "Sz", n2 - 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)) +function unit_cell_terms(::Model"heisenberg") + opsum = OpSum() + opsum += 0.5, "S+", 1, "S-", 2 + opsum += 0.5, "S-", 1, "S+", 2 + opsum += "Sz", 1, "Sz", 2 + return [opsum] end """ @@ -36,3 +25,22 @@ function reference(::Model"heisenberg", ::Observable"energy"; N=∞) correction = 1 + 0.375 / log(N)^3 return (E∞ - Eᶠⁱⁿⁱᵗᵉ * correction) / (2N) end + +# function ITensors.OpSum(::Model"heisenberg", n1, n2) +# opsum = OpSum() +# opsum += 0.5, "S+", n1, "S-", n2 +# opsum += 0.5, "S-", n1, "S+", n2 +# opsum += "Sz", n1, "Sz", n2 +# return opsum +# end +# +# 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 diff --git a/src/models/hubbard.jl b/src/models/hubbard.jl index 8528371..f14be4c 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) 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}, @@ -71,3 +42,43 @@ function reference(::Model"hubbard", ::Observable"energy"; U, Npoints=50) g(i) = -2 * (2π) / Npoints * cos(-π + (i - 1) / Npoints * 2π) * v[i] return sum(g, 1:Npoints) end + +## function ITensors.OpSum(::Model"hubbard", n1, n2; t, U, V) +## 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 +## 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 diff --git a/src/models/ising.jl b/src/models/ising.jl index fdaada6..c8c26de 100644 --- a/src/models/ising.jl +++ b/src/models/ising.jl @@ -1,38 +1,9 @@ # 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 - 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])) + opsum += -J, "X", 1, "X", 2 + opsum += -h, "Z", 1 + return [opsum] end # P. Pfeuty, The one-dimensional Ising model with a transverse field, Annals of Physics 57, p. 79 (1970) @@ -57,3 +28,39 @@ function reference(::Model"ising_classical", ::Observable"critical_inverse_tempe # log1p(x) = log(1 + x) return log1p(sqrt(2)) / 2 end + +# 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) +# 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 +# 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 diff --git a/src/models/models.jl b/src/models/models.jl index 802faa3..b2a6b1f 100644 --- a/src/models/models.jl +++ b/src/models/models.jl @@ -25,8 +25,15 @@ function InfiniteSum{MPO}(model::Model, s::CelledVector; kwargs...) 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] + return InfiniteSum{ITensor}(itensors, translator(s)) +end + # 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)...) @@ -126,7 +133,9 @@ function opsum_infinite(model::Model, nsites::Int; 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))") + 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 @@ -168,12 +177,30 @@ function opsum_finite(model::Model, n::Int; kwargs...) end # The ITensor of a single term `n` of the model. -function ITensors.ITensor(model::Model, n::Int, s::Index...; kwargs...) +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, 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 + # # Deprecated # @@ -185,12 +212,6 @@ end # return prod(MPO(model, s, n; kwargs...)) #modification to allow for more than two sites per term in the Hamiltonians # end -# function InfiniteSum{ITensor}(model::Model, s::CelledVector; kwargs...) -# N = length(s) -# 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; kwargs...) # opsum = OpSum(model, nsites(s); kwargs...) diff --git a/src/models/xx.jl b/src/models/xx.jl index 89a9407..099d287 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 += 1, "X", 1, "X", 2 + os += 1, "Y", 1, "Y", 2 + return os end function reference(::Model"xx", ::Observable"energy"; N=∞) @@ -36,3 +12,34 @@ function reference(::Model"xx", ::Observable"energy"; N=∞) λ(k) = cos(k * π / (N + 1)) return 4 * sum(k -> min(λ(k), 0.0), 1:N) / N end + +## function ITensors.MPO(::Model"xx", s) +## N = length(s) +## 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])) +## end diff --git a/src/models/xx_extended.jl b/src/models/xx_extended.jl index cb94c7f..f71c24b 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 ) @@ -68,3 +36,49 @@ function ITensorInfiniteMPS.reference( temp = sort(λ.(0:(N - 1))) return sum(temp[1:round(Int64, filling * N)]) / N - h * (2 * filling - 1) end + +## function ITensors.OpSum(::Model"xx_extended", n1, n2; 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 +## 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 diff --git a/test/test_vumps.jl b/test/test_vumps.jl index 6718438..22f3cf9 100644 --- a/test/test_vumps.jl +++ b/test/test_vumps.jl @@ -88,10 +88,10 @@ using Random function energy(ψ1, ψ2, h::ITensor) ϕ = ψ1 * ψ2 - return (noprime(ϕ * h) * dag(ϕ))[] + return inner(ϕ, apply(h, ϕ)) end - energy(ψ1, ψ2, h::MPO) = energy(ψ1, ψ2, prod(h)) + energy(ψ1, ψ2, h::MPO) = energy(ψ1, ψ2, contract(h)) function expect(ψ, o) return (noprime(ψ * op(o, filterinds(ψ, "Site")...)) * dag(ψ))[] diff --git a/test/test_vumps_extendedising.jl b/test/test_vumps_extendedising.jl index 9c5ece7..3d4791a 100644 --- a/test/test_vumps_extendedising.jl +++ b/test/test_vumps_extendedising.jl @@ -23,7 +23,9 @@ using Random sfinite = siteinds("S=1/2", Nfinite; conserve_szparity=true) Hfinite = MPO(model, sfinite; model_kwargs...) ψfinite = randomMPS(sfinite, initstate) - energy_finite_total, ψfinite = dmrg(Hfinite, ψfinite; outputlevel=0, nsweeps=30, maxdim=30, cutoff=1e-10) + 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) @@ -33,7 +35,12 @@ using Random nfinite = Nfinite ÷ 2 hnfinite = ITensor( - model, nfinite, 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) From 55512501557997f5c38107578ed84092a29708ce Mon Sep 17 00:00:00 2001 From: Matthew Fishman Date: Wed, 22 Jun 2022 19:29:08 -0400 Subject: [PATCH 06/32] Replace broken subtraction in OpSum with addition --- src/models/ising_extended.jl | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/src/models/ising_extended.jl b/src/models/ising_extended.jl index 91fc707..de9eac2 100644 --- a/src/models/ising_extended.jl +++ b/src/models/ising_extended.jl @@ -3,9 +3,9 @@ # H = -J Σⱼ XⱼXⱼ₊₁ - h Σⱼ Zⱼ - J₂ Σⱼ XⱼZⱼ₊₁Xⱼ₊₂ function unit_cell_terms(::Model"ising_extended"; J=1.0, h=1.0, J₂=0.0) opsum = OpSum() - opsum -= J, "X", 1, "X", 2 - opsum -= J₂, "X", 1, "Z", 2, "X", 3 - opsum -= h, "Z", 1 + opsum += -J, "X", 1, "X", 2 + opsum += -J₂, "X", 1, "Z", 2, "X", 3 + opsum += -h, "Z", 1 return [opsum] end From 0e409bf3b181b59e171172e5753d0c74c02481bf Mon Sep 17 00:00:00 2001 From: Matthew Fishman Date: Wed, 22 Jun 2022 20:19:29 -0400 Subject: [PATCH 07/32] Start testing examples --- examples/vumps/transfer_matrix_spectrum.jl | 15 ++++++++------- examples/vumps/vumps_heisenberg.jl | 6 +++--- src/ITensorInfiniteMPS.jl | 1 + src/abstractinfinitemps.jl | 12 ++++++++++-- src/celledvectors.jl | 5 +++-- src/models/hubbard.jl | 2 +- src/models/models.jl | 4 ++++ test/Project.toml | 1 + test/test_examples.jl | 11 +++++++++++ 9 files changed, 42 insertions(+), 15 deletions(-) create mode 100644 test/test_examples.jl diff --git a/examples/vumps/transfer_matrix_spectrum.jl b/examples/vumps/transfer_matrix_spectrum.jl index 9ff369a..44c5bfa 100644 --- a/examples/vumps/transfer_matrix_spectrum.jl +++ b/examples/vumps/transfer_matrix_spectrum.jl @@ -94,7 +94,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) @@ -115,7 +115,7 @@ Sz2_infinite = expect(ψ.AL[2] * ψ.C[2], "Sz") # using Arpack -using KrylovKit +using KrylovKit: eigsolve using LinearAlgebra T = TransferMatrix(ψ.AL) @@ -138,7 +138,7 @@ tol = 1e-10 neigs = length(v⃗ᴿ) # 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⃗) @@ -152,18 +152,19 @@ 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ⁿ) + [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ⁱᴿ - (translatecell(v⃗ᴸ[1], 1) * vⁱᴿ)[] / norm(v⃗ᴿ[1]) * v⃗ᴿ[1] +#vⁱᴿ² = vⁱᴿ - (translatecelltags(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 λ⃗ᴾᴿ +# 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) diff --git a/examples/vumps/vumps_heisenberg.jl b/examples/vumps/vumps_heisenberg.jl index d042862..9580a95 100644 --- a/examples/vumps/vumps_heisenberg.jl +++ b/examples/vumps/vumps_heisenberg.jl @@ -37,11 +37,11 @@ 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 +@time for _ in 1:outer_iters println("\nIncrease bond dimension") - global ψ = subspace_expansion(ψ, H; subspace_expansion_kwargs...) + global ψ = @time subspace_expansion(ψ, H; subspace_expansion_kwargs...) println("Run VUMPS with new bond dimension") - global ψ = vumps(H, ψ; vumps_kwargs...) + global ψ = @time vumps(H, ψ; vumps_kwargs...) end # Check translational invariance diff --git a/src/ITensorInfiniteMPS.jl b/src/ITensorInfiniteMPS.jl index a136574..738bc0e 100644 --- a/src/ITensorInfiniteMPS.jl +++ b/src/ITensorInfiniteMPS.jl @@ -73,6 +73,7 @@ export Cell, reference, subspace_expansion, translatecell, + translatecelltags, translator, tdvp, vumps, diff --git a/src/abstractinfinitemps.jl b/src/abstractinfinitemps.jl index 2b32f16..9b8ea45 100644 --- a/src/abstractinfinitemps.jl +++ b/src/abstractinfinitemps.jl @@ -206,13 +206,21 @@ 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)) +ITensors.siteinds(ψ::AbstractInfiniteMPS) = CelledVector(siteinds(ψ, Cell(1)), translator(ψ)) +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 c026510..1eddbea 100644 --- a/src/celledvectors.jl +++ b/src/celledvectors.jl @@ -67,8 +67,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) diff --git a/src/models/hubbard.jl b/src/models/hubbard.jl index f14be4c..3ed9cd0 100644 --- a/src/models/hubbard.jl +++ b/src/models/hubbard.jl @@ -6,7 +6,7 @@ function unit_cell_terms(::Model"hubbard"; t, U, V) opsum += -t, "Cdagdn", 2, "Cdn", 1 opsum += U, "Nupdn", 1 opsum += V, "Ntot", 1, "Ntot", 2 - return opsum + return [opsum] end """ diff --git a/src/models/models.jl b/src/models/models.jl index b2a6b1f..1acc7a7 100644 --- a/src/models/models.jl +++ b/src/models/models.jl @@ -183,6 +183,10 @@ function ITensors.ITensor(model::Model, s::Vector{<:Index}, n::Int; kwargs...) 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 diff --git a/test/Project.toml b/test/Project.toml index 340ddf6..69b1a82 100644 --- a/test/Project.toml +++ b/test/Project.toml @@ -1,4 +1,5 @@ [deps] ITensors = "9136182c-28ba-11e9-034c-db9fb085ebd5" 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..7985c7a --- /dev/null +++ b/test/test_examples.jl @@ -0,0 +1,11 @@ +using Test +using Suppressor + +@testset "examples" begin + @testset "$file" for file in joinpath(readdir(@__DIR__), "..", "examples", "vumps") + if endswith(file, ".jl") + println("Running $file") + @suppress include(file) + end + end +end From e111000a95921beafa789c2453269b7001e4812f Mon Sep 17 00:00:00 2001 From: Matt Fishman Date: Wed, 22 Jun 2022 20:20:57 -0400 Subject: [PATCH 08/32] Update examples/vumps/transfer_matrix_spectrum.jl Co-authored-by: github-actions[bot] <41898282+github-actions[bot]@users.noreply.github.com> --- examples/vumps/transfer_matrix_spectrum.jl | 4 +++- 1 file changed, 3 insertions(+), 1 deletion(-) diff --git a/examples/vumps/transfer_matrix_spectrum.jl b/examples/vumps/transfer_matrix_spectrum.jl index 44c5bfa..2759295 100644 --- a/examples/vumps/transfer_matrix_spectrum.jl +++ b/examples/vumps/transfer_matrix_spectrum.jl @@ -152,7 +152,9 @@ 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ⁿ) + [translatecelltags(Lⁿ, 1), translatecelltags(Rⁿ, -1)]; + input_inds=inds(Rⁿ), + output_inds=inds(Lⁿ), ) end From 73ea4492eac2e4a605080e80fb3e636b3ae71d75 Mon Sep 17 00:00:00 2001 From: Matt Fishman Date: Wed, 22 Jun 2022 20:21:05 -0400 Subject: [PATCH 09/32] Update src/abstractinfinitemps.jl Co-authored-by: github-actions[bot] <41898282+github-actions[bot]@users.noreply.github.com> --- src/abstractinfinitemps.jl | 4 +++- 1 file changed, 3 insertions(+), 1 deletion(-) diff --git a/src/abstractinfinitemps.jl b/src/abstractinfinitemps.jl index 9b8ea45..a14c7c0 100644 --- a/src/abstractinfinitemps.jl +++ b/src/abstractinfinitemps.jl @@ -218,7 +218,9 @@ end ITensors.siteinds(ψ::AbstractInfiniteMPS, c::Cell) = siteinds(ψ, siterange(ψ, c)) -ITensors.siteinds(ψ::AbstractInfiniteMPS) = CelledVector(siteinds(ψ, Cell(1)), translator(ψ)) +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]) From f80abc25e8d8fa6dcd1e59b30e06fe34a6fd72b5 Mon Sep 17 00:00:00 2001 From: Matthew Fishman Date: Wed, 22 Jun 2022 20:48:21 -0400 Subject: [PATCH 10/32] Fix path in examples --- test/test_examples.jl | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/test/test_examples.jl b/test/test_examples.jl index 7985c7a..a284669 100644 --- a/test/test_examples.jl +++ b/test/test_examples.jl @@ -1,8 +1,9 @@ using Test using Suppressor +examples_dir = joinpath(@__DIR__, "..", "examples", "vumps") @testset "examples" begin - @testset "$file" for file in joinpath(readdir(@__DIR__), "..", "examples", "vumps") + @testset "$file" for file in readdir(examples_dir; join=true) if endswith(file, ".jl") println("Running $file") @suppress include(file) From e1ec11ac80505e996edfdc11ae3ab43ab5934c49 Mon Sep 17 00:00:00 2001 From: Matthew Fishman Date: Thu, 23 Jun 2022 08:42:19 -0400 Subject: [PATCH 11/32] Move some old examples to development directory --- .../transfer_matrix_spectrum_arpack.jl | 200 ++++++++++++++++++ .../vumps_ising_noncontiguous.jl | 2 +- .../vumps/{ => development}/vumps_localham.jl | 2 - .../{ => development}/vumps_localham_qns.jl | 0 examples/vumps/transfer_matrix_spectrum.jl | 53 +---- 5 files changed, 207 insertions(+), 50 deletions(-) create mode 100644 examples/vumps/development/transfer_matrix_spectrum_arpack.jl rename examples/vumps/{ => development}/vumps_ising_noncontiguous.jl (97%) rename examples/vumps/{ => development}/vumps_localham.jl (97%) rename examples/vumps/{ => development}/vumps_localham_qns.jl (100%) diff --git a/examples/vumps/development/transfer_matrix_spectrum_arpack.jl b/examples/vumps/development/transfer_matrix_spectrum_arpack.jl new file mode 100644 index 0000000..2759295 --- /dev/null +++ b/examples/vumps/development/transfer_matrix_spectrum_arpack.jl @@ -0,0 +1,200 @@ +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 = true +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"() + +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) = "↑" +ψ = 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 = length(v⃗ᴿ) + +# 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_ising_noncontiguous.jl b/examples/vumps/development/vumps_ising_noncontiguous.jl similarity index 97% rename from examples/vumps/vumps_ising_noncontiguous.jl rename to examples/vumps/development/vumps_ising_noncontiguous.jl index 3fc08bf..e830d1d 100644 --- a/examples/vumps/vumps_ising_noncontiguous.jl +++ b/examples/vumps/development/vumps_ising_noncontiguous.jl @@ -70,7 +70,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/vumps/development/vumps_localham.jl similarity index 97% rename from examples/vumps/vumps_localham.jl rename to examples/vumps/development/vumps_localham.jl index 55082c8..14ff15c 100644 --- a/examples/vumps/vumps_localham.jl +++ b/examples/vumps/development/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/vumps/development/vumps_localham_qns.jl similarity index 100% rename from examples/vumps/vumps_localham_qns.jl rename to examples/vumps/development/vumps_localham_qns.jl diff --git a/examples/vumps/transfer_matrix_spectrum.jl b/examples/vumps/transfer_matrix_spectrum.jl index 2759295..05f8141 100644 --- a/examples/vumps/transfer_matrix_spectrum.jl +++ b/examples/vumps/transfer_matrix_spectrum.jl @@ -23,7 +23,7 @@ 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 @@ -114,7 +114,6 @@ Sz2_infinite = expect(ψ.AL[2] * ψ.C[2], "Sz") # Compute eigenspace of the transfer matrix # -using Arpack using KrylovKit: eigsolve using LinearAlgebra @@ -143,48 +142,11 @@ 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 - +# 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)) @@ -195,6 +157,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]] - λ⃗ᴿᴬ From 0637def80f1fd8d1e0cb97f0e5bb0658b651d104 Mon Sep 17 00:00:00 2001 From: Matthew Fishman Date: Thu, 23 Jun 2022 10:42:12 -0400 Subject: [PATCH 12/32] Start simplifying interface for defining site indices --- .../vumps}/transfer_matrix_spectrum_arpack.jl | 0 .../vumps}/vumps_from_dmrg.jl | 0 .../vumps}/vumps_ising_noncontiguous.jl | 0 .../vumps}/vumps_localham.jl | 0 .../vumps}/vumps_localham_qns.jl | 0 .../vumps}/vumps_mpo.jl | 0 examples/vumps/vumps_heisenberg.jl | 5 +- src/celledvectors.jl | 8 +- src/infinitecanonicalmps.jl | 86 ++++++++++++------- 9 files changed, 59 insertions(+), 40 deletions(-) rename examples/{vumps/development => development/vumps}/transfer_matrix_spectrum_arpack.jl (100%) rename examples/{vumps/development => development/vumps}/vumps_from_dmrg.jl (100%) rename examples/{vumps/development => development/vumps}/vumps_ising_noncontiguous.jl (100%) rename examples/{vumps/development => development/vumps}/vumps_localham.jl (100%) rename examples/{vumps/development => development/vumps}/vumps_localham_qns.jl (100%) rename examples/{vumps/development => development/vumps}/vumps_mpo.jl (100%) diff --git a/examples/vumps/development/transfer_matrix_spectrum_arpack.jl b/examples/development/vumps/transfer_matrix_spectrum_arpack.jl similarity index 100% rename from examples/vumps/development/transfer_matrix_spectrum_arpack.jl rename to examples/development/vumps/transfer_matrix_spectrum_arpack.jl diff --git a/examples/vumps/development/vumps_from_dmrg.jl b/examples/development/vumps/vumps_from_dmrg.jl similarity index 100% rename from examples/vumps/development/vumps_from_dmrg.jl rename to examples/development/vumps/vumps_from_dmrg.jl diff --git a/examples/vumps/development/vumps_ising_noncontiguous.jl b/examples/development/vumps/vumps_ising_noncontiguous.jl similarity index 100% rename from examples/vumps/development/vumps_ising_noncontiguous.jl rename to examples/development/vumps/vumps_ising_noncontiguous.jl diff --git a/examples/vumps/development/vumps_localham.jl b/examples/development/vumps/vumps_localham.jl similarity index 100% rename from examples/vumps/development/vumps_localham.jl rename to examples/development/vumps/vumps_localham.jl diff --git a/examples/vumps/development/vumps_localham_qns.jl b/examples/development/vumps/vumps_localham_qns.jl similarity index 100% rename from examples/vumps/development/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 100% rename from examples/vumps/development/vumps_mpo.jl rename to examples/development/vumps/vumps_mpo.jl diff --git a/examples/vumps/vumps_heisenberg.jl b/examples/vumps/vumps_heisenberg.jl index 9580a95..d96db62 100644 --- a/examples/vumps/vumps_heisenberg.jl +++ b/examples/vumps/vumps_heisenberg.jl @@ -17,11 +17,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=false, initstate) ψ = InfMPS(s, initstate) model = Model("heisenberg") diff --git a/src/celledvectors.jl b/src/celledvectors.jl index 1eddbea..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) diff --git a/src/infinitecanonicalmps.jl b/src/infinitecanonicalmps.jl index 3f1f04e..14f93f4 100644 --- a/src/infinitecanonicalmps.jl +++ b/src/infinitecanonicalmps.jl @@ -1,4 +1,18 @@ -# Get the promoted type of the Index objects in a collection +# TODO: Move to ITensors.jl +setval(qnval::ITensors.QNVal, val::Int) = ITensors.QNVal(ITensors.name(qnval), val, ITensors.modulus(qnval)) + +# TODO: Move to ITensors.jl +function Base.:/(qnval::ITensors.QNVal, n::Int) + return setval(qnval, Int(ITensors.val(qnval) / n)) +end + +# TODO: Move to ITensors.jl +function Base.:/(qn::QN, n::Int) + @show qn + @show n + 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 +23,34 @@ 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] -end +shift_flux(qnblock::Pair{QN,Int}, flux_density::QN) = ((ITensors.qn(qnblock) - flux_density) => ITensors.blockdim(qnblock)) +shift_flux(space::Vector{Pair{QN,Int}}, flux_density::QN) = map(qnblock -> shift_flux(qnblock, flux_density), space) +shift_flux(i::Index, flux_density::QN) = ITensors.setspace(i, shift_flux(space(i), flux_density)) -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_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; 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) - 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 +58,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 +66,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 @@ -162,3 +169,20 @@ end function ITensors.expect(ψ::InfiniteCanonicalMPS, h::InfiniteSum) return [expect(ψ, h[j]) for j in 1:nsites(ψ)] end + +## XXX: Delete +## # More general siteind that allows specifying +## # the space +## function _siteind(site_tag, n::Int; space) +## return addtags(Index(space, "Site,n=$n"), site_tag) +## end +## +## _siteinds(site_tag, N::Int; space) = __siteinds(site_tag, N, space) +## +## function __siteinds(site_tag, N::Int, space::Vector) +## return [_siteind(site_tag, n; space=space[n]) for n in 1:N] +## end +## +## function __siteinds(site_tag, N::Int, space) +## return [_siteind(site_tag, n; space=space) for n in 1:N] +## end From a974f025129131190cfe92740aff6f4d352b4070 Mon Sep 17 00:00:00 2001 From: Matthew Fishman Date: Thu, 23 Jun 2022 12:44:56 -0400 Subject: [PATCH 13/32] Update examples for simpler interface --- .../vumps/src/vumps_subspace_expansion.jl | 20 +++++++ examples/vumps/transfer_matrix_spectrum.jl | 33 +++++------- examples/vumps/vumps_heisenberg.jl | 50 ++++++++--------- examples/vumps/vumps_hubbard_extended.jl | 53 +++++++++++-------- examples/vumps/vumps_ising.jl | 30 +++-------- examples/vumps/vumps_ising_extended.jl | 29 ++++------ src/infinitecanonicalmps.jl | 8 +-- 7 files changed, 110 insertions(+), 113 deletions(-) create mode 100644 examples/vumps/src/vumps_subspace_expansion.jl diff --git a/examples/vumps/src/vumps_subspace_expansion.jl b/examples/vumps/src/vumps_subspace_expansion.jl new file mode 100644 index 0000000..70951d8 --- /dev/null +++ b/examples/vumps/src/vumps_subspace_expansion.jl @@ -0,0 +1,20 @@ +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 05f8141..63f817e 100644 --- a/examples/vumps/transfer_matrix_spectrum.jl +++ b/examples/vumps/transfer_matrix_spectrum.jl @@ -1,6 +1,8 @@ using ITensors using ITensorInfiniteMPS +include(joinpath(@__DIR__, "src", "vumps_subspace_expansion.jl")) + ############################################################################## # VUMPS parameters # @@ -25,17 +27,8 @@ model_params = (J=1.0, h=0.9) 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 +44,8 @@ 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 +55,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) @@ -127,6 +113,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]) @@ -134,7 +122,9 @@ 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⃗ = [(translatecelltags(v⃗ᴸ[n], 1) * v⃗ᴿ[n])[] for n in 1:neigs] @@ -150,6 +140,7 @@ v⃗ᴸ ./= sqrt.(N⃗) 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)) diff --git a/examples/vumps/vumps_heisenberg.jl b/examples/vumps/vumps_heisenberg.jl index d96db62..43df200 100644 --- a/examples/vumps/vumps_heisenberg.jl +++ b/examples/vumps/vumps_heisenberg.jl @@ -1,6 +1,8 @@ using ITensors using ITensorInfiniteMPS +include(joinpath(@__DIR__, "src", "vumps_subspace_expansion.jl")) + ############################################################################## # VUMPS parameters # @@ -8,8 +10,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 @@ -18,7 +21,7 @@ outer_iters = 8 # Number of times to increase the bond dimension N = 2 # Number of sites in the unit cell initstate(n) = isodd(n) ? "↑" : "↓" -s = infsiteinds("S=1/2", N; conserve_qns=false, initstate) +s = infsiteinds("S=1/2", N; conserve_qns, initstate) ψ = InfMPS(s, initstate) model = Model("heisenberg") @@ -27,21 +30,16 @@ 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 -@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 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) @@ -65,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(ψ))[] + return inner(ψ, apply(op(o, filterinds(ψ, "Site")...), ψ)) end nfinite = Nfinite ÷ 2 @@ -97,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..a1c32b9 100644 --- a/examples/vumps/vumps_hubbard_extended.jl +++ b/examples/vumps/vumps_hubbard_extended.jl @@ -1,6 +1,8 @@ using ITensors using ITensorInfiniteMPS +include(joinpath(@__DIR__, "src", "vumps_subspace_expansion.jl")) + ############################################################################## # VUMPS parameters # @@ -8,6 +10,7 @@ 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 eager = true @@ -23,18 +26,19 @@ 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 +## 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) -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) ψ = InfMPS(s, initstate) model = Model"hubbard"() @@ -48,7 +52,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=outputlevel, eager) subspace_expansion_kwargs = (cutoff=cutoff, maxdim=maxdim) # For now, to increase the bond dimension you must alternate @@ -57,14 +61,16 @@ 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...) +ψ = vumps_subspace_expansion(H, ψ; outer_iters, subspace_expansion_kwargs, 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 +# ψ = @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 # Check translational invariance println("\nCheck translational invariance of optimized infinite MPS") @@ -108,14 +114,17 @@ 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..31165a2 100644 --- a/examples/vumps/vumps_ising.jl +++ b/examples/vumps/vumps_ising.jl @@ -1,5 +1,7 @@ -using ITensorInfiniteMPS using ITensors +using ITensorInfiniteMPS + +include(joinpath(@__DIR__, "src", "vumps_subspace_expansion.jl")) ############################################################################## # VUMPS/TDVP parameters @@ -24,22 +26,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 +46,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 +56,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 2f9cf36..7119e22 100644 --- a/examples/vumps/vumps_ising_extended.jl +++ b/examples/vumps/vumps_ising_extended.jl @@ -1,6 +1,8 @@ using ITensors using ITensorInfiniteMPS +include(joinpath(@__DIR__, "src", "vumps_subspace_expansion.jl")) + ############################################################################## # VUMPS parameters # @@ -9,24 +11,20 @@ maxdim = 64 # Maximum bond dimension cutoff = 1e-6 # Singular value cutoff when increasing the 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,26 +34,19 @@ 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) dmrg_kwargs = (nsweeps=10, maxdim, cutoff) energy_finite_total, ψfinite = dmrg(Hfinite, ψfinite; dmrg_kwargs...) -energy_finite_total, ψfinite = dmrg(Hfinite, ψfinite, sweeps) nfinite = Nfinite ÷ 2 Sz_finite = expect(ψfinite, "Sz")[nfinite:(nfinite + N - 1)] @@ -66,5 +57,5 @@ println("\nEnergy") @show reference(model, Observable("energy"); J=J, h=h, J₂=J₂) println("\nSz") -@show Sz +@show Sz_infinite @show Sz_finite diff --git a/src/infinitecanonicalmps.jl b/src/infinitecanonicalmps.jl index 14f93f4..e73a7cf 100644 --- a/src/infinitecanonicalmps.jl +++ b/src/infinitecanonicalmps.jl @@ -3,13 +3,15 @@ setval(qnval::ITensors.QNVal, val::Int) = ITensors.QNVal(ITensors.name(qnval), v # TODO: Move to ITensors.jl function Base.:/(qnval::ITensors.QNVal, n::Int) - return setval(qnval, Int(ITensors.val(qnval) / n)) + 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) - @show qn - @show n return QN(map(qnval -> qnval / n, qn.data)) end From 6ae2e2370bb983cd350ba7e25da67130fcf111e7 Mon Sep 17 00:00:00 2001 From: Matthew Fishman Date: Thu, 23 Jun 2022 15:06:31 -0400 Subject: [PATCH 14/32] Continue updating tests for new interface --- .../vumps/src/vumps_subspace_expansion.jl | 16 +++-- examples/vumps/transfer_matrix_spectrum.jl | 4 +- src/infinitecanonicalmps.jl | 20 +++++-- test/test_readwrite.jl | 9 +-- test/test_vumps.jl | 60 +++++++------------ test/test_vumps_extendedising.jl | 15 +---- test/test_vumpsmpo.jl | 30 +++++----- 7 files changed, 71 insertions(+), 83 deletions(-) diff --git a/examples/vumps/src/vumps_subspace_expansion.jl b/examples/vumps/src/vumps_subspace_expansion.jl index 70951d8..bca3068 100644 --- a/examples/vumps/src/vumps_subspace_expansion.jl +++ b/examples/vumps/src/vumps_subspace_expansion.jl @@ -3,10 +3,14 @@ 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) +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])") + 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...) @@ -15,6 +19,10 @@ function tdvp_subspace_expansion(H, ψ; time_step, outer_iters, subspace_expansi 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) +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 63f817e..300324d 100644 --- a/examples/vumps/transfer_matrix_spectrum.jl +++ b/examples/vumps/transfer_matrix_spectrum.jl @@ -45,7 +45,9 @@ vumps_kwargs = ( ) subspace_expansion_kwargs = (cutoff=cutoff, maxdim=maxdim) -ψ = tdvp_subspace_expansion(H, ψ; time_step, outer_iters, subspace_expansion_kwargs, vumps_kwargs) +ψ = 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]...)) diff --git a/src/infinitecanonicalmps.jl b/src/infinitecanonicalmps.jl index e73a7cf..60b83ee 100644 --- a/src/infinitecanonicalmps.jl +++ b/src/infinitecanonicalmps.jl @@ -1,5 +1,7 @@ # TODO: Move to ITensors.jl -setval(qnval::ITensors.QNVal, val::Int) = ITensors.QNVal(ITensors.name(qnval), val, ITensors.modulus(qnval)) +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) @@ -36,9 +38,15 @@ function shift_flux_to_zero(s::Vector{<:Index}, initstate::Function) return shift_flux_to_zero(s, flux(MPS(s, initstate))) end -shift_flux(qnblock::Pair{QN,Int}, flux_density::QN) = ((ITensors.qn(qnblock) - flux_density) => ITensors.blockdim(qnblock)) -shift_flux(space::Vector{Pair{QN,Int}}, flux_density::QN) = map(qnblock -> shift_flux(qnblock, flux_density), space) -shift_flux(i::Index, flux_density::QN) = ITensors.setspace(i, shift_flux(space(i), flux_density)) +function shift_flux(qnblock::Pair{QN,Int}, flux_density::QN) + return ((ITensors.qn(qnblock) - flux_density) => ITensors.blockdim(qnblock)) +end +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 shift_flux_to_zero(s::Vector{<:Index}, flux::QN) if iszero(flux) @@ -49,7 +57,9 @@ function shift_flux_to_zero(s::Vector{<:Index}, flux::QN) return map(sₙ -> shift_flux(sₙ, flux_density), s) end -function infsiteinds(site_tag, n::Int; translator=translatecelltags, initstate=nothing, kwargs...) +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) 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 22f3cf9..063105d 100644 --- a/test/test_vumps.jl +++ b/test/test_vumps.jl @@ -6,17 +6,9 @@ using Random @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 +23,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 +37,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 @@ -94,7 +85,7 @@ using Random energy(ψ1, ψ2, h::MPO) = energy(ψ1, ψ2, contract(h)) function expect(ψ, o) - return (noprime(ψ * op(o, filterinds(ψ, "Site")...)) * dag(ψ))[] + return inner(ψ, apply(op(o, filterinds(ψ, "Site")...), ψ)) end nfinite = Nfinite ÷ 2 @@ -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-10outputlevel = 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 diff --git a/test/test_vumps_extendedising.jl b/test/test_vumps_extendedising.jl index 3d4791a..2202929 100644 --- a/test/test_vumps_extendedising.jl +++ b/test/test_vumps_extendedising.jl @@ -6,18 +6,10 @@ using Random @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) @@ -30,7 +22,7 @@ using Random function energy(ψ, h, n) ϕ = ψ[n] * ψ[n + 1] * ψ[n + 2] - return (noprime(ϕ * h) * dag(ϕ))[] + return inner(ϕ, apply(h, ϕ)) end nfinite = Nfinite ÷ 2 @@ -64,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) diff --git a/test/test_vumpsmpo.jl b/test/test_vumpsmpo.jl index 51b3dec..f8d6d75 100644 --- a/test/test_vumpsmpo.jl +++ b/test/test_vumpsmpo.jl @@ -7,16 +7,16 @@ using Random @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 + ## 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 @@ -32,6 +32,7 @@ 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) @@ -40,7 +41,7 @@ using Random function energy(ψ, h, n) ϕ = ψ[n] * ψ[n + 1] * ψ[n + 2] - return (noprime(ϕ * h) * dag(ϕ))[] + return inner(ϕ, (apply(h, ϕ))) end nfinite = Nfinite ÷ 2 @@ -65,8 +66,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_qns) ψ = InfMPS(s, initstate) Hmpo = InfiniteMPOMatrix(model, s; model_kwargs...) @@ -121,10 +121,10 @@ 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) From 1186278db6a96c885ab3f915178a4899921814b4 Mon Sep 17 00:00:00 2001 From: Matthew Fishman Date: Thu, 23 Jun 2022 15:15:01 -0400 Subject: [PATCH 15/32] Format --- test/test_vumpsmpo.jl | 42 +++++++++++------------------------------- 1 file changed, 11 insertions(+), 31 deletions(-) diff --git a/test/test_vumpsmpo.jl b/test/test_vumpsmpo.jl index f8d6d75..6912565 100644 --- a/test/test_vumpsmpo.jl +++ b/test/test_vumpsmpo.jl @@ -10,14 +10,6 @@ using Random 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 @@ -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 @@ -129,7 +113,7 @@ end function energy(ψ, h, n) ϕ = ψ[n] * ψ[n + 1] * ψ[n + 2] - return (noprime(ϕ * h) * dag(ϕ))[] + return inner(ϕ, apply(h, ϕ)) end nfinite = Nfinite ÷ 2 @@ -153,8 +137,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 +164,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 @@ -240,9 +215,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...) From 82b23c4045de8ed0b29b2741d9cf2aa1581a7059 Mon Sep 17 00:00:00 2001 From: Matthew Fishman Date: Thu, 23 Jun 2022 16:05:47 -0400 Subject: [PATCH 16/32] Fix more tests --- examples/vumps/transfer_matrix_spectrum.jl | 6 +++--- examples/vumps/vumps_heisenberg.jl | 10 +++------ examples/vumps/vumps_hubbard_extended.jl | 24 ---------------------- src/infinitecanonicalmps.jl | 6 ++++-- test/Project.toml | 1 + test/test_vumps.jl | 2 +- 6 files changed, 12 insertions(+), 37 deletions(-) diff --git a/examples/vumps/transfer_matrix_spectrum.jl b/examples/vumps/transfer_matrix_spectrum.jl index 300324d..f383906 100644 --- a/examples/vumps/transfer_matrix_spectrum.jl +++ b/examples/vumps/transfer_matrix_spectrum.jl @@ -69,11 +69,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 diff --git a/examples/vumps/vumps_heisenberg.jl b/examples/vumps/vumps_heisenberg.jl index 43df200..8133cd9 100644 --- a/examples/vumps/vumps_heisenberg.jl +++ b/examples/vumps/vumps_heisenberg.jl @@ -42,18 +42,14 @@ subspace_expansion_kwargs = (cutoff=cutoff, maxdim=maxdim) 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] @@ -79,7 +75,7 @@ energy_finite_total, ψfinite = dmrg(Hfinite, ψfinite; nsweeps, maxdim, cutoff) energy_exact_finite = reference(model, Observable("energy"); N=Nfinite) -function ITensors.expect(ψ, o) +function ITensors.expect(ψ::ITensor, o::String) return inner(ψ, apply(op(o, filterinds(ψ, "Site")...), ψ)) end diff --git a/examples/vumps/vumps_hubbard_extended.jl b/examples/vumps/vumps_hubbard_extended.jl index a1c32b9..34baefd 100644 --- a/examples/vumps/vumps_hubbard_extended.jl +++ b/examples/vumps/vumps_hubbard_extended.jl @@ -26,17 +26,6 @@ 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) - initstate(n) = isodd(n) ? "↑" : "↓" s = infsiteinds("Electron", N; initstate) ψ = InfMPS(s, initstate) @@ -63,23 +52,10 @@ subspace_expansion_kwargs = (cutoff=cutoff, maxdim=maxdim) println("\nRun VUMPS on initial product state, unit cell size $N") ψ = vumps_subspace_expansion(H, ψ; outer_iters, subspace_expansion_kwargs, vumps_kwargs) -# ψ = @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 - # 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] diff --git a/src/infinitecanonicalmps.jl b/src/infinitecanonicalmps.jl index 60b83ee..8983935 100644 --- a/src/infinitecanonicalmps.jl +++ b/src/infinitecanonicalmps.jl @@ -148,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/test/Project.toml b/test/Project.toml index 69b1a82..b022cab 100644 --- a/test/Project.toml +++ b/test/Project.toml @@ -1,5 +1,6 @@ [deps] ITensors = "9136182c-28ba-11e9-034c-db9fb085ebd5" +KrylovKit = "0b1a1467-8014-51b9-945f-bf0ae24f4b77" Random = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c" Suppressor = "fd094767-a336-5f1f-9728-57cf17d0bbfb" Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40" diff --git a/test/test_vumps.jl b/test/test_vumps.jl index 063105d..ef6ffe9 100644 --- a/test/test_vumps.jl +++ b/test/test_vumps.jl @@ -142,7 +142,7 @@ end ψfinite = randomMPS(sfinite, initstate) nsweeps = 20 energy_finite_total, ψfinite = dmrg( - Hfinite, ψfinite; nsweeps, maxdims=10, cutoff=1e-10outputlevel = 0 + Hfinite, ψfinite; nsweeps, maxdims=10, cutoff=1e-10, outputlevel=0 ) multisite_update_alg = "sequential" From 3f98fbc62d0e4c828fbb8d77a0ba99d3b4278795 Mon Sep 17 00:00:00 2001 From: Matthew Fishman Date: Thu, 23 Jun 2022 16:15:12 -0400 Subject: [PATCH 17/32] Clean up tests --- test/test_vumps.jl | 47 +++++++++++++------------------- test/test_vumps_extendedising.jl | 12 ++++---- test/test_vumpsmpo.jl | 32 ++++++++-------------- 3 files changed, 36 insertions(+), 55 deletions(-) diff --git a/test/test_vumps.jl b/test/test_vumps.jl index ef6ffe9..bed4ecf 100644 --- a/test/test_vumps.jl +++ b/test/test_vumps.jl @@ -3,6 +3,17 @@ 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(ψ::ITensor, o::String) + return inner(ψ, apply(op(o, filterinds(ψ, "Site")...), ψ)) +end + @testset "vumps" begin Random.seed!(1234) @@ -77,29 +88,18 @@ 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 inner(ϕ, apply(h, ϕ)) - end - - energy(ψ1, ψ2, h::MPO) = energy(ψ1, ψ2, contract(h)) - - function expect(ψ, o) - return inner(ψ, apply(op(o, filterinds(ψ, "Site")...), ψ)) - 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") @@ -179,27 +179,18 @@ 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") diff --git a/test/test_vumps_extendedising.jl b/test/test_vumps_extendedising.jl index 2202929..7842cf2 100644 --- a/test/test_vumps_extendedising.jl +++ b/test/test_vumps_extendedising.jl @@ -3,6 +3,11 @@ 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) @@ -20,11 +25,6 @@ using Random ) Szs_finite = expect(ψfinite, "Sz") - function energy(ψ, h, n) - ϕ = ψ[n] * ψ[n + 1] * ψ[n + 2] - return inner(ϕ, apply(h, ϕ)) - end - nfinite = Nfinite ÷ 2 hnfinite = ITensor( model, @@ -35,7 +35,7 @@ using Random model_kwargs..., ) orthogonalize!(ψfinite, nfinite) - energy_finite = energy(ψfinite, hnfinite, nfinite) + energy_finite = expect_three_site(ψfinite, hnfinite, nfinite) cutoff = 1e-8 maxdim = 100 diff --git a/test/test_vumpsmpo.jl b/test/test_vumpsmpo.jl index 6912565..d6cb7d9 100644 --- a/test/test_vumpsmpo.jl +++ b/test/test_vumpsmpo.jl @@ -3,6 +3,11 @@ 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) @@ -31,15 +36,10 @@ using Random 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 inner(ϕ, (apply(h, ϕ))) - 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 [ @@ -58,7 +58,7 @@ using Random ) subspace_expansion_kwargs = (cutoff=cutoff, maxdim=maxdim) - s = infsiteinds("S=1/2", nsites; initstate, conserve_qns) + s = infsiteinds("S=1/2", nsites; initstate, conserve_szparity=conserve_qns) ψ = InfMPS(s, initstate) Hmpo = InfiniteMPOMatrix(model, s; model_kwargs...) @@ -111,20 +111,15 @@ end ) Szs_finite = expect(ψfinite, "Sz") - function energy(ψ, h, n) - ϕ = ψ[n] * ψ[n + 1] * ψ[n + 2] - return inner(ϕ, apply(h, ϕ)) - 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] @@ -187,22 +182,17 @@ 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] From a145680e7bd1faa37481de41126402bd9d686e04 Mon Sep 17 00:00:00 2001 From: Matthew Fishman Date: Thu, 23 Jun 2022 16:36:44 -0400 Subject: [PATCH 18/32] Small test fix --- test/test_vumps.jl | 18 +++++++++--------- 1 file changed, 9 insertions(+), 9 deletions(-) diff --git a/test/test_vumps.jl b/test/test_vumps.jl index bed4ecf..556540c 100644 --- a/test/test_vumps.jl +++ b/test/test_vumps.jl @@ -10,7 +10,7 @@ end expect_two_site(ψ1::ITensor, ψ2::ITensor, h::MPO) = expect_two_site(ψ1, ψ2, contract(h)) -function expect(ψ::ITensor, o::String) +function expect_one_site(ψ::ITensor, o::String) return inner(ψ, apply(op(o, filterinds(ψ, "Site")...), ψ)) end @@ -102,13 +102,13 @@ end 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 @@ -193,13 +193,13 @@ end 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 From 847ac332a8d3676a6ce73d1ede3e4b4daf68a96a Mon Sep 17 00:00:00 2001 From: Matthew Fishman Date: Thu, 23 Jun 2022 17:10:19 -0400 Subject: [PATCH 19/32] Fix Hubbard example, fix more tests --- examples/vumps/vumps_hubbard_extended.jl | 7 ++++--- test/Project.toml | 1 + test/test_vumpsmpo.jl | 5 +++++ 3 files changed, 10 insertions(+), 3 deletions(-) diff --git a/examples/vumps/vumps_hubbard_extended.jl b/examples/vumps/vumps_hubbard_extended.jl index 34baefd..880a5e0 100644 --- a/examples/vumps/vumps_hubbard_extended.jl +++ b/examples/vumps/vumps_hubbard_extended.jl @@ -13,6 +13,7 @@ max_vumps_iters = 200 # Maximum number of iterations of the VUMPS algorithm at e 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) @@ -27,7 +28,7 @@ N = 2 # Unit cell size @show localham_type initstate(n) = isodd(n) ? "↑" : "↓" -s = infsiteinds("Electron", N; initstate) +s = infsiteinds("Electron", N; initstate, conserve_qns) ψ = InfMPS(s, initstate) model = Model"hubbard"() @@ -41,7 +42,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=vumps_tol, 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 @@ -85,7 +86,7 @@ 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") diff --git a/test/Project.toml b/test/Project.toml index b022cab..b9fb740 100644 --- a/test/Project.toml +++ b/test/Project.toml @@ -1,6 +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_vumpsmpo.jl b/test/test_vumpsmpo.jl index d6cb7d9..98b1185 100644 --- a/test/test_vumpsmpo.jl +++ b/test/test_vumpsmpo.jl @@ -196,6 +196,11 @@ end 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, From 9cfd0d3ec3ad9944d3446af03cc9caf2c7716d0e Mon Sep 17 00:00:00 2001 From: Matthew Fishman Date: Thu, 23 Jun 2022 18:05:34 -0400 Subject: [PATCH 20/32] Update FQHE example for new interface --- test/test_vumpsmpo_fqhe.jl | 49 +++++++++++++++++++------------------- 1 file changed, 24 insertions(+), 25 deletions(-) diff --git a/test/test_vumpsmpo_fqhe.jl b/test/test_vumpsmpo_fqhe.jl index 49b8b9a..1f227b8 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,14 @@ using Random nsites in [6], time_step in [-Inf] - vumps_kwargs = ( - multisite_update_alg=multisite_update_alg, - tol=tol, - maxiter=maxiter, + vumps_kwargs = (; + multisite_update_alg, + tol, + maxiter, outputlevel=0, - time_step=time_step, + time_step, ) - subspace_expansion_kwargs = (cutoff=cutoff, maxdim=maxdim) + subspace_expansion_kwargs = (; cutoff, maxdim) function fermion_momentum_translator(i::Index, n::Integer; N=nsites) ts = tags(i) @@ -75,13 +79,8 @@ 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) From f4451e744dcfeee4cbf30c30414988f398797fe7 Mon Sep 17 00:00:00 2001 From: Matt Fishman Date: Thu, 23 Jun 2022 18:18:06 -0400 Subject: [PATCH 21/32] Update test/test_vumpsmpo_fqhe.jl Co-authored-by: github-actions[bot] <41898282+github-actions[bot]@users.noreply.github.com> --- test/test_vumpsmpo_fqhe.jl | 8 +++++++- 1 file changed, 7 insertions(+), 1 deletion(-) diff --git a/test/test_vumpsmpo_fqhe.jl b/test/test_vumpsmpo_fqhe.jl index 1f227b8..626deac 100644 --- a/test/test_vumpsmpo_fqhe.jl +++ b/test/test_vumpsmpo_fqhe.jl @@ -80,7 +80,13 @@ end end s = infsiteinds( - "FermionK", nsites; initstate, translator=fermion_momentum_translator, p, q, conserve_momentum + "FermionK", + nsites; + initstate, + translator=fermion_momentum_translator, + p, + q, + conserve_momentum, ) ψ = InfMPS(s, initstate) From 24663c3f48ad6ad5779f2705990e600a5f02a275 Mon Sep 17 00:00:00 2001 From: Matthew Fishman Date: Thu, 23 Jun 2022 18:25:03 -0400 Subject: [PATCH 22/32] Format and clean up code --- src/infinitecanonicalmps.jl | 17 ------------- src/infinitempo.jl | 9 ------- src/models/fqhe13.jl | 30 ----------------------- src/models/heisenberg.jl | 19 --------------- src/models/hubbard.jl | 40 ------------------------------- src/models/ising.jl | 36 ---------------------------- src/models/ising_extended.jl | 29 ----------------------- src/models/models.jl | 30 ----------------------- src/models/xx.jl | 31 ------------------------ src/models/xx_extended.jl | 46 ------------------------------------ test/test_vumpsmpo_fqhe.jl | 16 ++++++------- 11 files changed, 8 insertions(+), 295 deletions(-) diff --git a/src/infinitecanonicalmps.jl b/src/infinitecanonicalmps.jl index 8983935..2831ba0 100644 --- a/src/infinitecanonicalmps.jl +++ b/src/infinitecanonicalmps.jl @@ -183,20 +183,3 @@ end function ITensors.expect(ψ::InfiniteCanonicalMPS, h::InfiniteSum) return [expect(ψ, h[j]) for j in 1:nsites(ψ)] end - -## XXX: Delete -## # More general siteind that allows specifying -## # the space -## function _siteind(site_tag, n::Int; space) -## return addtags(Index(space, "Site,n=$n"), site_tag) -## end -## -## _siteinds(site_tag, N::Int; space) = __siteinds(site_tag, N, space) -## -## function __siteinds(site_tag, N::Int, space::Vector) -## return [_siteind(site_tag, n; space=space[n]) for n in 1:N] -## end -## -## function __siteinds(site_tag, N::Int, space) -## return [_siteind(site_tag, n; space=space) for n in 1:N] -## end 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/models/fqhe13.jl b/src/models/fqhe13.jl index 10b792c..851aefb 100644 --- a/src/models/fqhe13.jl +++ b/src/models/fqhe13.jl @@ -303,33 +303,3 @@ function generate_Hamiltonian(coeff::Dict; global_factor=1, prec=1e-12) mpo = OpSum() return generate_Hamiltonian(mpo, coeff; global_factor=global_factor, prec=prec) end - -## function ITensors.OpSum( -## ::Model"fqhe_2b_pot", n1, n2; 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 diff --git a/src/models/heisenberg.jl b/src/models/heisenberg.jl index 2388b99..0d095f7 100644 --- a/src/models/heisenberg.jl +++ b/src/models/heisenberg.jl @@ -25,22 +25,3 @@ function reference(::Model"heisenberg", ::Observable"energy"; N=∞) correction = 1 + 0.375 / log(N)^3 return (E∞ - Eᶠⁱⁿⁱᵗᵉ * correction) / (2N) end - -# function ITensors.OpSum(::Model"heisenberg", n1, n2) -# opsum = OpSum() -# opsum += 0.5, "S+", n1, "S-", n2 -# opsum += 0.5, "S-", n1, "S+", n2 -# opsum += "Sz", n1, "Sz", n2 -# return opsum -# end -# -# 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 diff --git a/src/models/hubbard.jl b/src/models/hubbard.jl index 3ed9cd0..f9636f4 100644 --- a/src/models/hubbard.jl +++ b/src/models/hubbard.jl @@ -42,43 +42,3 @@ function reference(::Model"hubbard", ::Observable"energy"; U, Npoints=50) g(i) = -2 * (2π) / Npoints * cos(-π + (i - 1) / Npoints * 2π) * v[i] return sum(g, 1:Npoints) end - -## function ITensors.OpSum(::Model"hubbard", n1, n2; t, U, V) -## 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 -## 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 diff --git a/src/models/ising.jl b/src/models/ising.jl index c8c26de..20897ac 100644 --- a/src/models/ising.jl +++ b/src/models/ising.jl @@ -28,39 +28,3 @@ function reference(::Model"ising_classical", ::Observable"critical_inverse_tempe # log1p(x) = log(1 + x) return log1p(sqrt(2)) / 2 end - -# 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) -# 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 -# 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 diff --git a/src/models/ising_extended.jl b/src/models/ising_extended.jl index de9eac2..e7d4ddc 100644 --- a/src/models/ising_extended.jl +++ b/src/models/ising_extended.jl @@ -13,32 +13,3 @@ function reference(::Model"ising_extended", ::Observable"energy"; J=1.0, h=1.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), -π, π) end - -# -# Deprecated -# - -# nrange(::Model"ising_extended") = 3 - -# function ITensors.ITensor(::Model"ising_extended", s1::Index, s2::Index, s3::Index; 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 - -# # For finite Hamiltonian with open boundary conditions -# # H = -J Σⱼ XⱼXⱼ₊₁ - h Σⱼ Zⱼ - J₂ Σⱼ XⱼZⱼ₊₁Xⱼ₊₂ -# function opsum_finite(::Model"ising_extended", n; J=1.0, h=1.0, J₂=0.0) -# # TODO: define using opsum_infinite, filtering sites -# # that extend over the boundary. -# opsum = OpSum() -# for j in 1:n -# if j ≤ (n - 1) -# opsum += -J, "X", j, "X", j + 1 -# end -# if j ≤ (n - 2) -# opsum += -J₂, "X", j, "Z", j + 1, "X", j + 2 -# end -# opsum -= h, "Z", j -# end -# return opsum -# end diff --git a/src/models/models.jl b/src/models/models.jl index 1acc7a7..54a86fc 100644 --- a/src/models/models.jl +++ b/src/models/models.jl @@ -204,33 +204,3 @@ function ITensors.ITensor(model::Model, s::CelledVector, n::Int64; kwargs...) # Deprecated version # return contract(MPO(model, s, n; kwargs...)) end - -# -# Deprecated -# - -#ITensorInfiniteMPS.nrange(model::Model) = 2; #required to keep everything compatible with the current implementation for 2 band models - -# 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 -# end - -# # MPO building version -# function ITensors.MPO(model::Model, s::CelledVector; kwargs...) -# opsum = OpSum(model, nsites(s); kwargs...) -# -# # @show nrange(ops) -# # @show n:(n + nrange(opsum) - 1) -# # @show opsum -# -# # return MPO(opsum, [s[x] for x in n:(n + nrange(opsum) - 1)]) #modification to allow for more than two sites per term in the Hamiltonians -# return MPO(opsum, [s[x] for x in ITensors.sites(opsum)]) #modification to allow for more than two sites per term in the Hamiltonians -# end - -# function InfiniteSum{MPO}(model::Model, s::CelledVector; kwargs...) -# N = length(s) -# # mpos = [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 -# mpo = MPO(model, s; kwargs...) -# return InfiniteSum{MPO}([mpo], translator(s)) -# end diff --git a/src/models/xx.jl b/src/models/xx.jl index 099d287..0141024 100644 --- a/src/models/xx.jl +++ b/src/models/xx.jl @@ -12,34 +12,3 @@ function reference(::Model"xx", ::Observable"energy"; N=∞) λ(k) = cos(k * π / (N + 1)) return 4 * sum(k -> min(λ(k), 0.0), 1:N) / N end - -## function ITensors.MPO(::Model"xx", s) -## N = length(s) -## 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])) -## end diff --git a/src/models/xx_extended.jl b/src/models/xx_extended.jl index f71c24b..9e6939e 100644 --- a/src/models/xx_extended.jl +++ b/src/models/xx_extended.jl @@ -36,49 +36,3 @@ function ITensorInfiniteMPS.reference( temp = sort(λ.(0:(N - 1))) return sum(temp[1:round(Int64, filling * N)]) / N - h * (2 * filling - 1) end - -## function ITensors.OpSum(::Model"xx_extended", n1, n2; 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 -## 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 diff --git a/test/test_vumpsmpo_fqhe.jl b/test/test_vumpsmpo_fqhe.jl index 1f227b8..76a9d22 100644 --- a/test/test_vumpsmpo_fqhe.jl +++ b/test/test_vumpsmpo_fqhe.jl @@ -56,13 +56,7 @@ end nsites in [6], time_step in [-Inf] - vumps_kwargs = (; - multisite_update_alg, - tol, - maxiter, - outputlevel=0, - time_step, - ) + 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) @@ -80,7 +74,13 @@ end end s = infsiteinds( - "FermionK", nsites; initstate, translator=fermion_momentum_translator, p, q, conserve_momentum + "FermionK", + nsites; + initstate, + translator=fermion_momentum_translator, + p, + q, + conserve_momentum, ) ψ = InfMPS(s, initstate) From d30675d75b11ea5f2d79e77d5f8bd5f6414b2983 Mon Sep 17 00:00:00 2001 From: Matthew Fishman Date: Fri, 24 Jun 2022 10:23:32 -0400 Subject: [PATCH 23/32] Default to V=0 for Hubbard model --- src/models/hubbard.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/models/hubbard.jl b/src/models/hubbard.jl index f9636f4..f2b1dd1 100644 --- a/src/models/hubbard.jl +++ b/src/models/hubbard.jl @@ -1,4 +1,4 @@ -function unit_cell_terms(::Model"hubbard"; t, U, V) +function unit_cell_terms(::Model"hubbard"; t, U, V=0.0) opsum = OpSum() opsum += -t, "Cdagup", 1, "Cup", 2 opsum += -t, "Cdagup", 2, "Cup", 1 From baea0061d72d4add43f898f9c5d6fc4e9943cef6 Mon Sep 17 00:00:00 2001 From: Matthew Fishman Date: Mon, 27 Jun 2022 11:59:16 -0400 Subject: [PATCH 24/32] Start updating development examples --- .../vumps/transfer_matrix_spectrum_arpack.jl | 19 ++++++------------- .../vumps/vumps_ising_noncontiguous.jl | 7 +------ examples/development/vumps/vumps_mpo.jl | 19 +++++++------------ 3 files changed, 14 insertions(+), 31 deletions(-) diff --git a/examples/development/vumps/transfer_matrix_spectrum_arpack.jl b/examples/development/vumps/transfer_matrix_spectrum_arpack.jl index 2759295..b7c65da 100644 --- a/examples/development/vumps/transfer_matrix_spectrum_arpack.jl +++ b/examples/development/vumps/transfer_matrix_spectrum_arpack.jl @@ -13,7 +13,7 @@ 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 = true +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 @@ -23,19 +23,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 @@ -135,7 +126,9 @@ 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⃗ = [(translatecelltags(v⃗ᴸ[n], 1) * v⃗ᴿ[n])[] for n in 1:neigs] diff --git a/examples/development/vumps/vumps_ising_noncontiguous.jl b/examples/development/vumps/vumps_ising_noncontiguous.jl index e830d1d..032a9a2 100644 --- a/examples/development/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"); diff --git a/examples/development/vumps/vumps_mpo.jl b/examples/development/vumps/vumps_mpo.jl index 652a05c..acfa204 100644 --- a/examples/development/vumps/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])) From 742a5aa8e0cf1085d75ee1f171a81ae3a7b2e81c Mon Sep 17 00:00:00 2001 From: Matthew Fishman Date: Fri, 15 Jul 2022 18:28:07 -0400 Subject: [PATCH 25/32] Update directsum for latest ITensors interface --- Project.toml | 2 +- examples/vumps/transfer_matrix_spectrum.jl | 2 +- examples/vumps/vumps_heisenberg.jl | 2 +- examples/vumps/vumps_hubbard_extended.jl | 2 +- examples/vumps/vumps_ising.jl | 2 +- examples/vumps/vumps_ising_extended.jl | 2 +- src/subspace_expansion.jl | 4 ++-- 7 files changed, 8 insertions(+), 8 deletions(-) diff --git a/Project.toml b/Project.toml index d39631a..e9ed03a 100644 --- a/Project.toml +++ b/Project.toml @@ -18,7 +18,7 @@ Requires = "ae029012-a4dd-5104-9daa-d747884805df" [compat] Compat = "3, 4" HDF5 = "0.15, 0.16" -ITensors = "0.3" +ITensors = "0.3.18" Infinities = "0.1" IterTools = "1" KrylovKit = "0.5" diff --git a/examples/vumps/transfer_matrix_spectrum.jl b/examples/vumps/transfer_matrix_spectrum.jl index f383906..f70725a 100644 --- a/examples/vumps/transfer_matrix_spectrum.jl +++ b/examples/vumps/transfer_matrix_spectrum.jl @@ -1,7 +1,7 @@ using ITensors using ITensorInfiniteMPS -include(joinpath(@__DIR__, "src", "vumps_subspace_expansion.jl")) +include(joinpath(pkgdir(ITensorInfiniteMPS), "examples", "src", "vumps_subspace_expansion.jl")) ############################################################################## # VUMPS parameters diff --git a/examples/vumps/vumps_heisenberg.jl b/examples/vumps/vumps_heisenberg.jl index 8133cd9..45a1cd6 100644 --- a/examples/vumps/vumps_heisenberg.jl +++ b/examples/vumps/vumps_heisenberg.jl @@ -1,7 +1,7 @@ using ITensors using ITensorInfiniteMPS -include(joinpath(@__DIR__, "src", "vumps_subspace_expansion.jl")) +include(joinpath(pkgdir(ITensorInfiniteMPS), "examples", "src", "vumps_subspace_expansion.jl")) ############################################################################## # VUMPS parameters diff --git a/examples/vumps/vumps_hubbard_extended.jl b/examples/vumps/vumps_hubbard_extended.jl index 880a5e0..17aa357 100644 --- a/examples/vumps/vumps_hubbard_extended.jl +++ b/examples/vumps/vumps_hubbard_extended.jl @@ -1,7 +1,7 @@ using ITensors using ITensorInfiniteMPS -include(joinpath(@__DIR__, "src", "vumps_subspace_expansion.jl")) +include(joinpath(pkgdir(ITensorInfiniteMPS), "examples", "src", "vumps_subspace_expansion.jl")) ############################################################################## # VUMPS parameters diff --git a/examples/vumps/vumps_ising.jl b/examples/vumps/vumps_ising.jl index 31165a2..8efaf12 100644 --- a/examples/vumps/vumps_ising.jl +++ b/examples/vumps/vumps_ising.jl @@ -1,7 +1,7 @@ using ITensors using ITensorInfiniteMPS -include(joinpath(@__DIR__, "src", "vumps_subspace_expansion.jl")) +include(joinpath(pkgdir(ITensorInfiniteMPS), "examples", "src", "vumps_subspace_expansion.jl")) ############################################################################## # VUMPS/TDVP parameters diff --git a/examples/vumps/vumps_ising_extended.jl b/examples/vumps/vumps_ising_extended.jl index 7119e22..20dcf91 100644 --- a/examples/vumps/vumps_ising_extended.jl +++ b/examples/vumps/vumps_ising_extended.jl @@ -1,7 +1,7 @@ using ITensors using ITensorInfiniteMPS -include(joinpath(@__DIR__, "src", "vumps_subspace_expansion.jl")) +include(joinpath(pkgdir(ITensorInfiniteMPS), "examples", "src", "vumps_subspace_expansion.jl")) ############################################################################## # VUMPS parameters diff --git a/src/subspace_expansion.jl b/src/subspace_expansion.jl index a9786d1..1fb58ab 100644 --- a/src/subspace_expansion.jl +++ b/src/subspace_expansion.jl @@ -213,10 +213,10 @@ 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)...) From 275ea00f4558b5ee2d2feee4af572926e25b585a Mon Sep 17 00:00:00 2001 From: Matt Fishman Date: Fri, 15 Jul 2022 19:04:38 -0400 Subject: [PATCH 26/32] Update examples/vumps/transfer_matrix_spectrum.jl Co-authored-by: github-actions[bot] <41898282+github-actions[bot]@users.noreply.github.com> --- examples/vumps/transfer_matrix_spectrum.jl | 4 +++- 1 file changed, 3 insertions(+), 1 deletion(-) diff --git a/examples/vumps/transfer_matrix_spectrum.jl b/examples/vumps/transfer_matrix_spectrum.jl index f70725a..f7d326d 100644 --- a/examples/vumps/transfer_matrix_spectrum.jl +++ b/examples/vumps/transfer_matrix_spectrum.jl @@ -1,7 +1,9 @@ using ITensors using ITensorInfiniteMPS -include(joinpath(pkgdir(ITensorInfiniteMPS), "examples", "src", "vumps_subspace_expansion.jl")) +include( + joinpath(pkgdir(ITensorInfiniteMPS), "examples", "src", "vumps_subspace_expansion.jl") +) ############################################################################## # VUMPS parameters From 1ae9ed43f2518d67349a0fbc8ac942cf0a3d37dd Mon Sep 17 00:00:00 2001 From: Matt Fishman Date: Fri, 15 Jul 2022 19:04:46 -0400 Subject: [PATCH 27/32] Update examples/vumps/vumps_heisenberg.jl Co-authored-by: github-actions[bot] <41898282+github-actions[bot]@users.noreply.github.com> --- examples/vumps/vumps_heisenberg.jl | 4 +++- 1 file changed, 3 insertions(+), 1 deletion(-) diff --git a/examples/vumps/vumps_heisenberg.jl b/examples/vumps/vumps_heisenberg.jl index 45a1cd6..c0a807d 100644 --- a/examples/vumps/vumps_heisenberg.jl +++ b/examples/vumps/vumps_heisenberg.jl @@ -1,7 +1,9 @@ using ITensors using ITensorInfiniteMPS -include(joinpath(pkgdir(ITensorInfiniteMPS), "examples", "src", "vumps_subspace_expansion.jl")) +include( + joinpath(pkgdir(ITensorInfiniteMPS), "examples", "src", "vumps_subspace_expansion.jl") +) ############################################################################## # VUMPS parameters From b03215c4334bf2d201c1eb5866d7bb7439a1968b Mon Sep 17 00:00:00 2001 From: Matthew Fishman Date: Fri, 15 Jul 2022 19:08:05 -0400 Subject: [PATCH 28/32] Format --- examples/vumps/transfer_matrix_spectrum.jl | 4 +++- examples/vumps/vumps_heisenberg.jl | 4 +++- examples/vumps/vumps_hubbard_extended.jl | 4 +++- examples/vumps/vumps_ising.jl | 4 +++- examples/vumps/vumps_ising_extended.jl | 4 +++- src/subspace_expansion.jl | 8 ++++++-- 6 files changed, 21 insertions(+), 7 deletions(-) diff --git a/examples/vumps/transfer_matrix_spectrum.jl b/examples/vumps/transfer_matrix_spectrum.jl index f70725a..f7d326d 100644 --- a/examples/vumps/transfer_matrix_spectrum.jl +++ b/examples/vumps/transfer_matrix_spectrum.jl @@ -1,7 +1,9 @@ using ITensors using ITensorInfiniteMPS -include(joinpath(pkgdir(ITensorInfiniteMPS), "examples", "src", "vumps_subspace_expansion.jl")) +include( + joinpath(pkgdir(ITensorInfiniteMPS), "examples", "src", "vumps_subspace_expansion.jl") +) ############################################################################## # VUMPS parameters diff --git a/examples/vumps/vumps_heisenberg.jl b/examples/vumps/vumps_heisenberg.jl index 45a1cd6..c0a807d 100644 --- a/examples/vumps/vumps_heisenberg.jl +++ b/examples/vumps/vumps_heisenberg.jl @@ -1,7 +1,9 @@ using ITensors using ITensorInfiniteMPS -include(joinpath(pkgdir(ITensorInfiniteMPS), "examples", "src", "vumps_subspace_expansion.jl")) +include( + joinpath(pkgdir(ITensorInfiniteMPS), "examples", "src", "vumps_subspace_expansion.jl") +) ############################################################################## # VUMPS parameters diff --git a/examples/vumps/vumps_hubbard_extended.jl b/examples/vumps/vumps_hubbard_extended.jl index 17aa357..03f4f06 100644 --- a/examples/vumps/vumps_hubbard_extended.jl +++ b/examples/vumps/vumps_hubbard_extended.jl @@ -1,7 +1,9 @@ using ITensors using ITensorInfiniteMPS -include(joinpath(pkgdir(ITensorInfiniteMPS), "examples", "src", "vumps_subspace_expansion.jl")) +include( + joinpath(pkgdir(ITensorInfiniteMPS), "examples", "src", "vumps_subspace_expansion.jl") +) ############################################################################## # VUMPS parameters diff --git a/examples/vumps/vumps_ising.jl b/examples/vumps/vumps_ising.jl index 8efaf12..bc5a7d4 100644 --- a/examples/vumps/vumps_ising.jl +++ b/examples/vumps/vumps_ising.jl @@ -1,7 +1,9 @@ using ITensors using ITensorInfiniteMPS -include(joinpath(pkgdir(ITensorInfiniteMPS), "examples", "src", "vumps_subspace_expansion.jl")) +include( + joinpath(pkgdir(ITensorInfiniteMPS), "examples", "src", "vumps_subspace_expansion.jl") +) ############################################################################## # VUMPS/TDVP parameters diff --git a/examples/vumps/vumps_ising_extended.jl b/examples/vumps/vumps_ising_extended.jl index 20dcf91..5f5457b 100644 --- a/examples/vumps/vumps_ising_extended.jl +++ b/examples/vumps/vumps_ising_extended.jl @@ -1,7 +1,9 @@ using ITensors using ITensorInfiniteMPS -include(joinpath(pkgdir(ITensorInfiniteMPS), "examples", "src", "vumps_subspace_expansion.jl")) +include( + joinpath(pkgdir(ITensorInfiniteMPS), "examples", "src", "vumps_subspace_expansion.jl") +) ############################################################################## # VUMPS parameters diff --git a/src/subspace_expansion.jl b/src/subspace_expansion.jl index 1fb58ab..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] => uniqueinds(ψ.AL[n1], NL), dag(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] => uniqueinds(ψ.AR[n2], NR), dag(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)...) From ab7f0cd98b3e321afde6cac6a0cf242520f66c78 Mon Sep 17 00:00:00 2001 From: Matthew Fishman Date: Tue, 16 Aug 2022 11:26:34 -0400 Subject: [PATCH 29/32] Fix examples. Simplify overload requirements for new models --- examples/vumps/transfer_matrix_spectrum.jl | 2 +- examples/vumps/vumps_heisenberg.jl | 2 +- examples/vumps/vumps_hubbard_extended.jl | 2 +- examples/vumps/vumps_ising.jl | 2 +- examples/vumps/vumps_ising_extended.jl | 2 +- src/models/fqhe13.jl | 3 +-- src/models/heisenberg.jl | 2 +- src/models/ising.jl | 2 +- src/models/ising_extended.jl | 5 ++--- src/models/models.jl | 26 +++++++++++++++++++++- src/models/xx.jl | 4 ++-- 11 files changed, 37 insertions(+), 15 deletions(-) diff --git a/examples/vumps/transfer_matrix_spectrum.jl b/examples/vumps/transfer_matrix_spectrum.jl index f7d326d..e62e271 100644 --- a/examples/vumps/transfer_matrix_spectrum.jl +++ b/examples/vumps/transfer_matrix_spectrum.jl @@ -2,7 +2,7 @@ using ITensors using ITensorInfiniteMPS include( - joinpath(pkgdir(ITensorInfiniteMPS), "examples", "src", "vumps_subspace_expansion.jl") + joinpath(pkgdir(ITensorInfiniteMPS), "examples", "vumps", "src", "vumps_subspace_expansion.jl") ) ############################################################################## diff --git a/examples/vumps/vumps_heisenberg.jl b/examples/vumps/vumps_heisenberg.jl index c0a807d..c3da2e7 100644 --- a/examples/vumps/vumps_heisenberg.jl +++ b/examples/vumps/vumps_heisenberg.jl @@ -2,7 +2,7 @@ using ITensors using ITensorInfiniteMPS include( - joinpath(pkgdir(ITensorInfiniteMPS), "examples", "src", "vumps_subspace_expansion.jl") + joinpath(pkgdir(ITensorInfiniteMPS), "examples", "vumps", "src", "vumps_subspace_expansion.jl") ) ############################################################################## diff --git a/examples/vumps/vumps_hubbard_extended.jl b/examples/vumps/vumps_hubbard_extended.jl index 03f4f06..e2754be 100644 --- a/examples/vumps/vumps_hubbard_extended.jl +++ b/examples/vumps/vumps_hubbard_extended.jl @@ -2,7 +2,7 @@ using ITensors using ITensorInfiniteMPS include( - joinpath(pkgdir(ITensorInfiniteMPS), "examples", "src", "vumps_subspace_expansion.jl") + joinpath(pkgdir(ITensorInfiniteMPS), "examples", "vumps", "src", "vumps_subspace_expansion.jl") ) ############################################################################## diff --git a/examples/vumps/vumps_ising.jl b/examples/vumps/vumps_ising.jl index bc5a7d4..29facb0 100644 --- a/examples/vumps/vumps_ising.jl +++ b/examples/vumps/vumps_ising.jl @@ -2,7 +2,7 @@ using ITensors using ITensorInfiniteMPS include( - joinpath(pkgdir(ITensorInfiniteMPS), "examples", "src", "vumps_subspace_expansion.jl") + joinpath(pkgdir(ITensorInfiniteMPS), "examples", "vumps", "src", "vumps_subspace_expansion.jl") ) ############################################################################## diff --git a/examples/vumps/vumps_ising_extended.jl b/examples/vumps/vumps_ising_extended.jl index 5f5457b..6de5b69 100644 --- a/examples/vumps/vumps_ising_extended.jl +++ b/examples/vumps/vumps_ising_extended.jl @@ -2,7 +2,7 @@ using ITensors using ITensorInfiniteMPS include( - joinpath(pkgdir(ITensorInfiniteMPS), "examples", "src", "vumps_subspace_expansion.jl") + joinpath(pkgdir(ITensorInfiniteMPS), "examples", "vumps", "src", "vumps_subspace_expansion.jl") ) ############################################################################## diff --git a/src/models/fqhe13.jl b/src/models/fqhe13.jl index 851aefb..40d87e6 100644 --- a/src/models/fqhe13.jl +++ b/src/models/fqhe13.jl @@ -6,8 +6,7 @@ function unit_cell_terms( 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)] + return generate_Hamiltonian(opt) end #Please contact Loic Herviou before using this part of the code for production diff --git a/src/models/heisenberg.jl b/src/models/heisenberg.jl index 0d095f7..2d9a338 100644 --- a/src/models/heisenberg.jl +++ b/src/models/heisenberg.jl @@ -4,7 +4,7 @@ function unit_cell_terms(::Model"heisenberg") opsum += 0.5, "S+", 1, "S-", 2 opsum += 0.5, "S-", 1, "S+", 2 opsum += "Sz", 1, "Sz", 2 - return [opsum] + return opsum end """ diff --git a/src/models/ising.jl b/src/models/ising.jl index 20897ac..960638a 100644 --- a/src/models/ising.jl +++ b/src/models/ising.jl @@ -3,7 +3,7 @@ function unit_cell_terms(::Model"ising"; J=1.0, h=1.0) opsum = OpSum() opsum += -J, "X", 1, "X", 2 opsum += -h, "Z", 1 - return [opsum] + return opsum end # P. Pfeuty, The one-dimensional Ising model with a transverse field, Annals of Physics 57, p. 79 (1970) diff --git a/src/models/ising_extended.jl b/src/models/ising_extended.jl index e7d4ddc..e234efb 100644 --- a/src/models/ising_extended.jl +++ b/src/models/ising_extended.jl @@ -1,12 +1,11 @@ -# A vector of the terms associated with each site of the -# unit cell. +# The terms of the Hamiltonian in the first unit cell. # H = -J Σⱼ XⱼXⱼ₊₁ - h Σⱼ Zⱼ - J₂ Σⱼ XⱼZⱼ₊₁Xⱼ₊₂ function unit_cell_terms(::Model"ising_extended"; J=1.0, h=1.0, J₂=0.0) opsum = OpSum() opsum += -J, "X", 1, "X", 2 opsum += -J₂, "X", 1, "Z", 2, "X", 3 opsum += -h, "Z", 1 - return [opsum] + return opsum end function reference(::Model"ising_extended", ::Observable"energy"; J=1.0, h=1.0, J₂=0.0) diff --git a/src/models/models.jl b/src/models/models.jl index 54a86fc..9bb2d83 100644 --- a/src/models/models.jl +++ b/src/models/models.jl @@ -119,8 +119,32 @@ end translatecell(translator, opsum::OpSum, n::Integer) = translator(opsum, n) +function split_by_first_site(opsum::OpSum) + cell_terms = Ops.terms(opsum) + # The first site of each term of the `opsum` + term_first_site = [minimum(Ops.sites(o)) for o in cell_terms] + nsites_in_unit_cell = maximum(term_first_site) + opsums = [OpSum() for _ in 1:nsites_in_unit_cell] + for j in eachindex(cell_terms) + siteⱼ = term_first_site[j] + termⱼ = cell_terms[j] + opsums[siteⱼ] += termⱼ + end + return opsums +end + function infinite_terms(model::Model; kwargs...) - opsum_cell = unit_cell_terms(model; kwargs...) + # An `OpSum` storing all of the terms in the + # first unit cell. + # TODO: Allow specifying the unit cell size + # explicitly. + opsum = unit_cell_terms(model; 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 = split_by_first_site(opsum) nsites = length(opsum_cell) function _shift_cell(opsum::OpSum, cell::Int) return shift_sites(opsum, nsites * cell) diff --git a/src/models/xx.jl b/src/models/xx.jl index 0141024..061a49d 100644 --- a/src/models/xx.jl +++ b/src/models/xx.jl @@ -1,8 +1,8 @@ # H = Σⱼ XⱼXⱼ₊₁ + YⱼYⱼ₊₁ function unit_cell_terms(::Model"xx") os = OpSum() - os += 1, "X", 1, "X", 2 - os += 1, "Y", 1, "Y", 2 + os += "X", 1, "X", 2 + os += "Y", 1, "Y", 2 return os end From 6fa08e185669bc551d6a216e33c47bba78eed33d Mon Sep 17 00:00:00 2001 From: Matthew Fishman Date: Tue, 16 Aug 2022 11:30:03 -0400 Subject: [PATCH 30/32] Format --- examples/vumps/transfer_matrix_spectrum.jl | 4 +++- examples/vumps/vumps_heisenberg.jl | 4 +++- examples/vumps/vumps_hubbard_extended.jl | 4 +++- examples/vumps/vumps_ising.jl | 4 +++- examples/vumps/vumps_ising_extended.jl | 4 +++- 5 files changed, 15 insertions(+), 5 deletions(-) diff --git a/examples/vumps/transfer_matrix_spectrum.jl b/examples/vumps/transfer_matrix_spectrum.jl index e62e271..ed6a4e4 100644 --- a/examples/vumps/transfer_matrix_spectrum.jl +++ b/examples/vumps/transfer_matrix_spectrum.jl @@ -2,7 +2,9 @@ using ITensors using ITensorInfiniteMPS include( - joinpath(pkgdir(ITensorInfiniteMPS), "examples", "vumps", "src", "vumps_subspace_expansion.jl") + joinpath( + pkgdir(ITensorInfiniteMPS), "examples", "vumps", "src", "vumps_subspace_expansion.jl" + ), ) ############################################################################## diff --git a/examples/vumps/vumps_heisenberg.jl b/examples/vumps/vumps_heisenberg.jl index c3da2e7..d2842a6 100644 --- a/examples/vumps/vumps_heisenberg.jl +++ b/examples/vumps/vumps_heisenberg.jl @@ -2,7 +2,9 @@ using ITensors using ITensorInfiniteMPS include( - joinpath(pkgdir(ITensorInfiniteMPS), "examples", "vumps", "src", "vumps_subspace_expansion.jl") + joinpath( + pkgdir(ITensorInfiniteMPS), "examples", "vumps", "src", "vumps_subspace_expansion.jl" + ), ) ############################################################################## diff --git a/examples/vumps/vumps_hubbard_extended.jl b/examples/vumps/vumps_hubbard_extended.jl index e2754be..cba2dd9 100644 --- a/examples/vumps/vumps_hubbard_extended.jl +++ b/examples/vumps/vumps_hubbard_extended.jl @@ -2,7 +2,9 @@ using ITensors using ITensorInfiniteMPS include( - joinpath(pkgdir(ITensorInfiniteMPS), "examples", "vumps", "src", "vumps_subspace_expansion.jl") + joinpath( + pkgdir(ITensorInfiniteMPS), "examples", "vumps", "src", "vumps_subspace_expansion.jl" + ), ) ############################################################################## diff --git a/examples/vumps/vumps_ising.jl b/examples/vumps/vumps_ising.jl index 29facb0..6085045 100644 --- a/examples/vumps/vumps_ising.jl +++ b/examples/vumps/vumps_ising.jl @@ -2,7 +2,9 @@ using ITensors using ITensorInfiniteMPS include( - joinpath(pkgdir(ITensorInfiniteMPS), "examples", "vumps", "src", "vumps_subspace_expansion.jl") + joinpath( + pkgdir(ITensorInfiniteMPS), "examples", "vumps", "src", "vumps_subspace_expansion.jl" + ), ) ############################################################################## diff --git a/examples/vumps/vumps_ising_extended.jl b/examples/vumps/vumps_ising_extended.jl index 6de5b69..ed4af6f 100644 --- a/examples/vumps/vumps_ising_extended.jl +++ b/examples/vumps/vumps_ising_extended.jl @@ -2,7 +2,9 @@ using ITensors using ITensorInfiniteMPS include( - joinpath(pkgdir(ITensorInfiniteMPS), "examples", "vumps", "src", "vumps_subspace_expansion.jl") + joinpath( + pkgdir(ITensorInfiniteMPS), "examples", "vumps", "src", "vumps_subspace_expansion.jl" + ), ) ############################################################################## From a6b7279eee9d87f69b9aed1a1424291d495539ee Mon Sep 17 00:00:00 2001 From: Matthew Fishman Date: Tue, 16 Aug 2022 12:26:44 -0400 Subject: [PATCH 31/32] Fix test --- src/ITensorInfiniteMPS.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/ITensorInfiniteMPS.jl b/src/ITensorInfiniteMPS.jl index 738bc0e..ff3b9d9 100644 --- a/src/ITensorInfiniteMPS.jl +++ b/src/ITensorInfiniteMPS.jl @@ -22,7 +22,7 @@ using KrylovKit: eigsolve, linsolve, exponentiate import Base: getindex, length, setindex!, +, -, * -import ITensors: AbstractMPS +import ITensors: AbstractMPS, ⊕ include("ITensors.jl") include("ITensorNetworks.jl") From 11ce0f8eebb7c5580ffc3ea06c3f5b49805d1787 Mon Sep 17 00:00:00 2001 From: Matthew Fishman Date: Tue, 16 Aug 2022 13:11:25 -0400 Subject: [PATCH 32/32] Clean up code a bit with SplitApplyCombine --- Project.toml | 2 ++ src/ITensorInfiniteMPS.jl | 4 ++++ src/models/hubbard.jl | 2 +- src/models/models.jl | 32 +++++++++++++++----------------- 4 files changed, 22 insertions(+), 18 deletions(-) diff --git a/Project.toml b/Project.toml index e9ed03a..7d2104d 100644 --- a/Project.toml +++ b/Project.toml @@ -14,6 +14,7 @@ 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" @@ -25,6 +26,7 @@ KrylovKit = "0.5" OffsetArrays = "1" QuadGK = "2" Requires = "1" +SplitApplyCombine = "1.2.2" julia = "1.6" [extras] diff --git a/src/ITensorInfiniteMPS.jl b/src/ITensorInfiniteMPS.jl index ff3b9d9..996cfbb 100644 --- a/src/ITensorInfiniteMPS.jl +++ b/src/ITensorInfiniteMPS.jl @@ -16,6 +16,10 @@ 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 diff --git a/src/models/hubbard.jl b/src/models/hubbard.jl index f2b1dd1..77960e1 100644 --- a/src/models/hubbard.jl +++ b/src/models/hubbard.jl @@ -6,7 +6,7 @@ function unit_cell_terms(::Model"hubbard"; t, U, V=0.0) opsum += -t, "Cdagdn", 2, "Cdn", 1 opsum += U, "Nupdn", 1 opsum += V, "Ntot", 1, "Ntot", 2 - return [opsum] + return opsum end """ diff --git a/src/models/models.jl b/src/models/models.jl index 9bb2d83..75d5085 100644 --- a/src/models/models.jl +++ b/src/models/models.jl @@ -119,33 +119,31 @@ end translatecell(translator, opsum::OpSum, n::Integer) = translator(opsum, n) -function split_by_first_site(opsum::OpSum) - cell_terms = Ops.terms(opsum) - # The first site of each term of the `opsum` - term_first_site = [minimum(Ops.sites(o)) for o in cell_terms] - nsites_in_unit_cell = maximum(term_first_site) - opsums = [OpSum() for _ in 1:nsites_in_unit_cell] - for j in eachindex(cell_terms) - siteⱼ = term_first_site[j] - termⱼ = cell_terms[j] - opsums[siteⱼ] += termⱼ - end - return opsums -end - 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. - opsum = unit_cell_terms(model; kwargs...) + 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 = split_by_first_site(opsum) - nsites = length(opsum_cell) + 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