diff --git a/Project.toml b/Project.toml
index b8594b64..daec0556 100644
--- a/Project.toml
+++ b/Project.toml
@@ -1,7 +1,7 @@
 name = "ITensorNetworks"
 uuid = "2919e153-833c-4bdc-8836-1ea460a35fc7"
 authors = ["Matthew Fishman <mfishman@flatironinstitute.org>, Joseph Tindall <jtindall@flatironinstitute.org> and contributors"]
-version = "0.12.2"
+version = "0.13.0"
 
 [deps]
 AbstractTrees = "1520ce14-60c1-5f80-bbc7-55ef81b5835c"
@@ -13,7 +13,6 @@ Dictionaries = "85a47980-9c8c-11e8-2b9f-f7ca1fa99fb4"
 Distributions = "31c24e10-a181-5473-b8eb-7969acd0382f"
 DocStringExtensions = "ffbed154-4ef7-542d-bbb7-c09d3a79fcae"
 Graphs = "86223c79-3864-5bf0-83f7-82e725a168b6"
-ITensorMPS = "0d1a4710-d33b-49a5-8f18-73bdf49b47e2"
 ITensors = "9136182c-28ba-11e9-034c-db9fb085ebd5"
 IsApprox = "28f27b66-4bd8-47e7-9110-e2746eb8bed7"
 IterTools = "c8e1da08-722c-5040-9ed9-7db0dc04731e"
@@ -63,7 +62,6 @@ DocStringExtensions = "0.9"
 EinExprs = "0.6.4"
 Graphs = "1.8"
 GraphsFlows = "0.1.1"
-ITensorMPS = "0.3"
 ITensors = "0.7, 0.8"
 IsApprox = "0.1, 1, 2"
 IterTools = "1.4.0"
diff --git a/docs/Project.toml b/docs/Project.toml
index 41ba6ca9..55b577dd 100644
--- a/docs/Project.toml
+++ b/docs/Project.toml
@@ -5,5 +5,5 @@ Literate = "98b081ad-f1c9-55d3-8b20-4c87d4299306"
 
 [compat]
 Documenter = "1.10.0"
-ITensorNetworks = "0.12.0"
+ITensorNetworks = "0.13.0"
 Literate = "2.20.1"
diff --git a/examples/Project.toml b/examples/Project.toml
index c4ae6bbb..7c06c098 100644
--- a/examples/Project.toml
+++ b/examples/Project.toml
@@ -2,4 +2,4 @@
 ITensorNetworks = "2919e153-833c-4bdc-8836-1ea460a35fc7"
 
 [compat]
-ITensorNetworks = "0.12.0"
+ITensorNetworks = "0.13.0"
diff --git a/src/ITensorNetworks.jl b/src/ITensorNetworks.jl
index 2b16b030..56ab768f 100644
--- a/src/ITensorNetworks.jl
+++ b/src/ITensorNetworks.jl
@@ -23,7 +23,6 @@ include("contract_approx/partition.jl")
 include("contract_approx/binary_tree_partition.jl")
 include("contract.jl")
 include("specialitensornetworks.jl")
-include("boundarymps.jl")
 include("partitioneditensornetwork.jl")
 include("edge_sequences.jl")
 include("caches/abstractbeliefpropagationcache.jl")
diff --git a/src/abstractitensornetwork.jl b/src/abstractitensornetwork.jl
index 6c80e36e..e35fe063 100644
--- a/src/abstractitensornetwork.jl
+++ b/src/abstractitensornetwork.jl
@@ -38,7 +38,6 @@ using ITensors:
   settags,
   sim,
   swaptags
-using ITensorMPS: ITensorMPS, add, linkdim, linkinds, siteinds
 using .ITensorsExtensions: ITensorsExtensions, indtype, promote_indtype
 using LinearAlgebra: LinearAlgebra, factorize
 using MacroTools: @capture
@@ -255,7 +254,7 @@ indsnetwork(tn::AbstractITensorNetwork) = IndsNetwork(tn)
 
 # TODO: Output a `VertexDataGraph`? Unfortunately
 # `IndsNetwork` doesn't allow iterating over vertex data.
-function ITensorMPS.siteinds(tn::AbstractITensorNetwork)
+function siteinds(tn::AbstractITensorNetwork)
   is = IndsNetwork(underlying_graph(tn))
   for v in vertices(tn)
     is[v] = uniqueinds(tn, v)
@@ -268,7 +267,7 @@ function flatten_siteinds(tn::AbstractITensorNetwork)
   return identity.(flatten(map(v -> siteinds(tn, v), vertices(tn))))
 end
 
-function ITensorMPS.linkinds(tn::AbstractITensorNetwork)
+function linkinds(tn::AbstractITensorNetwork)
   is = IndsNetwork(underlying_graph(tn))
   for e in edges(tn)
     is[e] = commoninds(tn, e)
@@ -302,7 +301,11 @@ function ITensors.uniqueinds(tn::AbstractITensorNetwork, edge::Pair)
   return uniqueinds(tn, edgetype(tn)(edge))
 end
 
-function ITensors.siteinds(tn::AbstractITensorNetwork, vertex)
+function siteinds(tn::AbstractITensorNetwork, vertex)
+  return uniqueinds(tn, vertex)
+end
+# Fix ambiguity error with IndsNetwork constructor.
+function siteinds(tn::AbstractITensorNetwork, vertex::Int)
   return uniqueinds(tn, vertex)
 end
 
@@ -311,7 +314,7 @@ function ITensors.commoninds(tn::AbstractITensorNetwork, edge)
   return commoninds(tn[src(e)], tn[dst(e)])
 end
 
-function ITensorMPS.linkinds(tn::AbstractITensorNetwork, edge)
+function linkinds(tn::AbstractITensorNetwork, edge)
   return commoninds(tn, edge)
 end
 
@@ -807,7 +810,7 @@ end
 # Link dimensions
 # 
 
-function ITensorMPS.maxlinkdim(tn::AbstractITensorNetwork)
+function maxlinkdim(tn::AbstractITensorNetwork)
   md = 1
   for e in edges(tn)
     md = max(md, linkdim(tn, e))
@@ -815,16 +818,16 @@ function ITensorMPS.maxlinkdim(tn::AbstractITensorNetwork)
   return md
 end
 
