diff --git a/benchmark/bench_dmrg.jl b/benchmark/bench_dmrg.jl
index 3732634b82..7b94bbc32a 100644
--- a/benchmark/bench_dmrg.jl
+++ b/benchmark/bench_dmrg.jl
@@ -45,6 +45,12 @@ let N = 100
   dmrg(H, psi, sweeps; outputlevel=0)
 
   suite["1d_S=1_heisenberg_qn"] = @benchmarkable dmrg($H, $psi, $sweeps; outputlevel=0)
+
+  # Precompile
+  H_splitblocks = splitblocks(linkinds, H)
+  dmrg(H_splitblocks, psi, sweeps; outputlevel=0)
+
+  suite["1d_S=1_heisenberg_qn_splitblocks"] = @benchmarkable dmrg($H_splitblocks, $psi, $sweeps; outputlevel=0)
 end
 
 #let Ny = 6,
diff --git a/docs/src/MPSandMPO.md b/docs/src/MPSandMPO.md
index a5652ee0d6..cad2e747f9 100644
--- a/docs/src/MPSandMPO.md
+++ b/docs/src/MPSandMPO.md
@@ -15,10 +15,10 @@ MPS(::Type{<:Number}, ::Vector{<:Index})
 randomMPS(sites::Vector{<:Index})
 randomMPS(::Type{<:Number}, sites::Vector{<:Index})
 randomMPS(::Vector{<:Index}, ::Any)
-productMPS(::Vector{<:Index}, ::Any)
-productMPS(::Type{<:Number}, ::Vector{<:Index}, ::Any)
-productMPS(::Vector{<:IndexVal})
-productMPS(::Type{<:Number}, ::Vector{<:IndexVal})
+MPS(::Vector{<:Index}, ::Any)
+MPS(::Type{<:Number}, ::Vector{<:Index}, ::Any)
+MPS(::Vector{<:IndexVal})
+MPS(::Type{<:Number}, ::Vector{<:IndexVal})
 ```
 
 ## MPO Constructors
@@ -40,6 +40,7 @@ deepcopy(::ITensors.AbstractMPS)
 
 ```@docs
 eltype(::ITensors.AbstractMPS)
+ITensors.promote_itensor_eltype(::ITensors.AbstractMPS)
 flux(::ITensors.AbstractMPS)
 hasqns(::ITensors.AbstractMPS)
 length(::ITensors.AbstractMPS)
diff --git a/docs/src/examples/DMRG.md b/docs/src/examples/DMRG.md
index 90840f159d..ae5476e9d1 100644
--- a/docs/src/examples/DMRG.md
+++ b/docs/src/examples/DMRG.md
@@ -52,7 +52,7 @@ psi0 = randomMPS(sites,2)
 Here we have made a random MPS of bond dimension 2. We could have used a random product
 state instead, but choosing a slightly larger bond dimension can help DMRG avoid getting
 stuck in local minima. We could also set psi to some specific initial state using the 
-function `productMPS`, which is actually required if we were conserving QNs.
+function `MPS`, which is actually required if we were conserving QNs.
 
 Finally, we are ready to call DMRG:
 
diff --git a/docs/src/getting_started/Tutorials.md b/docs/src/getting_started/Tutorials.md
index db7efbe784..5f8b6aa015 100644
--- a/docs/src/getting_started/Tutorials.md
+++ b/docs/src/getting_started/Tutorials.md
@@ -216,7 +216,7 @@ by the lines
 
 ```julia
 state = [isodd(n) ? "Up" : "Dn" for n=1:N]
-psi0 = productMPS(sites,state)
+psi0 = MPS(sites,state)
 ```
 
 The first line of the new code above makes an array of strings which 
@@ -251,13 +251,13 @@ computing the quantum-number flux of `psi0`
     use the initial state:
     ```julia
     state = ["Up" for n=1:N]
-    psi0 = productMPS(sites,state)
+    psi0 = MPS(sites,state)
     ```
     Or to initialize this 10-site system to have a total "Sz" of +16
     in ITensor units (``S^z=8`` in physical units):
     ```julia
     state = ["Dn","Up","Up","Up","Up","Up","Up","Up","Up","Up"]
-    psi0 = productMPS(sites,state)
+    psi0 = MPS(sites,state)
     ```
     would work (as would any `state` with one "Dn" and nine "Up"'s
     in any order).
@@ -265,7 +265,7 @@ computing the quantum-number flux of `psi0`
     in ITensor units (``S^z=9`` in physical units) as
     ```julia
     state = ["Z0","Up","Up","Up","Up","Up","Up","Up","Up","Up"]
-    psi0 = productMPS(sites,state)
+    psi0 = MPS(sites,state)
     ```
     where "Z0" refers to the ``S^z=0`` state of a spin-one spin.
 
@@ -294,7 +294,7 @@ let
   H = MPO(ampo,sites)
 
   state = [isodd(n) ? "Up" : "Dn" for n=1:N]
-  psi0 = productMPS(sites,state)
+  psi0 = MPS(sites,state)
   @show flux(psi0)
 
   sweeps = Sweeps(5)
@@ -413,7 +413,7 @@ let
   end
 
   # Initialize psi to be a product state (alternating up and down)
-  psi = productMPS(s, n -> isodd(n) ? "Up" : "Dn")
+  psi = MPS(s, n -> isodd(n) ? "Up" : "Dn")
 
   c = div(N,2)
 
@@ -490,7 +490,7 @@ we will take to be near the center of the MPS. The details of this
 function are outside the scope of this tutorial, but are explained in the
 example code for measuring MPS.
 
-The line of code `psi = productMPS(s, n -> isodd(n) ? "Up" : "Dn")`
+The line of code `psi = MPS(s, n -> isodd(n) ? "Up" : "Dn")`
 initializes our MPS `psi` as a product state of alternating 
 up and down spins. We call `measure_Sz` before starting the
 time evolution.
