diff --git a/Project.toml b/Project.toml
index e28bc99..20c6724 100644
--- a/Project.toml
+++ b/Project.toml
@@ -1,20 +1,30 @@
 name = "ITensorMPS"
 uuid = "0d1a4710-d33b-49a5-8f18-73bdf49b47e2"
 authors = ["Matthew Fishman <mfishman@flatironinstitute.org>", "Miles Stoudenmire <mstoudenmire@flatironinstitute.org>"]
-version = "0.3.18"
+version = "0.4.0"
 
 [deps]
 Adapt = "79e6a3ab-5dfb-504d-930d-738a2a938a0e"
+BackendSelection = "680c2d7c-f67a-4cc9-ae9c-da132b1447a5"
 Compat = "34da2185-b29b-5c13-b0c7-acf172513d20"
+DiagonalArrays = "74fd4be6-21e2-4f6f-823a-4360d37c7a77"
+FillArrays = "1a297f60-69ca-5386-bcde-b61e274b549b"
+GradedArrays = "bc96ca6e-b7c8-4bb6-888e-c93f838762c2"
 ITensors = "9136182c-28ba-11e9-034c-db9fb085ebd5"
 IsApprox = "28f27b66-4bd8-47e7-9110-e2746eb8bed7"
 KrylovKit = "0b1a1467-8014-51b9-945f-bf0ae24f4b77"
 LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
-NDTensors = "23ae76d9-e61a-49c4-8f12-3f1a16adf9cf"
+NamedDimsArrays = "60cbd0c0-df58-4cb7-918c-6f5607b73fde"
 Printf = "de0858da-6303-5e67-8744-51eddeeeb8d7"
+QuantumOperatorAlgebra = "c7a6d0f7-daa6-4368-ba67-89ed64127c3b"
+QuantumOperatorDefinitions = "826dd319-6fd5-459a-a990-3a4f214664bf"
 Random = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c"
 SerializedElementArrays = "d3ce8812-9567-47e9-a7b5-65a6d70a3065"
+SparseArraysBase = "0d5efcca-f356-4864-8770-e1ed8d78f208"
+TensorAlgebra = "68bd88dc-f39d-4e12-b2ca-f046b68fcc6a"
 TupleTools = "9d95972d-f1c8-5527-a6e0-b4b365fa01f6"
+TypeParameterAccessors = "7e5a90cf-f82e-492e-a09b-e3e26432c138"
+VectorInterface = "409d34a3-91d5-4945-b6ec-7529ddf182d8"
 
 [weakdeps]
 ChainRulesCore = "d360d2e6-b24c-11e9-a2a3-2a2ae2dbcce4"
@@ -31,22 +41,32 @@ ITensorMPSPackageCompilerExt = "PackageCompiler"
 ITensorMPSZygoteRulesExt = ["ChainRulesCore", "ZygoteRules"]
 
 [compat]
-Adapt = "4.1.0"
+Adapt = "4.1"
+BackendSelection = "0.1"
 ChainRulesCore = "1.10"
 Compat = "4.16.0"
+DiagonalArrays = "0.3"
+FillArrays = "1.13"
+GradedArrays = "0.4"
 HDF5 = "0.14, 0.15, 0.16, 0.17"
-ITensors = "0.8, 0.9"
-IsApprox = "2.0.0"
+ITensors = "0.10"
+IsApprox = "2.0"
 KrylovKit = "0.8.1, 0.9"
 LinearAlgebra = "1.10"
-NDTensors = "0.4"
+NamedDimsArrays = "0.7"
 Observers = "0.2"
 PackageCompiler = "1, 2"
 Printf = "1.10"
+QuantumOperatorAlgebra = "0.1"
+QuantumOperatorDefinitions = "0.2"
 Random = "1.10"
-SerializedElementArrays = "0.1.0"
+SerializedElementArrays = "0.1"
+SparseArraysBase = "0.5"
+TensorAlgebra = "0.3"
 Test = "1.10"
-TupleTools = "1.5.0"
+TupleTools = "1.5"
+TypeParameterAccessors = "0.3"
+VectorInterface = "0.5"
 ZygoteRules = "0.2.2"
 julia = "1.10"
 
diff --git a/examples/dmrg/1d_heisenberg.jl b/examples/dmrg/1d_heisenberg.jl
index 4718721..8e7febe 100644
--- a/examples/dmrg/1d_heisenberg.jl
+++ b/examples/dmrg/1d_heisenberg.jl
@@ -2,37 +2,62 @@ using ITensors, ITensorMPS
 using Printf
 using Random
 
-Random.seed!(1234)
+using Adapt: adapt
+using JLArrays: JLArray
+using Metal: MtlArray
+using SerializedArrays: SerializedArray
 
-let
-  N = 100
+dev = adapt(SerializedArray)
 
-  # Create N spin-one degrees of freedom
-  sites = siteinds("S=1", N)
-  # Alternatively can make spin-half sites instead
-  #sites = siteinds("S=1/2", N)
+using ConstructionBase: constructorof
+using GPUArraysCore: AbstractGPUArray
+using TensorAlgebra: TensorAlgebra, factorize
+function TensorAlgebra.factorize(
+  a::AbstractGPUArray, labels, codomain_labels, domain_labels; kwargs...
+)
+  x, y = factorize(Array(a), labels, codomain_labels, domain_labels; kwargs...)
+  return constructorof(typeof(a)).((x, y))
+end
 
-  # Input operator terms which define a Hamiltonian
+function heisenberg(N)
   os = OpSum()
   for j in 1:(N - 1)
     os += "Sz", j, "Sz", j + 1
     os += 0.5, "S+", j, "S-", j + 1
     os += 0.5, "S-", j, "S+", j + 1
   end
-  # Convert these terms to an MPO tensor network
-  H = MPO(os, sites)
-
-  # Create an initial random matrix product state
-  psi0 = random_mps(sites; linkdims=10)
-
-  # Plan to do 5 DMRG sweeps:
-  nsweeps = 5
-  # Set maximum MPS bond dimensions for each sweep
-  maxdim = [10, 20, 100, 100, 200]
-  # Set maximum truncation error allowed when adapting bond dimensions
-  cutoff = [1E-11]
-
-  # Run the DMRG algorithm, returning energy and optimized MPS
-  energy, psi = dmrg(H, psi0; nsweeps, maxdim, cutoff)
-  @printf("Final energy = %.12f\n", energy)
+  return os
 end