-function ITensorMPS.linkdim(tn::AbstractITensorNetwork, edge::Pair)
+function linkdim(tn::AbstractITensorNetwork, edge::Pair)
   return linkdim(tn, edgetype(tn)(edge))
 end
 
-function ITensorMPS.linkdim(tn::AbstractITensorNetwork{V}, edge::AbstractEdge{V}) where {V}
+function linkdim(tn::AbstractITensorNetwork{V}, edge::AbstractEdge{V}) where {V}
   ls = linkinds(tn, edge)
   return prod([isnothing(l) ? 1 : dim(l) for l in ls])
 end
 
-function ITensorMPS.linkdims(tn::AbstractITensorNetwork{V}) where {V}
+function linkdims(tn::AbstractITensorNetwork{V}) where {V}
   ld = DataGraph{V}(
     copy(underlying_graph(tn)); vertex_data_eltype=Nothing, edge_data_eltype=Int
   )
@@ -882,7 +885,7 @@ is_multi_edge(tn::AbstractITensorNetwork, e) = length(linkinds(tn, e)) > 1
 is_multi_edge(tn::AbstractITensorNetwork) = Base.Fix1(is_multi_edge, tn)
 
 """Add two itensornetworks together by growing the bond dimension. The network structures need to be have the same vertex names, same site index on each vertex """
-function ITensorMPS.add(tn1::AbstractITensorNetwork, tn2::AbstractITensorNetwork)
+function add(tn1::AbstractITensorNetwork, tn2::AbstractITensorNetwork)
   @assert issetequal(vertices(tn1), vertices(tn2))
 
   tn1 = combine_linkinds(tn1; edges=filter(is_multi_edge(tn1), edges(tn1)))
diff --git a/src/apply.jl b/src/apply.jl
index 3d2cfc0e..b4946210 100644
--- a/src/apply.jl
+++ b/src/apply.jl
@@ -23,7 +23,6 @@ using ITensors:
   unioninds,
   uniqueinds
 using ITensors.ContractionSequenceOptimization: optimal_contraction_sequence
-using ITensorMPS: siteinds
 using KrylovKit: linsolve
 using LinearAlgebra: eigen, norm, svd
 using NamedGraphs: NamedEdge, has_edge
diff --git a/src/boundarymps.jl b/src/boundarymps.jl
deleted file mode 100644
index b37f4249..00000000
--- a/src/boundarymps.jl
+++ /dev/null
@@ -1,16 +0,0 @@
-using ITensors: inner
-using ITensorMPS: ITensorMPS
-
-#Given an ITensorNetwork on an Lx*Ly grid with sites indexed as (i,j) then perform contraction using a sequence of mps-mpo contractions
-function contract_boundary_mps(tn::ITensorNetwork; kwargs...)
-  dims = maximum(vertices(tn))
-  d1, d2 = dims
-  vL = ITensorMPS.MPS([tn[i1, 1] for i1 in 1:d1])
-  for i2 in 2:(d2 - 2)
-    T = ITensorMPS.MPO([tn[i1, i2] for i1 in 1:d1])
-    vL = contract(T, vL; kwargs...)
-  end
-  T = ITensorMPS.MPO([tn[i1, d2 - 1] for i1 in 1:d1])
-  vR = ITensorMPS.MPS([tn[i1, d2] for i1 in 1:d1])
-  return inner(dag(vL), T, vR)[]
-end
diff --git a/src/caches/abstractbeliefpropagationcache.jl b/src/caches/abstractbeliefpropagationcache.jl
index b5163d3b..63b7e43d 100644
--- a/src/caches/abstractbeliefpropagationcache.jl
+++ b/src/caches/abstractbeliefpropagationcache.jl
@@ -2,7 +2,6 @@ using Graphs: IsDirected
 using SplitApplyCombine: group
 using LinearAlgebra: diag, dot
 using ITensors: dir
-using ITensorMPS: ITensorMPS
 using NamedGraphs.PartitionedGraphs:
   PartitionedGraphs,
   PartitionedGraph,
@@ -140,7 +139,7 @@ for f in [
   :(PartitionedGraphs.partitionvertices),
   :(PartitionedGraphs.vertices),
   :(PartitionedGraphs.boundary_partitionedges),
-  :(ITensorMPS.linkinds),
+  :(linkinds),
 ]
   @eval begin
     function $f(bpc::AbstractBeliefPropagationCache, args...; kwargs...)
diff --git a/src/caches/beliefpropagationcache.jl b/src/caches/beliefpropagationcache.jl
index 31d43e9b..1e820107 100644
--- a/src/caches/beliefpropagationcache.jl
+++ b/src/caches/beliefpropagationcache.jl
@@ -2,7 +2,6 @@ using Graphs: IsDirected
 using SplitApplyCombine: group
 using LinearAlgebra: diag, dot
 using ITensors: dir
-using ITensorMPS: ITensorMPS
 using NamedGraphs.PartitionedGraphs:
   PartitionedGraphs,
   PartitionedGraph,
diff --git a/src/expect.jl b/src/expect.jl
index 69ce42ae..ac9752cf 100644
--- a/src/expect.jl
+++ b/src/expect.jl
@@ -1,10 +1,9 @@
 using Dictionaries: Dictionary, set!
-using ITensors: Op, op, contract, siteinds, which_op
-using ITensorMPS: ITensorMPS, expect
+using ITensors: Op, op, contract, which_op
 
 default_expect_alg() = "bp"
 