diff --git a/src/ITensors.jl b/src/ITensors.jl
index 2bef5b5c52..199c239444 100644
--- a/src/ITensors.jl
+++ b/src/ITensors.jl
@@ -7,6 +7,10 @@ The ITensor library also includes composable and extensible algorithms for optim
 """
 module ITensors
 
+#if isdefined(Base, :Experimental) && isdefined(Base.Experimental, Symbol("@compiler_options"))
+#  @eval Base.Experimental.@compiler_options compile=min optimize=1 infer=false
+#end
+
 #####################################
 # External packages
 #
diff --git a/src/indexset.jl b/src/indexset.jl
index 23072dc4d7..6a721b4b99 100644
--- a/src/indexset.jl
+++ b/src/indexset.jl
@@ -456,28 +456,45 @@ function swaptags(is::Indices, tags1, tags2, args...; kwargs...)
   return swaptags(fmatch(args...; kwargs...), is, tags1, tags2)
 end
 
-function replaceinds(is::Indices, rep_inds::Pair{<:Index,<:Index}...)
-  return replaceinds(is, zip(rep_inds...)...)
-end
+# Check that the QNs are all the same
+hassameflux(i1::Index, i2::Index) = (dim(i1) == dim(i2))
 
-function replaceinds(is::Indices, rep_inds::Vector{<:Pair{<:Index,<:Index}})
-  return replaceinds(is, rep_inds...)
+# For notation:
+# `replaceinds(is, [i, j] => [k, l])`
+# `replaceinds(is, (i, j) => (k, l))`
+function replaceinds(is::Indices, rep_inds::Pair)
+  return replaceinds(is, first(rep_inds), last(rep_inds))
 end
 
-function replaceinds(is::Indices, rep_inds::Tuple{Vararg{Pair{<:Index,<:Index}}})
-  return replaceinds(is, rep_inds...)
+# For notation:
+# `replaceinds(is, i => k, j => l)`
+function replaceinds(is::Indices, rep_inds::Pair{<:Index,<:Index}...)
+  return replaceinds(is, first.(rep_inds), last.(rep_inds))
 end
 
-function replaceinds(is::Indices, rep_inds::Pair)
-  return replaceinds(is, Tuple(first(rep_inds)) .=> Tuple(last(rep_inds)))
+# For notation:
+# `replaceinds(is, i, k)`
+function replaceinds(is::Indices, inds1::Index, inds2::Index)
+  return replaceinds(is, (inds1,), (inds2,))
 end
 
-# Check that the QNs are all the same
-hassameflux(i1::Index, i2::Index) = (dim(i1) == dim(i2))
+# For notation:
+# `replaceinds(is, [i => k, j => l])`
+# `replaceinds(is, (i => k, j => l))`
+function replaceinds(is::Indices, inds1_inds2)
+  @assert all(x -> x isa Pair, inds1_inds2)
+  return replaceinds(is, first.(inds1_inds2), last.(inds1_inds2))
+end
 
+# For notation:
+# `replaceinds(is, [i, j], [k, l])`
+# `replaceinds(is, (i, j), (k, l))`
 function replaceinds(is::Indices, inds1, inds2)
+  @assert length(inds1) == length(inds2)
+  isempty(inds1) && return is
   is1 = (inds1)
   poss = indexin(is1, is)
+  # TODO: don't convert to Tuple in the case of a Vector
   is_tuple = Tuple(is)
   for (j, pos) in enumerate(poss)
     isnothing(pos) && continue
diff --git a/src/mps/abstractmps.jl b/src/mps/abstractmps.jl
index 0eac94f5c7..0565a0c793 100644
--- a/src/mps/abstractmps.jl
+++ b/src/mps/abstractmps.jl
@@ -42,7 +42,7 @@ end
 The element type of the MPS/MPO. Always returns `ITensor`.
 
 For the element type of the ITensors of the MPS/MPO,
-use `common_itensor_eltype`.
+use [`ITensors.promote_itensor_eltype`](@ref).
 """
 eltype(::AbstractMPS) = ITensor
 
@@ -813,7 +813,7 @@ function hassamenuminds(::typeof(siteinds), M1::AbstractMPS, M2::AbstractMPS)
 end
 
 for fname in
-    (:sim, :prime, :setprime, :noprime, :addtags, :removetags, :replacetags, :settags)
+    (:sim, :prime, :setprime, :noprime, :replaceprime, :swapprime, :addtags, :removetags, :replacetags, :settags)
   fname! = Symbol(fname, :!)
   @eval begin
     """
@@ -1437,8 +1437,8 @@ end
 # TODO: add a version that determines the sites
 # from common site indices of ψ and A
 """
-    setindex!(ψ::Union{MPS, MPO}, A::ITensor, r::UnitRange{Int};
-              orthocenter::Int = last(r), perm = nothing, kwargs...)
+    setindex!(ψ::Union{MPS,MPO}, A::ITensor, r::UnitRange{Int};
+              orthocenter::Int=last(r), perm=nothing, kwargs...)
     replacesites!([...])
     replacesites([...])
 