+
+Random.seed!(1234)
+
+N = 10
+
+# Create N spin-one degrees of freedom
+sites = siteinds("S=1", N)
+# Alternatively can make spin-half sites instead
+# sites = siteinds("S=1/2", N)
+
+os = heisenberg(N)
+
+# Input operator terms which define a Hamiltonian
+# Convert these terms to an MPO tensor network
+H = dev(MPO(os, sites))
+
+# Create an initial random matrix product state
+psi0 = dev(random_mps(sites; maxdim=10))
+# psi0 = random_mps(j -> isodd(j) ? "↑" : "↓", sites; maxdim=10)
+# psi0 = MPS(j -> isodd(j) ? "↑" : "↓", sites)
+
+# Plan to do 5 DMRG sweeps:
+nsweeps = 5
+# Set maximum MPS bond dimensions for each sweep
+maxdim = [10]
+# Set maximum truncation error allowed when adapting bond dimensions
+cutoff = [1E-11]
+
+# Run the DMRG algorithm, returning energy and optimized MPS
+energy, psi = dmrg(H, psi0; nsweeps, maxdim, cutoff)
+@show inner(psi', H, psi)
+@show inner(psi, psi)
+@printf("Final energy = %.12f\n", energy)
diff --git a/examples/dmrg/1d_heisenberg_conserve_spin.jl b/examples/dmrg/1d_heisenberg_conserve_spin.jl
index f8ec1fb..6b2ec08 100644
--- a/examples/dmrg/1d_heisenberg_conserve_spin.jl
+++ b/examples/dmrg/1d_heisenberg_conserve_spin.jl
@@ -1,39 +1,44 @@
+using BlockSparseArrays
+using GradedArrays
 using ITensors, ITensorMPS
-using LinearAlgebra
 using Printf
 using Random
-using Strided
-
-Random.seed!(1234)
-BLAS.set_num_threads(1)
-Strided.set_num_threads(1)
-ITensors.enable_threaded_blocksparse()
-#ITensors.disable_threaded_blocksparse()
-
-let
-  N = 100
-
-  sites = siteinds("S=1", N; conserve_qns=true)
 
+function heisenberg(N)
   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
-  H = MPO(os, sites)
+  return os
+end
+
+Random.seed!(1234)
 
-  state = [isodd(n) ? "Up" : "Dn" for n in 1:N]
-  psi0 = random_mps(sites, state; linkdims=10)
+N = 10
 
-  # Plan to do 5 DMRG sweeps:
-  nsweeps = 5
-  # Set maximum MPS bond dimensions for each sweep
-  maxdim = [10, 20, 100, 100, 200]
-  # Set maximum truncation error allowed when adapting bond dimensions
-  cutoff = [1E-11]
+s = siteinds("S=1/2", N; gradings=("Sz",))
 
-  # Run the DMRG algorithm, returning energy and optimized MPS
-  energy, psi = dmrg(H, psi0; nsweeps, maxdim, cutoff)
-  @printf("Final energy = %.12f\n", energy)
-end
+os = heisenberg(N)
+
+# Create an initial random matrix product state
+psi0 = random_mps(j -> isodd(j) ? "↑" : "↓", s; maxdim=10)
+# psi0 = MPS(j -> isodd(j) ? "↑" : "↓", s)
+
+# Input operator terms which define a Hamiltonian
+# Convert these terms to an MPO tensor network
+H = MPO(os, s)
+
+# Plan to do 5 DMRG sweeps:
+nsweeps = 5
+# Set maximum MPS bond dimensions for each sweep
+maxdim = [10]
+# Set maximum truncation error allowed when adapting bond dimensions
+cutoff = [1E-11]
+
+# Run the DMRG algorithm, returning energy and optimized MPS
+energy, psi = dmrg(H, psi0; nsweeps, maxdim, cutoff)
+@show inner(psi', H, psi)
+@show inner(psi, psi)
+@printf("Final energy = %.12f\n", energy)
diff --git a/examples/dmrg/Project.toml b/examples/dmrg/Project.toml
index e6b9177..0ca1a4d 100644
--- a/examples/dmrg/Project.toml
+++ b/examples/dmrg/Project.toml
@@ -1,6 +1,18 @@
 [deps]
+Adapt = "79e6a3ab-5dfb-504d-930d-738a2a938a0e"
+BlockSparseArrays = "2c9a651f-6452-4ace-a6ac-809f4280fbb4"
+ConstructionBase = "187b0558-2788-49d3-abe0-74a17ed4e7c9"
+GPUArraysCore = "46192b85-c4d5-4398-a991-12ede77f4527"
+GradedArrays = "bc96ca6e-b7c8-4bb6-888e-c93f838762c2"
+ITensorBase = "4795dd04-0d67-49bb-8f44-b89c448a1dc7"
 ITensorMPS = "0d1a4710-d33b-49a5-8f18-73bdf49b47e2"
 ITensors = "9136182c-28ba-11e9-034c-db9fb085ebd5"
+JLArrays = "27aeb0d3-9eb9-45fb-866b-73c2ecf80fcb"
 LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
+MatrixAlgebraKit = "6c742aac-3347-4629-af66-fc926824e5e4"
+Metal = "dde4c033-4e86-420c-a63e-0dd931031962"
+NamedDimsArrays = "60cbd0c0-df58-4cb7-918c-6f5607b73fde"
 Random = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c"
+SerializedArrays = "621c0da3-e96e-4f80-bd06-5ae31cdfcb39"
 Strided = "5e0ebb24-38b0-5f93-81fe-25c709ecae67"
+TensorAlgebra = "68bd88dc-f39d-4e12-b2ca-f046b68fcc6a"
diff --git a/ext/ITensorMPSChainRulesCoreExt/abstractmps.jl b/ext/ITensorMPSChainRulesCoreExt/abstractmps.jl
index 09f396c..972cca5 100644
--- a/ext/ITensorMPSChainRulesCoreExt/abstractmps.jl
+++ b/ext/ITensorMPSChainRulesCoreExt/abstractmps.jl
@@ -3,7 +3,7 @@ using ChainRulesCore: ChainRulesCore, HasReverseMode, NoTangent, RuleConfig, rru
 using ITensors:
   ITensors, ITensor, dag, hassameinds, inds, itensor, mapprime, replaceprime, swapprime
 using ITensorMPS: ITensorMPS, MPO, MPS, apply, inner, siteinds
-using NDTensors: datatype
+## using NDTensors: datatype
 
 function ChainRulesCore.rrule(
   ::Type{T}, x::Vector{<:ITensor}; kwargs...
diff --git a/ext/ITensorMPSPackageCompilerExt/compile.jl b/ext/ITensorMPSPackageCompilerExt/compile.jl
index c14d0e9..58f3452 100644
--- a/ext/ITensorMPSPackageCompilerExt/compile.jl
+++ b/ext/ITensorMPSPackageCompilerExt/compile.jl
@@ -1,4 +1,4 @@
-using NDTensors: @Algorithm_str
+using BackendSelection: @Algorithm_str
 using ITensors: ITensors
 using PackageCompiler: PackageCompiler
 
diff --git a/src/ITensorMPS.jl b/src/ITensorMPS.jl
index 1c42b83..e044803 100644
--- a/src/ITensorMPS.jl
+++ b/src/ITensorMPS.jl
@@ -2,6 +2,7 @@ module ITensorMPS
 using ITensors
 include("imports.jl")
 include("exports.jl")
+include("siteinds.jl")
 include("abstractmps.jl")
 include("mps.jl")
 include("mpo.jl")
@@ -20,7 +21,6 @@ include("opsum_to_mpo/qnmatelem.jl")
 include("opsum_to_mpo/opsum_to_mpo_generic.jl")
 include("opsum_to_mpo/opsum_to_mpo.jl")
 include("opsum_to_mpo/opsum_to_mpo_qn.jl")
-include("deprecated.jl")
 include("defaults.jl")
 include("update_observer.jl")
 include("lattices/lattices.jl")
diff --git a/src/abstractmps.jl b/src/abstractmps.jl
index d838ee1..25b5780 100644
--- a/src/abstractmps.jl
+++ b/src/abstractmps.jl
@@ -1,10 +1,24 @@
+using BackendSelection: @Algorithm_str, Algorithm
+using GradedArrays: GradedArrays, dag
 using IsApprox: Approx, IsApprox
-using ITensors: ITensors
-using NDTensors: NDTensors, using_auto_fermion, scalartype, tensor
-using ITensors.Ops: Prod
-using ITensors.QuantumNumbers: QuantumNumbers, removeqn
-using ITensors.SiteTypes: SiteTypes, siteinds
-using ITensors.TagSets: TagSets
+using ITensors:
+  ITensors,
+  ITensor,
+  Index,
+  commonind,
+  commoninds,
+  factorize,
+  hasqns,
+  replaceinds,
+  sim,
+  tags,
+  uniqueind,
+  uniqueinds
+using QuantumOperatorAlgebra: Prod
+## using ITensors.QuantumNumbers: QuantumNumbers, removeqn
+## using ITensors.TagSets: TagSets
+using LinearAlgebra: LinearAlgebra
+using VectorInterface: VectorInterface, scalartype
 
 abstract type AbstractMPS end
 
@@ -53,9 +67,11 @@ have type `ComplexF64`, return `ComplexF64`.
 """
 promote_itensor_eltype(m::AbstractMPS) = LinearAlgebra.promote_leaf_eltypes(m)
 
-NDTensors.scalartype(m::AbstractMPS) = LinearAlgebra.promote_leaf_eltypes(m)
-NDTensors.scalartype(m::Array{ITensor}) = LinearAlgebra.promote_leaf_eltypes(m)
-NDTensors.scalartype(m::Array{<:Array{ITensor}}) = LinearAlgebra.promote_leaf_eltypes(m)
+VectorInterface.scalartype(m::AbstractMPS) = LinearAlgebra.promote_leaf_eltypes(m)
+
+## TODO: Are these definitions needed?
+## VectorInterface.scalartype(m::Array{ITensor}) = LinearAlgebra.promote_leaf_eltypes(m)
+## VectorInterface.scalartype(m::Array{<:Array{ITensor}}) = LinearAlgebra.promote_leaf_eltypes(m)
 
 """
     eltype(m::MPS)
@@ -74,7 +90,7 @@ Base.imag(ψ::AbstractMPS) = imag.(ψ)
 Base.conj(ψ::AbstractMPS) = conj.(ψ)
 
 function convert_leaf_eltype(eltype::Type, ψ::AbstractMPS)
-  return map(ψᵢ -> convert_leaf_eltype(eltype, ψᵢ), ψ; set_limits=false)
+  return map(ψᵢ -> eltype.(ψᵢ), ψ; set_limits=false)
 end
 
 """
@@ -351,7 +367,11 @@ If there is no link Index, return `nothing`.
 function linkind(M::AbstractMPS, j::Integer)
   N = length(M)
   (j ≥ length(M) || j < 1) && return nothing
-  return commonind(M[j], M[j + 1])
+  is = commoninds(M[j], M[j + 1])
+  if iszero(length(is))
+    return nothing
+  end
+  return only(is)
 end
 
 """
@@ -363,7 +383,7 @@ MPS or MPO tensor on site j to site j+1.
 """
 function linkinds(M::AbstractMPS, j::Integer)
   N = length(M)
-  (j ≥ length(M) || j < 1) && return IndexSet()
+  (j ≥ length(M) || j < 1) && return []
   return commoninds(M[j], M[j + 1])
 end
 
@@ -468,7 +488,7 @@ end
 
 Get the site index (or indices) of MPO `A` that is unique to `A` (not shared with MPS/MPO `B`).
 """
-function SiteTypes.siteinds(
+function siteinds(
   f::Union{typeof(uniqueinds),typeof(uniqueind)},
   A::AbstractMPS,
   B::AbstractMPS,
@@ -488,7 +508,7 @@ end
 
 Get the site indices of MPO `A` that are unique to `A` (not shared with MPS/MPO `B`), as a `Vector{<:Index}`.
 """
-function SiteTypes.siteinds(
+function siteinds(
   f::Union{typeof(uniqueinds),typeof(uniqueind)}, A::AbstractMPS, B::AbstractMPS; kwargs...
 )
   return [siteinds(f, A, B, j; kwargs...) for j in eachindex(A)]
@@ -500,7 +520,7 @@ end
 
 Get the site index (or indices) of  the `j`th MPO tensor of `A` that is shared with MPS/MPO `B`.
 """
-function SiteTypes.siteinds(
+function siteinds(
   f::Union{typeof(commoninds),typeof(commonind)},
   A::AbstractMPS,
   B::AbstractMPS,
@@ -516,7 +536,7 @@ end
 
 Get a vector of the site index (or indices) of MPO `A` that is shared with MPS/MPO `B`.
 """
-function SiteTypes.siteinds(
+function siteinds(
   f::Union{typeof(commoninds),typeof(commonind)}, A::AbstractMPS, B::AbstractMPS; kwargs...
 )
   return [siteinds(f, A, B, j) for j in eachindex(A)]
@@ -584,7 +604,7 @@ findsites(M, (s[4]', s[4])) == [4]
 findsites(M, (s[4]', s[3])) == [3, 4]
 ```
 """
-findsites(ψ::AbstractMPS, is) = findall(hascommoninds(is), ψ)
+findsites(ψ::AbstractMPS, is) = findall(ψᵢ -> !isdisjoint(inds(ψᵢ), inds(is)), ψ)
 
 findsites(ψ::AbstractMPS, s::Index) = findsites(ψ, IndexSet(s))
 
@@ -601,7 +621,7 @@ findsites(ψ::AbstractMPS, s::Index) = findsites(ψ, IndexSet(s))
 Return the first site of the MPS or MPO that has the
 site index `i`.
 """
-findfirstsiteind(ψ::AbstractMPS, s::Index) = findfirst(hasind(s), ψ)
+findfirstsiteind(ψ::AbstractMPS, s::Index) = findfirst(ψᵢ -> s ∈ inds(ψᵢ), ψ)
 
 # TODO: depracate in favor of findsite.
 """
@@ -611,7 +631,7 @@ findfirstsiteind(ψ::AbstractMPS, s::Index) = findfirst(hasind(s), ψ)
 Return the first site of the MPS or MPO that has the
 site indices `is`.
 """
-findfirstsiteinds(ψ::AbstractMPS, s) = findfirst(hasinds(s), ψ)
+findfirstsiteinds(ψ::AbstractMPS, s) = findfirst(ψᵢ -> s ⊆ inds(ψᵢ), ψ)
 
 """
     siteind(::typeof(first), M::Union{MPS,MPO}, j::Integer; kwargs...)
@@ -622,7 +642,7 @@ Return the first site Index found on the MPS or MPO
 You can choose different filters, like prime level
 and tags, with the `kwargs`.
 """
-function SiteTypes.siteind(::typeof(first), M::AbstractMPS, j::Integer; kwargs...)
+function siteind(::typeof(first), M::AbstractMPS, j::Integer; kwargs...)
   N = length(M)
   (N == 1) && return firstind(M[1]; kwargs...)
   if j == 1
@@ -644,7 +664,7 @@ at the site `j` as an IndexSet.
 Optionally filter prime tags and prime levels with
 keyword arguments like `plev` and `tags`.
 """
-function SiteTypes.siteinds(M::AbstractMPS, j::Integer; kwargs...)
+function siteinds(M::AbstractMPS, j::Integer; kwargs...)
   N = length(M)
   (N == 1) && return inds(M[1]; kwargs...)
   if j == 1
@@ -657,19 +677,19 @@ function SiteTypes.siteinds(M::AbstractMPS, j::Integer; kwargs...)
   return si
 end
 
-function SiteTypes.siteinds(::typeof(all), ψ::AbstractMPS, n::Integer; kwargs...)
+function siteinds(::typeof(all), ψ::AbstractMPS, n::Integer; kwargs...)
   return siteinds(ψ, n; kwargs...)
 end
 
-function SiteTypes.siteinds(::typeof(first), ψ::AbstractMPS; kwargs...)
+function siteinds(::typeof(first), ψ::AbstractMPS; kwargs...)
   return [siteind(first, ψ, j; kwargs...) for j in 1:length(ψ)]
 end
 
-function SiteTypes.siteinds(::typeof(only), ψ::AbstractMPS; kwargs...)
+function siteinds(::typeof(only), ψ::AbstractMPS; kwargs...)
   return [siteind(only, ψ, j; kwargs...) for j in 1:length(ψ)]
 end
 
-function SiteTypes.siteinds(::typeof(all), ψ::AbstractMPS; kwargs...)
+function siteinds(::typeof(all), ψ::AbstractMPS; kwargs...)
   return [siteinds(ψ, j; kwargs...) for j in 1:length(ψ)]
 end
 
@@ -687,12 +707,12 @@ function replaceinds!(::typeof(linkinds), M::AbstractMPS, l̃s::Vector{<:Index})
   return M
 end
 
-function replaceinds(::typeof(linkinds), M::AbstractMPS, l̃s::Vector{<:Index})
+function ITensors.replaceinds(::typeof(linkinds), M::AbstractMPS, l̃s::Vector{<:Index})
   return replaceinds!(linkinds, copy(M), l̃s)
 end
 
 # TODO: change kwarg from `set_limits` to `preserve_ortho`
-function map!(f::Function, M::AbstractMPS; set_limits::Bool=true)
+function Base.map!(f::Function, M::AbstractMPS; set_limits::Bool=true)
   for i in eachindex(M)
     M[i, set_limits = set_limits] = f(M[i])
   end
@@ -700,21 +720,22 @@ function map!(f::Function, M::AbstractMPS; set_limits::Bool=true)
 end
 
 # TODO: change kwarg from `set_limits` to `preserve_ortho`
-function map(f::Function, M::AbstractMPS; set_limits::Bool=true)
+function Base.map(f::Function, M::AbstractMPS; set_limits::Bool=true)
   return map!(f, copy(M); set_limits=set_limits)
 end
 
 for (fname, fname!) in [
-  (:(ITensors.dag), :(dag!)),
-  (:(ITensors.prime), :(ITensors.prime!)),
-  (:(ITensors.setprime), :(ITensors.setprime!)),
-  (:(ITensors.noprime), :(ITensors.noprime!)),
-  (:(ITensors.swapprime), :(ITensors.swapprime!)),
-  (:(ITensors.replaceprime), :(ITensors.replaceprime!)),
-  (:(TagSets.addtags), :(ITensors.addtags!)),
-  (:(TagSets.removetags), :(ITensors.removetags!)),
-  (:(TagSets.replacetags), :(ITensors.replacetags!)),
-  (:(ITensors.settags), :(ITensors.settags!)),
+  (:(GradedArrays.dag), :(dag!)),
+  (:(ITensors.prime), :(prime!)),
+  (:(ITensors.setprime), :(setprime!)),
+  (:(ITensors.noprime), :(noprime!)),
+  (:(ITensors.swapprime), :(swapprime!)),
+  (:(ITensors.replaceprime), :(replaceprime!)),
+  # TODO: Delete these, use `settag` instead.
+  # (:(ITensors.addtags), :(addtags!)),
+  # (:(ITensors.removetags), :(removetags!)),
+  # (:(ITensors.replacetags), :(replacetags!)),
+  # (:(ITensors.settags), :(settags!)),
 ]
   @eval begin
     """
@@ -755,61 +776,61 @@ function check_hascommoninds(::typeof(siteinds), A::AbstractMPS, B::AbstractMPS)
     )
   end
   for n in 1:N
-    !hascommoninds(siteinds(A, n), siteinds(B, n)) && error(
+    isdisjoint(siteinds(A, n), siteinds(B, n)) && error(
       "$(typeof(A)) A and $(typeof(B)) B must share site indices. On site $n, A has site indices $(siteinds(A, n)) while B has site indices $(siteinds(B, n)).",
     )
   end
   return nothing
 end
 
-function map!(f::Function, ::typeof(linkinds), M::AbstractMPS)
+function Base.map!(f::Function, ::typeof(linkinds), M::AbstractMPS)
   for i in eachindex(M)[1:(end - 1)]
     l = linkinds(M, i)
     if !isempty(l)
-      l̃ = f(l)
+      l̃ = f.(l)
       @preserve_ortho M begin
-        M[i] = replaceinds(M[i], l, l̃)
-        M[i + 1] = replaceinds(M[i + 1], l, l̃)
+        M[i] = replaceinds(M[i], (l .=> l̃)...)
+        M[i + 1] = replaceinds(M[i + 1], (l .=> l̃)...)
       end
     end
   end
   return M
 end
 
-map(f::Function, ::typeof(linkinds), M::AbstractMPS) = map!(f, linkinds, copy(M))
+Base.map(f::Function, ::typeof(linkinds), M::AbstractMPS) = map!(f, linkinds, copy(M))
 
-function map!(f::Function, ::typeof(siteinds), M::AbstractMPS)
+function Base.map!(f::Function, ::typeof(siteinds), M::AbstractMPS)
   for i in eachindex(M)
     s = siteinds(M, i)
     if !isempty(s)
       @preserve_ortho M begin
-        M[i] = replaceinds(M[i], s, f(s))
+        M[i] = replaceinds(M[i], s, f.(s))
       end
     end
   end
   return M
 end
 
-map(f::Function, ::typeof(siteinds), M::AbstractMPS) = map!(f, siteinds, copy(M))
+Base.map(f::Function, ::typeof(siteinds), M::AbstractMPS) = map!(f, siteinds, copy(M))
 
-function map!(
+function Base.map!(
   f::Function, ::typeof(siteinds), ::typeof(commoninds), M1::AbstractMPS, M2::AbstractMPS
 )
   length(M1) != length(M2) && error("MPOs/MPSs must be the same length")
   for i in eachindex(M1)
     s = siteinds(commoninds, M1, M2, i)
     if !isempty(s)
-      s̃ = f(s)
+      s̃ = f.(s)
       @preserve_ortho (M1, M2) begin
-        M1[i] = replaceinds(M1[i], s .=> s̃)
-        M2[i] = replaceinds(M2[i], s .=> s̃)
+        M1[i] = replaceinds(M1[i], (s .=> s̃)...)
+        M2[i] = replaceinds(M2[i], (s .=> s̃)...)
       end
     end
   end
   return M1, M2
 end
 
-function map!(
+function Base.map!(
   f::Function, ::typeof(siteinds), ::typeof(uniqueinds), M1::AbstractMPS, M2::AbstractMPS
 )
   length(M1) != length(M2) && error("MPOs/MPSs must be the same length")
@@ -817,14 +838,14 @@ function map!(
     s = siteinds(uniqueinds, M1, M2, i)
     if !isempty(s)
       @preserve_ortho M1 begin
-        M1[i] = replaceinds(M1[i], s .=> f(s))
+        M1[i] = replaceinds(M1[i], (s .=> f(s))...)
       end
     end
   end
   return M1
 end
 
-function map(
+function Base.map(
   f::Function,
   ffilter::typeof(siteinds),
   fsubset::Union{typeof(commoninds),typeof(uniqueinds)},
@@ -837,7 +858,7 @@ end
 function hassameinds(::typeof(siteinds), M1::AbstractMPS, M2::AbstractMPS)
   length(M1) ≠ length(M2) && return false
   for n in 1:length(M1)
-    !hassameinds(siteinds(all, M1, n), siteinds(all, M2, n)) && return false
+    !issetequal(siteinds(all, M1, n), siteinds(all, M2, n)) && return false
   end
   return true
 end
@@ -851,14 +872,15 @@ function hassamenuminds(::typeof(siteinds), M1::AbstractMPS, M2::AbstractMPS)
 end
 
 for (fname, fname!) in [
-  (:(NDTensors.sim), :(sim!)),
-  (:(ITensors.prime), :(ITensors.prime!)),
-  (:(ITensors.setprime), :(ITensors.setprime!)),
-  (:(ITensors.noprime), :(ITensors.noprime!)),
-  (:(TagSets.addtags), :(ITensors.addtags!)),
-  (:(TagSets.removetags), :(ITensors.removetags!)),
-  (:(TagSets.replacetags), :(ITensors.replacetags!)),
-  (:(ITensors.settags), :(ITensors.settags!)),
+  (:(ITensors.sim), :(sim!)),
+  (:(ITensors.prime), :(prime!)),
+  (:(ITensors.setprime), :(setprime!)),
+  (:(ITensors.noprime), :(noprime!)),
+  # TODO: Delete these, use `settag` instead.
+  # (:(ITensors.addtags), :(addtags!)),
+  # (:(ITensors.removetags), :(removetags!)),
+  # (:(ITensors.replacetags), :(replacetags!)),
+  # (:(ITensors.settags), :(settags!)),
 ]
   @eval begin
     """
@@ -938,7 +960,7 @@ for (fname, fname!) in [
       args...;
       kwargs...,
     )
-      return map(i -> $fname(i, args...; kwargs...), ffilter, fsubset, M1, M2)
+      return map(i -> $fname.(i, args...; kwargs...), ffilter, fsubset, M1, M2)
     end
 
     function $(fname!)(
@@ -1089,7 +1111,7 @@ function deprecate_make_inds_match!(
     make_inds_match = false
   end
   if !hassameinds(siteinds, M1dag, M2) && make_inds_match
-    ITensors.warn_once(inner_mps_mpo_mps_deprecation_warning(), :inner_mps_mps)
+    #ITensors.warn_once(inner_mps_mpo_mps_deprecation_warning(), :inner_mps_mps)
     replace_siteinds!(M1dag, siteindsM2)
   end
   return M1dag, M2
@@ -1229,7 +1251,7 @@ the full inner product of the MPS/MPO with itself.
 
 See also [`lognorm`](@ref).
 """
-function norm(M::AbstractMPS)
+function LinearAlgebra.norm(M::AbstractMPS)
   if isortho(M)
     return norm(M[orthocenter(M)])
   end
@@ -1305,7 +1327,7 @@ of the orthogonality center to avoid numerical overflow in the case of diverging
 
 See also [`normalize!`](@ref), [`norm`](@ref), [`lognorm`](@ref).
 """
-function normalize(M::AbstractMPS; (lognorm!)=[])
+function LinearAlgebra.normalize(M::AbstractMPS; (lognorm!)=[])
   return normalize!(deepcopy_ortho_center(M); (lognorm!)=(lognorm!))
 end
 
@@ -1577,7 +1599,7 @@ truncation.
 - `cutoff::Real`: singular value truncation cutoff
 - `maxdim::Int`: maximum MPS/MPO bond dimension
 """
-function sum(ψ⃗::Vector{T}; kwargs...) where {T<:AbstractMPS}
+function Base.sum(ψ⃗::Vector{T}; kwargs...) where {T<:AbstractMPS}
   iszero(length(ψ⃗)) && return T()
   isone(length(ψ⃗)) && return copy(only(ψ⃗))
   return +(ψ⃗...; kwargs...)
@@ -1602,11 +1624,11 @@ out-of-place with `orthogonalize`.
 """
 function orthogonalize!(M::AbstractMPS, j::Int; maxdim=nothing, normalize=nothing)
   # TODO: Delete `maxdim` and `normalize` keyword arguments.
-  @debug_check begin
-    if !(1 <= j <= length(M))
-      error("Input j=$j to `orthogonalize!` out of range (valid range = 1:$(length(M)))")
-    end
-  end
+  ## @debug_check begin
+  ##   if !(1 <= j <= length(M))
+  ##     error("Input j=$j to `orthogonalize!` out of range (valid range = 1:$(length(M)))")
+  ##   end
+  ## end
   while leftlim(M) < (j - 1)
     (leftlim(M) < 0) && setleftlim!(M, 0)
     b = leftlim(M) + 1
@@ -1615,7 +1637,7 @@ function orthogonalize!(M::AbstractMPS, j::Int; maxdim=nothing, normalize=nothin
     if !isnothing(lb)
       ltags = tags(lb)
     else
-      ltags = TagSet("Link,l=$b")
+      ltags = ["l" => "$b"]
     end
     L, R = factorize(M[b], linds; tags=ltags, maxdim)
     M[b] = L
@@ -1625,9 +1647,7 @@ function orthogonalize!(M::AbstractMPS, j::Int; maxdim=nothing, normalize=nothin
       setrightlim!(M, leftlim(M) + 2)
     end
   end
-
   N = length(M)
-
   while rightlim(M) > (j + 1)
     (rightlim(M) > (N + 1)) && setrightlim!(M, N + 1)
     b = rightlim(M) - 2
@@ -1636,12 +1656,11 @@ function orthogonalize!(M::AbstractMPS, j::Int; maxdim=nothing, normalize=nothin
     if !isnothing(lb)
       ltags = tags(lb)
     else
-      ltags = TagSet("Link,l=$b")
+      ltags = ["l" => "$b"]
     end
     L, R = factorize(M[b + 1], rinds; tags=ltags, maxdim)
     M[b + 1] = L
     M[b] *= R
-
     setrightlim!(M, b + 1)
     if leftlim(M) > rightlim(M) - 2
       setleftlim!(M, rightlim(M) - 2)
@@ -1787,12 +1806,12 @@ end
 
 _isodd_fermionic_parity(s::Index, ::Integer) = false
 
-function _isodd_fermionic_parity(s::QNIndex, n::Integer)
-  qn_n = qn(space(s)[n])
-  fermionic_qn_pos = findfirst(q -> isfermionic(q), qn_n)
-  isnothing(fermionic_qn_pos) && return false
-  return isodd(val(qn_n[fermionic_qn_pos]))
-end
+## function _isodd_fermionic_parity(s::QNIndex, n::Integer)
+##   qn_n = qn(space(s)[n])
+##   fermionic_qn_pos = findfirst(q -> isfermionic(q), qn_n)
+##   isnothing(fermionic_qn_pos) && return false
+##   return isodd(val(qn_n[fermionic_qn_pos]))
+## end
 
 function _fermionic_swap(s1::Index, s2::Index)
   T = ITensor(QN(), s1', s2', dag(s1), dag(s2))
@@ -1858,57 +1877,55 @@ function setindex!(
 
   # Check that A has the proper common
   # indices with ψ
-  lind = linkind(ψ, firstsite - 1)
-  rind = linkind(ψ, lastsite)
-
+  linds = linkinds(ψ, firstsite - 1)
+  rinds = linkinds(ψ, lastsite)
   sites = [siteinds(ψ, j) for j in firstsite:lastsite]
 
-  #s = collect(Iterators.flatten(sites))
-  indsA = filter(x -> !isnothing(x), [lind, Iterators.flatten(sites)..., rind])
-  @assert hassameinds(A, indsA)
+  @assert issetequal(inds(A), [linds; reduce(vcat, sites); rinds])
 
   # For MPO case, restrict to 0 prime level
   #sites = filter(hasplev(0), sites)
 
-  if !isnothing(perm)
-    sites0 = sites
-    sites = sites0[[perm...]]
-    # Check if the site indices
-    # are fermionic
-    if !using_auto_fermion() && any(ITensors.anyfermionic, sites)
-      if length(sites) == 2 && ψ isa MPS
-        if all(ITensors.allfermionic, sites)
-          s0 = Index.(sites0)
-
-          # TODO: the Fermionic swap is could be diagonal,
-          # if we combine the site indices
-          #C = combiner(s0[1], s0[2])
-          #c = combinedind(C)
-          #AC = A * C
-          #AC = noprime(AC * _fermionic_swap(c))
-          #A = AC * dag(C)
-
-          FSWAP = adapt(datatype(A), _fermionic_swap(s0[1], s0[2]))
-          A = noprime(A * FSWAP)
-        end
-      elseif ψ isa MPO
-        @warn "In setindex!(MPO, ::ITensor, ::UnitRange), " *
-          "fermionic signs are only not handled properly for non-trivial " *
-          "permutations of sites. Please inform the developers of ITensors " *
-          "if you require this feature (otherwise, fermionic signs can be " *
-          "put in manually with fermionic swap gates)."
-      else
-        @warn "In setindex!(::Union{MPS, MPO}, ::ITensor, ::UnitRange), " *
-          "fermionic signs are only handled properly for permutations involving 2 sites. " *
-          "The original sites are $sites0, with a permutation $perm. " *
-          "To have the fermion sign handled correctly, we recommend performing your permutation " *
-          "pairwise."
-      end
-    end
-  end
+  ## TODO: Bring this package to handle fermions.
+  ## if !isnothing(perm)
+  ##   sites0 = sites
+  ##   sites = sites0[[perm...]]
+  ##   # Check if the site indices
+  ##   # are fermionic
+  ##   if !using_auto_fermion() && any(ITensors.anyfermionic, sites)
+  ##     if length(sites) == 2 && ψ isa MPS
+  ##       if all(ITensors.allfermionic, sites)
+  ##         s0 = Index.(sites0)
+
+  ##         # TODO: the Fermionic swap is could be diagonal,
+  ##         # if we combine the site indices
+  ##         #C = combiner(s0[1], s0[2])
+  ##         #c = combinedind(C)
+  ##         #AC = A * C
+  ##         #AC = noprime(AC * _fermionic_swap(c))
+  ##         #A = AC * dag(C)
+
+  ##         FSWAP = adapt(datatype(A), _fermionic_swap(s0[1], s0[2]))
+  ##         A = noprime(A * FSWAP)
+  ##       end
+  ##     elseif ψ isa MPO
+  ##       @warn "In setindex!(MPO, ::ITensor, ::UnitRange), " *
+  ##         "fermionic signs are only not handled properly for non-trivial " *
+  ##         "permutations of sites. Please inform the developers of ITensors " *
+  ##         "if you require this feature (otherwise, fermionic signs can be " *
+  ##         "put in manually with fermionic swap gates)."
+  ##     else
+  ##       @warn "In setindex!(::Union{MPS, MPO}, ::ITensor, ::UnitRange), " *
+  ##         "fermionic signs are only handled properly for permutations involving 2 sites. " *
+  ##         "The original sites are $sites0, with a permutation $perm. " *
+  ##         "To have the fermion sign handled correctly, we recommend performing your permutation " *
+  ##         "pairwise."
+  ##     end
+  ##   end
+  ## end
 
   # use the first link index if present, otherwise use the default tag
-  linktags = TagSet[
+  linktags = [
     (b=linkind(ψ, i); isnothing(b) ? defaultlinkindtags(i) : tags(b)) for
     i in firstsite:(lastsite - 1)
   ]
@@ -1916,7 +1933,7 @@ function setindex!(
   ψA = MPST(
     A,
     sites;
-    leftinds=lind,
+    leftinds=linds,
     orthocenter=orthocenter - first(r) + 1,
     tags=linktags,
     kwargs...,
@@ -1939,7 +1956,7 @@ replacesites!(ψ::AbstractMPS, args...; kwargs...) = setindex!(ψ, args...; kwar
 replacesites(ψ::AbstractMPS, args...; kwargs...) = setindex!(copy(ψ), args...; kwargs...)
 
 _number_inds(s::Index) = 1
-_number_inds(s::IndexSet) = length(s)
+_number_inds(s::Union{Vector{<:Index},Tuple{Vararg{Index}}}) = length(s)
 _number_inds(sites) = sum(_number_inds(s) for s in sites)
 
 """
@@ -1951,7 +1968,7 @@ by site according to the site indices `sites`.
 
 # Keywords
 
-- `leftinds = nothing`: optional left dangling indices. Indices that are not
+- `leftinds = []`: optional left dangling indices. Indices that are not
    in `sites` and `leftinds` will be dangling off of the right side of the MPS/MPO.
 - `orthocenter::Integer = length(sites)`: the desired final orthogonality
    center of the output MPS/MPO.
@@ -1972,10 +1989,9 @@ function (::Type{MPST})(
 ) where {MPST<:AbstractMPS}
   N = length(sites)
   for s in sites
-    @assert hasinds(A, s)
+    @assert s ⊆ inds(A)
   end
-  @assert isnothing(leftinds) || hasinds(A, leftinds)
-
+  @assert leftinds ⊆ inds(A)
   @assert 1 ≤ orthocenter ≤ N
 
   ψ = Vector{ITensor}(undef, N)
@@ -1985,12 +2001,9 @@ function (::Type{MPST})(
   # 1:orthocenter and reverse(orthocenter:N)
   # so the orthogonality center is set correctly.
   for n in 1:(N - 1)
-    Lis = IndexSet(sites[n])
-    if !isnothing(l)
-      Lis = unioninds(Lis, l)
-    end
-    L, R = factorize(Ã, Lis; kwargs..., tags=tags[n], ortho="left")
-    l = commonind(L, R)
+    Lis = sites[n] ∪ l
+    L, R = factorize(Ã, Lis; kwargs..., tags=["site" => "$n"], ortho="left")
+    l = commoninds(L, R)
     ψ[n] = L
     Ã = R
   end
@@ -2003,7 +2016,7 @@ function (::Type{MPST})(
 end
 
 function (::Type{MPST})(A::AbstractArray, sites; kwargs...) where {MPST<:AbstractMPS}
-  return MPST(itensor(A, sites...), sites; kwargs...)
+  return MPST(ITensor(A, sites...), sites; kwargs...)
 end
 
 """
@@ -2116,6 +2129,70 @@ function movesites(ψ::AbstractMPS, ns, ns′; kwargs...)
   return ψ
 end
 
+# TODO: Fix this, move to `ITensorBase.jl`, rewrite in terms of
+# a new `ITensorMap` type to generalize beyond prime levels.
+using ITensors: setprime
+function apply(A::ITensor, B::ITensor; apply_dag::Bool=false)
+  commonindsAB = filter(i -> iszero(plev(i)), commoninds(A, B))
+  isempty(commonindsAB) && error("In product, must have common indices with prime level 0.")
+  common_paired_indsA = filter(
+    i -> ((i ∈ commonindsAB) && (setprime(i, 1) ∈ inds(A))), inds(A)
+  )
+  common_paired_indsB = filter(
+    i -> ((i ∈ commonindsAB) && (setprime(i, 1) ∈ inds(B))), inds(B)
+  )
+
+  if !isempty(common_paired_indsA)
+    ## Old version: `commoninds_pairs = unioninds(common_paired_indsA, common_paired_indsA')`
+    commoninds_pairs = common_paired_indsA ∪ prime.(common_paired_indsA)
+  elseif !isempty(common_paired_indsB)
+    ## Old version: `commoninds_pairs = unioninds(common_paired_indsB, common_paired_indsB')`
+    commoninds_pairs = common_paired_indsB ∪ prime.(common_paired_indsB)
+  else
+    # vector-vector product
+    apply_dag && error("apply_dag not supported for vector-vector product")
+    return A * B
+  end
+  danglings_indsA = setdiff(inds(A), commoninds_pairs)
+  danglings_indsB = setdiff(inds(B), commoninds_pairs)
+  danglings_inds = danglings_indsA ∪ danglings_indsB
+  if issetequal(common_paired_indsA, common_paired_indsB)
+    # matrix-matrix product
+    ## Old version: `prime(A; inds=!danglings_inds)`
+    A′ = replaceinds(i -> i ∉ danglings_inds ? i' : i, A)
+    A′B = A′ * B
+    ## Old version: `AB = mapprime(A′B, 2 => 1; inds=!danglings_inds)`
+    AB = replaceinds(A′B) do i
+      if i ∉ danglings_inds && plev(i) == 2
+        return setprime(i, 1)
+      end
+      return i
+    end
+    if apply_dag
+      AB′ = prime(AB; inds=(!danglings_inds))
+      Adag = swapprime(dag(A), 0 => 1; inds=(!danglings_inds))
+      return mapprime(AB′ * Adag, 2 => 1; inds=(!danglings_inds))
+    end
+    return AB
+  elseif isempty(common_paired_indsA) && !isempty(common_paired_indsB)
+    # vector-matrix product
+    apply_dag && error("apply_dag not supported for matrix-vector product")
+    A′ = prime(A; inds=(!danglings_inds))
+    return A′ * B
+  elseif !isempty(common_paired_indsA) && isempty(common_paired_indsB)
+    # matrix-vector product
+    apply_dag && error("apply_dag not supported for vector-matrix product")
+    ## Old version: `replaceprime(A * B, 1 => 0; inds=!danglings_inds)`
+    AB = A * B
+    return replaceinds(AB) do i
+      if i ∉ danglings_inds && plev(i) == 1
+        return setprime(i, 0)
+      end
+      return i
+    end
+  end
+end
+
 """
     apply(o::ITensor, ψ::Union{MPS, MPO}, [ns::Vector{Int}]; kwargs...)
     product([...])
@@ -2140,7 +2217,7 @@ to false.
 - `move_sites_back::Bool = true`: after the ITensors are applied to the MPS or
    MPO, move the sites of the MPS or MPO back to their original locations.
 """
-function product(
+function apply(
   o::ITensor,
   ψ::AbstractMPS,
   ns=findsites(ψ, o);
@@ -2166,7 +2243,7 @@ function product(
   for n in 2:N
     ϕ *= ψ[ns′[n]]
   end
-  ϕ = product(o, ϕ; apply_dag=apply_dag)
+  ϕ = apply(o, ϕ; apply_dag=apply_dag)
   ψ[ns′[1]:ns′[end], kwargs...] = ϕ
   if move_sites_back
     # Move the sites back to their original positions
@@ -2237,7 +2314,7 @@ function ITensors.op(::OpName"CX", ::SiteType"S=1/2", s1::Index, s2::Index)
          0 1 0 0
          0 0 0 1
          0 0 1 0]
-  return itensor(mat, s2', s1', s2, s1)
+  return ITensor(mat, s2', s1', s2, s1)
 end
 
 os = [("CX", 1, 3), ("σz", 3)]
@@ -2273,7 +2350,7 @@ expτH = ops(os, s)
 ψτ = apply(expτH, ψ0)
 ```
 """
-function product(
+function apply(
   As::Vector{ITensor},
   ψ::AbstractMPS;
   move_sites_back_between_gates::Bool=true,
@@ -2282,7 +2359,7 @@ function product(
 )
   Aψ = ψ
   for A in As
-    Aψ = product(A, Aψ; move_sites_back=move_sites_back_between_gates, kwargs...)
+    Aψ = apply(A, Aψ; move_sites_back=move_sites_back_between_gates, kwargs...)
   end
   if !move_sites_back_between_gates && move_sites_back
     s = siteinds(Aψ)
@@ -2306,8 +2383,8 @@ end
 #
 # # U|ψ⟩ = Z₁X₁|ψ⟩
 # apply(U,
-function product(o::Prod{ITensor}, ψ::AbstractMPS; kwargs...)
-  return product(reverse(terms(o)), ψ; kwargs...)
+function apply(o::Prod{ITensor}, ψ::AbstractMPS; kwargs...)
+  return apply(reverse(terms(o)), ψ; kwargs...)
 end
 
 function (o::Prod{ITensor})(ψ::AbstractMPS; kwargs...)
@@ -2326,7 +2403,7 @@ end
 Return true if the MPS or MPO has
 tensors which carry quantum numbers.
 """
-hasqns(M::AbstractMPS) = any(hasqns, data(M))
+ITensors.hasqns(M::AbstractMPS) = any(hasqns, data(M))
 
 # Trait type version of hasqns
 # Note this is not inferrable, so hasqns would be preferred
@@ -2393,9 +2470,9 @@ function splitblocks(::typeof(linkinds), M::AbstractMPS; tol=0)
 end
 
 removeqns(M::AbstractMPS) = map(removeqns, M; set_limits=false)
-function QuantumNumbers.removeqn(M::AbstractMPS, qn_name::String)
-  return map(m -> removeqn(m, qn_name), M; set_limits=false)
-end
+## function QuantumNumbers.removeqn(M::AbstractMPS, qn_name::String)
+##   return map(m -> removeqn(m, qn_name), M; set_limits=false)
+## end
 
 #
 # Broadcasting
@@ -2442,7 +2519,7 @@ function Base.show(io::IO, M::AbstractMPS)
       println(io, "#undef")
     else
       A = M[i]
-      if order(A) != 0
+      if ndims(A) != 0
         println(io, "[$i] $(inds(A))")
       else
         println(io, "[$i] ITensor()")
diff --git a/src/abstractprojmpo/abstractprojmpo.jl b/src/abstractprojmpo/abstractprojmpo.jl
index 05aae2d..8dd5a87 100644
--- a/src/abstractprojmpo/abstractprojmpo.jl
+++ b/src/abstractprojmpo/abstractprojmpo.jl
@@ -1,3 +1,7 @@
+using FillArrays: Ones
+using GradedArrays: dag
+using ITensors: noprime
+
 abstract type AbstractProjMPO end
 
 copy(::AbstractProjMPO) = error("Not implemented")
@@ -23,18 +27,18 @@ the length of the MPO used to construct it
 """
 Base.length(P::AbstractProjMPO) = length(P.H)
 
-function lproj(P::AbstractProjMPO)::Union{ITensor,OneITensor}
-  (P.lpos <= 0) && return OneITensor()
+function lproj(P::AbstractProjMPO)
+  (P.lpos <= 0) && return ITensor(Ones{Bool}())
   return P.LR[P.lpos]
 end
 
-function rproj(P::AbstractProjMPO)::Union{ITensor,OneITensor}
-  (P.rpos >= length(P) + 1) && return OneITensor()
+function rproj(P::AbstractProjMPO)
+  (P.rpos >= length(P) + 1) && return ITensor(Ones{Bool}())
   return P.LR[P.rpos]
 end
 
-function ITensors.contract(P::AbstractProjMPO, v::ITensor)::ITensor
-  itensor_map = Union{ITensor,OneITensor}[lproj(P)]
+function ITensors.contract(P::AbstractProjMPO, v::ITensor)
+  itensor_map = [lproj(P)]
   append!(itensor_map, P.H[site_range(P)])
   push!(itensor_map, rproj(P))
 
@@ -69,7 +73,7 @@ shorthand for `product(P,v)`.
 """
 function product(P::AbstractProjMPO, v::ITensor)::ITensor
   Pv = contract(P, v)
-  if order(Pv) != order(v)
+  if ndims(Pv) != ndims(v)
     error(
       string(
         "The order of the ProjMPO-ITensor product P*v is not equal to the order of the ITensor v, ",
diff --git a/src/abstractprojmpo/diskprojmpo.jl b/src/abstractprojmpo/diskprojmpo.jl
index dd9d1df..f70d3b2 100644
--- a/src/abstractprojmpo/diskprojmpo.jl
+++ b/src/abstractprojmpo/diskprojmpo.jl
@@ -1,3 +1,6 @@
+using FillArrays: Ones
+using ITensors: ITensor
+using SerializedElementArrays: SerializedElementVector
 
 """
 A DiskProjMPO computes and stores the projection of an
@@ -26,10 +29,10 @@ mutable struct DiskProjMPO <: AbstractProjMPO
   rpos::Int
   nsite::Int
   H::MPO
-  LR::DiskVector{ITensor}
-  Lcache::Union{ITensor,OneITensor}
+  LR::SerializedElementVector{ITensor}
+  Lcache::ITensor
   lposcache::Union{Int,Nothing}
-  Rcache::Union{ITensor,OneITensor}
+  Rcache::ITensor
   rposcache::Union{Int,Nothing}
 end
 
@@ -59,9 +62,9 @@ function DiskProjMPO(H::MPO)
     2,
     H,
     disk(Vector{ITensor}(undef, length(H))),
-    OneITensor,
+    ITensor,
     nothing,
-    OneITensor,
+    ITensor,
     nothing,
   )
 end
@@ -84,8 +87,8 @@ disk(pm::DiskProjMPO; kwargs...) = pm
 # Special overload of lproj which uses the cached
 # version of the left projected MPO, and if the
 # cache doesn't exist it loads it from disk.
-function lproj(P::DiskProjMPO)::Union{ITensor,OneITensor}
-  (P.lpos <= 0) && return OneITensor()
+function lproj(P::DiskProjMPO)
+  (P.lpos <= 0) && return ITensor(Ones{Bool}())
   if (P.lpos ≠ P.lposcache) || (P.lpos == 1)
     # Need to update the cache
     P.Lcache = P.LR[P.lpos]
@@ -97,8 +100,8 @@ end
 # Special overload of rproj which uses the cached
 # version of the right projected MPO, and if the
 # cache doesn't exist it loads it from disk.
-function rproj(P::DiskProjMPO)::Union{ITensor,OneITensor}
-  (P.rpos >= length(P) + 1) && return OneITensor()
+function rproj(P::DiskProjMPO)
+  (P.rpos >= length(P) + 1) && return ITensor(Ones{Bool}())
   if (P.rpos ≠ P.rposcache) || (P.rpos == length(P))
     # Need to update the cache
     P.Rcache = P.LR[P.rpos]
diff --git a/src/abstractprojmpo/projmposum.jl b/src/abstractprojmpo/projmposum.jl
index 3d5d5d9..2b3372b 100644
--- a/src/abstractprojmpo/projmposum.jl
+++ b/src/abstractprojmpo/projmposum.jl
@@ -1,8 +1,8 @@
-using Compat: allequal
+using QuantumOperatorAlgebra: QuantumOperatorAlgebra
 
 abstract type AbstractSum end
 
-terms(sum::AbstractSum) = sum.terms
+QuantumOperatorAlgebra.terms(sum::AbstractSum) = sum.terms
 
 function set_terms(sum::AbstractSum, terms)
   return error("Please implement `set_terms` for the `AbstractSum` type `$(typeof(sum))`.")
diff --git a/src/abstractprojmpo/projmps.jl b/src/abstractprojmpo/projmps.jl
index 0ad60a5..5ff2785 100644
--- a/src/abstractprojmpo/projmps.jl
+++ b/src/abstractprojmpo/projmps.jl
@@ -1,3 +1,4 @@
+using GradedArrays: dag
 
 mutable struct ProjMPS
   lpos::Int
diff --git a/src/dmrg.jl b/src/dmrg.jl
index 19c3e6f..501e38b 100644
--- a/src/dmrg.jl
+++ b/src/dmrg.jl
@@ -1,10 +1,15 @@
 using Adapt: adapt
+using ITensors: ITensors, plev
 using KrylovKit: eigsolve
-using NDTensors: scalartype, timer
+using NamedDimsArrays: NamedDimsArrays, aligndims
 using Printf: @printf
 using TupleTools: TupleTools
+using VectorInterface: scalartype
 
-function permute(
+## TODO: Add this back?
+## using NDTensors: timer
+
+function NamedDimsArrays.aligndims(
   M::AbstractMPS, ::Tuple{typeof(linkind),typeof(siteinds),typeof(linkind)}
 )::typeof(M)
   M̃ = typeof(M)(length(M))
@@ -12,7 +17,7 @@ function permute(
     lₙ₋₁ = linkind(M, n - 1)
     lₙ = linkind(M, n)
     s⃗ₙ = TupleTools.sort(Tuple(siteinds(M, n)); by=plev)
-    M̃[n] = ITensors.permute(M[n], filter(!isnothing, (lₙ₋₁, s⃗ₙ..., lₙ)))
+    M̃[n] = aligndims(M[n], filter(!isnothing, (lₙ₋₁, s⃗ₙ..., lₙ)))
   end
   set_ortho_lims!(M̃, ortho_lims(M))
   return M̃
@@ -23,7 +28,7 @@ function dmrg(H::MPO, psi0::MPS, sweeps::Sweeps; kwargs...)
   check_hascommoninds(siteinds, H, psi0')
   # Permute the indices to have a better memory layout
   # and minimize permutations
-  H = permute(H, (linkind, siteinds, linkind))
+  H = aligndims(H, (linkind, siteinds, linkind))
   PH = ProjMPO(H)
   return dmrg(PH, psi0, sweeps; kwargs...)
 end
@@ -33,7 +38,7 @@ function dmrg(Hs::Vector{MPO}, psi0::MPS, sweeps::Sweeps; kwargs...)
     check_hascommoninds(siteinds, H, psi0)
     check_hascommoninds(siteinds, H, psi0')
   end
-  Hs .= permute.(Hs, Ref((linkind, siteinds, linkind)))
+  Hs .= aligndims.(Hs, Ref((linkind, siteinds, linkind)))
   PHS = ProjMPOSum(Hs)
   return dmrg(PHS, psi0, sweeps; kwargs...)
 end
@@ -44,8 +49,8 @@ function dmrg(H::MPO, Ms::Vector{MPS}, psi0::MPS, sweeps::Sweeps; weight=true, k
   for M in Ms
     check_hascommoninds(siteinds, M, psi0)
   end
-  H = permute(H, (linkind, siteinds, linkind))
-  Ms .= permute.(Ms, Ref((linkind, siteinds, linkind)))
+  H = aligndims(H, (linkind, siteinds, linkind))
+  Ms .= aligndims.(Ms, Ref((linkind, siteinds, linkind)))
   if weight <= 0
     error(
       "weight parameter should be > 0.0 in call to excited-state dmrg (value passed was weight=$weight)",
@@ -55,7 +60,8 @@ function dmrg(H::MPO, Ms::Vector{MPS}, psi0::MPS, sweeps::Sweeps; weight=true, k
   return dmrg(PMM, psi0, sweeps; kwargs...)
 end
 
-using NDTensors.TypeParameterAccessors: unwrap_array_type
+using TypeParameterAccessors: unwrap_array_type
+
 """
     dmrg(H::MPO, psi0::MPS; kwargs...)
     dmrg(H::MPO, psi0::MPS, sweeps::Sweeps; kwargs...)
@@ -179,12 +185,12 @@ function dmrg(
     )
   end
 
-  @debug_check begin
-    # Debug level checks
-    # Enable with ITensors.enable_debug_checks()
-    checkflux(psi0)
-    checkflux(PH)
-  end
+  ## @debug_check begin
+  ##   # Debug level checks
+  ##   # Enable with ITensors.enable_debug_checks()
+  ##   checkflux(psi0)
+  ##   checkflux(PH)
+  ## end
 
   psi = copy(psi0)
   N = length(psi)
@@ -217,87 +223,85 @@ function dmrg(
       end
 
       for (b, ha) in sweepnext(N)
-        @debug_check begin
-          checkflux(psi)
-          checkflux(PH)
-        end
-
-        @timeit_debug timer "dmrg: position!" begin
-          PH = position!(PH, psi, b)
-        end
-
-        @debug_check begin
-          checkflux(psi)
-          checkflux(PH)
-        end
-
-        @timeit_debug timer "dmrg: psi[b]*psi[b+1]" begin
-          phi = psi[b] * psi[b + 1]
-        end
-
-        @timeit_debug timer "dmrg: eigsolve" begin
-          vals, vecs = eigsolve(
-            PH,
-            phi,
-            1,
-            eigsolve_which_eigenvalue;
-            ishermitian,
-            tol=eigsolve_tol,
-            krylovdim=eigsolve_krylovdim,
-            maxiter=eigsolve_maxiter,
-            verbosity=eigsolve_verbosity,
-          )
-        end
+        ## TODO: Add back `@debug_check`, `checkflux`.
+        ## @debug_check begin
+        ##   checkflux(psi)
+        ##   checkflux(PH)
+        ## end
+
+        PH = position!(PH, psi, b)
+
+        ## TODO: Add back `@debug_check`, `checkflux`.
+        ## @debug_check begin
+        ##   checkflux(psi)
+        ##   checkflux(PH)
+        ## end
+
+        phi = psi[b] * psi[b + 1]
+
+        vals, vecs = eigsolve(
+          PH,
+          phi,
+          1,
+          eigsolve_which_eigenvalue;
+          ishermitian,
+          tol=eigsolve_tol,
+          krylovdim=eigsolve_krylovdim,
+          maxiter=eigsolve_maxiter,
+          verbosity=eigsolve_verbosity,
+        )
 
         energy = vals[1]
+
+        ## TODO: Bring back some version of this logic.
         ## Right now there is a conversion problem in CUDA.jl where `UnifiedMemory` Arrays are being converted
         ## into `DeviceMemory`. This conversion line is here temporarily to fix that problem when it arises
         ## Adapt is only called when using CUDA backend. CPU will work as implemented previously.
         ## TODO this might be the only place we really need iscu if its not fixed.
-        phi = if NDTensors.iscu(phi) && NDTensors.iscu(vecs[1])
-          adapt(ITensors.set_eltype(unwrap_array_type(phi), eltype(vecs[1])), vecs[1])
-        else
-          vecs[1]
-        end
+        ## phi = if NDTensors.iscu(phi) && NDTensors.iscu(vecs[1])
+        ##   adapt(ITensors.set_eltype(unwrap_array_type(phi), eltype(vecs[1])), vecs[1])
+        ## else
+        ##  vecs[1]
+        ## end
+
+        phi = vecs[1]
 
         ortho = ha == 1 ? "left" : "right"
 
         drho = nothing
         if noise(sweeps, sw) > 0
-          @timeit_debug timer "dmrg: noiseterm" begin
-            # Use noise term when determining new MPS basis.
-            # This is used to preserve the element type of the MPS.
-            elt = real(scalartype(psi))
-            drho = elt(noise(sweeps, sw)) * noiseterm(PH, phi, ortho)
-          end
+          # Use noise term when determining new MPS basis.
+          # This is used to preserve the element type of the MPS.
+          elt = real(scalartype(psi))
+          drho = elt(noise(sweeps, sw)) * noiseterm(PH, phi, ortho)
         end
 
-        @debug_check begin
-          checkflux(phi)
-        end
+        ## TODO: Add back `@debug_check`, `checkflux`.
+        ## @debug_check begin
+        ##   checkflux(phi)
+        ## end
 
-        @timeit_debug timer "dmrg: replacebond!" begin
-          spec = replacebond!(
-            PH,
-            psi,
-            b,
-            phi;
-            maxdim=maxdim(sweeps, sw),
-            mindim=mindim(sweeps, sw),
-            cutoff=cutoff(sweeps, sw),
-            eigen_perturbation=drho,
-            ortho,
-            normalize=true,
-            which_decomp,
-            svd_alg,
-          )
-        end
-        maxtruncerr = max(maxtruncerr, spec.truncerr)
+        spec = replacebond!(
+          PH,
+          psi,
+          b,
+          phi;
+          maxdim=maxdim(sweeps, sw),
+          mindim=mindim(sweeps, sw),
+          cutoff=cutoff(sweeps, sw),
+          eigen_perturbation=drho,
+          ortho,
+          normalize=true,
+          which_decomp,
+          svd_alg,
+        )
+        # maxtruncerr = max(maxtruncerr, spec.truncerr)
 
-        @debug_check begin
-          checkflux(psi)
-          checkflux(PH)
-        end
+        ## TODO: Add back `@debug_check`, `checkflux`.
+        ## @debug_check begin
+        ##   checkflux(psi)
+        ##   checkflux(PH)
+        ## end
 
         if outputlevel >= 2
           @printf("Sweep %d, half %d, bond (%d,%d) energy=%s\n", sw, ha, b, b + 1, energy)
diff --git a/src/exports.jl b/src/exports.jl
index 749a73c..417a90b 100644
--- a/src/exports.jl
+++ b/src/exports.jl
@@ -7,10 +7,9 @@ export
   @StateName_str,
   @TagType_str,
   @ValName_str,
-  Apply,
+  ## Apply,
   Op,
   OpName,
-  Ops,
   Prod,
   Scaled,
   SiteType,
diff --git a/src/imports.jl b/src/imports.jl
index d3a4b93..1083896 100644
--- a/src/imports.jl
+++ b/src/imports.jl
@@ -1,20 +1,13 @@
 # Primarily used to import names into the `ITensorMPS`
 # module from submodules or from `ITensors` so they can
 # be reexported.
-using ITensors.SiteTypes:
-  @OpName_str,
-  @SiteType_str,
-  @StateName_str,
-  @TagType_str,
-  @ValName_str,
-  OpName,
-  SiteType,
-  StateName,
-  TagType,
-  ValName,
-  ops
-using ITensors.Ops: Trotter
+# TODO: Delete these and expect users to load
+# `QuantumOperatorDefinitions`?
+using QuantumOperatorDefinitions:
+  @OpName_str, @SiteType_str, @StateName_str, OpName, SiteType, StateName
+using QuantumOperatorAlgebra: Trotter
 
+# TODO: Delete this, overload explicitly.
 import Base:
   # types
   Array,
@@ -85,6 +78,7 @@ import Base:
   # macros
   @propagate_inbounds
 
+# TODO: Delete this, overload explicitly.
 import Base.Broadcast:
   # types
   AbstractArrayStyle,
@@ -98,77 +92,8 @@ import Base.Broadcast:
   broadcastable,
   instantiate
 
-import ..ITensors.NDTensors:
-  Algorithm,
-  @Algorithm_str,
-  EmptyNumber,
-  _Tuple,
-  _NTuple,
-  blas_get_num_threads,
-  datatype,
-  dense,
-  diagind,
-  disable_auto_fermion,
-  double_precision,
-  eachblock,
-  eachdiagblock,
-  enable_auto_fermion,
-  fill!!,
-  randn!!,
-  permutedims,
-  permutedims!
-
-import ..ITensors:
-  AbstractRNG,
-  Apply,
-  apply,
-  argument,
-  Broadcasted,
-  @Algorithm_str,
-  checkflux,
-  convert_leaf_eltype,
-  commontags,
-  @debug_check,
-  dag,
-  data,
-  DefaultArrayStyle,
-  DiskVector,
-  flux,
-  hascommoninds,
-  hasqns,
-  hassameinds,
-  inner,
-  isfermionic,
-  maxdim,
-  mindim,
-  noprime,
-  noprime!,
-  norm,
-  normalize,
-  outer,
-  OneITensor,
-  permute,
-  prime,
-  prime!,
-  product,
-  QNIndex,
-  replaceinds,
-  replaceprime,
-  replacetags,
-  setprime,
-  sim,
-  site,
-  splitblocks,
-  store,
-  Style,
-  sum,
-  swapprime,
-  symmetrystyle,
-  terms,
-  @timeit_debug,
-  truncate!,
-  which_op
-
-import ..ITensors.Ops: params
+# TODO: Delete this, overload explicitly.
+import QuantumOperatorAlgebra: params
 
+# TODO: Delete this, overload explicitly.
 import SerializedElementArrays: disk
diff --git a/src/mpo.jl b/src/mpo.jl
index 9acd561..fb87e65 100644
--- a/src/mpo.jl
+++ b/src/mpo.jl
@@ -1,8 +1,12 @@
 using Adapt: adapt
+using BackendSelection: @Algorithm_str, Algorithm
+using DiagonalArrays: δ
+using GradedArrays: dag
 using LinearAlgebra: dot
-using Random: Random
-using ITensors.Ops: OpSum
-using ITensors.SiteTypes: SiteTypes, siteind, siteinds
+# TODO: Add this back?
+# using ITensors: outer
+using QuantumOperatorAlgebra: OpSum
+using Random: Random, AbstractRNG
 
 """
     MPO
@@ -42,7 +46,7 @@ function MPO(::Type{ElT}, sites::Vector{<:Index}) where {ElT<:Number}
     return MPO(v)
   end
   space_ii = all(hasqns, sites) ? [QN() => 1] : 1
-  l = [Index(space_ii, "Link,l=$ii") for ii in 1:(N - 1)]
+  l = [settag(Index(space_ii), "l", "$ii") for ii in 1:(N - 1)]
   for ii in eachindex(sites)
     s = sites[ii]
     if ii == 1
@@ -146,6 +150,8 @@ function outer_mps_mps_deprecation_warning()
   return "Calling `outer(ψ::MPS, ϕ::MPS)` for MPS `ψ` and `ϕ` with shared indices is deprecated. Currently, we automatically prime `ψ` to make sure the site indices don't clash, but that will no longer be the case in ITensors v0.4. To upgrade your code, call `outer(ψ', ϕ)`. Although the new interface seems less convenient, it will allow `outer` to accept more general outer products going forward, such as outer products where some indices are shared (a batched outer product) or outer products of MPS between site indices that aren't just related by a single prime level."
 end
 
+function outer end
+
 function deprecate_make_inds_unmatch(::typeof(outer), ψ::MPS, ϕ::MPS; kw...)
   if hassameinds(siteinds, ψ, ϕ)
     ITensors.warn_once(outer_mps_mps_deprecation_warning(), :outer_mps_mps)
@@ -241,7 +247,7 @@ end
 Get the first site Index of the MPO found, by
 default with prime level 0.
 """
-SiteTypes.siteind(M::MPO, j::Int; kwargs...) = siteind(first, M, j; plev=0, kwargs...)
+siteind(M::MPO, j::Int; kwargs...) = siteind(first, M, j; plev=0, kwargs...)
 
 # TODO: make this return the site indices that would have
 # been used to create the MPO? I.e.:
@@ -251,9 +257,9 @@ SiteTypes.siteind(M::MPO, j::Int; kwargs...) = siteind(first, M, j; plev=0, kwar
 
 Get a Vector of IndexSets of all the site indices of M.
 """
-SiteTypes.siteinds(M::MPO; kwargs...) = siteinds(all, M; kwargs...)
+siteinds(M::MPO; kwargs...) = siteinds(all, M; kwargs...)
 
-function SiteTypes.siteinds(Mψ::Tuple{MPO,MPS}, n::Int; kwargs...)
+function siteinds(Mψ::Tuple{MPO,MPS}, n::Int; kwargs...)
   return siteinds(uniqueinds, Mψ[1], Mψ[2], n; kwargs...)
 end
 
@@ -264,7 +270,7 @@ function nsites(Mψ::Tuple{MPO,MPS})
   return N
 end
 
-function SiteTypes.siteinds(Mψ::Tuple{MPO,MPS}; kwargs...)
+function siteinds(Mψ::Tuple{MPO,MPS}; kwargs...)
   return [siteinds(Mψ, n; kwargs...) for n in 1:nsites(Mψ)]
 end
 
@@ -282,7 +288,7 @@ function hassameinds(::typeof(siteinds), ψ::MPS, Hϕ::Tuple{MPO,MPS})
   N = length(ψ)
   @assert N == length(Hϕ[1]) == length(Hϕ[1])
   for n in 1:N
-    !hassameinds(siteinds(Hϕ, n), siteinds(ψ, n)) && return false
+    !issetequal(siteinds(Hϕ, n), siteinds(ψ, n)) && return false
   end
   return true
 end
@@ -449,9 +455,10 @@ Same as [`dot`](@ref).
 """
 inner(y::MPS, A::MPO, x::MPS; kwargs...) = dot(y, A, x; kwargs...)
 
-function inner(y::MPS, Ax::Apply{Tuple{MPO,MPS}})
-  return inner(y', Ax.args[1], Ax.args[2])
-end
+## TODO: Add this back?
+## function inner(y::MPS, Ax::Apply{Tuple{MPO,MPS}})
+##   return inner(y', Ax.args[1], Ax.args[2])
+## end
 
 """
     loginner(y::MPS, A::MPO, x::MPS)
@@ -540,7 +547,7 @@ function logdot(M1::MPO, M2::MPO; make_inds_match::Bool=false, kwargs...)
   return _log_or_not_dot(M1, M2, true; make_inds_match=make_inds_match)
 end
 
-function LinearAlgebra.tr(M::MPO; plev::Pair{Int,Int}=0 => 1, tags::Pair=ts"" => ts"")
+function LinearAlgebra.tr(M::MPO; plev::Pair{Int,Int}=0 => 1, tags::Pair="" => "")
   N = length(M)
   #
   # TODO: choose whether to contract or trace
@@ -551,10 +558,10 @@ function LinearAlgebra.tr(M::MPO; plev::Pair{Int,Int}=0 => 1, tags::Pair=ts"" =>
   #
   # So tracing first is better if d > √χ.
   #
-  L = tr(M[1]; plev=plev, tags=tags)
+  L = tr(M[1]; plev, tags)
   for j in 2:N
     L *= M[j]
-    L = tr(L; plev=plev, tags=tags)
+    L = tr(L; plev, tags)
   end
   return L
 end
@@ -608,13 +615,14 @@ end
 
 (A::MPO)(ψ::MPS; kwargs...) = apply(A, ψ; kwargs...)
 
-function Apply(A::MPO, ψ::MPS; kwargs...)
-  return ITensors.LazyApply.Applied(apply, (A, ψ), NamedTuple(kwargs))
-end
+## TODO: Decide what to do about this.
+## function Apply(A::MPO, ψ::MPS; kwargs...)
+##   return ITensors.LazyApply.Applied(apply, (A, ψ), NamedTuple(kwargs))
+## end
 
 function ITensors.contract(A::MPO, ψ::MPS; alg=nothing, method=alg, kwargs...)
   # TODO: Delete `method` since it is deprecated.
-  alg = NDTensors.replace_nothing(method, "densitymatrix")
+  alg = replace_nothing(method, "densitymatrix")
 
   # Keyword argument deprecations
   # TODO: Delete these.
diff --git a/src/mps.jl b/src/mps.jl
index de71318..ac7c7d1 100644
--- a/src/mps.jl
+++ b/src/mps.jl
@@ -1,7 +1,13 @@
 using Adapt: adapt
-using NDTensors: using_auto_fermion
-using Random: Random
-using ITensors.SiteTypes: SiteTypes, siteind, siteinds, state
+using GradedArrays: dual
+using ITensors: hasqns
+using LinearAlgebra: qr
+using QuantumOperatorDefinitions: QuantumOperatorDefinitions, state
+using Random: Random, AbstractRNG
+using SparseArraysBase: oneelement
+
+## TODO: Add this back.
+## using NDTensors: using_auto_fermion
 
 """
     MPS
@@ -44,334 +50,42 @@ function MPS(N::Int; ortho_lims::UnitRange=1:N)
   return MPS(Vector{ITensor}(undef, N); ortho_lims=ortho_lims)
 end
 
-"""
-    MPS([::Type{ElT} = Float64, ]sites; linkdims=1)
-
-Construct an MPS filled with Empty ITensors of type `ElT` from a collection of indices.
-
-Optionally specify the link dimension with the keyword argument `linkdims`, which by default is 1.
-
-In the future we may generalize `linkdims` to allow specifying each individual link dimension as a vector,
-and additionally allow specifying quantum numbers.
-"""
-function MPS(
-  ::Type{T}, sites::Vector{<:Index}; linkdims::Union{Integer,Vector{<:Integer}}=1
-) where {T<:Number}
-  _linkdims = _fill_linkdims(linkdims, sites)
-  N = length(sites)
-  v = Vector{ITensor}(undef, N)
-  if N == 1
-    v[1] = ITensor(T, sites[1])
-    return MPS(v)
-  end
-
-  spaces = if hasqns(sites)
-    [[QN() => _linkdims[j]] for j in 1:(N - 1)]
-  else
-    [_linkdims[j] for j in 1:(N - 1)]
-  end
-
-  l = [Index(spaces[ii], "Link,l=$ii") for ii in 1:(N - 1)]
-  for ii in eachindex(sites)
-    s = sites[ii]
-    if ii == 1
-      v[ii] = ITensor(T, l[ii], s)
-    elseif ii == N
-      v[ii] = ITensor(T, dag(l[ii - 1]), s)
-    else
-      v[ii] = ITensor(T, dag(l[ii - 1]), s, l[ii])
-    end
-  end
-  return MPS(v)
-end
-
-MPS(sites::Vector{<:Index}, args...; kwargs...) = MPS(Float64, sites, args...; kwargs...)
-
-function randomU(eltype::Type{<:Number}, s1::Index, s2::Index)
-  return randomU(Random.default_rng(), eltype, s1, s2)
-end
-
-function randomU(rng::AbstractRNG, eltype::Type{<:Number}, s1::Index, s2::Index)
-  if !hasqns(s1) && !hasqns(s2)
-    mdim = dim(s1) * dim(s2)
-    RM = randn(rng, eltype, mdim, mdim)
-    Q, _ = NDTensors.qr_positive(RM)
-    G = itensor(Q, dag(s1), dag(s2), s1', s2')
-  else
-    M = random_itensor(rng, eltype, QN(), s1', s2', dag(s1), dag(s2))
-    U, S, V = svd(M, (s1', s2'))
-    u = commonind(U, S)
-    v = commonind(S, V)
-    replaceind!(U, u, v)
-    G = U * V
-  end
-  return G
-end
-
-function randomizeMPS!(eltype::Type{<:Number}, M::MPS, sites::Vector{<:Index}, linkdims=1)
-  return randomizeMPS!(Random.default_rng(), eltype, M, sites, linkdims)
-end
-
-function randomizeMPS!(
-  rng::AbstractRNG, eltype::Type{<:Number}, M::MPS, sites::Vector{<:Index}, linkdims=1
+function apply_random_staircase_circuit(
+  rng::AbstractRNG, elt::Type, state::MPS; depth, kwargs...
 )
-  _linkdims = _fill_linkdims(linkdims, sites)
-  if isone(length(sites))
-    randn!(rng, M[1])
-    normalize!(M)
-    return M
-  end
-  N = length(sites)
-  c = div(N, 2)
-  max_pass = 100
-  for pass in 1:max_pass, half in 1:2
-    if half == 1
-      (db, brange) = (+1, 1:1:(N - 1))
-    else
-      (db, brange) = (-1, N:-1:2)
-    end
-    for b in brange
-      s1 = sites[b]
-      s2 = sites[b + db]
-      G = randomU(rng, eltype, s1, s2)
-      T = noprime(G * M[b] * M[b + db])
-      rinds = uniqueinds(M[b], M[b + db])
-
-      b_dim = half == 1 ? b : b + db
-      U, S, V = svd(T, rinds; maxdim=_linkdims[b_dim], utags="Link,l=$(b-1)")
-      M[b] = U
-      M[b + db] = S * V
-      M[b + db] /= norm(M[b + db])
-    end
-    if half == 2 && dim(commonind(M[c], M[c + 1])) >= _linkdims[c]
-      break
-    end
+  n = length(state)
+  layer = [(j - 1, j) for j in reverse(2:n)]
+  s = siteinds(state)
+  gate_layers = mapreduce(vcat, 1:depth) do _
+    # TODO: Pass `elt` and `rng` as kwargs to `op`, to be
+    # used as parameters in `OpName` by `QuantumOperaterDefinitions.jl`.
+    return map(((i, j),) -> op("RandomUnitary", (s[i], s[j]); rng, eltype=elt), layer)
   end
-  setleftlim!(M, 0)
-  setrightlim!(M, 2)
-  if dim(commonind(M[c], M[c + 1])) < _linkdims[c]
-    @warn "MPS center bond dimension is less than requested (you requested $(_linkdims[c]), but in practice it is $(dim(commonind(M[c], M[c + 1]))). This is likely due to technicalities of truncating quantum number sectors."
-  end
-end
-
-function randomCircuitMPS(
-  eltype::Type{<:Number}, sites::Vector{<:Index}, linkdims::Vector{<:Integer}; kwargs...
-)
-  return randomCircuitMPS(Random.default_rng(), eltype, sites, linkdims; kwargs...)
-end
-
-function randomCircuitMPS(
-  rng::AbstractRNG,
-  eltype::Type{<:Number},
-  sites::Vector{<:Index},
-  linkdims::Vector{<:Integer};
-  kwargs...,
-)
-  N = length(sites)
-  M = MPS(N)
-
-  if N == 1
-    M[1] = ITensor(randn(rng, eltype, dim(sites[1])), sites[1])
-    M[1] /= norm(M[1])
-    return M
-  end
-
-  l = Vector{Index}(undef, N)
-
-  d = dim(sites[N])
-  chi = min(linkdims[N - 1], d)
-  l[N - 1] = Index(chi, "Link,l=$(N-1)")
-  O = NDTensors.random_unitary(rng, eltype, chi, d)
-  M[N] = itensor(O, l[N - 1], sites[N])
-
-  for j in (N - 1):-1:2
-    chi *= dim(sites[j])
-    chi = min(linkdims[j - 1], chi)
-    l[j - 1] = Index(chi, "Link,l=$(j-1)")
-    O = NDTensors.random_unitary(rng, eltype, chi, dim(sites[j]) * dim(l[j]))
-    T = reshape(O, (chi, dim(sites[j]), dim(l[j])))
-    M[j] = itensor(T, l[j - 1], sites[j], l[j])
-  end
-
-  O = NDTensors.random_unitary(rng, eltype, 1, dim(sites[1]) * dim(l[1]))
-  l0 = Index(1, "Link,l=0")
-  T = reshape(O, (1, dim(sites[1]), dim(l[1])))
-  M[1] = itensor(T, l0, sites[1], l[1])
-  M[1] *= onehot(eltype, l0 => 1)
-
-  M.llim = 0
-  M.rlim = 2
-
-  return M
+  # TODO: Use `apply`, for some reason that can't be found right now.
+  return apply(gate_layers, state; kwargs...)
 end
 
-function randomCircuitMPS(sites::Vector{<:Index}, linkdims::Vector{<:Integer}; kwargs...)
-  return randomCircuitMPS(Random.default_rng(), sites, linkdims; kwargs...)
-end
-
-function randomCircuitMPS(
-  rng::AbstractRNG, sites::Vector{<:Index}, linkdims::Vector{<:Integer}; kwargs...
-)
-  return randomCircuitMPS(rng, Float64, sites, linkdims; kwargs...)
-end
-
-function _fill_linkdims(linkdims::Vector{<:Integer}, sites::Vector{<:Index})
-  @assert length(linkdims) == length(sites) - 1
-  return linkdims
-end
-
-function _fill_linkdims(linkdims::Integer, sites::Vector{<:Index})
-  return fill(linkdims, length(sites) - 1)
-end
-
-"""
-    random_mps(eltype::Type{<:Number}, sites::Vector{<:Index}; linkdims=1)
-
-Construct a random MPS with link dimension `linkdims` of
-type `eltype`.
-
-`linkdims` can also accept a `Vector{Int}` with
-`length(linkdims) == length(sites) - 1` for constructing an
-MPS with non-uniform bond dimension.
-"""
 function random_mps(
-  ::Type{ElT}, sites::Vector{<:Index}; linkdims::Union{Integer,Vector{<:Integer}}=1
-) where {ElT<:Number}
-  return random_mps(Random.default_rng(), ElT, sites; linkdims)
-end
-
-function random_mps(
-  rng::AbstractRNG,
-  ::Type{ElT},
-  sites::Vector{<:Index};
-  linkdims::Union{Integer,Vector{<:Integer}}=1,
-) where {ElT<:Number}
-  _linkdims = _fill_linkdims(linkdims, sites)
-  if any(hasqns, sites)
-    error("initial state required to use random_mps with QNs")
-  end
-
-  # For non-QN-conserving MPS, instantiate
-  # the random MPS directly as a circuit:
-  return randomCircuitMPS(rng, ElT, sites, _linkdims)
-end
-
-"""
-    random_mps(sites::Vector{<:Index}; linkdims=1)
-    random_mps(eltype::Type{<:Number}, sites::Vector{<:Index}; linkdims=1)
-
-Construct a random MPS with link dimension `linkdims` which by
-default has element type `Float64`.
-
-`linkdims` can also accept a `Vector{Int}` with
-`length(linkdims) == length(sites) - 1` for constructing an
-MPS with non-uniform bond dimension.
-"""
-function random_mps(sites::Vector{<:Index}; linkdims::Union{Integer,Vector{<:Integer}}=1)
-  return random_mps(Random.default_rng(), sites; linkdims)
-end
-
-function random_mps(
-  rng::AbstractRNG, sites::Vector{<:Index}; linkdims::Union{Integer,Vector{<:Integer}}=1
+  rng::AbstractRNG, elt::Type, state, sites::Vector{<:Index}; maxdim, kwargs...
 )
-  return random_mps(rng, Float64, sites; linkdims)
+  x = MPS(elt, state, sites)
+  depth = ceil(Int, log(minimum(Int ∘ length, sites), maxdim))
+  return apply_random_staircase_circuit(rng, elt, x; depth, maxdim, kwargs...)
 end
-
-function random_mps(
-  sites::Vector{<:Index}, state; linkdims::Union{Integer,Vector{<:Integer}}=1
-)
-  return random_mps(Random.default_rng(), sites, state; linkdims)
+function random_mps(rng::AbstractRNG, state, sites::Vector{<:Index}; kwargs...)
+  return random_mps(rng, Float64, state, sites; kwargs...)
 end
-
-function random_mps(
-  rng::AbstractRNG,
-  sites::Vector{<:Index},
-  state;
-  linkdims::Union{Integer,Vector{<:Integer}}=1,
-)
-  return random_mps(rng, Float64, sites, state; linkdims)
-end
-
-function random_mps(
-  eltype::Type{<:Number},
-  sites::Vector{<:Index},
-  state;
-  linkdims::Union{Integer,Vector{<:Integer}}=1,
-)
-  return random_mps(Random.default_rng(), eltype, sites, state; linkdims)
+function random_mps(elt::Type, state, sites::Vector{<:Index}; kwargs...)
+  return random_mps(Random.default_rng(), elt, state, sites; kwargs...)
 end
-
-function random_mps(
-  rng::AbstractRNG,
-  eltype::Type{<:Number},
-  sites::Vector{<:Index},
-  state;
-  linkdims::Union{Integer,Vector{<:Integer}}=1,
-)::MPS
-  M = MPS(eltype, sites, state)
-  if any(>(1), linkdims)
-    randomizeMPS!(rng, eltype, M, sites, linkdims)
-  end
-  return M
+function random_mps(state, sites::Vector{<:Index}; kwargs...)
+  return random_mps(Random.default_rng(), Float64, state, sites; kwargs...)
 end
-
-@doc """
-    random_mps(sites::Vector{<:Index}, state; linkdims=1)
-
-Construct a real, random MPS with link dimension `linkdims`,
-made by randomizing an initial product state specified by
-`state`. This version of `random_mps` is necessary when creating
-QN-conserving random MPS (consisting of QNITensors). The initial
-`state` array provided determines the total QN of the resulting
-random MPS.
-""" random_mps(::Vector{<:Index}, ::Any)
-
-"""
-    MPS(::Type{T<:Number}, ivals::Vector{<:Pair{<:Index}})
-
-Construct a product state MPS with element type `T` and
-nonzero values determined from the input IndexVals.
-"""
-function MPS(::Type{T}, ivals::Vector{<:Pair{<:Index}}) where {T<:Number}
-  N = length(ivals)
-  M = MPS(N)
-
-  if N == 1
-    M[1] = ITensor(T, ind(ivals[1]))
-    M[1][ivals[1]] = one(T)
-    return M
-  end
-
-  if hasqns(ind(ivals[1]))
-    lflux = QN()
-    for j in 1:(N - 1)
-      lflux += qn(ivals[j])
-    end
-    links = Vector{QNIndex}(undef, N - 1)
-    for j in (N - 1):-1:1
-      links[j] = dag(Index(lflux => 1; tags="Link,l=$j"))
-      lflux -= qn(ivals[j])
-    end
-  else
-    links = [Index(1, "Link,l=$n") for n in 1:(N - 1)]
-  end
-
-  M[1] = ITensor(T, ind(ivals[1]), links[1])
-  M[1][ivals[1], links[1] => 1] = one(T)
-  for n in 2:(N - 1)
-    s = ind(ivals[n])
-    M[n] = ITensor(T, dag(links[n - 1]), s, links[n])
-    M[n][links[n - 1] => 1, ivals[n], links[n] => 1] = one(T)
-  end
-  M[N] = ITensor(T, dag(links[N - 1]), ind(ivals[N]))
-  M[N][links[N - 1] => 1, ivals[N]] = one(T)
-
-  return M
+function random_mps(sites; kwargs...)
+  state = fill("0", length(sites))
+  return random_mps(Random.default_rng(), Float64, state, sites; kwargs...)
 end
 
-# For backwards compatibility
-const productMPS = MPS
-
 """
     MPS(ivals::Vector{<:Pair{<:Index}})
 
@@ -406,59 +120,23 @@ psi = MPS(ComplexF64, sites, states)
 phi = MPS(sites, "Up")
 ```
 """
-function MPS(eltype::Type{<:Number}, sites::Vector{<:Index}, states_)
-  if length(sites) != length(states_)
-    throw(DimensionMismatch("Number of sites and and initial vals don't match"))
-  end
-  N = length(states_)
-  M = MPS(N)
-
-  if N == 1
-    M[1] = state(sites[1], states_[1])
-    return convert_leaf_eltype(eltype, M)
-  end
-
-  states = [state(sites[j], states_[j]) for j in 1:N]
-
-  if hasqns(states[1])
-    lflux = QN()
-    for j in 1:(N - 1)
-      lflux += flux(states[j])
-    end
-    links = Vector{QNIndex}(undef, N - 1)
-    for j in (N - 1):-1:1
-      links[j] = dag(Index(lflux => 1; tags="Link,l=$j"))
-      lflux -= flux(states[j])
-    end
-  else
-    links = [Index(1; tags="Link,l=$n") for n in 1:N]
-  end
-
-  M[1] = ITensor(sites[1], links[1])
-  M[1] += states[1] * state(links[1], 1)
-  for n in 2:(N - 1)
-    M[n] = ITensor(dag(links[n - 1]), sites[n], links[n])
-    M[n] += state(dag(links[n - 1]), 1) * states[n] * state(links[n], 1)
-  end
-  M[N] = ITensor(dag(links[N - 1]), sites[N])
-  M[N] += state(dag(links[N - 1]), 1) * states[N]
-
-  return convert_leaf_eltype(eltype, M)
+function MPS(elt::Type{<:Number}, states_, sites::Vector{<:Index})
+  return MPS([state(elt, states_[j], sites[j]) for j in 1:length(sites)])
 end
 
 function MPS(
-  ::Type{T}, sites::Vector{<:Index}, state::Union{String,Integer}
+  ::Type{T}, state::Union{String,Integer}, sites::Vector{<:Index}
 ) where {T<:Number}
-  return MPS(T, sites, fill(state, length(sites)))
+  return MPS(T, fill(state, length(sites)), sites)
 end
 
-function MPS(::Type{T}, sites::Vector{<:Index}, states::Function) where {T<:Number}
+function MPS(::Type{T}, states::Function, sites::Vector{<:Index}) where {T<:Number}
   states_vec = [states(n) for n in 1:length(sites)]
-  return MPS(T, sites, states_vec)
+  return MPS(T, states_vec, sites)
 end
 
 """
-    MPS(sites::Vector{<:Index},states)
+    MPS(states, sites::Vector{<:Index})
 
 Construct a product state MPS having
 site indices `sites`, and which corresponds to the initial
@@ -476,21 +154,21 @@ states = [isodd(n) ? "Up" : "Dn" for n in 1:N]
 psi = MPS(sites, states)
 ```
 """
-MPS(sites::Vector{<:Index}, states) = MPS(Float64, sites, states)
+MPS(state, sites::Vector{<:Index}) = MPS(Float64, state, sites)
 
 """
     siteind(M::MPS, j::Int; kwargs...)
 
 Get the first site Index of the MPS. Return `nothing` if none is found.
 """
-SiteTypes.siteind(M::MPS, j::Int; kwargs...) = siteind(first, M, j; kwargs...)
+siteind(M::MPS, j::Int; kwargs...) = siteind(first, M, j; kwargs...)
 
 """
     siteind(::typeof(only), M::MPS, j::Int; kwargs...)
 
 Get the only site Index of the MPS. Return `nothing` if none is found.
 """
-function SiteTypes.siteind(::typeof(only), M::MPS, j::Int; kwargs...)
+function siteind(::typeof(only), M::MPS, j::Int; kwargs...)
   is = siteinds(M, j; kwargs...)
   if isempty(is)
     return nothing
@@ -512,7 +190,7 @@ Get a vector of the only site Index found on each tensor of the MPS. Errors if m
 
 Get a vector of the all site Indices found on each tensor of the MPS. Returns a Vector of IndexSets.
 """
-SiteTypes.siteinds(M::MPS; kwargs...) = siteinds(first, M; kwargs...)
+siteinds(M::MPS; kwargs...) = siteinds(first, M; kwargs...)
 
 function replace_siteinds!(M::MPS, sites)
   for j in eachindex(M)
@@ -550,9 +228,9 @@ function replacebond!(
   use_relative_cutoff=nothing,
   min_blockdim=nothing,
 )
-  normalize = NDTensors.replace_nothing(normalize, false)
-  swapsites = NDTensors.replace_nothing(swapsites, false)
-  ortho = NDTensors.replace_nothing(ortho, "left")
+  normalize = replace_nothing(normalize, false)
+  swapsites = replace_nothing(swapsites, false)
+  ortho = replace_nothing(ortho, "left")
 
   indsMb = inds(M[b])
   if swapsites
@@ -560,7 +238,8 @@ function replacebond!(
     sbp1 = siteind(M, b + 1)
     indsMb = replaceind(indsMb, sb, sbp1)
   end
-  L, R, spec = factorize(
+  indsMb = indsMb ∩ inds(phi)
+  L, R = factorize(
     phi,
     indsMb;
     mindim,
@@ -590,7 +269,7 @@ function replacebond!(
       "In replacebond!, got ortho = $ortho, only currently supports `left` and `right`."
     )
   end
-  return spec
+  return nothing
 end
 
 """
@@ -802,7 +481,7 @@ function correlation_matrix(
     L = ITensor(1.0)
   else
     lind = commonind(psi[start_site], psi[start_site - 1])
-    L = delta(dag(lind), lind')
+    L = delta(dual(lind), lind')
   end
   pL = start_site - 1
 
@@ -896,7 +575,7 @@ function correlation_matrix(
             Li21 *= oᵢ * dag(psi[pL21])'
           else
             sᵢ = siteind(psi, pL21)
-            Li21 *= prime(dag(si[pL21]), !sᵢ)
+            Li21 *= prime(dual(si[pL21]), !sᵢ)
           end
           Li21 *= psi[pL21]
         end
diff --git a/src/opsum_to_mpo/opsum_to_mpo.jl b/src/opsum_to_mpo/opsum_to_mpo.jl
index 114eec3..916d3e7 100644
--- a/src/opsum_to_mpo/opsum_to_mpo.jl
+++ b/src/opsum_to_mpo/opsum_to_mpo.jl
@@ -1,4 +1,96 @@
-using NDTensors: using_auto_fermion
+## using NDTensors: using_auto_fermion
+using ITensors: dim, inds, settag
+
+replace_nothing(::Nothing, y) = y
+replace_nothing(x, y) = x
+
+default_mindim(a) = true
+default_use_absolute_cutoff(a) = false
+default_use_relative_cutoff(a) = true
+
+function truncate!(
+  P::AbstractVector;
+  mindim=nothing,
+  maxdim=nothing,
+  cutoff=nothing,
+  use_absolute_cutoff=nothing,
+  use_relative_cutoff=nothing,
+)
+  mindim = replace_nothing(mindim, default_mindim(P))
+  maxdim = replace_nothing(maxdim, length(P))
+  cutoff = replace_nothing(cutoff, typemin(eltype(P)))
+  use_absolute_cutoff = replace_nothing(use_absolute_cutoff, default_use_absolute_cutoff(P))
+  use_relative_cutoff = replace_nothing(use_relative_cutoff, default_use_relative_cutoff(P))
+
+  origm = length(P)
+  docut = zero(eltype(P))
+
+  #if P[1] <= 0.0
+  #  P[1] = 0.0
+  #  resize!(P, 1)
+  #  return 0.0, 0.0
+  #end
+
+  if origm == 1
+    docut = abs(P[1]) / 2
+    return zero(eltype(P)), docut
+  end
+
+  s = sign(P[1])
+  s < 0 && (P .*= s)
+
+  #Zero out any negative weight
+  for n in origm:-1:1
+    (P[n] >= zero(eltype(P))) && break
+    P[n] = zero(eltype(P))
+  end
+
+  n = origm
+  truncerr = zero(eltype(P))
+  while n > maxdim
+    truncerr += P[n]
+    n -= 1
+  end
+
+  if use_absolute_cutoff
+    #Test if individual prob. weights fall below cutoff
+    #rather than using *sum* of discarded weights
+    while P[n] <= cutoff && n > mindim
+      truncerr += P[n]
+      n -= 1
+    end
+  else
+    scale = one(eltype(P))
+    if use_relative_cutoff
+      scale = sum(P)
+      (scale == zero(eltype(P))) && (scale = one(eltype(P)))
+    end
+
+    #Continue truncating until *sum* of discarded probability
+    #weight reaches cutoff reached (or m==mindim)
+    while (truncerr + P[n] <= cutoff * scale) && (n > mindim)
+      truncerr += P[n]
+      n -= 1
+    end
+
+    truncerr /= scale
+  end
+
+  if n < 1
+    n = 1
+  end
+
+  if n < origm
+    docut = (P[n] + P[n + 1]) / 2
+    if abs(P[n] - P[n + 1]) < eltype(P)(1e-3) * P[n]
+      docut += eltype(P)(1e-3) * P[n]
+    end
+  end
+
+  s < 0 && (P .*= s)
+  resize!(P, n)
+  return truncerr, docut
+end
 
 # `ValType::Type{<:Number}` is used instead of `ValType::Type` for efficiency, possibly due to increased method specialization.
 # See https://github.com/ITensor/ITensors.jl/pull/1183.
@@ -67,7 +159,7 @@ function svdMPO(
   end
 
   llinks = Vector{Index{Int}}(undef, N + 1)
-  llinks[1] = Index(2, "Link,l=0")
+  llinks[1] = settag(Index(2), "l", "0")
 
   H = MPO(sites)
 
@@ -79,7 +171,7 @@ function svdMPO(
     VR = Vs[n]
     tdim = isempty(rightmaps[n]) ? 0 : size(VR, 2)
 
-    llinks[n + 1] = Index(2 + tdim, "Link,l=$n")
+    llinks[n + 1] = settag(Index(2 + tdim), "l", "$n")
 
     ll = llinks[n]
     rl = llinks[n + 1]
@@ -114,7 +206,8 @@ function svdMPO(
         end
       end
 
-      T = itensor(M, ll, rl)
+      T = ITensor(M, ll, rl)
+
       H[n] += T * computeSiteProd(sites, argument(t))
     end
 
@@ -125,7 +218,7 @@ function svdMPO(
     idM = zeros(ValType, dim(ll), dim(rl))
     idM[1, 1] = 1.0
     idM[end, end] = 1.0
-    T = itensor(idM, ll, rl)
+    T = ITensor(idM, ll, rl)
     H[n] += T * computeSiteProd(sites, Prod([Op("Id", n)]))
   end
 
diff --git a/src/opsum_to_mpo/opsum_to_mpo_generic.jl b/src/opsum_to_mpo/opsum_to_mpo_generic.jl
index 915fd0e..97c94ed 100644
--- a/src/opsum_to_mpo/opsum_to_mpo_generic.jl
+++ b/src/opsum_to_mpo/opsum_to_mpo_generic.jl
@@ -1,6 +1,37 @@
-using NDTensors: using_auto_fermion
-using ITensors.Ops: Ops, Op, OpSum, Scaled, Sum, coefficient
-using ITensors.SiteTypes: has_fermion_string, op
+## using NDTensors: using_auto_fermion
+using QuantumOperatorAlgebra:
+  QuantumOperatorAlgebra,
+  Op,
+  OpSum,
+  Scaled,
+  Sum,
+  argument,
+  coefficient,
+  site,
+  sites,
+  terms,
+  which_op
+using QuantumOperatorDefinitions: has_fermion_string, op
+
+"""
+    parity_sign(P)
+
+Given an array or tuple of integers representing
+a permutation or a subset of a permutation,
+compute the parity sign defined as -1 for a
+permutation consisting of an odd number of swaps
+and +1 for an even number of swaps. This
+implementation uses an O(n^2) algorithm and is
+intended for small permutations only.
+"""
+function parity_sign(P)::Int
+  L = length(P)
+  s = +1
+  for i in 1:L, j in (i + 1):L
+    s *= sign(P[j] - P[i])
+  end
+  return s
+end
 
 # TODO: Deprecate.
 const AutoMPO = OpSum
@@ -64,10 +95,10 @@ end
 add!(os::OpSum, o::Op) = add!(os, Prod{Op}() * o)
 add!(os::OpSum, o::Scaled{C,Op}) where {C} = add!(os, Prod{Op}() * o)
 add!(os::OpSum, o::Prod{Op}) = add!(os, one(Float64) * o)
-add!(os::OpSum, o::Tuple) = add!(os, Ops.op_term(o))
+add!(os::OpSum, o::Tuple) = add!(os, QuantumOperatorAlgebra.op_term(o))
 add!(os::OpSum, a1::String, args...) = add!(os, (a1, args...))
 add!(os::OpSum, a1::Number, args...) = add!(os, (a1, args...))
-subtract!(os::OpSum, o::Tuple) = add!(os, -Ops.op_term(o))
+subtract!(os::OpSum, o::Tuple) = add!(os, -QuantumOperatorAlgebra.op_term(o))
 
 function isfermionic(t::Vector{Op}, sites)
   p = +1
@@ -110,7 +141,7 @@ function Base.copyto!(os, bc::Broadcast.Broadcasted{OpSumAddTermStyle,<:Any,type
 end
 
 # XXX: Create a new function name for this.
-isempty(op_qn::Pair{Vector{Op},QN}) = isempty(op_qn.first)
+## isempty(op_qn::Pair{Vector{Op},QN}) = isempty(op_qn.first)
 
 # the key type is Prod{Op} for the dense case
 # and is Pair{Prod{Op},QN} for the QN conserving case
@@ -134,10 +165,10 @@ end
 
 function computeSiteProd(sites, ops::Prod{Op})::ITensor
   i = only(site(ops[1]))
-  T = op(sites[i], which_op(ops[1]); params(ops[1])...)
+  T = op(which_op(ops[1]), sites[i]; params(ops[1])...)
   for j in 2:length(ops)
     (only(site(ops[j])) != i) && error("Mismatch of site number in computeSiteProd")
-    opj = op(sites[i], which_op(ops[j]); params(ops[j])...)
+    opj = op(which_op(ops[j]), sites[i]; params(ops[j])...)
     T = product(T, opj)
   end
   return T
@@ -167,7 +198,7 @@ function sorteachterm(os::OpSum, sites)
   os = copy(os)
 
   for (j, t) in enumerate(os)
-    if maximum(ITensors.sites(t)) > length(sites)
+    if maximum(QuantumOperatorAlgebra.sites(t)) > length(sites)
       error(
         "The OpSum contains a term $t that extends beyond the number of sites $(length(sites)).",
       )
@@ -188,7 +219,8 @@ function sorteachterm(os::OpSum, sites)
     t_parity = +1
     for n in reverse(1:Nt)
       site_n = only(site(t[n]))
-      if !using_auto_fermion() && (t_parity == -1) && (site_n < prevsite)
+      ## if !using_auto_fermion() && (t_parity == -1) && (site_n < prevsite)
+      if (t_parity == -1) && (site_n < prevsite)
         # Insert local piece of Jordan-Wigner string emanating
         # from fermionic operators to the right
         # (Remaining F operators will be put in by svdMPO)
@@ -212,7 +244,7 @@ function sorteachterm(os::OpSum, sites)
     filter!(!iszero, perm)
     # and account for anti-commuting, fermionic operators
     # during above sort; put resulting sign into coef
-    t *= ITensors.parity_sign(perm)
+    t *= parity_sign(perm)
     terms(os)[j] = t
   end
 
@@ -298,9 +330,10 @@ function MPO(os::OpSum, sites::Vector{<:Index}; splitblocks=true, kwargs...)::MP
     return qn_svdMPO(os, sites; kwargs...)
   end
   M = svdMPO(os, sites; kwargs...)
-  if splitblocks
-    M = ITensors.splitblocks(linkinds, M)
-  end
+  ## TODO: Add this back.
+  ## if splitblocks
+  ##   M = ITensors.splitblocks(linkinds, M)
+  ## end
   return M
 end
 
@@ -339,7 +372,7 @@ function MPO(
   return MPO(eltype, OpSum{C}() + o, s; kwargs...)
 end
 
-# Like `Ops.OpSumLike` but without `OpSum` included.
+# Like `QuantumOperatorAlgebra.OpSumLike` but without `OpSum` included.
 const OpSumLikeWithoutOpSum{C} = Union{
   Op,Scaled{C,Op},Sum{Op},Prod{Op},Scaled{C,Prod{Op}},Sum{Scaled{C,Op}}
 }
diff --git a/src/opsum_to_mpo/opsum_to_mpo_qn.jl b/src/opsum_to_mpo/opsum_to_mpo_qn.jl
index 25f75b3..1f9b98c 100644
--- a/src/opsum_to_mpo/opsum_to_mpo_qn.jl
+++ b/src/opsum_to_mpo/opsum_to_mpo_qn.jl
@@ -1,261 +1,261 @@
-using NDTensors: using_auto_fermion
-
-# `ValType::Type{<:Number}` is used instead of `ValType::Type` for efficiency, possibly due to increased method specialization.
-# See https://github.com/ITensor/ITensors.jl/pull/1183.
-function qn_svdMPO(
-  ValType::Type{<:Number}, os::OpSum{C}, sites; mindim=1, maxdim=typemax(Int), cutoff=1e-15
-)::MPO where {C}
-  N = length(sites)
-
-  # Specifying the element type with `Dict{QN,Matrix{ValType}}[...]` improves type inference and therefore efficiency.
-  # See https://github.com/ITensor/ITensors.jl/pull/1183.
-  Vs = Dict{QN,Matrix{ValType}}[Dict{QN,Matrix{ValType}}() for n in 1:(N + 1)]
-  sparse_MPO = [QNMatElem{Scaled{C,Prod{Op}}}[] for n in 1:N]
-
-  function crosses_bond(t::Scaled{C,Prod{Op}}, n::Int)
-    return (only(site(t[1])) <= n <= only(site(t[end])))
-  end
-
-  # A cache of the ITensor operators on a certain site
-  # of a certain type
-  op_cache = Dict{Pair{String,Int},ITensor}()
-  function calcQN(term::Vector{Op})
-    q = QN()
-    for st in term
-      op_tensor = get(op_cache, which_op(st) => only(site(st)), nothing)
-      if op_tensor === nothing
-        op_tensor = op(sites[only(site(st))], which_op(st); params(st)...)
-        op_cache[which_op(st) => only(site(st))] = op_tensor
-      end
-      q -= flux(op_tensor)
-    end
-    return q
-  end
-
-  Hflux = -calcQN(terms(first(terms(os))))
-
-  rightmap = Dict{Pair{Vector{Op},QN},Int}()
-  next_rightmap = Dict{Pair{Vector{Op},QN},Int}()
-
-  for n in 1:N
-    h_sparse = Dict{QN,Vector{MatElem{ValType}}}()
-
-    leftmap = Dict{Pair{Vector{Op},QN},Int}()
-    for term in os
-      crosses_bond(term, n) || continue
-
-      left = filter(t -> (only(site(t)) < n), terms(term))
-      onsite = filter(t -> (only(site(t)) == n), terms(term))
-      right = filter(t -> (only(site(t)) > n), terms(term))
-
-      lqn = calcQN(left)
-      sqn = calcQN(onsite)
-
-      bond_row = -1
-      bond_col = -1
-      if !isempty(left)
-        bond_row = posInLink!(leftmap, left => lqn)
-        bond_col = posInLink!(rightmap, vcat(onsite, right) => lqn)
-        bond_coef = convert(ValType, coefficient(term))
-        q_h_sparse = get!(h_sparse, lqn, MatElem{ValType}[])
-        push!(q_h_sparse, MatElem(bond_row, bond_col, bond_coef))
-      end
-
-      rqn = sqn + lqn
-      A_row = bond_col
-      A_col = posInLink!(next_rightmap, right => rqn)
-      site_coef = one(C)
-      if A_row == -1
-        site_coef = coefficient(term)
-      end
-      if isempty(onsite)
-        if !using_auto_fermion() && isfermionic(right, sites)
-          push!(onsite, Op("F", n))
-        else
-          push!(onsite, Op("Id", n))
-        end
-      end
-      el = QNMatElem(lqn, rqn, A_row, A_col, site_coef * Prod(onsite))
-      push!(sparse_MPO[n], el)
-    end
-    remove_dups!(sparse_MPO[n])
-
-    if n > 1 && !isempty(h_sparse)
-      for (q, mat) in h_sparse
-        h = toMatrix(mat)
-        U, S, V = svd(h)
-        P = S .^ 2
-        truncate!(P; maxdim, cutoff, mindim)
-        tdim = length(P)
-        Vs[n][q] = Matrix{ValType}(V[:, 1:tdim])
-      end
-    end
-
-    rightmap = next_rightmap
-    next_rightmap = Dict{Pair{Vector{Op},QN},Int}()
-  end
-
-  #
-  # Make MPO link indices
-  #
-  llinks = Vector{QNIndex}(undef, N + 1)
-  # Set dir=In for fermionic ordering, avoid arrow sign
-  # <fermions>:
-  linkdir = using_auto_fermion() ? ITensors.In : ITensors.Out
-  llinks[1] = Index([QN() => 1, Hflux => 1]; tags="Link,l=0", dir=linkdir)
-  for n in 1:N
-    qi = Vector{Pair{QN,Int}}()
-    push!(qi, QN() => 1)
-    for (q, Vq) in Vs[n + 1]
-      cols = size(Vq, 2)
-      if using_auto_fermion() # <fermions>
-        push!(qi, (-q) => cols)
-      else
-        push!(qi, q => cols)
-      end
-    end
-    push!(qi, Hflux => 1)
-    llinks[n + 1] = Index(qi...; tags="Link,l=$n", dir=linkdir)
-  end
-
-  H = MPO(N)
-
-  # Find location where block of Index i
-  # matches QN q, but *not* 1 or dim(i)
-  # which are special ending/starting states
-  function qnblock(i::Index, q::QN)
-    for b in 2:(nblocks(i) - 1)
-      flux(i, Block(b)) == q && return b
-    end
-    return error("Could not find block of QNIndex with matching QN")
-  end
-  qnblockdim(i::Index, q::QN) = blockdim(i, qnblock(i, q))
-
-  for n in 1:N
-    ll = llinks[n]
-    rl = llinks[n + 1]
-
-    begin_block = Dict{Tuple{QN,Vector{Op}},Matrix{ValType}}()
-    cont_block = Dict{Tuple{QN,Vector{Op}},Matrix{ValType}}()
-    end_block = Dict{Tuple{QN,Vector{Op}},Matrix{ValType}}()
-    onsite_block = Dict{Tuple{QN,Vector{Op}},Matrix{ValType}}()
-
-    for el in sparse_MPO[n]
-      t = el.val
-      (abs(coefficient(t)) > eps()) || continue
-      A_row = el.row
-      A_col = el.col
-      ct = convert(ValType, coefficient(t))
-
-      ldim = (A_row == -1) ? 1 : qnblockdim(ll, el.rowqn)
-      rdim = (A_col == -1) ? 1 : qnblockdim(rl, el.colqn)
-      zero_mat() = zeros(ValType, ldim, rdim)
-
-      if A_row == -1 && A_col == -1
-        # Onsite term
-        M = get!(onsite_block, (el.rowqn, terms(t)), zeros(ValType, 1, 1))
-        M[1, 1] += ct
-      elseif A_row == -1
-        # Operator beginning a term on site n
-        M = get!(begin_block, (el.rowqn, terms(t)), zero_mat())
-        VR = Vs[n + 1][el.colqn]
-        for c in 1:size(VR, 2)
-          M[1, c] += ct * VR[A_col, c]
-        end
-      elseif A_col == -1
-        # Operator ending a term on site n
-        M = get!(end_block, (el.rowqn, terms(t)), zero_mat())
-        VL = Vs[n][el.rowqn]
-        for r in 1:size(VL, 2)
-          M[r, 1] += ct * conj(VL[A_row, r])
-        end
-      else
-        # Operator continuing a term on site n
-        M = get!(cont_block, (el.rowqn, terms(t)), zero_mat())
-        VL = Vs[n][el.rowqn]
-        VR = Vs[n + 1][el.colqn]
-        for r in 1:size(VL, 2), c in 1:size(VR, 2)
-          M[r, c] += ct * conj(VL[A_row, r]) * VR[A_col, c]
-        end
-      end
-    end
-
-    H[n] = ITensor()
-
-    # Helper functions to compute block locations
-    # of various blocks within the onsite blocks,
-    # begin blocks, etc.
-    loc_onsite(rq, cq) = Block(nblocks(ll), 1)
-    loc_begin(rq, cq) = Block(nblocks(ll), qnblock(rl, cq))
-    loc_cont(rq, cq) = Block(qnblock(ll, rq), qnblock(rl, cq))
-    loc_end(rq, cq) = Block(qnblock(ll, rq), 1)
-
-    for (loc, block) in (
-      (loc_onsite, onsite_block),
-      (loc_begin, begin_block),
-      (loc_end, end_block),
-      (loc_cont, cont_block),
-    )
-      for (q_op, M) in block
-        op_prod = q_op[2]
-        Op = computeSiteProd(sites, Prod(op_prod))
-        (nnzblocks(Op) == 0) && continue
-
-        rq = q_op[1]
-        sq = flux(Op)
-        cq = rq
-        if !isnothing(sq)
-          # By convention, if `Op` has no blocks it has a flux
-          # of `nothing`, catch this case
-          cq -= sq
-
-          if using_auto_fermion()
-            # <fermions>:
-            # MPO is defined with Index order
-            # of (rl,s[n]',s[n],cl) where rl = row link, cl = col link
-            # so compute sign that would result by permuting cl from
-            # second position to last position:
-            if fparity(sq) == 1 && fparity(cq) == 1
-              Op .*= -1
-            end
-          end
-        end
-
-        b = loc(rq, cq)
-        T = ITensors.NDTensors.BlockSparseTensor(ValType, [b], (dag(ll), rl))
-        T[b] .= M
-
-        H[n] += (itensor(T) * Op)
-      end
-    end
-
-    # Put in ending identity operator
-    Id = op("Id", sites[n])
-    b = Block(1, 1)
-    T = ITensors.NDTensors.BlockSparseTensor(ValType, [b], (dag(ll), rl))
-    T[b] = 1
-    H[n] += (itensor(T) * Id)
-
-    # Put in starting identity operator
-    b = Block(nblocks(ll), nblocks(rl))
-    T = ITensors.NDTensors.BlockSparseTensor(ValType, [b], (dag(ll), rl))
-    T[b] = 1
-    H[n] += (itensor(T) * Id)
-  end # for n in 1:N
-
-  L = ITensor(llinks[1])
-  L[llinks[1] => end] = 1.0
-  H[1] *= L
-
-  R = ITensor(dag(llinks[N + 1]))
-  R[dag(llinks[N + 1]) => 1] = 1.0
-  H[N] *= R
-
-  return H
-end #qn_svdMPO
-
-function qn_svdMPO(os::OpSum{C}, sites; kwargs...)::MPO where {C}
-  # Function barrier to improve type stability
-  ValType = determineValType(terms(os))
-  return qn_svdMPO(ValType, os, sites; kwargs...)
-end
+## ## using NDTensors: using_auto_fermion
+## 
+## # `ValType::Type{<:Number}` is used instead of `ValType::Type` for efficiency, possibly due to increased method specialization.
+## # See https://github.com/ITensor/ITensors.jl/pull/1183.
+## function qn_svdMPO(
+##   ValType::Type{<:Number}, os::OpSum{C}, sites; mindim=1, maxdim=typemax(Int), cutoff=1e-15
+## )::MPO where {C}
+##   N = length(sites)
+## 
+##   # Specifying the element type with `Dict{QN,Matrix{ValType}}[...]` improves type inference and therefore efficiency.
+##   # See https://github.com/ITensor/ITensors.jl/pull/1183.
+##   Vs = Dict{QN,Matrix{ValType}}[Dict{QN,Matrix{ValType}}() for n in 1:(N + 1)]
+##   sparse_MPO = [QNMatElem{Scaled{C,Prod{Op}}}[] for n in 1:N]
+## 
+##   function crosses_bond(t::Scaled{C,Prod{Op}}, n::Int)
+##     return (only(site(t[1])) <= n <= only(site(t[end])))
+##   end
+## 
+##   # A cache of the ITensor operators on a certain site
+##   # of a certain type
+##   op_cache = Dict{Pair{String,Int},ITensor}()
+##   function calcQN(term::Vector{Op})
+##     q = QN()
+##     for st in term
+##       op_tensor = get(op_cache, which_op(st) => only(site(st)), nothing)
+##       if op_tensor === nothing
+##         op_tensor = op(sites[only(site(st))], which_op(st); params(st)...)
+##         op_cache[which_op(st) => only(site(st))] = op_tensor
+##       end
+##       q -= flux(op_tensor)
+##     end
+##     return q
+##   end
+## 
+##   Hflux = -calcQN(terms(first(terms(os))))
+## 
+##   rightmap = Dict{Pair{Vector{Op},QN},Int}()
+##   next_rightmap = Dict{Pair{Vector{Op},QN},Int}()
+## 
+##   for n in 1:N
+##     h_sparse = Dict{QN,Vector{MatElem{ValType}}}()
+## 
+##     leftmap = Dict{Pair{Vector{Op},QN},Int}()
+##     for term in os
+##       crosses_bond(term, n) || continue
+## 
+##       left = filter(t -> (only(site(t)) < n), terms(term))
+##       onsite = filter(t -> (only(site(t)) == n), terms(term))
+##       right = filter(t -> (only(site(t)) > n), terms(term))
+## 
+##       lqn = calcQN(left)
+##       sqn = calcQN(onsite)
+## 
+##       bond_row = -1
+##       bond_col = -1
+##       if !isempty(left)
+##         bond_row = posInLink!(leftmap, left => lqn)
+##         bond_col = posInLink!(rightmap, vcat(onsite, right) => lqn)
+##         bond_coef = convert(ValType, coefficient(term))
+##         q_h_sparse = get!(h_sparse, lqn, MatElem{ValType}[])
+##         push!(q_h_sparse, MatElem(bond_row, bond_col, bond_coef))
+##       end
+## 
+##       rqn = sqn + lqn
+##       A_row = bond_col
+##       A_col = posInLink!(next_rightmap, right => rqn)
+##       site_coef = one(C)
+##       if A_row == -1
+##         site_coef = coefficient(term)
+##       end
+##       if isempty(onsite)
+##         if !using_auto_fermion() && isfermionic(right, sites)
+##           push!(onsite, Op("F", n))
+##         else
+##           push!(onsite, Op("Id", n))
+##         end
+##       end
+##       el = QNMatElem(lqn, rqn, A_row, A_col, site_coef * Prod(onsite))
+##       push!(sparse_MPO[n], el)
+##     end
+##     remove_dups!(sparse_MPO[n])
+## 
+##     if n > 1 && !isempty(h_sparse)
+##       for (q, mat) in h_sparse
+##         h = toMatrix(mat)
+##         U, S, V = svd(h)
+##         P = S .^ 2
+##         truncate!(P; maxdim, cutoff, mindim)
+##         tdim = length(P)
+##         Vs[n][q] = Matrix{ValType}(V[:, 1:tdim])
+##       end
+##     end
+## 
+##     rightmap = next_rightmap
+##     next_rightmap = Dict{Pair{Vector{Op},QN},Int}()
+##   end
+## 
+##   #
+##   # Make MPO link indices
+##   #
+##   llinks = Vector{QNIndex}(undef, N + 1)
+##   # Set dir=In for fermionic ordering, avoid arrow sign
+##   # <fermions>:
+##   linkdir = using_auto_fermion() ? ITensors.In : ITensors.Out
+##   llinks[1] = Index([QN() => 1, Hflux => 1]; tags="Link,l=0", dir=linkdir)
+##   for n in 1:N
+##     qi = Vector{Pair{QN,Int}}()
+##     push!(qi, QN() => 1)
+##     for (q, Vq) in Vs[n + 1]
+##       cols = size(Vq, 2)
+##       if using_auto_fermion() # <fermions>
+##         push!(qi, (-q) => cols)
+##       else
+##         push!(qi, q => cols)
+##       end
+##     end
+##     push!(qi, Hflux => 1)
+##     llinks[n + 1] = Index(qi...; tags="Link,l=$n", dir=linkdir)
+##   end
+## 
+##   H = MPO(N)
+## 
+##   # Find location where block of Index i
+##   # matches QN q, but *not* 1 or dim(i)
+##   # which are special ending/starting states
+##   function qnblock(i::Index, q::QN)
+##     for b in 2:(nblocks(i) - 1)
+##       flux(i, Block(b)) == q && return b
+##     end
+##     return error("Could not find block of QNIndex with matching QN")
+##   end
+##   qnblockdim(i::Index, q::QN) = blockdim(i, qnblock(i, q))
+## 
+##   for n in 1:N
+##     ll = llinks[n]
+##     rl = llinks[n + 1]
+## 
+##     begin_block = Dict{Tuple{QN,Vector{Op}},Matrix{ValType}}()
+##     cont_block = Dict{Tuple{QN,Vector{Op}},Matrix{ValType}}()
+##     end_block = Dict{Tuple{QN,Vector{Op}},Matrix{ValType}}()
+##     onsite_block = Dict{Tuple{QN,Vector{Op}},Matrix{ValType}}()
+## 
+##     for el in sparse_MPO[n]
+##       t = el.val
+##       (abs(coefficient(t)) > eps()) || continue
+##       A_row = el.row
+##       A_col = el.col
+##       ct = convert(ValType, coefficient(t))
+## 
+##       ldim = (A_row == -1) ? 1 : qnblockdim(ll, el.rowqn)
+##       rdim = (A_col == -1) ? 1 : qnblockdim(rl, el.colqn)
+##       zero_mat() = zeros(ValType, ldim, rdim)
+## 
+##       if A_row == -1 && A_col == -1
+##         # Onsite term
+##         M = get!(onsite_block, (el.rowqn, terms(t)), zeros(ValType, 1, 1))
+##         M[1, 1] += ct
+##       elseif A_row == -1
+##         # Operator beginning a term on site n
+##         M = get!(begin_block, (el.rowqn, terms(t)), zero_mat())
+##         VR = Vs[n + 1][el.colqn]
+##         for c in 1:size(VR, 2)
+##           M[1, c] += ct * VR[A_col, c]
+##         end
+##       elseif A_col == -1
+##         # Operator ending a term on site n
+##         M = get!(end_block, (el.rowqn, terms(t)), zero_mat())
+##         VL = Vs[n][el.rowqn]
+##         for r in 1:size(VL, 2)
+##           M[r, 1] += ct * conj(VL[A_row, r])
+##         end
+##       else
+##         # Operator continuing a term on site n
+##         M = get!(cont_block, (el.rowqn, terms(t)), zero_mat())
+##         VL = Vs[n][el.rowqn]
+##         VR = Vs[n + 1][el.colqn]
+##         for r in 1:size(VL, 2), c in 1:size(VR, 2)
+##           M[r, c] += ct * conj(VL[A_row, r]) * VR[A_col, c]
+##         end
+##       end
+##     end
+## 
+##     H[n] = ITensor()
+## 
+##     # Helper functions to compute block locations
+##     # of various blocks within the onsite blocks,
+##     # begin blocks, etc.
+##     loc_onsite(rq, cq) = Block(nblocks(ll), 1)
+##     loc_begin(rq, cq) = Block(nblocks(ll), qnblock(rl, cq))
+##     loc_cont(rq, cq) = Block(qnblock(ll, rq), qnblock(rl, cq))
+##     loc_end(rq, cq) = Block(qnblock(ll, rq), 1)
+## 
+##     for (loc, block) in (
+##       (loc_onsite, onsite_block),
+##       (loc_begin, begin_block),
+##       (loc_end, end_block),
+##       (loc_cont, cont_block),
+##     )
+##       for (q_op, M) in block
+##         op_prod = q_op[2]
+##         Op = computeSiteProd(sites, Prod(op_prod))
+##         (nnzblocks(Op) == 0) && continue
+## 
+##         rq = q_op[1]
+##         sq = flux(Op)
+##         cq = rq
+##         if !isnothing(sq)
+##           # By convention, if `Op` has no blocks it has a flux
+##           # of `nothing`, catch this case
+##           cq -= sq
+## 
+##           if using_auto_fermion()
+##             # <fermions>:
+##             # MPO is defined with Index order
+##             # of (rl,s[n]',s[n],cl) where rl = row link, cl = col link
+##             # so compute sign that would result by permuting cl from
+##             # second position to last position:
+##             if fparity(sq) == 1 && fparity(cq) == 1
+##               Op .*= -1
+##             end
+##           end
+##         end
+## 
+##         b = loc(rq, cq)
+##         T = ITensors.NDTensors.BlockSparseTensor(ValType, [b], (dag(ll), rl))
+##         T[b] .= M
+## 
+##         H[n] += (ITensor(T) * Op)
+##       end
+##     end
+## 
+##     # Put in ending identity operator
+##     Id = op("Id", sites[n])
+##     b = Block(1, 1)
+##     T = ITensors.NDTensors.BlockSparseTensor(ValType, [b], (dag(ll), rl))
+##     T[b] = 1
+##     H[n] += (ITensor(T) * Id)
+## 
+##     # Put in starting identity operator
+##     b = Block(nblocks(ll), nblocks(rl))
+##     T = ITensors.NDTensors.BlockSparseTensor(ValType, [b], (dag(ll), rl))
+##     T[b] = 1
+##     H[n] += (ITensor(T) * Id)
+##   end # for n in 1:N
+## 
+##   L = ITensor(llinks[1])
+##   L[llinks[1] => end] = 1.0
+##   H[1] *= L
+## 
+##   R = ITensor(dag(llinks[N + 1]))
+##   R[dag(llinks[N + 1]) => 1] = 1.0
+##   H[N] *= R
+## 
+##   return H
+## end #qn_svdMPO
+## 
+## function qn_svdMPO(os::OpSum{C}, sites; kwargs...)::MPO where {C}
+##   # Function barrier to improve type stability
+##   ValType = determineValType(terms(os))
+##   return qn_svdMPO(ValType, os, sites; kwargs...)
+## end
diff --git a/src/opsum_to_mpo/qnmatelem.jl b/src/opsum_to_mpo/qnmatelem.jl
index 7ec55c4..9245187 100644
--- a/src/opsum_to_mpo/qnmatelem.jl
+++ b/src/opsum_to_mpo/qnmatelem.jl
@@ -1,30 +1,30 @@
-struct QNMatElem{T}
-  rowqn::QN
-  colqn::QN
-  row::Int
-  col::Int
-  val::T
-end
-
-function Base.:(==)(m1::QNMatElem{T}, m2::QNMatElem{T})::Bool where {T}
-  return (
-    m1.row == m2.row &&
-    m1.col == m2.col &&
-    m1.val == m2.val &&
-    m1.rowqn == m2.rowqn &&
-    m1.colqn == m2.colqn
-  )
-end
-
-function Base.isless(m1::QNMatElem{T}, m2::QNMatElem{T})::Bool where {T}
-  if m1.rowqn != m2.rowqn
-    return m1.rowqn < m2.rowqn
-  elseif m1.colqn != m2.colqn
-    return m1.colqn < m2.colqn
-  elseif m1.row != m2.row
-    return m1.row < m2.row
-  elseif m1.col != m2.col
-    return m1.col < m2.col
-  end
-  return m1.val < m2.val
-end
+## struct QNMatElem{T}
+##   rowqn::QN
+##   colqn::QN
+##   row::Int
+##   col::Int
+##   val::T
+## end
+## 
+## function Base.:(==)(m1::QNMatElem{T}, m2::QNMatElem{T})::Bool where {T}
+##   return (
+##     m1.row == m2.row &&
+##     m1.col == m2.col &&
+##     m1.val == m2.val &&
+##     m1.rowqn == m2.rowqn &&
+##     m1.colqn == m2.colqn
+##   )
+## end
+## 
+## function Base.isless(m1::QNMatElem{T}, m2::QNMatElem{T})::Bool where {T}
+##   if m1.rowqn != m2.rowqn
+##     return m1.rowqn < m2.rowqn
+##   elseif m1.colqn != m2.colqn
+##     return m1.colqn < m2.colqn
+##   elseif m1.row != m2.row
+##     return m1.row < m2.row
+##   elseif m1.col != m2.col
+##     return m1.col < m2.col
+##   end
+##   return m1.val < m2.val
+## end
diff --git a/src/siteinds.jl b/src/siteinds.jl
new file mode 100644
index 0000000..819946b
--- /dev/null
+++ b/src/siteinds.jl
@@ -0,0 +1,9 @@
+using ITensors: ITensors, Index, settag
+using QuantumOperatorDefinitions: QuantumOperatorDefinitions
+
+# TODO: Move to `ITensorSites.jl`?
+function siteinds(t::String, n::Int; kwargs...)
+  # TODO: QuantumOperatorAlgebra.jl defines its own  `QuantumOperatorAlgebra.sites`
+  # with a different meaning, decide if we should keep both.
+  return QuantumOperatorDefinitions.sites(Index, t, n; kwargs...)
+end
diff --git a/src/solvers/ITensorsExtensions.jl b/src/solvers/ITensorsExtensions.jl
index 90ed8b1..2782697 100644
--- a/src/solvers/ITensorsExtensions.jl
+++ b/src/solvers/ITensorsExtensions.jl
@@ -1,9 +1,13 @@
 module ITensorsExtensions
-using ITensors: ITensor, array, inds, itensor
+
+using ITensors: ITensor, inds
+using NamedDimsArrays: unname
+
 function to_vec(x::ITensor)
   function to_itensor(x_vec)
-    return itensor(x_vec, inds(x))
+    return ITensor(x_vec, inds(x))
   end
-  return vec(array(x)), to_itensor
+  return vec(unname(x)), to_itensor
 end
+
 end
diff --git a/src/solvers/alternating_update.jl b/src/solvers/alternating_update.jl
index 6206123..79c607f 100644
--- a/src/solvers/alternating_update.jl
+++ b/src/solvers/alternating_update.jl
@@ -1,4 +1,5 @@
-using ITensors: ITensors, permute
+using ITensors: ITensors
+using VectorInterface: scalartype
 
 function _extend_sweeps_param(param, nsweeps)
   if param isa Number
@@ -39,7 +40,7 @@ function alternating_update(
   normalize=default_normalize(),
   maxdim=default_maxdim(),
   mindim=default_mindim(),
-  cutoff=default_cutoff(ITensors.scalartype(init)),
+  cutoff=default_cutoff(scalartype(init)),
   noise=default_noise(),
 )
   reduced_operator = ITensorMPS.reduced_operator(operator)
diff --git a/src/solvers/contract.jl b/src/solvers/contract.jl
index 1919fdd..664f307 100644
--- a/src/solvers/contract.jl
+++ b/src/solvers/contract.jl
@@ -1,4 +1,6 @@
-using ITensors: ITensors, Index, ITensor, @Algorithm_str, commoninds, contract, hasind, sim
+using BackendSelection: @Algorithm_str
+using ITensors: ITensors, ITensor, Index
+using TensorAlgebra: contract
 
 function contract_operator_state_updater(operator, init; internal_kwargs)
   # TODO: Use `contract(operator)`.
diff --git a/src/solvers/dmrg.jl b/src/solvers/dmrg.jl
index 90fc9df..d6be0ff 100644
--- a/src/solvers/dmrg.jl
+++ b/src/solvers/dmrg.jl
@@ -1,5 +1,6 @@
 using ITensors: ITensors
 using KrylovKit: eigsolve
+using VectorInterface: scalartype
 
 function eigsolve_updater(
   operator,
@@ -7,7 +8,7 @@ function eigsolve_updater(
   internal_kwargs,
   which_eigval=:SR,
   ishermitian=true,
-  tol=10^2 * eps(real(ITensors.scalartype(init))),
+  tol=10^2 * eps(real(scalartype(init))),
   krylovdim=3,
   maxiter=1,
   verbosity=0,
diff --git a/src/solvers/dmrg_x.jl b/src/solvers/dmrg_x.jl
index 48201e6..8d9972a 100644
--- a/src/solvers/dmrg_x.jl
+++ b/src/solvers/dmrg_x.jl
@@ -1,13 +1,17 @@
-using ITensors: array, contract, dag, uniqueind, onehot
+using GradedArrays: dag
+using ITensors: uniqueind
 using LinearAlgebra: eigen
+using NamedDimsArrays: unname
+using SparseArraysBase: oneelement
+using TensorAlgebra: contract
 
 function eigen_updater(operator, state; internal_kwargs)
   contracted_operator = contract(operator, ITensor(true))
   d, u = eigen(contracted_operator; ishermitian=true)
   u_ind = uniqueind(u, contracted_operator)
   u′_ind = uniqueind(d, u)
-  max_overlap, max_index = findmax(abs, array(state * dag(u)))
-  u_max = u * dag(onehot(eltype(u), u_ind => max_index))
+  max_overlap, max_index = findmax(abs, unname(state * dag(u)))
+  u_max = u * dag(oneelement(eltype(u), u_ind => max_index))
   d_max = d[u′_ind => max_index, u_ind => max_index]
   return u_max, (; eigval=d_max)
 end
diff --git a/src/solvers/expand.jl b/src/solvers/expand.jl
index c970e45..85d5bc6 100644
--- a/src/solvers/expand.jl
+++ b/src/solvers/expand.jl
@@ -1,21 +1,12 @@
 using Adapt: adapt
+using BackendSelection: @Algorithm_str, Algorithm
+using DiagonalArrays: δ
+using GradedArrays: dag
 using ITensors:
-  ITensors,
-  Algorithm,
-  Index,
-  ITensor,
-  @Algorithm_str,
-  δ,
-  commonind,
-  dag,
-  denseblocks,
-  directsum,
-  hasqns,
-  prime,
-  scalartype,
-  uniqueinds
+  ITensors, Index, ITensor, commonind, denseblocks, directsum, hasqns, prime, uniqueinds
 using LinearAlgebra: normalize, svd, tr
-using NDTensors: unwrap_array_type
+using TypeParameterAccessors: unwrap_array_type
+using VectorInterface: scalartype
 
 # Possible improvements:
 #  - Allow a maxdim argument to be passed to `expand`.
diff --git a/src/solvers/reducedconstantterm.jl b/src/solvers/reducedconstantterm.jl
index 52c882b..9b8916b 100644
--- a/src/solvers/reducedconstantterm.jl
+++ b/src/solvers/reducedconstantterm.jl
@@ -1,4 +1,5 @@
-using ITensors: ITensors, ITensor, dag, dim, prime
+using GradedArrays: dag
+using ITensors: ITensors, ITensor, dim, prime
 
 """
 Holds the following data where basis
@@ -105,7 +106,7 @@ function ITensorMPS.makeR!(reduced_state::ReducedConstantTerm, basis::MPS, posit
 end
 
 function ITensors.contract(reduced_state::ReducedConstantTerm, v::ITensor)
-  reduced_state_tensors = Union{ITensor,OneITensor}[lproj(reduced_state)]
+  reduced_state_tensors = [lproj(reduced_state)]
   append!(
     reduced_state_tensors,
     [prime(t, "Link") for t in reduced_state.state[site_range(reduced_state)]],
@@ -129,7 +130,7 @@ end
 
 # Contract the reduced constant term down to a since ITensor.
 function ITensors.contract(reduced_state::ReducedConstantTerm)
-  reduced_state_tensors = Union{ITensor,OneITensor}[lproj(reduced_state)]
+  reduced_state_tensors = [lproj(reduced_state)]
   append!(
     reduced_state_tensors,
     [dag(prime(t, "Link")) for t in reduced_state.state[site_range(reduced_state)]],
diff --git a/src/solvers/sweep_update.jl b/src/solvers/sweep_update.jl
index ac6cc24..f4d33b4 100644
--- a/src/solvers/sweep_update.jl
+++ b/src/solvers/sweep_update.jl
@@ -1,6 +1,7 @@
 using ITensors: ITensors, uniqueinds
 using LinearAlgebra: norm, normalize!, svd
 using Printf: @printf
+using VectorInterface: scalartype
 
 function sweep_update(
   order::TDVPOrder,
@@ -76,7 +77,7 @@ function sub_sweep_update(
   outputlevel=default_outputlevel(),
   maxdim=default_maxdim(),
   mindim=default_mindim(),
-  cutoff=default_cutoff(ITensors.scalartype(state)),
+  cutoff=default_cutoff(scalartype(state)),
   noise=default_noise(),
 )
   reduced_operator = copy(reduced_operator)
diff --git a/src/solvers/tdvp.jl b/src/solvers/tdvp.jl
index 02d8c5c..a02e9b1 100644
--- a/src/solvers/tdvp.jl
+++ b/src/solvers/tdvp.jl
@@ -1,4 +1,4 @@
-using ITensors: Algorithm, @Algorithm_str
+using BackendSelection: @Algorithm_str, Algorithm
 using KrylovKit: exponentiate
 
 function exponentiate_updater(operator, init; internal_kwargs, kwargs...)
diff --git a/src/solvers/timedependentsum.jl b/src/solvers/timedependentsum.jl
index f2998de..ed8dd12 100644
--- a/src/solvers/timedependentsum.jl
+++ b/src/solvers/timedependentsum.jl
@@ -1,4 +1,6 @@
-using ITensors: ITensor, inds, permute
+using ITensors: ITensor, inds
+using NamedDimsArrays: aligndims
+using QuantumOperatorAlgebra: QuantumOperatorAlgebra
 
 # Represents a time-dependent sum of terms:
 #
@@ -10,7 +12,7 @@ struct TimeDependentSum{Coefficients,Terms}
 end
 
 coefficients(expr::TimeDependentSum) = expr.coefficients
-terms(expr::TimeDependentSum) = expr.terms
+QuantumOperatorAlgebra.terms(expr::TimeDependentSum) = expr.terms
 function Base.copy(expr::TimeDependentSum)
   return TimeDependentSum(coefficients(expr), copy.(terms(expr)))
 end
@@ -51,7 +53,7 @@ struct ScaledSum{Coefficients,Terms}
 end
 
 coefficients(expr::ScaledSum) = expr.coefficients
-terms(expr::ScaledSum) = expr.terms
+QuantumOperatorAlgebra.terms(expr::ScaledSum) = expr.terms
 
 # Apply the scaled sum of terms:
 #
@@ -65,4 +67,4 @@ function scaledsum_apply(expr, x)
   end
 end
 (expr::ScaledSum)(x) = scaledsum_apply(expr, x)
-(expr::ScaledSum)(x::ITensor) = permute(scaledsum_apply(expr, x), inds(x))
+(expr::ScaledSum)(x::ITensor) = aligndims(scaledsum_apply(expr, x), inds(x))
diff --git a/test/Project.toml b/test/Project.toml
index d6fd7a4..c4dc1b5 100644
--- a/test/Project.toml
+++ b/test/Project.toml
@@ -9,7 +9,6 @@ ITensors = "9136182c-28ba-11e9-034c-db9fb085ebd5"
 JLD2 = "033835bb-8acc-5ee8-8aae-3f567f8a3819"
 KrylovKit = "0b1a1467-8014-51b9-945f-bf0ae24f4b77"
 LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
-NDTensors = "23ae76d9-e61a-49c4-8f12-3f1a16adf9cf"
 Observers = "338f10d5-c7f1-4033-a7d1-f9dec39bcaa0"
 OptimKit = "77e91f04-9b3b-57a6-a776-40b61faaebe0"
 OrdinaryDiffEqTsit5 = "b1df2697-797e-41e3-8120-5422d3b24e4a"
diff --git a/test/base/Project.toml b/test/base/Project.toml
index e84dc5a..e42d08c 100644
--- a/test/base/Project.toml
+++ b/test/base/Project.toml
@@ -1,9 +1,8 @@
 [deps]
 Combinatorics = "861a8166-3701-5b0c-9a16-15d98fcdc6aa"
 HDF5 = "f67ccb44-e63f-5c2f-98bd-6dc0ccc4ba2f"
+ITensors = "4795dd04-0d67-49bb-8f44-b89c448a1dc7"
 ITensorMPS = "0d1a4710-d33b-49a5-8f18-73bdf49b47e2"
-ITensors = "9136182c-28ba-11e9-034c-db9fb085ebd5"
 JLD2 = "033835bb-8acc-5ee8-8aae-3f567f8a3819"
-NDTensors = "23ae76d9-e61a-49c4-8f12-3f1a16adf9cf"
 StableRNGs = "860ef19b-820b-49d6-a774-d7a799459cd3"
 Suppressor = "fd094767-a336-5f1f-9728-57cf17d0bbfb"
diff --git a/test/base/test_autompo.jl b/test/base/test_autompo.jl
index c4ebeec..cacdc39 100644
--- a/test/base/test_autompo.jl
+++ b/test/base/test_autompo.jl
@@ -1,6 +1,7 @@
 @eval module $(gensym())
 using ITensorMPS, ITensors, Test, Random, JLD2
 using NDTensors: scalartype
+using SparseArraysBase: oneelement
 
 include(joinpath(@__DIR__, "utils", "util.jl"))
 
@@ -98,26 +99,26 @@ function NNheisenbergMPO(sites, J1::Float64, J2::Float64)::MPO
     ll = dag(link[n])
     rl = link[n + 1]
     H[n] = ITensor(ll, dag(s), s', rl)
-    H[n] += onehot(ll => 1) * onehot(rl => 1) * op(sites, "Id", n)
-    H[n] += onehot(ll => 8) * onehot(rl => 8) * op(sites, "Id", n)
-
-    H[n] += onehot(ll => 2) * onehot(rl => 1) * op(sites, "S-", n)
-    H[n] += onehot(ll => 5) * onehot(rl => 2) * op(sites, "Id", n)
-    H[n] += onehot(ll => 8) * onehot(rl => 2) * op(sites, "S+", n) * J1 / 2
-    H[n] += onehot(ll => 8) * onehot(rl => 5) * op(sites, "S+", n) * J2 / 2
-
-    H[n] += onehot(ll => 3) * onehot(rl => 1) * op(sites, "S+", n)
-    H[n] += onehot(ll => 6) * onehot(rl => 3) * op(sites, "Id", n)
-    H[n] += onehot(ll => 8) * onehot(rl => 3) * op(sites, "S-", n) * J1 / 2
-    H[n] += onehot(ll => 8) * onehot(rl => 6) * op(sites, "S-", n) * J2 / 2
-
-    H[n] += onehot(ll => 4) * onehot(rl => 1) * op(sites, "Sz", n)
-    H[n] += onehot(ll => 7) * onehot(rl => 4) * op(sites, "Id", n)
-    H[n] += onehot(ll => 8) * onehot(rl => 4) * op(sites, "Sz", n) * J1
-    H[n] += onehot(ll => 8) * onehot(rl => 7) * op(sites, "Sz", n) * J2
+    H[n] += oneelement(ll => 1) * oneelement(rl => 1) * op(sites, "Id", n)
+    H[n] += oneelement(ll => 8) * oneelement(rl => 8) * op(sites, "Id", n)
+
+    H[n] += oneelement(ll => 2) * oneelement(rl => 1) * op(sites, "S-", n)
+    H[n] += oneelement(ll => 5) * oneelement(rl => 2) * op(sites, "Id", n)
+    H[n] += oneelement(ll => 8) * oneelement(rl => 2) * op(sites, "S+", n) * J1 / 2
+    H[n] += oneelement(ll => 8) * oneelement(rl => 5) * op(sites, "S+", n) * J2 / 2
+
+    H[n] += oneelement(ll => 3) * oneelement(rl => 1) * op(sites, "S+", n)
+    H[n] += oneelement(ll => 6) * oneelement(rl => 3) * op(sites, "Id", n)
+    H[n] += oneelement(ll => 8) * oneelement(rl => 3) * op(sites, "S-", n) * J1 / 2
+    H[n] += oneelement(ll => 8) * oneelement(rl => 6) * op(sites, "S-", n) * J2 / 2
+
+    H[n] += oneelement(ll => 4) * oneelement(rl => 1) * op(sites, "Sz", n)
+    H[n] += oneelement(ll => 7) * oneelement(rl => 4) * op(sites, "Id", n)
+    H[n] += oneelement(ll => 8) * oneelement(rl => 4) * op(sites, "Sz", n) * J1
+    H[n] += oneelement(ll => 8) * oneelement(rl => 7) * op(sites, "Sz", n) * J2
   end
-  H[1] *= onehot(link[1] => 8)
-  H[N] *= onehot(dag(link[N + 1]) => 1)
+  H[1] *= oneelement(link[1] => 8)
+  H[N] *= oneelement(dag(link[N + 1]) => 1)
   return H
 end
 
@@ -1062,8 +1063,8 @@ end
       q = flux(op(which_op, sites[j]))
       links = [Index([n < j ? q => 1 : QN() => 1], "Link,l=$n") for n in 1:N]
       for n in 1:(N - 1)
-        M[n] *= onehot(links[n] => 1)
-        M[n + 1] *= onehot(dag(links[n]) => 1)
+        M[n] *= oneelement(links[n] => 1)
+        M[n + 1] *= oneelement(dag(links[n]) => 1)
       end
       return M
     end
diff --git a/test/base/test_mps.jl b/test/base/test_mps.jl
index 47b45c2..1875980 100644
--- a/test/base/test_mps.jl
+++ b/test/base/test_mps.jl
@@ -4,6 +4,7 @@ using ITensorMPS
 using ITensors
 using LinearAlgebra: diag
 using Random
+using SparseArraysBase: oneelement
 using Test
 
 Random.seed!(1234)
@@ -1452,8 +1453,8 @@ end
     CSWAP = [op("CSWAP", s, n, m, k) for n in 1:N, m in 1:N, k in 1:N]
     CCCNOT = [op("CCCNOT", s, n, m, k, l) for n in 1:N, m in 1:N, k in 1:N, l in 1:N]
 
-    v0 = [onehot(s[n] => "0") for n in 1:N]
-    v1 = [onehot(s[n] => "1") for n in 1:N]
+    v0 = [oneelement(s[n] => "0") for n in 1:N]
+    v1 = [oneelement(s[n] => "1") for n in 1:N]
 
     # Single qubit
     @test product(I[1], v0[1]) ≈ v0[1]