-function ITensorMPS.expect(
+function expect(
   ψIψ::AbstractFormNetwork, op::Op; contract_kwargs=(; sequence="automatic"), kwargs...
 )
   v = only(op.sites)
@@ -18,7 +17,7 @@ function ITensorMPS.expect(
   return numerator / denominator
 end
 
-function ITensorMPS.expect(
+function expect(
   alg::Algorithm,
   ψ::AbstractITensorNetwork,
   ops;
@@ -40,25 +39,21 @@ function ITensorMPS.expect(
   return map(op -> expect(ψIψ, op; alg, cache!, update_cache=false, kwargs...), ops)
 end
 
-function ITensorMPS.expect(alg::Algorithm"exact", ψ::AbstractITensorNetwork, ops; kwargs...)
+function expect(alg::Algorithm"exact", ψ::AbstractITensorNetwork, ops; kwargs...)
   ψIψ = inner_network(ψ, ψ)
   return map(op -> expect(ψIψ, op; alg, kwargs...), ops)
 end
 
-function ITensorMPS.expect(
-  ψ::AbstractITensorNetwork, op::Op; alg=default_expect_alg(), kwargs...
-)
+function expect(ψ::AbstractITensorNetwork, op::Op; alg=default_expect_alg(), kwargs...)
   return expect(Algorithm(alg), ψ, [op]; kwargs...)
 end
 
-function ITensorMPS.expect(
+function expect(
   ψ::AbstractITensorNetwork, op::String, vertices; alg=default_expect_alg(), kwargs...
 )
   return expect(Algorithm(alg), ψ, [Op(op, vertex) for vertex in vertices]; kwargs...)
 end
 
-function ITensorMPS.expect(
-  ψ::AbstractITensorNetwork, op::String; alg=default_expect_alg(), kwargs...
-)
+function expect(ψ::AbstractITensorNetwork, op::String; alg=default_expect_alg(), kwargs...)
   return expect(ψ, op, vertices(ψ); alg, kwargs...)
 end
diff --git a/src/inner.jl b/src/inner.jl
index 6f1cdc86..70bced01 100644
--- a/src/inner.jl
+++ b/src/inner.jl
@@ -1,5 +1,4 @@
 using ITensors: inner, scalar
-using ITensorMPS: ITensorMPS, loginner
 using LinearAlgebra: norm, norm_sqr
 
 default_contract_alg(tns::Tuple) = "bp"
@@ -54,7 +53,7 @@ function ITensors.inner(
   return scalar(tn; sequence)
 end
 
-function ITensorMPS.loginner(
+function loginner(
   ϕ::AbstractITensorNetwork,
   ψ::AbstractITensorNetwork;
   alg=default_contract_alg((ϕ, ψ)),
@@ -63,7 +62,7 @@ function ITensorMPS.loginner(
   return loginner(Algorithm(alg), ϕ, ψ; kwargs...)
 end
 
-function ITensorMPS.loginner(
+function loginner(
   ϕ::AbstractITensorNetwork,
   A::AbstractITensorNetwork,
   ψ::AbstractITensorNetwork;
@@ -73,13 +72,13 @@ function ITensorMPS.loginner(
   return loginner(Algorithm(alg), ϕ, A, ψ; kwargs...)
 end
 
-function ITensorMPS.loginner(
+function loginner(
   alg::Algorithm"exact", ϕ::AbstractITensorNetwork, ψ::AbstractITensorNetwork; kwargs...
 )
   return log(inner(alg, ϕ, ψ); kwargs...)
 end
 
-function ITensorMPS.loginner(
+function loginner(
   alg::Algorithm"exact",
   ϕ::AbstractITensorNetwork,
   A::AbstractITensorNetwork,
@@ -89,7 +88,7 @@ function ITensorMPS.loginner(
   return log(inner(alg, ϕ, A, ψ); kwargs...)
 end
 
-function ITensorMPS.loginner(
+function loginner(
   alg::Algorithm,
   ϕ::AbstractITensorNetwork,
   ψ::AbstractITensorNetwork;
@@ -100,7 +99,7 @@ function ITensorMPS.loginner(
   return logscalar(alg, tn; kwargs...)
 end
 
-function ITensorMPS.loginner(
+function loginner(
   alg::Algorithm,
   ϕ::AbstractITensorNetwork,
   A::AbstractITensorNetwork,
diff --git a/src/mpo_mps_compatibility.jl b/src/mpo_mps_compatibility.jl
deleted file mode 100644
index ac6d43f5..00000000
--- a/src/mpo_mps_compatibility.jl
+++ /dev/null
@@ -1,20 +0,0 @@
-using ITensorMPS: ITensorMPS
-
-function ITensorMPS.MPO(opsum::OpSum, s::IndsNetwork)
-  s_linear = [only(s[v]) for v in 1:nv(s)]
-  return ITensorMPS.MPO(opsum, s_linear)
-end
-
-function ITensorMPS.MPO(opsum_sum::Sum{<:OpSum}, s::IndsNetwork)
-  return ITensorMPS.MPO(sum(Ops.terms(opsum_sum)), s)
-end
-
-function ITensorMPS.random_mps(s::IndsNetwork, args...; kwargs...)
-  s_linear = [only(s[v]) for v in 1:nv(s)]
-  return ITensorMPS.random_mps(s_linear, args...; kwargs...)
-end
-
-function ITensorMPS.MPS(s::IndsNetwork, args...; kwargs...)
-  s_linear = [only(s[v]) for v in 1:nv(s)]
-  return ITensorMPS.MPS(s_linear, args...; kwargs...)
-end
diff --git a/src/partitioneditensornetwork.jl b/src/partitioneditensornetwork.jl
index 7a023adc..eab99de1 100644
--- a/src/partitioneditensornetwork.jl
+++ b/src/partitioneditensornetwork.jl
@@ -1,10 +1,9 @@
 using Graphs: dst, src
 using ITensors: commoninds
-using ITensorMPS: ITensorMPS
 using NamedGraphs.GraphsExtensions: subgraph
 using NamedGraphs.PartitionedGraphs: PartitionedGraph, PartitionEdge
 
-function ITensorMPS.linkinds(pitn::PartitionedGraph, edge::PartitionEdge)
+function linkinds(pitn::PartitionedGraph, edge::PartitionEdge)
   src_e_itn = subgraph(pitn, src(edge))
   dst_e_itn = subgraph(pitn, dst(edge))
   return commoninds(src_e_itn, dst_e_itn)
diff --git a/src/sitetype.jl b/src/sitetype.jl
index 570f711c..0b893a82 100644
--- a/src/sitetype.jl
+++ b/src/sitetype.jl
@@ -1,6 +1,6 @@
 using Dictionaries: Dictionary
 using Graphs: AbstractGraph, nv, vertices
-using ITensors: ITensors, Index, siteind, siteinds
+using ITensors: ITensors, Index, siteind
 
 function ITensors.siteind(sitetype::String, v::Tuple; kwargs...)
   return addtags(siteind(sitetype; kwargs...), vertex_tag(v))
@@ -16,17 +16,22 @@ function to_siteinds_callable(x::AbstractDictionary)
   return Base.Fix1(getindex, x) ∘ keytype(x)
 end
 
-function ITensors.siteinds(x, g::AbstractGraph; kwargs...)
+function siteinds(x, g::AbstractGraph; kwargs...)
   return siteinds(to_siteinds_callable(x), g; kwargs...)
 end
 
+# Convenient syntax for path graph.
+function siteinds(x, nv::Int; kwargs...)
+  return siteinds(x, path_graph(nv); kwargs...)
+end
+
 function to_siteind(x, vertex; kwargs...)
   return [siteind(x, vertex_tag(vertex); kwargs...)]
 end
 
 to_siteind(x::Index, vertex; kwargs...) = [x]
 
-function ITensors.siteinds(f::Function, g::AbstractGraph; kwargs...)
+function siteinds(f::Function, g::AbstractGraph; kwargs...)
   is = IndsNetwork(g)
   for v in vertices(g)
     is[v] = to_siteind(f(v), v; kwargs...)
diff --git a/src/solvers/alternating_update/alternating_update.jl b/src/solvers/alternating_update/alternating_update.jl
index 750f3f36..69f965f0 100644
--- a/src/solvers/alternating_update/alternating_update.jl
+++ b/src/solvers/alternating_update/alternating_update.jl
@@ -1,5 +1,4 @@
 using ITensors: state
-using ITensorMPS: linkind
 using NamedGraphs.GraphsExtensions: GraphsExtensions
 
 function alternating_update(
diff --git a/src/solvers/contract.jl b/src/solvers/contract.jl
index 0fff3894..00e5c4d6 100644
--- a/src/solvers/contract.jl
+++ b/src/solvers/contract.jl
@@ -1,6 +1,5 @@
 using Graphs: nv, vertices
 using ITensors: ITensors, sim
-using ITensorMPS: linkinds
 using ITensors.NDTensors: Algorithm, @Algorithm_str, contract
 using NamedGraphs: vertextype
 
diff --git a/src/solvers/defaults.jl b/src/solvers/defaults.jl
index 09c2ae2f..5f25e087 100644
--- a/src/solvers/defaults.jl
+++ b/src/solvers/defaults.jl
@@ -1,5 +1,4 @@
 using Printf: @printf, @sprintf
-using ITensorMPS: maxlinkdim
 default_outputlevel() = 0
 default_nsites() = 2
 default_nsweeps() = 1 #? or nothing?
diff --git a/src/solvers/dmrg.jl b/src/solvers/dmrg.jl
index 19511505..1acbde35 100644
--- a/src/solvers/dmrg.jl
+++ b/src/solvers/dmrg.jl
@@ -1,10 +1,6 @@
-using ITensorMPS: ITensorMPS, dmrg
 using KrylovKit: KrylovKit
 
-"""
-Overload of `ITensors.ITensorMPS.dmrg`.
-"""
-function ITensorMPS.dmrg(
+function dmrg(
   operator,
   init_state;
   nsweeps,
@@ -24,7 +20,4 @@ function ITensorMPS.dmrg(
   return eigval, state
 end
 
-"""
-Overload of `KrylovKit.eigsolve`.
-"""
 KrylovKit.eigsolve(H, init::AbstractTTN; kwargs...) = dmrg(H, init; kwargs...)
diff --git a/src/treetensornetworks/abstracttreetensornetwork.jl b/src/treetensornetworks/abstracttreetensornetwork.jl
index f2f82ef1..5d57bb46 100644
--- a/src/treetensornetworks/abstracttreetensornetwork.jl
+++ b/src/treetensornetworks/abstracttreetensornetwork.jl
@@ -9,7 +9,6 @@ using NamedGraphs.GraphsExtensions:
 using NamedGraphs: namedgraph_a_star, steiner_tree
 using IsApprox: IsApprox, Approx
 using ITensors: ITensors, Algorithm, @Algorithm_str, directsum, hasinds, permute, plev
-using ITensorMPS: ITensorMPS, linkind, loginner, lognorm, orthogonalize
 using TupleTools: TupleTools
 
 abstract type AbstractTreeTensorNetwork{V} <: AbstractITensorNetwork{V} end
@@ -44,7 +43,7 @@ function gauge(alg::Algorithm, ttn::AbstractTTN, region; kwargs...)
   return gauge(alg, ttn, [region]; kwargs...)
 end
 
-function ITensorMPS.orthogonalize(ttn::AbstractTTN, region; kwargs...)
+function orthogonalize(ttn::AbstractTTN, region; kwargs...)
   return gauge(Algorithm("orthogonalize"), ttn, region; kwargs...)
 end
 
@@ -163,7 +162,7 @@ function LinearAlgebra.rmul!(tn::AbstractTTN, α::Number)
   return _apply_to_ortho_region!(*, tn, α)
 end
 
-function ITensorMPS.lognorm(tn::AbstractTTN)
+function lognorm(tn::AbstractTTN)
   if isone(length(ortho_region(tn)))
     return log(norm(tn[only(ortho_region(tn))]))
   end
@@ -181,7 +180,7 @@ function logdot(tn1::AbstractTTN, tn2::AbstractTTN; kwargs...)
 end
 
 # TODO: stick with this traversal or find optimal contraction sequence?
-function ITensorMPS.loginner(
+function loginner(
   tn1::AbstractTTN, tn2::AbstractTTN; root_vertex=GraphsExtensions.default_root_vertex(tn1)
 )
   N = nv(tn1)
@@ -270,13 +269,13 @@ end
 
 Base.:+(tn::AbstractTTN) = tn
 
-ITensorMPS.add(tns::AbstractTTN...; kwargs...) = +(tns...; kwargs...)
+add(tns::AbstractTTN...; kwargs...) = +(tns...; kwargs...)
 
 function Base.:-(tn1::AbstractTTN, tn2::AbstractTTN; kwargs...)
   return +(tn1, -tn2; kwargs...)
 end
 
-function ITensorMPS.add(tn1::AbstractTTN, tn2::AbstractTTN; kwargs...)
+function add(tn1::AbstractTTN, tn2::AbstractTTN; kwargs...)
   return +(tn1, tn2; kwargs...)
 end
 
@@ -346,7 +345,7 @@ function ITensors.inner(
   return O[]
 end
 
-function ITensorMPS.expect(
+function expect(
   operator::String,
   state::AbstractTTN;
   vertices=vertices(state),
diff --git a/src/treetensornetworks/opsum_to_ttn/opsum_to_ttn.jl b/src/treetensornetworks/opsum_to_ttn/opsum_to_ttn.jl
index 09098b68..796a1aa5 100644
--- a/src/treetensornetworks/opsum_to_ttn/opsum_to_ttn.jl
+++ b/src/treetensornetworks/opsum_to_ttn/opsum_to_ttn.jl
@@ -1,10 +1,9 @@
 #using FillArrays: OneElement
 #using DataGraphs: DataGraph
 using Graphs: degree, is_tree
-using ITensorMPS: ITensorMPS, cutoff, linkdims, ops, truncate!, val
 using ITensors: flux, has_fermion_string, itensor, removeqns, space
 using ITensors.LazyApply: Prod, Sum, coefficient
-using ITensors.NDTensors: Block, blockdim, maxdim, nblocks, nnzblocks
+using ITensors.NDTensors: Block, blockdim, maxdim, nblocks, nnzblocks, truncate!
 using ITensors.Ops: argument, coefficient, Op, OpSum, name, params, site, terms, which_op
 using NamedGraphs.GraphsExtensions:
   GraphsExtensions, boundary_edges, degrees, is_leaf_vertex, vertex_path
@@ -650,6 +649,10 @@ end
 function mpo(os::OpSum, external_inds::Vector; kwargs...)
   return ttn(os, path_indsnetwork(external_inds); kwargs...)
 end
+function mpo(os::OpSum, s::IndsNetwork; kwargs...)
+  # TODO: Check it is a path graph.
+  return ttn(os, s; kwargs...)
+end
 
 # Conversion from other formats
 function ttn(o::Op, s::IndsNetwork; kwargs...)
diff --git a/src/treetensornetworks/projttns/abstractprojttn.jl b/src/treetensornetworks/projttns/abstractprojttn.jl
index 55852677..faaaf053 100644
--- a/src/treetensornetworks/projttns/abstractprojttn.jl
+++ b/src/treetensornetworks/projttns/abstractprojttn.jl
@@ -1,7 +1,6 @@
 using DataGraphs: DataGraphs, underlying_graph
 using Graphs: neighbors
 using ITensors: ITensor, contract, order, product
-using ITensorMPS: ITensorMPS, nsite
 using NamedGraphs: NamedGraphs, NamedEdge, vertextype
 using NamedGraphs.GraphsExtensions: incident_edges
 
@@ -29,7 +28,7 @@ Graphs.edgetype(P::AbstractProjTTN) = edgetype(underlying_graph(P))
 
 on_edge(P::AbstractProjTTN) = isa(pos(P), edgetype(P))
 
-ITensorMPS.nsite(P::AbstractProjTTN) = on_edge(P) ? 0 : length(pos(P))
+nsite(P::AbstractProjTTN) = on_edge(P) ? 0 : length(pos(P))
 
 function sites(P::AbstractProjTTN{V}) where {V}
   on_edge(P) && return V[]
diff --git a/src/treetensornetworks/projttns/projttnsum.jl b/src/treetensornetworks/projttns/projttnsum.jl
index 42ae6a05..ceff7d6d 100644
--- a/src/treetensornetworks/projttns/projttnsum.jl
+++ b/src/treetensornetworks/projttns/projttnsum.jl
@@ -28,7 +28,7 @@ end
 
 on_edge(P::ProjTTNSum) = on_edge(terms(P)[1])
 
-ITensorMPS.nsite(P::ProjTTNSum) = nsite(terms(P)[1])
+nsite(P::ProjTTNSum) = nsite(terms(P)[1])
 
 function set_nsite(Ps::ProjTTNSum, nsite)
   return ProjTTNSum(map(p -> set_nsite(p, nsite), terms(Ps)), factors(Ps))
diff --git a/test/Project.toml b/test/Project.toml
index 09126944..39a36dbb 100644
--- a/test/Project.toml
+++ b/test/Project.toml
@@ -10,7 +10,6 @@ EinExprs = "b1794770-133b-4de1-afb4-526377e9f4c5"
 Glob = "c27321d9-0574-5035-807b-f59d2c89b15c"
 Graphs = "86223c79-3864-5bf0-83f7-82e725a168b6"
 GraphsFlows = "06909019-6f44-4949-96fc-b9d9aaa02889"
-ITensorGaussianMPS = "2be41995-7c9f-4653-b682-bfa4e7cebb93"
 ITensorMPS = "0d1a4710-d33b-49a5-8f18-73bdf49b47e2"
 ITensorNetworks = "2919e153-833c-4bdc-8836-1ea460a35fc7"
 ITensorUnicodePlots = "73163f41-4a9e-479f-8353-73bf94dbd758"
@@ -45,25 +44,24 @@ EinExprs = "0.6.8"
 Glob = "1.3.1"
 Graphs = "1.12.0"
 GraphsFlows = "0.1.1"
-ITensorGaussianMPS = "0.1.12"
 ITensorMPS = "0.3.6"
-ITensorNetworks = "0.12.2"
+ITensorNetworks = "0.13.0"
 ITensorUnicodePlots = "0.1.6"
 ITensors = "0.7, 0.8, 0.9"
 KrylovKit = "0.8, 0.9"
+LinearAlgebra = "1.10.0"
 Metis = "1.5.0"
 NDTensors = "0.3, 0.4"
 NamedGraphs = "0.6.6"
 OMEinsumContractionOrders = "0.9.5"
 Observers = "0.2.5"
 OrdinaryDiffEqTsit5 = "1.1.0"
+Pkg = "1.10.0"
+Random = "1.10.0"
 SafeTestsets = "0.1.0"
 SplitApplyCombine = "1.2.3"
 StableRNGs = "1.0.2"
 Suppressor = "0.2.8"
 TensorOperations = "5.1.4"
-UnicodePlots = "3.7.2"
-LinearAlgebra = "1.10.0"
-Pkg = "1.10.0"
-Random = "1.10.0"
 Test = "1.10.0"
+UnicodePlots = "3.7.2"
diff --git a/test/test_belief_propagation.jl b/test/test_belief_propagation.jl
index 66224386..76e49ebd 100644
--- a/test/test_belief_propagation.jl
+++ b/test/test_belief_propagation.jl
@@ -8,7 +8,6 @@ using ITensorNetworks:
   ⊗,
   combine_linkinds,
   contract,
-  contract_boundary_mps,
   contraction_sequence,
   eachtensor,
   environment,
diff --git a/test/test_binary_tree_partition.jl b/test/test_binary_tree_partition.jl
index 239c0c14..f44536e7 100644
--- a/test/test_binary_tree_partition.jl
+++ b/test/test_binary_tree_partition.jl
@@ -4,7 +4,6 @@ using Graphs: add_vertex!, vertices
 # Trigger package extension.
 using GraphsFlows: GraphsFlows
 using ITensors: Index, ITensor, contract, noncommoninds, random_itensor
-using ITensorMPS: MPS
 using ITensorNetworks:
   _DensityMartrixAlgGraph,
   _contract_deltas_ignore_leaf_partitions,
@@ -27,38 +26,6 @@ using StableRNGs: StableRNG
 using TensorOperations: TensorOperations
 using Test: @test, @testset
 
-@testset "test mincut functions on top of MPS" begin
-  i = Index(2, "i")
-  j = Index(2, "j")
-  k = Index(2, "k")
-  l = Index(2, "l")
-  m = Index(2, "m")
-  n = Index(2, "n")
-  o = Index(2, "o")
-  p = Index(2, "p")
-
-  T = random_itensor(i, j, k, l, m, n, o, p)
-  M = MPS(T, (i, j, k, l, m, n, o, p); cutoff=1e-5, maxdim=500)
-  tn = ITensorNetwork(M[:])
-  for out in [binary_tree_structure(tn), path_graph_structure(tn)]
-    @test out isa DataGraph
-    @test is_binary_arborescence(out)
-    @test length(vertex_data(out).values) == 8
-  end
-  out = _mps_partition_inds_order(tn, [o, p, i, j, k, l, m, n])
-  @test out in [[i, j, k, l, m, n, o, p], [p, o, n, m, l, k, j, i]]
-  p1, p2 = _mincut_partitions(tn, [k, l], [m, n])
-  # When MPS bond dimensions are large, the partition will not across internal inds
-  @test (length(p1) == 0) || (length(p2) == 0)
-
-  M = MPS(T, (i, j, k, l, m, n, o, p); cutoff=1e-5, maxdim=2)
-  tn = ITensorNetwork(M[:])
-  p1, p2 = _mincut_partitions(tn, [k, l], [m, n])
-  # When MPS bond dimensions are small, the partition will across internal inds
-  @test sort(p1) == [1, 2, 3, 4]
-  @test sort(p2) == [5, 6, 7, 8]
-end
-
 @testset "test _binary_tree_partition_inds of a 2D network" begin
   N = (3, 3, 3)
   linkdim = 2
@@ -83,77 +50,4 @@ end
   network = mapreduce(v -> collect(eachtensor(par[v])), vcat, vertices(par))
   @test isapprox(contract(tn), contract(network))
 end
-
-@testset "test partition with mincut_recursive_bisection alg and approx_itensornetwork" begin
-  i = Index(2, "i")
-  j = Index(2, "j")
-  k = Index(2, "k")
-  l = Index(2, "l")
-  m = Index(2, "m")
-  for dtype in (Float64, Complex{Float64})
-    T = random_itensor(dtype, i, j, k, l, m)
-    M = MPS(T, (i, j, k, l, m); cutoff=1e-5, maxdim=5)
-    network = M[:]
-    out1 = contract(network...)
-    tn = ITensorNetwork(network)
-    inds_btree = binary_tree_structure(tn)
-    par = _partition(tn, inds_btree; alg="mincut_recursive_bisection")
-    par = _contract_deltas_ignore_leaf_partitions(par; root=root_vertex(inds_btree))
-    networks = map(v -> par[v], vertices(par))
-    network2 = reduce(union, networks)
-    out2 = contract(network2)
-    @test isapprox(out1, out2)
-    # test approx_itensornetwork (here we call `contract` to test the interface)
-    for structure in [path_graph_structure, binary_tree_structure]
-      for alg in ["density_matrix", "ttn_svd"]
-        approx_tn, lognorm = contract(
-          tn;
-          alg,
-          output_structure=structure,
-          contraction_sequence_kwargs=(; alg="sa_bipartite"),
-        )
-        network3 = collect(eachtensor(approx_tn))
-        out3 = contract(network3...) * exp(lognorm)
-        i1 = noncommoninds(network...)
-        i3 = noncommoninds(network3...)
-        @test (length(i1) == length(i3))
-        @test isapprox(out1, out3)
-      end
-    end
-  end
-end
-
-@testset "test caching in approx_itensornetwork" begin
-  i = Index(2, "i")
-  j = Index(2, "j")
-  k = Index(2, "k")
-  l = Index(2, "l")
-  m = Index(2, "m")
-  T = random_itensor(i, j, k, l, m)
-  M = MPS(T, (i, j, k, l, m); cutoff=1e-5, maxdim=5)
-  tn = ITensorNetwork(M[:])
-  out_tree = path_graph_structure(tn)
-  input_partition = _partition(tn, out_tree; alg="mincut_recursive_bisection")
-  underlying_tree = underlying_graph(input_partition)
-  # Change type of each partition[v] since they will be updated
-  # with potential data type chage.
-  p = DataGraph(NamedGraph())
-  for v in vertices(input_partition)
-    add_vertex!(p, v)
-    p[v] = ITensorNetwork{Any}(input_partition[v])
-  end
-  alg_graph = _DensityMartrixAlgGraph(p, underlying_tree, root_vertex(out_tree))
-  path = post_order_dfs_vertices(underlying_tree, root_vertex(out_tree))
-  for v in path[1:2]
-    _rem_vertex!(
-      alg_graph,
-      v;
-      cutoff=1e-15,
-      maxdim=10000,
-      contraction_sequence_kwargs=(; alg="optimal"),
-    )
-  end
-  # Check that a specific density matrix info has been cached
-  @test haskey(alg_graph.caches.es_to_pdm, Set([NamedEdge(nothing, path[3])]))
-end
 end
diff --git a/test/test_expect.jl b/test/test_expect.jl
index 35322cc6..093b7df2 100644
--- a/test/test_expect.jl
+++ b/test/test_expect.jl
@@ -2,12 +2,12 @@
 using Graphs: SimpleGraph, uniform_tree
 using NamedGraphs: NamedGraph, vertices
 using NamedGraphs.NamedGraphGenerators: named_grid
-using ITensors: siteinds
 using ITensorNetworks:
   BeliefPropagationCache,
   ITensorNetwork,
   expect,
   random_tensornetwork,
+  siteinds,
   original_state_vertex
 using SplitApplyCombine: group
 using StableRNGs: StableRNG
diff --git a/test/test_inner.jl b/test/test_inner.jl
index f55581a9..1b91753f 100644
--- a/test/test_inner.jl
+++ b/test/test_inner.jl
@@ -7,10 +7,11 @@ using ITensorNetworks:
   logscalar,
   random_tensornetwork,
   scalar,
+  siteinds,
   ttn,
   underlying_graph
 using ITensorNetworks.ModelHamiltonians: heisenberg
-using ITensors: dag, siteinds
+using ITensors: dag
 using SplitApplyCombine: group
 using Graphs: SimpleGraph, uniform_tree
 using NamedGraphs: NamedGraph
diff --git a/test/test_opsum_to_ttn.jl b/test/test_opsum_to_ttn.jl
index d8b64310..602dab73 100644
--- a/test/test_opsum_to_ttn.jl
+++ b/test/test_opsum_to_ttn.jl
@@ -15,7 +15,6 @@ using ITensors:
   removeqns
 using ITensorMPS: ITensorMPS
 using ITensors.NDTensors: matrix
-using ITensorGaussianMPS: ITensorGaussianMPS
 using ITensorNetworks: ITensorNetworks, OpSum, ttn, siteinds
 using ITensorNetworks.ITensorsExtensions: replace_vertices
 using ITensorNetworks.ModelHamiltonians: ModelHamiltonians
@@ -189,8 +188,6 @@ end
       sites = [only(is[v]) for v in reverse(post_order_dfs_vertices(c, root_vertex))]
       vmap = Dictionary(reverse(post_order_dfs_vertices(c, root_vertex)), 1:length(sites))
       Hline = ITensorMPS.MPO(replace_vertices(v -> vmap[v], H), sites)
-      # compare resulting sparse Hamiltonians
-      Hmat_sp = ITensorGaussianMPS.hopping_hamiltonian(replace_vertices(v -> vmap[v], H))
       @disable_warn_order begin
         Tmpo = prod(Hline)
         Tttno = contract(Hsvd)
@@ -207,11 +204,6 @@ end
       dTmm = to_matrix(Tmpo)
       dTtm = to_matrix(Tttno)
       @test any(>(1e-14), dTmm - dTtm)
-
-      # also compare with energies obtained from single-particle Hamiltonian
-      GS_mb, _, _ = eigsolve(dTtm, 1, :SR, eltype(dTtm))
-      spectrum_sp = eigvals(Hmat_sp)
-      @test minimum(cumsum(spectrum_sp)) ≈ GS_mb[1] atol = 1e-8
     end
     if !auto_fermion_enabled
       ITensors.disable_auto_fermion()
diff --git a/test/test_tebd.jl b/test/test_tebd.jl
index cdec4a1b..71acd1ed 100644
--- a/test/test_tebd.jl
+++ b/test/test_tebd.jl
@@ -1,7 +1,6 @@
 @eval module $(gensym())
 using Graphs: vertices
 using ITensors: ITensors
-using ITensorMPS: ITensorMPS
 using ITensorNetworks: ITensorNetwork, cartesian_to_linear, dmrg, expect, siteinds, tebd
 using ITensorNetworks.ITensorsExtensions: group_terms
 using ITensorNetworks.ModelHamiltonians: ModelHamiltonians
@@ -20,18 +19,6 @@ ITensors.disable_warn_order()
 
   s = siteinds("S=1/2", g)
 
-  #
-  # DMRG comparison
-  #
-  g_dmrg = rename_vertices(v -> cartesian_to_linear(dims)[v], g)
-  ℋ_dmrg = ModelHamiltonians.ising(g_dmrg; h)
-  s_dmrg = [only(s[v]) for v in vertices(s)]
-  H_dmrg = ITensorMPS.MPO(ℋ_dmrg, s_dmrg)
-  ψ_dmrg_init = ITensorMPS.MPS(s_dmrg, j -> "↑")
-  E_dmrg, ψ_dmrg = dmrg(
-    H_dmrg, ψ_dmrg_init; nsweeps=20, maxdim=[fill(10, 10); 20], cutoff=1e-8, outputlevel=0
-  )
-
   #
   # PEPS TEBD optimization
   #
@@ -66,6 +53,5 @@ ITensors.disable_warn_order()
   #E2 = expect(ℋ, ψ)
   #@show E0, E1, E2, E_dmrg
   @test_broken (((abs((E2 - E1) / E2) < 1e-3) && (E1 < E0)) || (E2 < E1 < E0))
-  @test_broken E2 ≈ E_dmrg rtol = 1e-3
 end
 end
diff --git a/test/test_ttn_contract.jl b/test/test_ttn_contract.jl
index 8255eaec..500826c3 100644
--- a/test/test_ttn_contract.jl
+++ b/test/test_ttn_contract.jl
@@ -17,7 +17,6 @@ using ITensorNetworks:
   siteinds
 using ITensorNetworks.ModelHamiltonians: ModelHamiltonians
 using ITensors: prime, replaceinds, replaceprime
-using ITensorMPS: ITensorMPS
 using LinearAlgebra: norm, normalize
 using NamedGraphs.NamedGraphGenerators: named_comb_tree
 using StableRNGs: StableRNG
@@ -150,23 +149,4 @@ end
   Hpsi = contract(H, psi; alg="fit", nsites=1, nsweeps=10, init=Hpsi_guess)
   @test inner(psit, Hpsi) ≈ inner(psit, H, psi) rtol = 1e-2
 end
-
-@testset "Contract TTN with dangling inds" begin
-  nbit = 3
-  sites = siteinds("Qubit", nbit)
-
-  # random_mpo does not support linkdims keyword.
-  M1 = replaceprime(
-    ITensorMPS.random_mpo(sites) + ITensorMPS.random_mpo(sites), 1 => 2, 0 => 1
-  )
-  M2 = ITensorMPS.random_mpo(sites) + ITensorMPS.random_mpo(sites)
-  M12_ref = contract(M1, M2; alg="naive")
-  t12_ref = ttn([M12_ref[v] for v in eachindex(M12_ref)])
-
-  t1 = ttn([M1[v] for v in eachindex(M1)])
-  t2 = ttn([M2[v] for v in eachindex(M2)])
-
-  # Test with good initial guess
-  @test contract(t1, t2; alg="fit", init=t12_ref, nsweeps=1) ≈ t12_ref rtol = 1e-7
-end
 end
diff --git a/test/test_ttn_dmrg.jl b/test/test_ttn_dmrg.jl
index 004ec561..b8a8cdb8 100644
--- a/test/test_ttn_dmrg.jl
+++ b/test/test_ttn_dmrg.jl
@@ -55,7 +55,7 @@ ITensors.disable_auto_fermion()
   # Compare to `ITensors.MPO` version of `dmrg`
   H_mpo = ITensorMPS.MPO([H[v] for v in 1:nv(H)])
   psi_mps = ITensorMPS.MPS([psi[v] for v in 1:nv(psi)])
-  e2, psi2 = dmrg(H_mpo, psi_mps; nsweeps, maxdim, outputlevel=0)
+  e2, psi2 = ITensorMPS.dmrg(H_mpo, psi_mps; nsweeps, maxdim, outputlevel=0)
 
   e, psi = dmrg(
     H, psi; nsweeps, maxdim, cutoff, nsites, updater_kwargs=(; krylovdim=3, maxiter=1)
@@ -242,9 +242,9 @@ end
     Hline = ITensorMPS.MPO(replace_vertices(v -> vmap[v], os), sline)
     rng = StableRNG(1234)
     psiline = ITensorMPS.random_mps(rng, sline, i -> isodd(i) ? "Up" : "Dn"; linkdims=20)
-    e2, psi2 = dmrg(Hline, psiline; nsweeps, maxdim, cutoff, outputlevel=0)
+    e2, psi2 = ITensorMPS.dmrg(Hline, psiline; nsweeps, maxdim, cutoff, outputlevel=0)
 
-    @test inner(psi', H, psi) ≈ inner(psi2', Hline, psi2) atol = 1e-5
+    @test inner(psi', H, psi) ≈ ITensorMPS.inner(psi2', Hline, psi2) atol = 1e-5
 
     if !auto_fermion_enabled
       ITensors.disable_auto_fermion()
@@ -278,7 +278,7 @@ end
   Hline = ITensorMPS.MPO(replace_vertices(v -> vmap[v], os), sline)
   rng = StableRNG(1234)
   psiline = ITensorMPS.random_mps(rng, sline, i -> isodd(i) ? "Up" : "Dn"; linkdims=20)
-  e_jw, psi_jw = dmrg(Hline, psiline; nsweeps, maxdim, cutoff, outputlevel=0)
+  e_jw, psi_jw = ITensorMPS.dmrg(Hline, psiline; nsweeps, maxdim, cutoff, outputlevel=0)
   ITensors.enable_auto_fermion()
 
   # now get auto-fermion results 
@@ -298,9 +298,9 @@ end
   Hline = ITensorMPS.MPO(replace_vertices(v -> vmap[v], os), sline)
   rng = StableRNG(1234)
   psiline = ITensorMPS.random_mps(rng, sline, i -> isodd(i) ? "Up" : "Dn"; linkdims=20)
-  e2, psi2 = dmrg(Hline, psiline; nsweeps, maxdim, cutoff, outputlevel=0)
+  e2, psi2 = ITensorMPS.dmrg(Hline, psiline; nsweeps, maxdim, cutoff, outputlevel=0)
 
-  @test inner(psi', H, psi) ≈ inner(psi2', Hline, psi2) atol = 1e-5
+  @test inner(psi', H, psi) ≈ ITensorMPS.inner(psi2', Hline, psi2) atol = 1e-5
   @test e2 ≈ e_jw atol = 1e-5
   @test inner(psi2', Hline, psi2) ≈ e_jw atol = 1e-5
 
diff --git a/test/test_ttn_dmrg_x.jl b/test/test_ttn_dmrg_x.jl
index f718ba7d..4f2583ac 100644
--- a/test/test_ttn_dmrg_x.jl
+++ b/test/test_ttn_dmrg_x.jl
@@ -1,7 +1,6 @@
 @eval module $(gensym())
 using Dictionaries: Dictionary
 using Graphs: nv, vertices
-using ITensorMPS: linkdims
 using ITensorNetworks:
   OpSum, ttn, apply, contract, dmrg_x, inner, mpo, mps, random_mps, siteinds
 using ITensorNetworks.ModelHamiltonians: ModelHamiltonians
diff --git a/test/test_ttn_expect.jl b/test/test_ttn_expect.jl
index 9bf55248..f819f109 100644
--- a/test/test_ttn_expect.jl
+++ b/test/test_ttn_expect.jl
@@ -1,23 +1,10 @@
 @eval module $(gensym())
 using Graphs: vertices
-using ITensorMPS: MPS
 using ITensorNetworks: ttn, expect, random_mps, siteinds
 using LinearAlgebra: norm
 using NamedGraphs.NamedGraphGenerators: named_comb_tree
 using StableRNGs: StableRNG
 using Test: @test, @testset
-@testset "MPS expect comparison with ITensors" begin
-  N = 4
-  s = siteinds("S=1/2", N)
-  rng = StableRNG(1234)
-  a = random_mps(rng, s; link_space=100)
-  @test norm(a) ≈ 1
-  b = MPS([a[v] for v in vertices(a)])
-  res_a = expect("Sz", a)
-  res_b = expect(b, "Sz")
-  res_a = [res_a[v] for v in vertices(a)]
-  @test res_a ≈ res_b rtol = 1e-6
-end
 @testset "TTN expect" begin
   tooth_lengths = fill(2, 2)
   c = named_comb_tree(tooth_lengths)