@@ -1490,7 +1490,7 @@ function setindex!(
   # For MPO case, restrict to 0 prime level
   #sites = filter(hasplev(0), sites)
 
-  if !isnothing(perm)
+  if !isnothing(perm) # TODO: also check if it is a trivial permutation like [1,2,3]
     sites0 = sites
     sites = sites0[[perm...]]
     # Check if the site indices
@@ -1526,12 +1526,8 @@ function setindex!(
       end
     end
   end
-
   ψA = MPST(A, sites; leftinds=lind, orthocenter=orthocenter - first(r) + 1, kwargs...)
-  #@assert prod(ψA) ≈ A
-
   ψ[firstsite:lastsite] = ψA
-
   return ψ
 end
 
@@ -1570,30 +1566,34 @@ function (::Type{MPST})(
     @assert hasinds(A, s)
   end
   @assert isnothing(leftinds) || hasinds(A, leftinds)
-
   @assert 1 ≤ orthocenter ≤ N
 
+  l = isnothing(leftinds) ? Index[] : leftinds
+  # Need to flatten with reduce(vcat, sites)
+  r = uniqueinds(A, reduce(vcat, sites), l)
   ψ = Vector{ITensor}(undef, N)
   Ã = A
-  l = leftinds
-  # TODO: To minimize work, loop from
-  # 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="Link,n=$n", ortho="left")
-    l = commonind(L, R)
+
+  # Sweep left to orthocenter
+  for n in 1:(orthocenter - 1)
+    Lis = unioninds(sites[n], l)
+    L, R = factorize(Ã, Lis; kwargs..., tags="Link,l=$n", ortho="left")
+    l = commoninds(L, R)
     ψ[n] = L
     Ã = R
   end
-  ψ[N] = Ã
-  M = MPST(ψ)
-  setleftlim!(M, N - 1)
-  setrightlim!(M, N + 1)
-  orthogonalize!(M, orthocenter)
+
+  # Sweep right to orthocenter
+  for n in reverse((orthocenter + 1):N)
+    Ris = unioninds(sites[n], r)
+    R, L = factorize(Ã, Ris; kwargs..., tags="Link,l=$n", ortho="left")
+    r = commonind(R, L)
+    ψ[n] = R
+    Ã = L
+  end
+  ψ[orthocenter] = Ã
+  M = MPST(ψ; ortho_lims=orthocenter:orthocenter)
+  @assert ortho_lims(M) == orthocenter:orthocenter
   return M
 end
 
@@ -1666,13 +1666,19 @@ function _movesite(ns::Vector{Int}, n1n2::Pair{Int,Int})
   return ns
 end
 
-function _movesites(ψ::AbstractMPS, ns::Vector{Int}, ns′::Vector{Int}; kwargs...)
+function _movesites(ψ::AbstractMPS, ns::Vector{Int}, ns′::Vector{Int}; (nswaps!)=nothing, kwargs...)
   ψ = copy(ψ)
   N = length(ns)
   @assert N == length(ns′)
+  nswaps = 0
   for i in 1:N
-    ψ = movesite(ψ, ns[i] => ns′[i]; kwargs...)
-    ns = _movesite(ns, ns[i] => ns′[i])
+    n1, n2 = ns[i], ns′[i]
+    nswaps += n1 > n2 ? length(n2:(n1 - 1)) : length(n1:(n2 - 1))
+    ψ = movesite(ψ, n1 => n2; kwargs...)
+    ns = _movesite(ns, n1 => n2)
+  end
+  if !isnothing(nswaps!)
+    push!(nswaps!, nswaps)
   end
   return ψ, ns
 end
@@ -1680,7 +1686,7 @@ end
 # TODO: make a permutesites(::MPS/MPO, perm)
 # function that takes a permutation of the sites
 # p(1:N) for N sites
-function movesites(ψ::AbstractMPS, nsns′::Vector{Pair{Int,Int}}; kwargs...)
+function movesites(ψ::AbstractMPS, nsns′::Vector{Pair{Int,Int}}; (nswaps!)=nothing, kwargs...)
   ns = first.(nsns′)
   ns′ = last.(nsns′)
   ψ = copy(ψ)
@@ -1691,29 +1697,121 @@ function movesites(ψ::AbstractMPS, nsns′::Vector{Pair{Int,Int}}; kwargs...)
   ns′ = ns′[p]
   ns = collect(ns)
   while ns ≠ ns′
-    ψ, ns = _movesites(ψ, ns, ns′; kwargs...)
+    ψ, ns = _movesites(ψ, ns, ns′; (nswaps!)=nswaps!, kwargs...)
   end
   return ψ
 end
 
 # TODO: call the Vector{Pair{Int, Int}} version
-function movesites(ψ::AbstractMPS, ns, ns′; kwargs...)
-  ψ = copy(ψ)
+## function movesites(ψ::AbstractMPS, ns, ns′; kwargs...)
+##   ψ = copy(ψ)
+##   N = length(ns)
+##   @assert N == length(ns′)
+##   p = sortperm(ns′)
+##   ns = ns[p]
+##   ns′ = ns′[p]
+##   ns = collect(ns)
+##   for i in 1:N
+##     ψ = movesite(ψ, ns[i] => ns′[i]; kwargs...)
+##     ns = _movesite(ns, ns[i] => ns′[i])
+##   end
+##   return ψ
+## end
+
+eachfront(itr, n = 1) = Iterators.take(itr, length(itr) - n)
+
+function areconsecutive(v)
+  for n in eachfront(eachindex(v))
+    if v[n] + 1 ≠ v[n + 1]
+      return false
+    end
+  end
+  return true
+end
+
+abstract type SwapAlgorithm end
+struct SwapMinimal <: SwapAlgorithm end
+struct SwapMinimalSimple <: SwapAlgorithm end
+struct SwapLeft <: SwapAlgorithm end
+struct SwapRight <: SwapAlgorithm end
+struct SwapMiddle <: SwapAlgorithm end
+
+function consecutive_range(::SwapLeft, ns, ns_next)
   N = length(ns)
-  @assert N == length(ns′)
-  p = sortperm(ns′)
-  ns = ns[p]
-  ns′ = ns′[p]
-  ns = collect(ns)
-  for i in 1:N
-    ψ = movesite(ψ, ns[i] => ns′[i]; kwargs...)
-    ns = _movesite(ns, ns[i] => ns′[i])
+  new_start = first(ns)
+  return new_start:(new_start + N - 1)
+end
+
+function consecutive_range(::SwapRight, ns, ns_next)
+  N = length(ns)
+  new_stop = last(ns)
+  return (new_stop - N + 1):new_stop
+end
+
+function consecutive_range(::SwapMiddle, ns, ns_next)
+  N = length(ns)
+  new_start = (first(ns) + last(ns)) ÷ 2 - (length(ns) - 1) ÷ 2
+  return new_start:(new_start + N - 1)
+end
+
+function consecutive_range(::SwapMinimalSimple, ns, ns_next)
+  N = length(ns)
+  N_next = length(ns_next)
+  new_start = if last(ns) ≤ first(ns_next)
+    # Next gate is to the right of the current gate
+    # Move right towards the current gate
+    last(ns) - N + 1
+  elseif first(ns) ≥ last(ns_next)
+    # Next gate is to the left of the current gate
+    # Move left towards the current gate
+    first(ns)
+  else
+    # If neither is true, just swap left
+    first(ns)
   end
-  return ψ
+  return new_start:(new_start + N - 1)
+end
+
+consecutive_range(::SwapMinimalSimple, ns, ns_next::Nothing) = consecutive_range(SwapRight(), ns, ns_next)
+
+# Compute a contiguous range of sites 
+# based on the current range of sites the gate
+# is being applied to as well as the range of
+# sites the next gate will be applied to
+function consecutive_range(::SwapMinimal, ns, ns_next)
+  N = length(ns)
+  N_next = length(ns_next)
+  new_start = if last(ns) ≤ first(ns_next)
+    # Next gate is to the right of the current gate
+    # Move right towards the current gate
+    last(ns) - N + 1
+  elseif first(ns) ≥ last(ns_next)
+    # Next gate is to the left of the current gate
+    # Move left towards the current gate
+    first(ns)
+  elseif first(ns) ≤ first(ns_next) && last(ns) ≥ last(ns_next)
+    # Next gate is completely within the current gate,
+    # move inwards to the start of the current gate
+    first(ns_next)
+  elseif first(ns) ≥ first(ns_next) && last(ns) ≤ last(ns_next)
+    # Current gate is completely within the next gate
+    first(ns)
+  elseif first(ns) ≤ first(ns_next) && last(ns) ≤ last(ns_next)
+    # Next gate is overlapping to the right of the current gate
+    first(ns_next) - 1
+  elseif first(ns) ≥ first(ns_next) && last(ns) ≥ last(ns_next)
+    # Next gate is overlapping to the left of the current gate
+    last(ns_next) + 2 - N
+  else
+    error("The current gate is being applied to the sites $ns and the next gate is being applied to the sites $ns_next. Could not determine where to move the current sites to make the contiguous.")
+  end
+  return new_start:(new_start + N - 1)
 end
 
+consecutive_range(::SwapMinimal, ns, ns_next::Nothing) = consecutive_range(SwapRight(), ns, ns_next)
+
 """
-    product(o::ITensor, ψ::Union{MPS, MPO}, [ns::Vector{Int}]; <keyword argument>)
+    product(o::ITensor, ψ::Union{MPS,MPO}, [ns::Vector{Int}]; <keyword argument>)
     apply([...])
 
 Get the product of the operator `o` with the MPS/MPO `ψ`,
@@ -1736,28 +1834,65 @@ function product(
   ns=findsites(ψ, o);
   move_sites_back::Bool=true,
   apply_dag::Bool=false,
+  next_gate::Union{ITensor,Nothing}=nothing,
+  swap_alg="minimize",
+  (nswaps!)=nothing,
   kwargs...,
 )
   N = length(ns)
   ns = sort(ns)
 
-  # TODO: make this smarter by minimizing
-  # distance to orthogonalization.
-  # For example, if ITensors.orthocenter(ψ) > ns[end],
-  # set to ns[end].
-  ψ = orthogonalize(ψ, ns[1])
-  diff_ns = diff(ns)
+  ns_next = isnothing(next_gate) ? nothing : sort(findsites(ψ, next_gate))
   ns′ = ns
-  if any(!=(1), diff_ns)
-    ns′ = [ns[1] + n - 1 for n in 1:N]
-    ψ = movesites(ψ, ns .=> ns′; kwargs...)
+
+  if !areconsecutive(ns)
+    ns′ = if swap_alg == "minimize"
+      consecutive_range(SwapMinimal(), ns, ns_next)
+    elseif swap_alg == "simple_minimize"
+      consecutive_range(SwapMinimalSimple(), ns, ns_next)
+    elseif swap_alg == "left"
+      consecutive_range(SwapLeft(), ns, ns_next)
+    elseif swap_alg == "right"
+      consecutive_range(SwapRight(), ns, ns_next)
+    elseif swap_alg == "middle"
+      consecutive_range(SwapMiddle(), ns, ns_next)
+    end
+    ψ = movesites(ψ, ns .=> ns′; (nswaps!)=nswaps!, kwargs...)
+  else
+    if !isnothing(nswaps!)
+      push!(nswaps!, 0)
+    end
   end
-  ϕ = ψ[ns′[1]]
-  for n in 2:N
-    ϕ *= ψ[ns′[n]]
+  ns_range = first(ns′):last(ns′)
+  # Minimize distance to orthogonalize
+  if first(ortho_lims(ψ)) > last(ns_range)
+    ψ = orthogonalize(ψ, last(ns_range))
+  elseif last(ortho_lims(ψ)) < first(ns_range)
+    ψ = orthogonalize(ψ, first(ns_range))
+  elseif first(ortho_lims(ψ)) < first(ns_range)
+    ψ = orthogonalize(ψ, first(ns_range))
+  elseif last(ortho_lims(ψ)) > last(ns_range)
+    ψ = orthogonalize(ψ, last(ns_range))
   end
+
+  ϕ = prod(ψ[ns_range])
   ϕ = product(o, ϕ; apply_dag=apply_dag)
-  ψ[ns′[1]:ns′[end], kwargs...] = ϕ
+
+  # Anticipate where the next orthogonality
+  # center should be based on the position
+  # of the next gate being applied
+  orthocenter = if !isnothing(ns_next)
+    if first(ns_next) > last(ns_range)
+      last(ns_range)
+    elseif last(ns_next) < first(ns_range)
+      first(ns_range)
+    else
+      last(ns_range)
+    end
+  else
+    last(ns_range)
+  end
+  ψ[ns_range, orthocenter=orthocenter, kwargs...] = ϕ
   if move_sites_back
     # Move the sites back to their original positions
     ψ = movesites(ψ, ns′ .=> ns; kwargs...)
@@ -1852,8 +1987,10 @@ function product(
   As::Vector{<:ITensor}, ψ::AbstractMPS; move_sites_back::Bool=true, kwargs...
 )
   Aψ = ψ
-  for A in As
-    Aψ = product(A, Aψ; move_sites_back=false, kwargs...)
+  for n in eachindex(As)
+    An = As[n]
+    next_gate = n == lastindex(As) ? nothing : As[n + 1]
+    Aψ = product(An, Aψ; move_sites_back=false, next_gate=next_gate, kwargs...)
   end
   if move_sites_back
     s = siteinds(Aψ)
diff --git a/src/mps/mpo.jl b/src/mps/mpo.jl
index 0ab0cce1ca..0f12311b8b 100644
--- a/src/mps/mpo.jl
+++ b/src/mps/mpo.jl
@@ -457,8 +457,15 @@ function _contract_densitymatrix(A::MPO, ψ::MPS; kwargs...)::MPS
     error("In `contract(A::MPO, x::MPS)`, `A` and `x` must share a set of site indices")
 
   # In case A and ψ have the same link indices
+  # For example, ψ could be A projected onto a product state,
+  # like `|ψ⟩ = A|0⟩`.
   A = sim(linkinds, A)
 
+  # Dangling indices of ψ. For now, only support dangling
+  # indices on the edges of the system.
+  sψ = siteinds(uniqueinds, ψ, A)
+  any(n -> n ≠ firstindex(sψ) && n ≠ lastindex(sψ) && !isempty(sψ[n]), eachindex(sψ)) && error("In `contract(A::MPO, x::MPS)`, `x` can only have dangling indices on the first and last tensor.")
+
   ψ_c = dag(ψ)
   A_c = dag(A)
 
@@ -467,13 +474,17 @@ function _contract_densitymatrix(A::MPO, ψ::MPS; kwargs...)::MPS
   sim!(linkinds, ψ_c)
   sim!(siteinds, commoninds, A_c, ψ_c)
 
+  # Prime any dangling indices on the ends of ψ
+  prime!(siteinds, uniqueinds, ψ_c, A_c)
+
   # A version helpful for making the density matrix
   simA_c = sim(siteinds, uniqueinds, A_c, ψ_c)
 
   # Store the left environment tensors
   E = Vector{ITensor}(undef, n - 1)
-
-  E[1] = ψ[1] * A[1] * A_c[1] * ψ_c[1]
+  sψ1 = sψ[1]
+  sψ1_c = sψ1'
+  E[1] = ψ[1] * A[1] * A_c[1] * replaceinds(ψ_c[1], sψ1_c => sψ1)
   for j in 2:(n - 1)
     E[j] = E[j - 1] * ψ[j] * A[j] * A_c[j] * ψ_c[j]
   end
@@ -481,10 +492,13 @@ function _contract_densitymatrix(A::MPO, ψ::MPS; kwargs...)::MPS
   simR_c = ψ_c[n] * simA_c[n]
   ρ = E[n - 1] * R * simR_c
   l = linkind(ψ, n - 1)
-  ts = isnothing(l) ? "" : tags(l)
-  Lis = siteinds(uniqueinds, A, ψ, n)
-  Ris = siteinds(uniqueinds, simA_c, ψ_c, n)
-  F = eigen(ρ, Lis, Ris; ishermitian=true, tags=ts, kwargs...)
+  tsₙ = isnothing(l) ? "" : tags(l)
+  plₙ = isnothing(l) ? 0 : plev(l)
+  sψn = sψ[n]
+  sψn_c = sψn'
+  Lis = unioninds(siteinds(uniqueinds, A, ψ, n), sψn)
+  Ris = unioninds(siteinds(uniqueinds, simA_c, ψ_c, n), sψn_c)
+  F = eigen(ρ, Lis, Ris; ishermitian=true, tags=tsₙ, plev=plₙ, kwargs...)
   D, U, Ut = F.D, F.V, F.Vt
   l_renorm, r_renorm = F.l, F.r
   ψ_out[n] = Ut
@@ -495,10 +509,11 @@ function _contract_densitymatrix(A::MPO, ψ::MPS; kwargs...)::MPS
     s̃ = siteinds(uniqueinds, simA_c, ψ_c, j)
     ρ = E[j - 1] * R * simR_c
     l = linkind(ψ, j - 1)
-    ts = isnothing(l) ? "" : tags(l)
-    Lis = IndexSet(s..., l_renorm)
-    Ris = IndexSet(s̃..., r_renorm)
-    F = eigen(ρ, Lis, Ris; ishermitian=true, tags=ts, kwargs...)
+    tsⱼ = isnothing(l) ? "" : tags(l)
+    plⱼ = isnothing(l) ? 0 : plev(l)
+    Lis = [s..., l_renorm]
+    Ris = [s̃..., r_renorm]
+    F = eigen(ρ, Lis, Ris; ishermitian=true, tags=tsⱼ, plev=plⱼ, kwargs...)
     D, U, Ut = F.D, F.V, F.Vt
     l_renorm, r_renorm = F.l, F.r
     ψ_out[j] = Ut
@@ -511,6 +526,9 @@ function _contract_densitymatrix(A::MPO, ψ::MPS; kwargs...)::MPS
   ψ_out[1] = R
   setleftlim!(ψ_out, 0)
   setrightlim!(ψ_out, 2)
+  # TODO: set the prime levels to those
+  # of the link indices of ψ
+  setprime!(linkinds, ψ_out, 0)
   return ψ_out
 end
 
diff --git a/src/tagset.jl b/src/tagset.jl
index cd14725dbd..6618df38ce 100644
--- a/src/tagset.jl
+++ b/src/tagset.jl
@@ -56,11 +56,7 @@ end
 function _addtag!(ts::MTagSetStorage, ntags::Int, tag::IntTag)
   t = Tag(tag)
   if !isnull(t)
-    if isint(t)
-      error("Cannot use a bare integer as a tag.")
-    else
-      ntags = _addtag_ordered!(ts, ntags, tag)
-    end
+    ntags = _addtag_ordered!(ts, ntags, tag)
   end
   return ntags
 end
diff --git a/test/apply.jl b/test/apply.jl
new file mode 100644
index 0000000000..a91119bc07
--- /dev/null
+++ b/test/apply.jl
@@ -0,0 +1,137 @@
+using ITensors
+using Test
+
+@testset "apply" begin
+
+  @testset "ITensors.minimal_swap_range" begin
+    @test ITensors.consecutive_range(ITensors.SwapMinimal(), [1, 2], [1, 2]) == 1:2
+    @test ITensors.consecutive_range(ITensors.SwapMinimal(), [1, 2], [3, 4]) == 1:2
+    @test ITensors.consecutive_range(ITensors.SwapMinimal(), [1, 3], [5, 6]) == 2:3
+    @test ITensors.consecutive_range(ITensors.SwapMinimal(), [1, 4], [5, 6]) == 3:4
+    @test ITensors.consecutive_range(ITensors.SwapMinimal(), [3, 6], [1, 4]) == 4:5
+    @test ITensors.consecutive_range(ITensors.SwapMinimal(), [1, 5], [5, 6]) == 4:5
+    @test ITensors.consecutive_range(ITensors.SwapMinimal(), [1, 6], [5, 6]) == 5:6
+    @test ITensors.consecutive_range(ITensors.SwapMinimal(), [1, 6], [4, 5, 6]) == 4:5
+    @test ITensors.consecutive_range(ITensors.SwapMinimal(), [1, 3, 5], [7, 8]) == 3:5
+    @test ITensors.consecutive_range(ITensors.SwapMinimal(), [1, 5], [2, 3]) == 2:3
+    @test ITensors.consecutive_range(ITensors.SwapMinimal(), [1, 5], [2, 4]) == 2:3
+    @test ITensors.consecutive_range(ITensors.SwapMinimal(), [1, 3, 5], [2, 4]) == 2:4
+    @test ITensors.consecutive_range(ITensors.SwapMinimal(), [1, 3, 5], [5, 8]) == 3:5
+    @test ITensors.consecutive_range(ITensors.SwapMinimal(), [3, 4], [1, 2]) == 3:4
+    @test ITensors.consecutive_range(ITensors.SwapMinimal(), [4, 6], [1, 2]) == 4:5
+    @test ITensors.consecutive_range(ITensors.SwapMinimal(), [4, 6], [1, 3]) == 4:5
+    @test ITensors.consecutive_range(ITensors.SwapMinimal(), [4, 6], [1, 4]) == 4:5
+    @test ITensors.consecutive_range(ITensors.SwapMinimal(), [2, 6], [1, 6]) == 2:3
+    @test ITensors.consecutive_range(ITensors.SwapMinimal(), [2, 6], [1, 7]) == 2:3
+    @test ITensors.consecutive_range(ITensors.SwapMinimal(), [2, 6], [1, 5]) == 5:6
+    @test ITensors.consecutive_range(ITensors.SwapMinimal(), [1, 5], [2, 6]) == 1:2
+    @test ITensors.consecutive_range(ITensors.SwapMinimal(), [1, 3, 5], [2, 6]) == 1:3
+    @test ITensors.consecutive_range(ITensors.SwapMinimal(), [2, 5, 6], [1, 5]) == 4:6
+    @test ITensors.consecutive_range(ITensors.SwapMinimal(), [1, 5], [2, 6]) == 1:2
+    @test ITensors.consecutive_range(ITensors.SwapMinimal(), [4, 8], [3, 8]) == 4:5
+    @test ITensors.consecutive_range(ITensors.SwapMinimal(), [4, 6, 8], [3, 5]) == 4:6
+    @test ITensors.consecutive_range(ITensors.SwapMinimal(), [4, 6, 8], [3, 4]) == 4:6
+    @test ITensors.consecutive_range(ITensors.SwapMinimal(), [4, 6, 8], [3, 8]) == 4:6
+  end
+
+  @testset "apply swap" begin
+    N = 6
+    s = siteinds("Qubit", N)
+
+    Xop = [("X", n) for n in 1:N]
+    CXop = [("CX", i, j) for i in 1:N, j in 1:N]
+    CCXop = [("CCX", i, j, k) for i in 1:N, j in 1:N, k in 1:N]
+
+    X = op.(Xop, (s,))
+    CX = op.(CXop, (s,))
+    CCX = op.(CCXop, (s,))
+
+    ψ0 = MPS(s, "0")
+    orthogonalize!(ψ0, 1)
+    @test ortho_lims(ψ0) == 1:1
+
+    move_sites_back = false
+    kwargs = (; move_sites_back=move_sites_back)
+
+    nswaps = Int[]
+    ψ = apply(X[1], ψ0; (nswaps!)=nswaps, kwargs...)
+    @test ortho_lims(ψ) == 1:1
+    @test nswaps == [0]
+
+    nswaps = Int[]
+    ψ = apply(X[3], ψ0; (nswaps!)=nswaps, kwargs...)
+    @test ortho_lims(ψ) == 3:3
+    @test nswaps == [0]
+
+    nswaps = Int[]
+    ψ = apply(CX[1, 2], ψ0; (nswaps!)=nswaps, kwargs...)
+    @test ortho_lims(ψ) == 2:2
+    @test nswaps == [0]
+
+    nswaps = Int[]
+    ψ = apply(CX[2, 3], ψ0; (nswaps!)=nswaps, kwargs...)
+    @test ortho_lims(ψ) == 3:3
+    @test nswaps == [0]
+
+    nswaps = Int[]
+    ψ = apply(CX[3, 4], ψ0; (nswaps!)=nswaps, kwargs...)
+    @test ortho_lims(ψ) == 4:4
+    @test nswaps == [0]
+
+    nswaps = Int[]
+    ψ = apply([CX[3, 4], X[4]], ψ0; (nswaps!)=nswaps, kwargs...)
+    @test ortho_lims(ψ) == 4:4
+    @test nswaps == [0, 0]
+
+    nswaps = Int[]
+    ψ = apply(CX[1, 3], ψ0; (nswaps!)=nswaps, kwargs...)
+    @test ortho_lims(ψ) == 3:3
+    @test nswaps == [1]
+
+    nswaps = Int[]
+    ψ = apply(CX[1, 4], ψ0; (nswaps!)=nswaps, kwargs...)
+    @test ortho_lims(ψ) == 4:4
+    @test nswaps == [2]
+
+    nswaps = Int[]
+    ψ = apply([CCX[2, 3, 4], X[2]], ψ0; (nswaps!)=nswaps, kwargs...)
+    @test ortho_lims(ψ) == 2:2
+    @test nswaps == [0, 0]
+
+    nswaps = Int[]
+    ψ = apply([CX[3, 5], CX[1, 3]], ψ0; (nswaps!)=nswaps, kwargs...)
+
+    nswaps = Int[]
+    ψ = apply(CX[1, 4], ψ0; next_gate=CX[2, 3], (nswaps!)=nswaps, kwargs...)
+    @test nswaps == [2]
+
+    nswaps = Int[]
+    ψ = apply([CX[3, 6], CX[1, 3]], ψ0; (nswaps!)=nswaps, kwargs...)
+    @test nswaps == [2, 1]
+
+    @test ITensors.consecutive_range(ITensors.SwapMinimal(), [3, 6], [1, 4]) == 4:5
+    nswaps = Int[]
+    ψ = apply([CX[3, 6], CX[1, 4]], ψ0; (nswaps!)=nswaps, kwargs...)
+    @test nswaps == [2, 1]
+
+    # XXX broken
+    nswaps = Int[]
+    ψ = apply([CX[1, 3], CX[3, 5]], ψ0; next_gate=CX[3, 5], (nswaps!)=nswaps, kwargs...)
+    @test nswaps == [1, 1]
+
+    # XXX broken
+    nswaps = Int[]
+    ψ = apply([CX[1, 3], CX[3, 5]], ψ0; next_gate=CX[3, 5], (nswaps!)=nswaps, kwargs...)
+    @test nswaps == [1, 1]
+
+  end
+
+  @testset "Debugging" begin
+    N = 4
+    s = siteinds("Qubit", N)
+    X = [op("X", s, n) for n in 1:N]
+    ψ = MPS(s, "0")
+    @test prod(product(X[1], ψ)) ≈ prod(MPS(s, n -> n == 1 ? "1" : "0"))
+  end
+end
+
diff --git a/test/autompo.jl b/test/autompo.jl
index ff65919aea..ff6901f806 100644
--- a/test/autompo.jl
+++ b/test/autompo.jl
@@ -723,15 +723,15 @@ end
     a3 += "Cdag", 1, "N", 2, "C", 3
     M3 = MPO(a3, s)
 
-    p011 = productMPS(s, [1, 2, 2, 1, 1])
-    p110 = productMPS(s, [2, 2, 1, 1, 1])
+    p011 = MPS(s, [1, 2, 2, 1, 1])
+    p110 = MPS(s, [2, 2, 1, 1, 1])
 
     @test inner(p110, M1, p011) ≈ -1.0
     @test inner(p110, M2, p011) ≈ -1.0
     @test inner(p110, M3, p011) ≈ -1.0
 
-    p001 = productMPS(s, [1, 1, 2, 1, 1])
-    p100 = productMPS(s, [2, 1, 1, 1, 1])
+    p001 = MPS(s, [1, 1, 2, 1, 1])
+    p100 = MPS(s, [2, 1, 1, 1, 1])
 
     @test inner(p100, M1, p001) ≈ +1.0
     @test inner(p100, M2, p001) ≈ +1.0
@@ -752,14 +752,14 @@ end
     a2 += -1, "Cdn", 3, "Cdagdn", 1
     M2 = MPO(a2, s)
 
-    p0uu = productMPS(s, [1, 2, 2, 1, 1])
-    puu0 = productMPS(s, [2, 2, 1, 1, 1])
-    p0ud = productMPS(s, [1, 2, 3, 1, 1])
-    pdu0 = productMPS(s, [3, 2, 1, 1, 1])
-    p00u = productMPS(s, [1, 1, 2, 1, 1])
-    pu00 = productMPS(s, [2, 1, 1, 1, 1])
-    p00d = productMPS(s, [1, 1, 3, 1, 1])
-    pd00 = productMPS(s, [3, 1, 1, 1, 1])
+    p0uu = MPS(s, [1, 2, 2, 1, 1])
+    puu0 = MPS(s, [2, 2, 1, 1, 1])
+    p0ud = MPS(s, [1, 2, 3, 1, 1])
+    pdu0 = MPS(s, [3, 2, 1, 1, 1])
+    p00u = MPS(s, [1, 1, 2, 1, 1])
+    pu00 = MPS(s, [2, 1, 1, 1, 1])
+    p00d = MPS(s, [1, 1, 3, 1, 1])
+    pd00 = MPS(s, [3, 1, 1, 1, 1])
 
     @test inner(puu0, M1, p0uu) ≈ -1.0
     @test inner(pdu0, M2, p0ud) ≈ -1.0
@@ -778,8 +778,8 @@ end
         ampo += -1im, "S-", i, "S+", i + 1
       end
       H = MPO(ampo, sites)
-      psiud = productMPS(sites, [1, 2, 1, 2])
-      psidu = productMPS(sites, [2, 1, 1, 2])
+      psiud = MPS(sites, [1, 2, 1, 2])
+      psidu = MPS(sites, [2, 1, 1, 2])
       @test inner(psiud, H, psidu) ≈ +1im
       @test inner(psidu, H, psiud) ≈ -1im
     end
diff --git a/test/dmrg.jl b/test/dmrg.jl
index 25c9185898..ca3f6971a1 100644
--- a/test/dmrg.jl
+++ b/test/dmrg.jl
@@ -40,7 +40,7 @@ using ITensors, Test, Random
     H = MPO(ampo, sites)
 
     state = [isodd(n) ? "Up" : "Dn" for n in 1:N]
-    psi = randomMPS(sites, state, 4)
+    psi = randomMPS(sites, state; linkdims=4)
 
     sweeps = Sweeps(3)
     @test length(sweeps) == 3
@@ -67,7 +67,7 @@ using ITensors, Test, Random
     H = MPO(ampo, sites)
 
     state = [isodd(n) ? "Up" : "Dn" for n in 1:N]
-    psi = randomMPS(sites, state, 4)
+    psi = randomMPS(sites, state; linkdims=4)
 
     sweeps = Sweeps(3)
     @test length(sweeps) == 3
@@ -94,7 +94,7 @@ using ITensors, Test, Random
     H = MPO(ampo, sites)
 
     state = [isodd(n) ? "Up" : "Dn" for n in 1:N]
-    psi = randomMPS(sites, state, 4)
+    psi = randomMPS(sites, state; linkdims=4)
     PH = ProjMPO(H)
 
     n = 4
@@ -254,7 +254,7 @@ using ITensors, Test, Random
     end
     H = MPO(ampo, sites)
 
-    psi0i = randomMPS(sites, 10)
+    psi0i = randomMPS(sites; linkdims=10)
 
     sweeps = Sweeps(4)
     maxdim!(sweeps, 10, 20, 100, 100)
@@ -264,7 +264,7 @@ using ITensors, Test, Random
     energy0, psi0 = dmrg(H, psi0i, sweeps; outputlevel=0)
     @test energy0 < -11.5
 
-    psi1i = randomMPS(sites, 10)
+    psi1i = randomMPS(sites; linkdims=10)
     energy1, psi1 = dmrg(H, [psi0], psi1i, sweeps; outputlevel=0, weight=weight)
 
     @test energy1 > energy0
@@ -285,7 +285,7 @@ using ITensors, Test, Random
     state[3] = 2
     state[5] = 2
     state[7] = 2
-    psi0 = productMPS(s, state)
+    psi0 = MPS(s, state)
 
     ampo = OpSum()
     for j in 1:(N - 1)
@@ -331,7 +331,7 @@ using ITensors, Test, Random
     maxdim!(sweeps, 50, 100, 200, 400, 800, 800)
     cutoff!(sweeps, 1E-10)
     state = ["Up", "Dn", "Dn", "Up", "Emp", "Up", "Up", "Emp", "Dn", "Dn"]
-    psi0 = randomMPS(sites, state, 10)
+    psi0 = randomMPS(sites, state; linkdims=10)
     energy, psi = dmrg(H, psi0, sweeps; outputlevel=0)
     @test (-8.02 < energy < -8.01)
   end
@@ -352,7 +352,7 @@ using ITensors, Test, Random
     maxdim!(sweeps, 10)
     cutoff!(sweeps, 1E-11)
 
-    psi0 = randomMPS(sites, 4)
+    psi0 = randomMPS(sites; linkdims=4)
 
     # Test that input works with wrong ortho center:
     orthogonalize!(psi0, 5)
@@ -395,7 +395,7 @@ using ITensors, Test, Random
     H = MPO(ampo, sites)
     Hsplit = splitblocks(linkinds, H)
     state = [isodd(n) ? "↑" : "↓" for n in 1:N]
-    ψ0 = productMPS(sites, state)
+    ψ0 = MPS(sites, state)
     using_threaded_blocksparse = ITensors.enable_threaded_blocksparse()
     energy, _ = dmrg(H, ψ0, sweeps; outputlevel=outputlevel)
     energy_split, _ = dmrg(Hsplit, ψ0, sweeps; outputlevel=outputlevel)
diff --git a/test/mps.jl b/test/mps.jl
index a1b0cc0bb2..b05f0292d1 100644
--- a/test/mps.jl
+++ b/test/mps.jl
@@ -60,7 +60,7 @@ include("util.jl")
     @test psi ⋅ psi ≈ *(dag(psi)..., psi...)[]
   end
 
-  @testset "productMPS" begin
+  @testset "MPS" begin
     @testset "vector of string input" begin
       sites = siteinds("S=1/2", N)
       state = fill("", N)
@@ -72,13 +72,13 @@ include("util.jl")
         sign = isodd(j) ? +1.0 : -1.0
         @test (psi[j] * op(sites, "Sz", j) * dag(prime(psi[j], "Site")))[] ≈ sign / 2
       end
-      psi = productMPS(sites, state)
+      psi = MPS(sites, state)
       for j in 1:N
         sign = isodd(j) ? +1.0 : -1.0
         @test (psi[j] * op(sites, "Sz", j) * dag(prime(psi[j], "Site")))[] ≈ sign / 2
       end
       @test_throws DimensionMismatch MPS(sites, fill("", N - 1))
-      @test_throws DimensionMismatch productMPS(sites, fill("", N - 1))
+      @test_throws DimensionMismatch MPS(sites, fill("", N - 1))
     end
 
     @testset "String input" begin
@@ -88,7 +88,7 @@ include("util.jl")
         sign = -1.0
         @test (psi[j] * op(sites, "Sz", j) * dag(prime(psi[j], "Site")))[] ≈ sign / 2
       end
-      psi = productMPS(sites, "Dn")
+      psi = MPS(sites, "Dn")
       for j in 1:N
         sign = -1.0
         @test (psi[j] * op(sites, "Sz", j) * dag(prime(psi[j], "Site")))[] ≈ sign / 2
@@ -102,7 +102,7 @@ include("util.jl")
         sign = -1.0
         @test (psi[j] * op(sites, "Sz", j) * dag(prime(psi[j], "Site")))[] ≈ sign / 2
       end
-      psi = productMPS(sites, 2)
+      psi = MPS(sites, 2)
       for j in 1:N
         sign = -1.0
         @test (psi[j] * op(sites, "Sz", j) * dag(prime(psi[j], "Site")))[] ≈ sign / 2
@@ -120,7 +120,7 @@ include("util.jl")
         sign = isodd(j) ? +1.0 : -1.0
         @test (psi[j] * op(sites, "Sz", j) * dag(prime(psi[j], "Site")))[] ≈ sign / 2
       end
-      psi = productMPS(sites, state)
+      psi = MPS(sites, state)
       for j in 1:N
         sign = isodd(j) ? +1.0 : -1.0
         @test (psi[j] * op(sites, "Sz", j) * dag(prime(psi[j], "Site")))[] ≈ sign / 2
@@ -139,7 +139,7 @@ include("util.jl")
         sign = isodd(j) ? +1.0 : -1.0
         @test (psi[j] * op(sites, "Sz", j) * dag(prime(psi[j], "Site")))[] ≈ sign / 2
       end
-      psi = productMPS(ivals)
+      psi = MPS(ivals)
       for j in 1:N
         sign = isodd(j) ? +1.0 : -1.0
         @test (psi[j] * op(sites, "Sz", j) * dag(prime(psi[j], "Site")))[] ≈ sign / 2
@@ -151,7 +151,7 @@ include("util.jl")
         for j in 1:N
           @test eltype(psi[j]) == ComplexF64
         end
-        psi = productMPS(ComplexF64, sites, fill(1, N))
+        psi = MPS(ComplexF64, sites, fill(1, N))
         for j in 1:N
           @test eltype(psi[j]) == ComplexF64
         end
@@ -165,13 +165,13 @@ include("util.jl")
       psi = MPS([site], [1])
       @test psi[1][1] ≈ 1.0
       @test psi[1][2] ≈ 0.0
-      psi = productMPS([site], [1])
+      psi = MPS([site], [1])
       @test psi[1][1] ≈ 1.0
       @test psi[1][2] ≈ 0.0
       psi = MPS([site], [2])
       @test psi[1][1] ≈ 0.0
       @test psi[1][2] ≈ 1.0
-      psi = productMPS([site], [2])
+      psi = MPS([site], [2])
       @test psi[1][1] ≈ 0.0
       @test psi[1][2] ≈ 1.0
     end
@@ -260,7 +260,7 @@ include("util.jl")
   end
 
   @testset "norm MPS" begin
-    psi = randomMPS(sites, 10)
+    psi = randomMPS(sites; linkdims=10)
     psidag = sim(linkinds, dag(psi))
     psi² = ITensor(1)
     for j in 1:N
@@ -269,11 +269,11 @@ include("util.jl")
     @test psi²[] ≈ psi ⋅ psi
     @test sqrt(psi²[]) ≈ norm(psi)
 
-    psi = randomMPS(sites, 10)
+    psi = randomMPS(sites; linkdims=10)
     psi .*= 1:N
     @test norm(psi) ≈ factorial(N)
 
-    psi = randomMPS(sites, 10)
+    psi = randomMPS(sites; linkdims=10)
     for j in 1:N
       psi[j] .*= j
     end
@@ -285,7 +285,7 @@ include("util.jl")
   end
 
   @testset "lognorm MPS" begin
-    psi = randomMPS(sites, 10)
+    psi = randomMPS(sites; linkdims=10)
     for j in 1:N
       psi[j] .*= j
     end
@@ -351,9 +351,9 @@ include("util.jl")
     s = siteinds("S=1/2", N; conserve_qns=conserve_qns)
     state = n -> isodd(n) ? "↑" : "↓"
 
-    ψ₁ = randomMPS(s, state, 4)
-    ψ₂ = randomMPS(s, state, 4)
-    ψ₃ = randomMPS(s, state, 4)
+    ψ₁ = randomMPS(s, state; linkdims=4)
+    ψ₂ = randomMPS(s, state; linkdims=4)
+    ψ₃ = randomMPS(s, state; linkdims=4)
 
     ψ = ψ₁ + ψ₂
 
@@ -535,7 +535,7 @@ end
   @testset "sample! method" begin
     N = 10
     sites = [Index(3, "Site,n=$n") for n in 1:N]
-    psi = randomMPS(sites, 3)
+    psi = randomMPS(sites; linkdims=3)
     nrm2 = inner(psi, psi)
     psi[1] *= (1.0 / sqrt(nrm2))
 
@@ -566,7 +566,7 @@ end
     N = 20
     chi = 8
     sites = siteinds(2, N)
-    M = randomMPS(sites, chi)
+    M = randomMPS(sites; linkdims=chi)
 
     @test ITensors.leftlim(M) == 0
     @test ITensors.rightlim(M) == 2
@@ -586,7 +586,7 @@ end
     end
 
     # Complex case
-    Mc = randomMPS(sites, chi)
+    Mc = randomMPS(sites; linkdims=chi)
     @test inner(Mc, Mc) ≈ 1.0 + 0.0im
   end
 
@@ -597,7 +597,7 @@ end
 
     # Make flux-zero random MPS
     state = [isodd(n) ? 1 : 2 for n in 1:N]
-    M = randomMPS(sites, state, chi)
+    M = randomMPS(sites, state; linkdims=chi)
     @test flux(M) == QN("Sz", 0)
 
     @test ITensors.leftlim(M) == 0
@@ -610,10 +610,10 @@ end
 
     # Test making random MPS with different flux
     state[1] = 2
-    M = randomMPS(sites, state, chi)
+    M = randomMPS(sites, state; linkdims=chi)
     @test flux(M) == QN("Sz", -2)
     state[3] = 2
-    M = randomMPS(sites, state, chi)
+    M = randomMPS(sites, state; linkdims=chi)
     @test flux(M) == QN("Sz", -4)
   end
 
@@ -623,7 +623,7 @@ end
 
     # Non-fermionic case - spin system
     s = siteinds("S=1/2", N; conserve_qns=true)
-    psi = randomMPS(s, n -> isodd(n) ? "Up" : "Dn", m)
+    psi = randomMPS(s, n -> isodd(n) ? "Up" : "Dn"; linkdims=m)
     Cpm = correlation_matrix(psi, "S+", "S-")
     # Check using OpSum:
     for i in 1:N, j in i:N
@@ -639,7 +639,7 @@ end
 
     # With start_site, end_site arguments:
     s = siteinds("S=1/2", N)
-    psi = randomMPS(ComplexF64, s, m)
+    psi = randomMPS(ComplexF64, s; linkdims=m)
     ss, es = 3, 6
     Nb = es - ss + 1
     Cpm = correlation_matrix(psi, "S+", "S-"; site_range=ss:es)
@@ -654,7 +654,7 @@ end
 
     # Fermionic case
     s = siteinds("Electron", N)
-    psi = randomMPS(s, m)
+    psi = randomMPS(s; linkdims=m)
     Cuu = correlation_matrix(psi, "Cdagup", "Cup")
     # Check using OpSum:
     for i in 1:N, j in i:N
@@ -874,7 +874,7 @@ end
     @test ITensors.orthocenter(ψ) == N
     @test maxlinkdim(ψ) == 1
 
-    ψ0 = randomMPS(s, 2)
+    ψ0 = randomMPS(s; linkdims=2)
     A = prod(ψ0)
     ψ = MPS(A, s; cutoff=1e-15, orthocenter=2)
     @test prod(ψ) ≈ A
@@ -914,7 +914,7 @@ end
   @testset "Set range of MPS tensors" begin
     N = 5
     s = siteinds("S=1/2", N)
-    ψ0 = randomMPS(s, 3)
+    ψ0 = randomMPS(s; linkdims=3)
 
     ψ = orthogonalize(ψ0, 2)
     A = prod(ITensors.data(ψ)[2:(N - 1)])
@@ -1545,7 +1545,7 @@ end
   @testset "dense conversion of MPS" begin
     N = 4
     s = siteinds("S=1/2", N; conserve_qns=true)
-    QM = randomMPS(s, ["Up", "Dn", "Up", "Dn"], 4)
+    QM = randomMPS(s, ["Up", "Dn", "Up", "Dn"]; linkdims=4)
     qsz1 = scalar(QM[1] * op("Sz", s[1]) * dag(prime(QM[1], "Site")))
 
     M = dense(QM)
@@ -1579,7 +1579,7 @@ end
       a .+= "Sz", j, "Sz", j + 1
     end
     H = MPO(a, s)
-    ψ = randomMPS(s, n -> isodd(n) ? "↑" : "↓", 10)
+    ψ = randomMPS(s, n -> isodd(n) ? "↑" : "↓"; linkdims=10)
     # Create MPO/MPS with pairs of sites merged
     H2 = MPO([H[b] * H[b + 1] for b in 1:2:N])
     ψ2 = MPS([ψ[b] * ψ[b + 1] for b in 1:2:N])
diff --git a/test/runtests.jl b/test/runtests.jl
index acc03f4c1f..f2324ff862 100644
--- a/test/runtests.jl
+++ b/test/runtests.jl
@@ -32,6 +32,7 @@ using Test
     "lattices.jl",
     "mps.jl",
     "mpo.jl",
+    "apply.jl",
     "sweeps.jl",
     "sweepnext.jl",
     "autompo.jl",