diff --git a/ITensorGPU/test/runtests.jl b/ITensorGPU/test/runtests.jl
index 9a2a90d73e..71a592139e 100644
--- a/ITensorGPU/test/runtests.jl
+++ b/ITensorGPU/test/runtests.jl
@@ -21,7 +21,6 @@ CUDA.allowscalar(false)
     "test_cucontract.jl",
     "test_cumpo.jl",
     "test_cumps.jl",
-    "test_cuiterativesolvers.jl",
     "test_cutruncate.jl",
     #"test_pastaq.jl",
   )
diff --git a/ITensorGPU/test/test_cuiterativesolvers.jl b/ITensorGPU/test/test_cuiterativesolvers.jl
deleted file mode 100644
index 0bb4d32683..0000000000
--- a/ITensorGPU/test/test_cuiterativesolvers.jl
+++ /dev/null
@@ -1,27 +0,0 @@
-using ITensorGPU, ITensors
-using Test
-
-# Wrap an ITensor with pairs of primed and
-# unprimed indices to pass to davidson
-struct ITensorMap
-  A::ITensor
-end
-Base.eltype(M::ITensorMap) = eltype(M.A)
-Base.size(M::ITensorMap) = dim(IndexSet(inds(M.A; plev=0)...))
-(M::ITensorMap)(v::ITensor) = noprime(M.A * v)
-
-@testset "Complex davidson" begin
-  d = 10
-  i = Index(d, "i")
-  A = cuITensor(randomITensor(ComplexF64, i, prime(i)))
-  A = mapprime(A * mapprime(dag(A), 0, 2), 2, 1)
-  M = ITensorMap(A)
-
-  v = randomCuITensor(i)
-  λ, v = davidson(M, v; maxiter=10)
-  @test M(v) ≈ λ * v
-
-  v = cuITensor(randomITensor(ComplexF64, i))
-  λ, v = davidson(M, v; maxiter=10)
-  @test M(v) ≈ λ * v
-end
diff --git a/Project.toml b/Project.toml
index 6b82e6be7a..255e178f4a 100644
--- a/Project.toml
+++ b/Project.toml
@@ -21,6 +21,7 @@ PackageCompiler = "9b87118b-4619-50d2-8e1e-99f35a4d4d9d"
 Pkg = "44cfe95a-1eb2-52ea-b672-e2afdf69b78f"
 Printf = "de0858da-6303-5e67-8744-51eddeeeb8d7"
 Random = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c"
+Reexport = "189a3867-3050-52da-a836-e630ba90ab69"
 Requires = "ae029012-a4dd-5104-9daa-d747884805df"
 SerializedElementArrays = "d3ce8812-9567-47e9-a7b5-65a6d70a3065"
 SimpleTraits = "699a6c99-e7fa-54fc-8d76-47d257e15c1d"
diff --git a/src/ITensorMPS/ITensorMPS.jl b/src/ITensorMPS/ITensorMPS.jl
new file mode 100644
index 0000000000..f8fdb9dc07
--- /dev/null
+++ b/src/ITensorMPS/ITensorMPS.jl
@@ -0,0 +1,26 @@
+module ITensorMPS
+
+include("usings.jl")
+include("imports.jl")
+include("exports.jl")
+include("abstractmps.jl")
+include("mps.jl")
+include("mpo.jl")
+include("sweeps.jl")
+include("abstractprojmpo/abstractprojmpo.jl")
+include("abstractprojmpo/projmpo.jl")
+include("abstractprojmpo/diskprojmpo.jl")
+include("abstractprojmpo/projmposum.jl")
+include("abstractprojmpo/projmps.jl")
+include("abstractprojmpo/projmpo_mps.jl")
+include("observer.jl")
+include("dmrg.jl")
+include("adapt.jl")
+include("autompo/matelem.jl")
+include("autompo/qnmatelem.jl")
+include("autompo/opsum_to_mpo_generic.jl")
+include("autompo/opsum_to_mpo.jl")
+include("autompo/opsum_to_mpo_qn.jl")
+include("deprecated.jl")
+
+end # module ITensorMPS
diff --git a/src/mps/abstractmps.jl b/src/ITensorMPS/abstractmps.jl
similarity index 99%
rename from src/mps/abstractmps.jl
rename to src/ITensorMPS/abstractmps.jl
index 1fdf589b07..d5d38531d0 100644
--- a/src/mps/abstractmps.jl
+++ b/src/ITensorMPS/abstractmps.jl
@@ -79,7 +79,7 @@ know if there is functionality not available for MPS/MPO you would like.
 """
 data(m::AbstractMPS) = m.data
 
-contract(ψ::AbstractMPS) = contract(data(ψ))
+ITensors.contract(ψ::AbstractMPS) = contract(data(ψ))
 
 leftlim(m::AbstractMPS) = m.llim
 
@@ -1076,7 +1076,7 @@ function deprecate_make_inds_match!(
     make_inds_match = false
   end
   if !hassameinds(siteinds, M1dag, M2) && make_inds_match
-    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
@@ -1776,7 +1776,7 @@ function _fermionic_swap(s1::Index, s2::Index)
     if _isodd_fermionic_parity(s1, n1) && _isodd_fermionic_parity(s2, n2)
       dval = -1.0
     end
-    Tb = blockview(tensor(T), b)
+    Tb = ITensors.blockview(tensor(T), b)
     mat_dim = prod(dims(Tb)[1:2])
     Tbr = reshape(Tb, mat_dim, mat_dim)
     for i in diagind(Tbr)
@@ -1847,9 +1847,9 @@ function setindex!(
     sites = sites0[[perm...]]
     # Check if the site indices
     # are fermionic
-    if !using_auto_fermion() && any(anyfermionic, sites)
+    if !using_auto_fermion() && any(ITensors.anyfermionic, sites)
       if length(sites) == 2 && ψ isa MPS
-        if all(allfermionic, sites)
+        if all(ITensors.allfermionic, sites)
           s0 = Index.(sites0)
 
           # TODO: the Fermionic swap is could be diagonal,
diff --git a/src/mps/abstractprojmpo/abstractprojmpo.jl b/src/ITensorMPS/abstractprojmpo/abstractprojmpo.jl
similarity index 98%
rename from src/mps/abstractprojmpo/abstractprojmpo.jl
rename to src/ITensorMPS/abstractprojmpo/abstractprojmpo.jl
index 7cb07a0252..05aae2d152 100644
--- a/src/mps/abstractprojmpo/abstractprojmpo.jl
+++ b/src/ITensorMPS/abstractprojmpo/abstractprojmpo.jl
@@ -33,7 +33,7 @@ function rproj(P::AbstractProjMPO)::Union{ITensor,OneITensor}
   return P.LR[P.rpos]
 end
 
-function contract(P::AbstractProjMPO, v::ITensor)::ITensor
+function ITensors.contract(P::AbstractProjMPO, v::ITensor)::ITensor
   itensor_map = Union{ITensor,OneITensor}[lproj(P)]
   append!(itensor_map, P.H[site_range(P)])
   push!(itensor_map, rproj(P))
diff --git a/src/mps/abstractprojmpo/diskprojmpo.jl b/src/ITensorMPS/abstractprojmpo/diskprojmpo.jl
similarity index 100%
rename from src/mps/abstractprojmpo/diskprojmpo.jl
rename to src/ITensorMPS/abstractprojmpo/diskprojmpo.jl
diff --git a/src/mps/abstractprojmpo/projmpo.jl b/src/ITensorMPS/abstractprojmpo/projmpo.jl
similarity index 100%
rename from src/mps/abstractprojmpo/projmpo.jl
rename to src/ITensorMPS/abstractprojmpo/projmpo.jl
diff --git a/src/mps/abstractprojmpo/projmpo_mps.jl b/src/ITensorMPS/abstractprojmpo/projmpo_mps.jl
similarity index 100%
rename from src/mps/abstractprojmpo/projmpo_mps.jl
rename to src/ITensorMPS/abstractprojmpo/projmpo_mps.jl
diff --git a/src/mps/abstractprojmpo/projmposum.jl b/src/ITensorMPS/abstractprojmpo/projmposum.jl
similarity index 98%
rename from src/mps/abstractprojmpo/projmposum.jl
rename to src/ITensorMPS/abstractprojmpo/projmposum.jl
index 55e7f6503a..c47ffc368f 100644
--- a/src/mps/abstractprojmpo/projmposum.jl
+++ b/src/ITensorMPS/abstractprojmpo/projmposum.jl
@@ -44,7 +44,7 @@ shorthand for `product(P,v)`.
 """
 product(A::AbstractSum, v::ITensor) = sum(t -> product(t, v), terms(A))
 
-contract(A::AbstractSum, v::ITensor) = sum(t -> contract(t, v), terms(A))
+ITensors.contract(A::AbstractSum, v::ITensor) = sum(t -> contract(t, v), terms(A))
 
 """
     eltype(P::ProjMPOSum)
diff --git a/src/mps/abstractprojmpo/projmps.jl b/src/ITensorMPS/abstractprojmpo/projmps.jl
similarity index 100%
rename from src/mps/abstractprojmpo/projmps.jl
rename to src/ITensorMPS/abstractprojmpo/projmps.jl
diff --git a/src/mps/adapt.jl b/src/ITensorMPS/adapt.jl
similarity index 100%
rename from src/mps/adapt.jl
rename to src/ITensorMPS/adapt.jl
diff --git a/src/physics/autompo/matelem.jl b/src/ITensorMPS/autompo/matelem.jl
similarity index 100%
rename from src/physics/autompo/matelem.jl
rename to src/ITensorMPS/autompo/matelem.jl
diff --git a/src/physics/autompo/opsum_to_mpo.jl b/src/ITensorMPS/autompo/opsum_to_mpo.jl
similarity index 100%
rename from src/physics/autompo/opsum_to_mpo.jl
rename to src/ITensorMPS/autompo/opsum_to_mpo.jl
diff --git a/src/physics/autompo/opsum_to_mpo_generic.jl b/src/ITensorMPS/autompo/opsum_to_mpo_generic.jl
similarity index 97%
rename from src/physics/autompo/opsum_to_mpo_generic.jl
rename to src/ITensorMPS/autompo/opsum_to_mpo_generic.jl
index 05081b49ac..0cd310993c 100644
--- a/src/physics/autompo/opsum_to_mpo_generic.jl
+++ b/src/ITensorMPS/autompo/opsum_to_mpo_generic.jl
@@ -68,7 +68,7 @@ subtract!(os::OpSum, o::Tuple) = add!(os, -Ops.op_term(o))
 function isfermionic(t::Vector{Op}, sites)
   p = +1
   for op in t
-    if has_fermion_string(name(op), sites[site(op)])
+    if has_fermion_string(ITensors.name(op), sites[site(op)])
       p *= -1
     end
   end
@@ -207,7 +207,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 *= parity_sign(perm)
+    t *= ITensors.parity_sign(perm)
     terms(os)[j] = t
   end
 
@@ -299,8 +299,8 @@ function MPO(os::OpSum, sites::Vector{<:Index}; splitblocks=true, kwargs...)::MP
   return M
 end
 
-function MPO(eltype::Type{<:Number}, os::OpSum, sites::Vector{<:Index}; kwargs...)
-  return NDTensors.convert_scalartype(eltype, MPO(os, sites; kwargs...))
+function MPO(elt::Type{<:Number}, os::OpSum, sites::Vector{<:Index}; kwargs...)
+  return NDTensors.convert_scalartype(elt, MPO(os, sites; kwargs...))
 end
 
 # Conversion from other formats
diff --git a/src/physics/autompo/opsum_to_mpo_qn.jl b/src/ITensorMPS/autompo/opsum_to_mpo_qn.jl
similarity index 96%
rename from src/physics/autompo/opsum_to_mpo_qn.jl
rename to src/ITensorMPS/autompo/opsum_to_mpo_qn.jl
index 316a875491..0149ea7634 100644
--- a/src/physics/autompo/opsum_to_mpo_qn.jl
+++ b/src/ITensorMPS/autompo/opsum_to_mpo_qn.jl
@@ -99,7 +99,7 @@ function qn_svdMPO(
   llinks = Vector{QNIndex}(undef, N + 1)
   # Set dir=In for fermionic ordering, avoid arrow sign
   # <fermions>:
-  linkdir = using_auto_fermion() ? In : Out
+  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}}()
@@ -220,7 +220,7 @@ function qn_svdMPO(
         end
 
         b = loc(rq, cq)
-        T = BlockSparseTensor(ValType, [b], (dag(ll), rl))
+        T = ITensors.NDTensors.BlockSparseTensor(ValType, [b], (dag(ll), rl))
         T[b] .= M
 
         H[n] += (itensor(T) * Op)
@@ -230,13 +230,13 @@ function qn_svdMPO(
     # Put in ending identity operator
     Id = op("Id", sites[n])
     b = Block(1, 1)
-    T = BlockSparseTensor(ValType, [b], (dag(ll), rl))
+    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 = BlockSparseTensor(ValType, [b], (dag(ll), rl))
+    T = ITensors.NDTensors.BlockSparseTensor(ValType, [b], (dag(ll), rl))
     T[b] = 1
     H[n] += (itensor(T) * Id)
   end # for n in 1:N
diff --git a/src/physics/autompo/qnmatelem.jl b/src/ITensorMPS/autompo/qnmatelem.jl
similarity index 100%
rename from src/physics/autompo/qnmatelem.jl
rename to src/ITensorMPS/autompo/qnmatelem.jl
diff --git a/src/mps/deprecated.jl b/src/ITensorMPS/deprecated.jl
similarity index 58%
rename from src/mps/deprecated.jl
rename to src/ITensorMPS/deprecated.jl
index 97b88e3b7c..e94a6bfd7c 100644
--- a/src/mps/deprecated.jl
+++ b/src/ITensorMPS/deprecated.jl
@@ -1,3 +1,37 @@
+# mps/abstractmps.jl
+@deprecate orthoCenter(args...; kwargs...) orthocenter(args...; kwargs...)
+@deprecate store(m::AbstractMPS) data(m)
+@deprecate replacesites!(args...; kwargs...) ITensors.replace_siteinds!(args...; kwargs...)
+@deprecate applyMPO(args...; kwargs...) contract(args...; kwargs...)
+@deprecate applympo(args...; kwargs...) contract(args...; kwargs...)
+@deprecate errorMPOprod(args...; kwargs...) error_contract(args...; kwargs...)
+@deprecate error_mpoprod(args...; kwargs...) error_contract(args...; kwargs...)
+@deprecate error_mul(args...; kwargs...) error_contract(args...; kwargs...)
+@deprecate multMPO(args...; kwargs...) contract(args...; kwargs...)
+@deprecate sum(A::AbstractMPS, B::AbstractMPS; kwargs...) add(A, B; kwargs...)
+@deprecate multmpo(args...; kwargs...) contract(args...; kwargs...)
+@deprecate set_leftlim!(args...; kwargs...) ITensors.setleftlim!(args...; kwargs...)
+@deprecate set_rightlim!(args...; kwargs...) ITensors.setrightlim!(args...; kwargs...)
+@deprecate tensors(args...; kwargs...) ITensors.data(args...; kwargs...)
+@deprecate primelinks!(args...; kwargs...) ITensors.prime_linkinds!(args...; kwargs...)
+@deprecate simlinks!(args...; kwargs...) ITensors.sim_linkinds!(args...; kwargs...)
+@deprecate mul(A::AbstractMPS, B::AbstractMPS; kwargs...) contract(A, B; kwargs...)
+
+# mps/mpo.jl
+@deprecate MPO(A::MPS; kwargs...) outer(A', A; kwargs...)
+
+# mps/mps.jl
+@deprecate randomMPS(sites::Vector{<:Index}, linkdims::Integer) randomMPS(
+  sites; linkdims=linkdims
+)
+@deprecate randomMPS(ElType::Type, sites::Vector{<:Index}, linkdims::Integer) randomMPS(
+  ElType, sites; linkdims=linkdims
+)
+@deprecate randomMPS(sites::Vector{<:Index}, state, linkdims::Integer) randomMPS(
+  sites, state; linkdims=linkdims
+)
+
+# Pass throughs of old name to new name:
 
 unique_siteind(A::AbstractMPS, B::AbstractMPS, j::Integer) = siteinds(uniqueind, A, B, j)
 
diff --git a/src/mps/dmrg.jl b/src/ITensorMPS/dmrg.jl
similarity index 99%
rename from src/mps/dmrg.jl
rename to src/ITensorMPS/dmrg.jl
index 80e8af50d8..8bd94f6511 100644
--- a/src/mps/dmrg.jl
+++ b/src/ITensorMPS/dmrg.jl
@@ -1,7 +1,4 @@
-
-#function IndexSet_ignore_missing(is::Union{Index,Nothing}...)
-#  return IndexSet(filter(i -> i isa Index, is))
-#end
+using KrylovKit: eigsolve
 
 function permute(
   M::AbstractMPS, ::Tuple{typeof(linkind),typeof(siteinds),typeof(linkind)}
diff --git a/src/ITensorMPS/exports.jl b/src/ITensorMPS/exports.jl
new file mode 100644
index 0000000000..e3a253dba5
--- /dev/null
+++ b/src/ITensorMPS/exports.jl
@@ -0,0 +1,128 @@
+export
+  # dmrg.jl
+  dmrg,
+
+  # abstractmps.jl
+  # Macros
+  @preserve_ortho,
+  # Methods
+  AbstractMPS,
+  add,
+  common_siteind,
+  common_siteinds,
+  findfirstsiteind,
+  findfirstsiteinds,
+  findsite,
+  findsites,
+  firstsiteind,
+  firstsiteinds,
+  logdot,
+  loginner,
+  lognorm,
+  movesite,
+  movesites,
+  normalize,
+  ortho_lims,
+  orthocenter,
+  promote_itensor_eltype,
+  reset_ortho_lims!,
+  set_ortho_lims!,
+  siteinds,
+  sim!,
+
+  # autompo/
+  AutoMPO,
+  add!,
+
+  # mpo.jl
+  # Types
+  MPO,
+  # Methods
+  error_contract,
+  maxlinkdim,
+  orthogonalize,
+  orthogonalize!,
+  outer,
+  projector,
+  randomMPO,
+  truncate,
+  truncate!,
+  unique_siteind,
+  unique_siteinds,
+
+  # mps.jl
+  # Types
+  MPS,
+  # Methods
+  ⋅,
+  correlation_matrix,
+  expect,
+  inner,
+  isortho,
+  linkdim,
+  linkdims,
+  linkind,
+  linkinds,
+  productMPS,
+  randomMPS,
+  replacebond,
+  replacebond!,
+  sample,
+  sample!,
+  siteind,
+  siteinds,
+  replace_siteinds!,
+  replace_siteinds,
+  swapbondsites,
+  totalqn,
+
+  # observer.jl
+  # Types
+  AbstractObserver,
+  DMRGObserver,
+  DMRGMeasurement,
+  NoObserver,
+  # Methods
+  checkdone!,
+  energies,
+  measure!,
+  measurements,
+  truncerrors,
+
+  # projmpo.jl
+  checkflux,
+  disk,
+  ProjMPO,
+  lproj,
+  product,
+  rproj,
+  noiseterm,
+  nsite,
+  position!,
+
+  # projmposum.jl
+  ProjMPOSum,
+
+  # projmpo_mps.jl
+  ProjMPO_MPS,
+
+  # sweeps.jl
+  Sweeps,
+  cutoff,
+  cutoff!,
+  get_cutoffs,
+  get_maxdims,
+  get_mindims,
+  get_noises,
+  maxdim,
+  maxdim!,
+  mindim,
+  mindim!,
+  noise,
+  noise!,
+  nsweep,
+  setmaxdim!,
+  setmindim!,
+  setcutoff!,
+  setnoise!,
+  sweepnext
diff --git a/src/ITensorMPS/imports.jl b/src/ITensorMPS/imports.jl
new file mode 100644
index 0000000000..937f431ba5
--- /dev/null
+++ b/src/ITensorMPS/imports.jl
@@ -0,0 +1,171 @@
+import Base:
+  # types
+  Array,
+  CartesianIndices,
+  Vector,
+  NTuple,
+  Tuple,
+  # symbols
+  +,
+  -,
+  *,
+  ^,
+  /,
+  ==,
+  <,
+  >,
+  !,
+  # functions
+  adjoint,
+  allunique,
+  axes,
+  complex,
+  conj,
+  convert,
+  copy,
+  copyto!,
+  deepcopy,
+  deleteat!,
+  eachindex,
+  eltype,
+  fill!,
+  filter,
+  filter!,
+  findall,
+  findfirst,
+  getindex,
+  hash,
+  imag,
+  intersect,
+  intersect!,
+  isapprox,
+  isassigned,
+  isempty,
+  isless,
+  isreal,
+  iszero,
+  iterate,
+  keys,
+  lastindex,
+  length,
+  map,
+  map!,
+  print,
+  promote_rule,
+  push!,
+  real,
+  resize!,
+  setdiff,
+  setdiff!,
+  setindex!,
+  show,
+  similar,
+  size,
+  summary,
+  truncate,
+  zero,
+  # macros
+  @propagate_inbounds
+
+import Base.Broadcast:
+  # types
+  AbstractArrayStyle,
+  Broadcasted,
+  BroadcastStyle,
+  DefaultArrayStyle,
+  Style,
+  # functions
+  _broadcast_getindex,
+  broadcasted,
+  broadcastable,
+  instantiate
+
+import Adapt: adapt_structure, adapt_storage
+
+import LinearAlgebra: dot, normalize!, tr
+
+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!,
+  scalartype,
+  single_precision,
+  tensor,
+  timer,
+  using_auto_fermion
+
+import ..ITensors:
+  AbstractRNG,
+  addtags,
+  Apply,
+  apply,
+  argument,
+  Broadcasted,
+  @Algorithm_str,
+  checkflux,
+  convert_leaf_eltype,
+  commontags,
+  @debug_check,
+  dag,
+  data,
+  DefaultArrayStyle,
+  DiskVector,
+  flux,
+  hascommoninds,
+  hasqns,
+  hassameinds,
+  HDF5,
+  inner,
+  isfermionic,
+  maxdim,
+  mindim,
+  ndims,
+  noprime,
+  noprime!,
+  norm,
+  normalize,
+  outer,
+  OneITensor,
+  permute,
+  prime,
+  prime!,
+  product,
+  QNIndex,
+  replaceinds,
+  replaceprime,
+  replacetags,
+  replacetags!,
+  setprime,
+  sim,
+  site,
+  siteind,
+  siteinds,
+  splitblocks,
+  store,
+  Style,
+  sum,
+  swapprime,
+  symmetrystyle,
+  terms,
+  @timeit_debug,
+  truncate!,
+  which_op
+
+import ..ITensors.Ops: params
+
+import SerializedElementArrays: disk
diff --git a/src/mps/mpo.jl b/src/ITensorMPS/mpo.jl
similarity index 97%
rename from src/mps/mpo.jl
rename to src/ITensorMPS/mpo.jl
index ac052dd7ad..9be5974cfe 100644
--- a/src/mps/mpo.jl
+++ b/src/ITensorMPS/mpo.jl
@@ -143,7 +143,7 @@ end
 
 function deprecate_make_inds_unmatch(::typeof(outer), ψ::MPS, ϕ::MPS; kw...)
   if hassameinds(siteinds, ψ, ϕ)
-    warn_once(outer_mps_mps_deprecation_warning(), :outer_mps_mps)
+    ITensors.warn_once(outer_mps_mps_deprecation_warning(), :outer_mps_mps)
     ψ = ψ'
   end
   return ψ, ϕ
@@ -353,7 +353,7 @@ function deprecate_make_inds_match!(
       )
     end
     if !hassameinds(siteinds, ydag, (A, x)) && make_inds_match
-      warn_once(inner_mps_mpo_mps_deprecation_warning(), :inner_mps_mpo_mps)
+      ITensors.warn_once(inner_mps_mpo_mps_deprecation_warning(), :inner_mps_mpo_mps)
       replace_siteinds!(ydag, sAx)
     end
   end
@@ -595,9 +595,11 @@ end
 
 (A::MPO)(ψ::MPS; kwargs...) = apply(A, ψ; kwargs...)
 
-Apply(A::MPO, ψ::MPS; kwargs...) = Applied(apply, (A, ψ), NamedTuple(kwargs))
+function Apply(A::MPO, ψ::MPS; kwargs...)
+  return ITensors.LazyApply.Applied(apply, (A, ψ), NamedTuple(kwargs))
+end
 
-function contract(A::MPO, ψ::MPS; alg=nothing, method=alg, kwargs...)
+function ITensors.contract(A::MPO, ψ::MPS; alg=nothing, method=alg, kwargs...)
   # TODO: Delete `method` since it is deprecated.
   alg = NDTensors.replace_nothing(method, "densitymatrix")
 
@@ -659,9 +661,9 @@ See also [`apply`](@ref).
 
 @doc """
 $contract_mpo_mps_doc
-""" contract(::MPO, ::MPS)
+""" ITensors.contract(::MPO, ::MPS)
 
-contract(ψ::MPS, A::MPO; kwargs...) = contract(A, ψ; kwargs...)
+ITensors.contract(ψ::MPS, A::MPO; kwargs...) = contract(A, ψ; kwargs...)
 
 *(A::MPO, B::MPS; kwargs...) = contract(A, B; kwargs...)
 *(A::MPS, B::MPO; kwargs...) = contract(A, B; kwargs...)
@@ -674,7 +676,7 @@ contract(ψ::MPS, A::MPO; kwargs...) = contract(A, ψ; kwargs...)
 
 #@doc (@doc contract(::MPO, ::MPS)) *(::MPO, ::MPS)
 
-function contract(
+function ITensors.contract(
   ::Algorithm"densitymatrix",
   A::MPO,
   ψ::MPS;
@@ -793,19 +795,19 @@ function _contract(::Algorithm"naive", A, ψ; truncate=true, kwargs...)
   return ψ_out
 end
 
-function contract(alg::Algorithm"naive", A::MPO, ψ::MPS; kwargs...)
+function ITensors.contract(alg::Algorithm"naive", A::MPO, ψ::MPS; kwargs...)
   return _contract(alg, A, ψ; kwargs...)
 end
 
-function contract(A::MPO, B::MPO; alg="zipup", kwargs...)
+function ITensors.contract(A::MPO, B::MPO; alg="zipup", kwargs...)
   return contract(Algorithm(alg), A, B; kwargs...)
 end
 
-function contract(alg::Algorithm"naive", A::MPO, B::MPO; kwargs...)
+function ITensors.contract(alg::Algorithm"naive", A::MPO, B::MPO; kwargs...)
   return _contract(alg, A, B; kwargs...)
 end
 
-function contract(
+function ITensors.contract(
   ::Algorithm"zipup",
   A::MPO,
   B::MPO;
@@ -945,7 +947,7 @@ See also [`apply`](@ref) for details about the arguments available.
 
 @doc """
 $contract_mpo_mpo_doc
-""" contract(::MPO, ::MPO)
+""" ITensors.contract(::MPO, ::MPO)
 
 *(A::MPO, B::MPO; kwargs...) = contract(A, B; kwargs...)
 
diff --git a/src/mps/mps.jl b/src/ITensorMPS/mps.jl
similarity index 100%
rename from src/mps/mps.jl
rename to src/ITensorMPS/mps.jl
diff --git a/src/mps/observer.jl b/src/ITensorMPS/observer.jl
similarity index 100%
rename from src/mps/observer.jl
rename to src/ITensorMPS/observer.jl
diff --git a/src/mps/sweeps.jl b/src/ITensorMPS/sweeps.jl
similarity index 100%
rename from src/mps/sweeps.jl
rename to src/ITensorMPS/sweeps.jl
diff --git a/src/ITensorMPS/usings.jl b/src/ITensorMPS/usings.jl
new file mode 100644
index 0000000000..697e8d6870
--- /dev/null
+++ b/src/ITensorMPS/usings.jl
@@ -0,0 +1,7 @@
+using Adapt
+using ..ITensors
+using HDF5: attributes, create_group, open_group
+using IsApprox
+using Printf
+using Random
+using TupleTools
diff --git a/src/ITensorNetworkMaps/ITensorNetworkMaps.jl b/src/ITensorNetworkMaps/ITensorNetworkMaps.jl
index 9864b11a41..bdeaf84f5c 100644
--- a/src/ITensorNetworkMaps/ITensorNetworkMaps.jl
+++ b/src/ITensorNetworkMaps/ITensorNetworkMaps.jl
@@ -3,7 +3,7 @@ module ITensorNetworkMaps
 using ..ITensors
 using LinearMaps
 
-using ITensors: promote_itensor_eltype
+import ITensors.ITensorMPS: promote_itensor_eltype
 
 import Base: *
 
diff --git a/src/ITensors.jl b/src/ITensors.jl
index c820798060..37cef8adae 100644
--- a/src/ITensors.jl
+++ b/src/ITensors.jl
@@ -61,13 +61,13 @@ using DocStringExtensions
 using Functors
 using HDF5
 using IsApprox
-using KrylovKit
 using LinearAlgebra
 using NDTensors
 using PackageCompiler
 using Pkg
 using Printf
 using Random
+using Reexport
 using SerializedElementArrays
 using StaticArrays
 using TimerOutputs
@@ -137,7 +137,6 @@ include("tensor_operations/matrix_algebra.jl")
 include("tensor_operations/permutations.jl")
 include("broadcast.jl")
 include("tensor_operations/matrix_decomposition.jl")
-include("iterativesolvers.jl")
 include("adapt.jl")
 include("set_types.jl")
 
@@ -162,30 +161,6 @@ include("nullspace.jl")
 #
 include("Ops/ops_itensor.jl")
 
-#####################################
-# MPS/MPO
-#
-include("mps/abstractmps.jl")
-include("mps/deprecated.jl")
-include("mps/mps.jl")
-include("mps/mpo.jl")
-include("mps/sweeps.jl")
-include("mps/abstractprojmpo/abstractprojmpo.jl")
-include("mps/abstractprojmpo/projmpo.jl")
-include("mps/abstractprojmpo/diskprojmpo.jl")
-include("mps/abstractprojmpo/projmposum.jl")
-include("mps/abstractprojmpo/projmps.jl")
-include("mps/abstractprojmpo/projmpo_mps.jl")
-include("mps/observer.jl")
-include("mps/dmrg.jl")
-include("mps/adapt.jl")
-
-#####################################
-# ITensorsNamedDimsArraysExt
-# Requires `AbstractMPS`.
-include("ITensorsNamedDimsArraysExt/src/ITensorsNamedDimsArraysExt.jl")
-using .ITensorsNamedDimsArraysExt: ITensorsNamedDimsArraysExt
-
 #####################################
 # Physics
 #
@@ -202,11 +177,18 @@ include("physics/site_types/tj.jl")
 include("physics/site_types/qudit.jl") # EXPERIMENTAL
 include("physics/site_types/boson.jl") # EXPERIMENTAL
 include("physics/fermions.jl")
-include("physics/autompo/matelem.jl")
-include("physics/autompo/qnmatelem.jl")
-include("physics/autompo/opsum_to_mpo_generic.jl")
-include("physics/autompo/opsum_to_mpo.jl")
-include("physics/autompo/opsum_to_mpo_qn.jl")
+
+#####################################
+# MPS/MPO
+#
+include("ITensorMPS/ITensorMPS.jl")
+@reexport using .ITensorMPS
+
+#####################################
+# ITensorsNamedDimsArraysExt
+# Requires `AbstractMPS`.
+include("ITensorsNamedDimsArraysExt/src/ITensorsNamedDimsArraysExt.jl")
+using .ITensorsNamedDimsArraysExt: ITensorsNamedDimsArraysExt
 
 #####################################
 # Trotter-Suzuki decomposition
diff --git a/src/ITensorsNamedDimsArraysExt/src/to_nameddimsarray.jl b/src/ITensorsNamedDimsArraysExt/src/to_nameddimsarray.jl
index bb943bd0b6..3644b90c7b 100644
--- a/src/ITensorsNamedDimsArraysExt/src/to_nameddimsarray.jl
+++ b/src/ITensorsNamedDimsArraysExt/src/to_nameddimsarray.jl
@@ -61,6 +61,6 @@ to_nameddimsarray(t::ITensor) = ITensor(to_nameddimsarray(t.tensor))
 using ITensors: ITensors
 # TODO: Delete once `TensorStorage` is removed.
 # TODO: Move to `ITensors`?
-function to_nameddimsarray(x::ITensors.AbstractMPS)
+function to_nameddimsarray(x::ITensors.ITensorMPS.AbstractMPS)
   return to_nameddimsarray.(x)
 end
diff --git a/src/deprecated.jl b/src/deprecated.jl
index 465b38cacb..7cb108c5dc 100644
--- a/src/deprecated.jl
+++ b/src/deprecated.jl
@@ -32,39 +32,6 @@
 @deprecate setstore(T::ITensor, st) setstorage(T, st) false
 @deprecate uniqueindex(args...; kwargs...) uniqueind(args...; kwargs...)
 
-# mps/abstractmps.jl
-@deprecate orthoCenter(args...; kwargs...) orthocenter(args...; kwargs...)
-@deprecate store(m::AbstractMPS) data(m)
-@deprecate replacesites!(args...; kwargs...) ITensors.replace_siteinds!(args...; kwargs...)
-@deprecate applyMPO(args...; kwargs...) contract(args...; kwargs...)
-@deprecate applympo(args...; kwargs...) contract(args...; kwargs...)
-@deprecate errorMPOprod(args...; kwargs...) error_contract(args...; kwargs...)
-@deprecate error_mpoprod(args...; kwargs...) error_contract(args...; kwargs...)
-@deprecate error_mul(args...; kwargs...) error_contract(args...; kwargs...)
-@deprecate multMPO(args...; kwargs...) contract(args...; kwargs...)
-@deprecate sum(A::AbstractMPS, B::AbstractMPS; kwargs...) add(A, B; kwargs...)
-@deprecate multmpo(args...; kwargs...) contract(args...; kwargs...)
-@deprecate set_leftlim!(args...; kwargs...) ITensors.setleftlim!(args...; kwargs...)
-@deprecate set_rightlim!(args...; kwargs...) ITensors.setrightlim!(args...; kwargs...)
-@deprecate tensors(args...; kwargs...) ITensors.data(args...; kwargs...)
-@deprecate primelinks!(args...; kwargs...) ITensors.prime_linkinds!(args...; kwargs...)
-@deprecate simlinks!(args...; kwargs...) ITensors.sim_linkinds!(args...; kwargs...)
-@deprecate mul(A::AbstractMPS, B::AbstractMPS; kwargs...) contract(A, B; kwargs...)
-
-# mps/mpo.jl
-@deprecate MPO(A::MPS; kwargs...) outer(A', A; kwargs...)
-
-# mps/mps.jl
-@deprecate randomMPS(sites::Vector{<:Index}, linkdims::Integer) randomMPS(
-  sites; linkdims=linkdims
-)
-@deprecate randomMPS(ElType::Type, sites::Vector{<:Index}, linkdims::Integer) randomMPS(
-  ElType, sites; linkdims=linkdims
-)
-@deprecate randomMPS(sites::Vector{<:Index}, state, linkdims::Integer) randomMPS(
-  sites, state; linkdims=linkdims
-)
-
 # physics/autompo.jl
 @deprecate toMPO(args...; kwargs...) MPO(args...; kwargs...)
 
diff --git a/src/exports.jl b/src/exports.jl
index 83af4065b6..659331b788 100644
--- a/src/exports.jl
+++ b/src/exports.jl
@@ -48,6 +48,9 @@ export
   svd,
   diag,
 
+  # tensor_operations/tensor_algebra.jl
+  contract,
+
   # global_variables.jl
   # Methods
 
@@ -199,9 +202,6 @@ export
   vector,
   emptyITensor,
 
-  # iterativesolvers.jl
-  davidson,
-
   # not.jl
   not,
 
@@ -223,126 +223,6 @@ export
   Prod,
   coefficient,
 
-  # mps/dmrg.jl
-  dmrg,
-
-  # mps/abstractmps.jl
-  # Macros
-  @preserve_ortho,
-  # Methods
-  add,
-  common_siteind,
-  common_siteinds,
-  contract,
-  convert_eltype,
-  dag!,
-  findfirstsiteind,
-  findfirstsiteinds,
-  findsite,
-  findsites,
-  firstsiteind,
-  firstsiteinds,
-  logdot,
-  loginner,
-  lognorm,
-  movesite,
-  movesites,
-  normalize,
-  ortho_lims,
-  reset_ortho_lims!,
-  set_ortho_lims!,
-  siteinds,
-
-  # mps/mpo.jl
-  # Types
-  MPO,
-  # Methods
-  error_contract,
-  maxlinkdim,
-  orthogonalize,
-  orthogonalize!,
-  outer,
-  projector,
-  randomMPO,
-  truncate,
-  truncate!,
-  unique_siteind,
-  unique_siteinds,
-
-  # mps/mps.jl
-  # Types
-  MPS,
-  # Methods
-  ⋅,
-  correlation_matrix,
-  expect,
-  inner,
-  isortho,
-  linkdim,
-  linkdims,
-  linkind,
-  linkinds,
-  productMPS,
-  randomMPS,
-  replacebond,
-  replacebond!,
-  sample,
-  sample!,
-  siteind,
-  siteinds,
-  replace_siteinds!,
-  replace_siteinds,
-  swapbondsites,
-  totalqn,
-
-  # mps/observer.jl
-  # Types
-  AbstractObserver,
-  DMRGObserver,
-  DMRGMeasurement,
-  NoObserver,
-  # Methods
-  checkdone!,
-  energies,
-  measure!,
-  measurements,
-  truncerrors,
-
-  # mps/projmpo.jl
-  ProjMPO,
-  lproj,
-  product,
-  rproj,
-  noiseterm,
-  position!,
-
-  # mps/projmposum.jl
-  ProjMPOSum,
-
-  # mps/projmpo_mps.jl
-  ProjMPO_MPS,
-
-  # mps/sweeps.jl
-  Sweeps,
-  cutoff,
-  cutoff!,
-  get_cutoffs,
-  get_maxdims,
-  get_mindims,
-  get_noises,
-  maxdim,
-  maxdim!,
-  mindim,
-  mindim!,
-  noise,
-  noise!,
-  nsweep,
-  setmaxdim!,
-  setmindim!,
-  setcutoff!,
-  setnoise!,
-  sweepnext,
-
   # Ops/Ops.jl
   Ops,
   Op,
diff --git a/src/iterativesolvers.jl b/src/iterativesolvers.jl
deleted file mode 100644
index 83f8e7661f..0000000000
--- a/src/iterativesolvers.jl
+++ /dev/null
@@ -1,119 +0,0 @@
-
-function get_vecs!((phi, q), M, V, AV, ni)
-  F = eigen(Hermitian(M))
-  lambda = F.values[1]
-  u = F.vectors[:, 1]
-  phi .= u[1] .* V[1]
-  q .= u[1] .* AV[1]
-  for n in 2:ni
-    phi .+= u[n] .* V[n]
-    q .+= u[n] .* AV[n]
-  end
-  q .-= lambda .* phi
-  #Fix sign
-  if real(u[1]) < 0.0
-    phi .*= -1
-    q .*= -1
-  end
-  return lambda
-end
-
-function orthogonalize!(q::ITensor, V, ni)
-  q0 = copy(q)
-  for k in 1:ni
-    Vq0k = dot(V[k], q0)
-    q .+= -Vq0k .* V[k]
-  end
-  qnrm = norm(q)
-  if qnrm < 1E-10 #orthog failure, try randomizing
-    randn!(q)
-    qnrm = norm(q)
-  end
-  q .*= 1.0 / qnrm
-  return nothing
-end
-
-function expand_krylov_space(M::Matrix{elT}, V, AV, ni) where {elT}
-  newM = fill(zero(elT), (ni + 1, ni + 1))
-  newM[1:ni, 1:ni] = M
-  for k in 1:(ni + 1)
-    newM[k, ni + 1] = dot(V[k], AV[ni + 1])
-    newM[ni + 1, k] = conj(newM[k, ni + 1])
-  end
-  return newM
-end
-
-function davidson(
-  A, phi0::ITensorT; maxiter=2, miniter=1, errgoal=1e-14, Northo_pass=1
-) where {ITensorT<:ITensor}
-  elTA = eltype(A)
-  elTphi = eltype(phi0)
-
-  # if the matrix is complex and the starting vector is real,
-  # that's not going to last long
-  # and we're going to need phi to be complex
-  # down when we reuse the storage in get_vecs!
-  if !(elTA <: Real) && (elTphi <: Real)
-    phi = complex(phi0)
-  else
-    phi = copy(phi0)
-  end
-  approx0 = 1E-12
-
-  nrm = norm(phi)
-  if nrm < 1E-18
-    phi_ = similar(phi)
-    randn!(phi_)
-    phi = phi_
-    nrm = norm(phi)
-  end
-  phi .*= 1.0 / nrm
-
-  maxsize = size(A)[1]
-  actual_maxiter = min(maxiter, maxsize - 1)
-
-  if dim(inds(phi)) != maxsize
-    error("linear size of A and dimension of phi should match in davidson")
-  end
-
-  Aphi = A(phi)
-
-  V = ITensorT[copy(phi)]
-  AV = ITensorT[Aphi]
-
-  last_lambda = NaN
-  lambda::Float64 = real(dot(V[1], AV[1]))
-  q = AV[1] - lambda * V[1]
-
-  M = fill(elTA(lambda), (1, 1))
-
-  for ni in 1:actual_maxiter
-    qnorm = norm(q)
-
-    errgoal_reached = (qnorm < errgoal && abs(lambda - last_lambda) < errgoal)
-    small_qnorm = (qnorm < max(approx0, errgoal * 1E-3))
-    converged = errgoal_reached || small_qnorm
-
-    if (qnorm < 1E-20) || (converged && ni > miniter) #|| (ni >= actual_maxiter)
-      #@printf "  done with davidson, ni=%d, qnorm=%.3E\n" ni qnorm
-      break
-    end
-
-    last_lambda = lambda
-
-    for pass in 1:Northo_pass
-      orthogonalize!(q, V, ni)
-    end
-
-    Aq = A(q)
-
-    push!(V, copy(q))
-    push!(AV, Aq)
-
-    M = expand_krylov_space(M, V, AV, ni)
-
-    lambda = get_vecs!((phi, q), M, V, AV, ni + 1)
-  end #for ni=1:actual_maxiter+1
-
-  return lambda, phi
-end
diff --git a/test/ITensorLegacyMPS/ITensorChainRules/Project.toml b/test/ITensorMPS/ITensorChainRules/Project.toml
similarity index 100%
rename from test/ITensorLegacyMPS/ITensorChainRules/Project.toml
rename to test/ITensorMPS/ITensorChainRules/Project.toml
diff --git a/test/ITensorLegacyMPS/ITensorChainRules/runtests.jl b/test/ITensorMPS/ITensorChainRules/runtests.jl
similarity index 100%
rename from test/ITensorLegacyMPS/ITensorChainRules/runtests.jl
rename to test/ITensorMPS/ITensorChainRules/runtests.jl
diff --git a/test/ITensorLegacyMPS/ITensorChainRules/test_chainrules.jl b/test/ITensorMPS/ITensorChainRules/test_chainrules.jl
similarity index 100%
rename from test/ITensorLegacyMPS/ITensorChainRules/test_chainrules.jl
rename to test/ITensorMPS/ITensorChainRules/test_chainrules.jl
diff --git a/test/ITensorLegacyMPS/ITensorChainRules/test_optimization.jl b/test/ITensorMPS/ITensorChainRules/test_optimization.jl
similarity index 100%
rename from test/ITensorLegacyMPS/ITensorChainRules/test_optimization.jl
rename to test/ITensorMPS/ITensorChainRules/test_optimization.jl
diff --git a/test/ITensorLegacyMPS/ITensorChainRules/utils/circuit.jl b/test/ITensorMPS/ITensorChainRules/utils/circuit.jl
similarity index 100%
rename from test/ITensorLegacyMPS/ITensorChainRules/utils/circuit.jl
rename to test/ITensorMPS/ITensorChainRules/utils/circuit.jl
diff --git a/test/ITensorLegacyMPS/Ops/runtests.jl b/test/ITensorMPS/Ops/runtests.jl
similarity index 100%
rename from test/ITensorLegacyMPS/Ops/runtests.jl
rename to test/ITensorMPS/Ops/runtests.jl
diff --git a/test/ITensorLegacyMPS/Ops/test_ops_mpo.jl b/test/ITensorMPS/Ops/test_ops_mpo.jl
similarity index 100%
rename from test/ITensorLegacyMPS/Ops/test_ops_mpo.jl
rename to test/ITensorMPS/Ops/test_ops_mpo.jl
diff --git a/test/ITensorLegacyMPS/Ops/test_trotter.jl b/test/ITensorMPS/Ops/test_trotter.jl
similarity index 100%
rename from test/ITensorLegacyMPS/Ops/test_trotter.jl
rename to test/ITensorMPS/Ops/test_trotter.jl
diff --git a/test/ITensorLegacyMPS/base/Project.toml b/test/ITensorMPS/base/Project.toml
similarity index 100%
rename from test/ITensorLegacyMPS/base/Project.toml
rename to test/ITensorMPS/base/Project.toml
diff --git a/test/ITensorLegacyMPS/base/backup/test_arraystorage.jl b/test/ITensorMPS/base/backup/test_arraystorage.jl
similarity index 100%
rename from test/ITensorLegacyMPS/base/backup/test_arraystorage.jl
rename to test/ITensorMPS/base/backup/test_arraystorage.jl
diff --git a/test/ITensorLegacyMPS/base/runtests.jl b/test/ITensorMPS/base/runtests.jl
similarity index 100%
rename from test/ITensorLegacyMPS/base/runtests.jl
rename to test/ITensorMPS/base/runtests.jl
diff --git a/test/ITensorLegacyMPS/base/test_abstractprojmpo.jl b/test/ITensorMPS/base/test_abstractprojmpo.jl
similarity index 71%
rename from test/ITensorLegacyMPS/base/test_abstractprojmpo.jl
rename to test/ITensorMPS/base/test_abstractprojmpo.jl
index 4c18f13195..993aee5d81 100644
--- a/test/ITensorLegacyMPS/base/test_abstractprojmpo.jl
+++ b/test/ITensorMPS/base/test_abstractprojmpo.jl
@@ -1,6 +1,7 @@
 using ITensors
 using Random
 using Test
+using ITensors: ITensorMPS
 
 @testset "AbstractProjMPO (eltype=$elt, conserve_qns=$conserve_qns)" for elt in (
     Float32, Float64, Complex{Float32}, Complex{Float64}
@@ -17,10 +18,10 @@ using Test
     po = storage(pmpo)
 
     # `AbstractProjMPO` interface.
-    @test ITensors.nsite(po) == 2
-    @test ITensors.site_range(po) == 2:3
+    @test ITensorMPS.nsite(po) == 2
+    @test ITensorMPS.site_range(po) == 2:3
     @test eltype(po) == elt
-    @test isnothing(ITensors.checkflux(po))
+    @test isnothing(ITensorMPS.checkflux(po))
     po_contracted = contract(po, ITensor(one(Bool)))
     @test po_contracted isa ITensor
     @test ndims(po_contracted) == 8
@@ -38,42 +39,42 @@ using Test
     po = storage(ProjMPOSum([pmpo, pmpo]))
 
     # `AbstractProjMPO` interface.
-    @test ITensors.nsite(po) == 2
-    @test ITensors.site_range(po) == 2:3
+    @test ITensorMPS.nsite(po) == 2
+    @test ITensorMPS.site_range(po) == 2:3
     @test eltype(po) == elt
-    @test isnothing(ITensors.checkflux(po))
+    @test isnothing(ITensorMPS.checkflux(po))
     po_contracted = contract(po, ITensor(one(Bool)))
     @test po_contracted isa ITensor
     @test ndims(po_contracted) == 8
     @test eltype(po_contracted) == elt
 
     # Specific to `ProjMPOSum`.
-    @test length(ITensors.terms(po)) == 2
+    @test length(ITensorMPS.terms(po)) == 2
   end
-  @testset "ITensors.ProjMPS" begin
+  @testset "ITensorMPS.ProjMPS" begin
     # TODO: Replace with `ProjOuter`, make it into
     # a proper `AbstractProjMPO`.
-    px = ITensors.ProjMPS(x)
+    px = ITensorMPS.ProjMPS(x)
     position!(px, x, 2)
 
     # `AbstractProjMPO` interface.
-    @test ITensors.nsite(px) == 2
-    @test ITensors.site_range(px) == 2:3
+    @test ITensorMPS.nsite(px) == 2
+    @test ITensorMPS.site_range(px) == 2:3
     @test_broken eltype(px) == elt
-    @test isnothing(ITensors.checkflux(px))
+    @test isnothing(ITensorMPS.checkflux(px))
     @test_broken contract(px, ITensor(one(Bool)))
   end
-  @testset "ITensors.ProjMPO_MPS" begin
+  @testset "ITensorMPS.ProjMPO_MPS" begin
     # TODO: Replace with `ProjOuter`, make it into
     # a proper `AbstractProjMPO`.
-    po = ITensors.ProjMPO_MPS(o, [x])
+    po = ITensorMPS.ProjMPO_MPS(o, [x])
     position!(po, x, 2)
 
     # `AbstractProjMPO` interface.
-    @test ITensors.nsite(po) == 2
-    @test ITensors.site_range(po) == 2:3
+    @test ITensorMPS.nsite(po) == 2
+    @test ITensorMPS.site_range(po) == 2:3
     @test_broken eltype(po) == elt
-    @test isnothing(ITensors.checkflux(po))
+    @test isnothing(ITensorMPS.checkflux(po))
     @test_broken contract(po, ITensor(one(Bool)))
   end
 end
diff --git a/test/ITensorLegacyMPS/base/test_algorithm.jl b/test/ITensorMPS/base/test_algorithm.jl
similarity index 100%
rename from test/ITensorLegacyMPS/base/test_algorithm.jl
rename to test/ITensorMPS/base/test_algorithm.jl
diff --git a/test/ITensorLegacyMPS/base/test_autompo.jl b/test/ITensorMPS/base/test_autompo.jl
similarity index 99%
rename from test/ITensorLegacyMPS/base/test_autompo.jl
rename to test/ITensorMPS/base/test_autompo.jl
index d44c80e688..b44ef903ef 100644
--- a/test/ITensorLegacyMPS/base/test_autompo.jl
+++ b/test/ITensorMPS/base/test_autompo.jl
@@ -1248,5 +1248,3 @@ end
     @test_nowarn H = MPO(os, sites)
   end
 end
-
-nothing
diff --git a/test/ITensorLegacyMPS/base/test_dmrg.jl b/test/ITensorMPS/base/test_dmrg.jl
similarity index 99%
rename from test/ITensorLegacyMPS/base/test_dmrg.jl
rename to test/ITensorMPS/base/test_dmrg.jl
index 4854587cf9..a13a67e93b 100644
--- a/test/ITensorLegacyMPS/base/test_dmrg.jl
+++ b/test/ITensorMPS/base/test_dmrg.jl
@@ -1,6 +1,5 @@
 using ITensors, Test, Random
-
-using ITensors: nsite, set_nsite!
+using ITensors.ITensorMPS: nsite, set_nsite!, site_range
 
 @testset "Basic DMRG" begin
   @testset "Spin-one Heisenberg" begin
@@ -108,7 +107,7 @@ using ITensors: nsite, set_nsite!
 
     @test length(PH) == N
     @test length(PHdisk) == N
-    @test ITensors.site_range(PH) == n:(n + 1)
+    @test site_range(PH) == n:(n + 1)
     @test eltype(PH) == Float64
     ## TODO sometimes randomMPS gives a linkdim value of 3 
     ## which causes an error in `calculated_dim = 3^2 * 4^2`
@@ -467,5 +466,3 @@ using ITensors: nsite, set_nsite!
     energy, psi = dmrg(H, psi0, sweeps; outputlevel=0)
   end
 end
-
-nothing
diff --git a/test/ITensorLegacyMPS/base/test_examples.jl b/test/ITensorMPS/base/test_examples.jl
similarity index 98%
rename from test/ITensorLegacyMPS/base/test_examples.jl
rename to test/ITensorMPS/base/test_examples.jl
index d1a6292d54..6c14e7a650 100644
--- a/test/ITensorLegacyMPS/base/test_examples.jl
+++ b/test/ITensorMPS/base/test_examples.jl
@@ -21,5 +21,3 @@ using ITensors, Test, Suppressor
     end
   end
 end
-
-nothing
diff --git a/test/ITensorLegacyMPS/base/test_fermions.jl b/test/ITensorMPS/base/test_fermions.jl
similarity index 99%
rename from test/ITensorLegacyMPS/base/test_fermions.jl
rename to test/ITensorMPS/base/test_fermions.jl
index 6bfbcb056f..da31384d26 100644
--- a/test/ITensorLegacyMPS/base/test_fermions.jl
+++ b/test/ITensorMPS/base/test_fermions.jl
@@ -285,5 +285,3 @@ import ITensors: Out, In
 
   ITensors.disable_auto_fermion()
 end
-
-nothing
diff --git a/test/ITensorLegacyMPS/base/test_inference.jl b/test/ITensorMPS/base/test_inference.jl
similarity index 100%
rename from test/ITensorLegacyMPS/base/test_inference.jl
rename to test/ITensorMPS/base/test_inference.jl
diff --git a/test/ITensorLegacyMPS/base/test_mpo.jl b/test/ITensorMPS/base/test_mpo.jl
similarity index 99%
rename from test/ITensorLegacyMPS/base/test_mpo.jl
rename to test/ITensorMPS/base/test_mpo.jl
index d426f0d27b..cbe96e205b 100644
--- a/test/ITensorLegacyMPS/base/test_mpo.jl
+++ b/test/ITensorMPS/base/test_mpo.jl
@@ -852,5 +852,3 @@ end
     @test maxlinkdim(H) ≤ maxlinkdim(H₁) + maxlinkdim(H₂)
   end
 end
-
-nothing
diff --git a/test/ITensorLegacyMPS/base/test_mps.jl b/test/ITensorMPS/base/test_mps.jl
similarity index 95%
rename from test/ITensorLegacyMPS/base/test_mps.jl
rename to test/ITensorMPS/base/test_mps.jl
index 251916d345..d44e0edf1e 100644
--- a/test/ITensorLegacyMPS/base/test_mps.jl
+++ b/test/ITensorMPS/base/test_mps.jl
@@ -1,5 +1,6 @@
 using Combinatorics
 using ITensors
+using ITensors: ITensorMPS
 using Random
 using LinearAlgebra: diag
 using Test
@@ -160,7 +161,7 @@ include(joinpath(@__DIR__, "utils", "util.jl"))
           @test eltype(psi[j]) == ComplexF64
         end
         @test eltype(psi) == ITensor
-        @test ITensors.promote_itensor_eltype(psi) == ComplexF64
+        @test ITensorMPS.promote_itensor_eltype(psi) == ComplexF64
       end
     end
 
@@ -585,17 +586,17 @@ include(joinpath(@__DIR__, "utils", "util.jl"))
   for (l0, l) in zip(l0s, ls)
     @test tags(l0) == tags(l)
   end
-  @test ITensors.leftlim(psi) == length(psi) - 2
-  @test ITensors.rightlim(psi) == length(psi)
+  @test ITensorMPS.leftlim(psi) == length(psi) - 2
+  @test ITensorMPS.rightlim(psi) == length(psi)
   orthogonalize!(psi, 2)
-  @test ITensors.leftlim(psi) == 1
-  @test ITensors.rightlim(psi) == 3
+  @test ITensorMPS.leftlim(psi) == 1
+  @test ITensorMPS.rightlim(psi) == 3
   psi = randomMPS(sites)
-  ITensors.setrightlim!(psi, length(psi) + 1) # do this to test qr
+  ITensorMPS.setrightlim!(psi, length(psi) + 1) # do this to test qr
   # from rightmost tensor
   orthogonalize!(psi, div(length(psi), 2))
-  @test ITensors.leftlim(psi) == div(length(psi), 2) - 1
-  @test ITensors.rightlim(psi) == div(length(psi), 2) + 1
+  @test ITensorMPS.leftlim(psi) == div(length(psi), 2) - 1
+  @test ITensorMPS.rightlim(psi) == div(length(psi), 2) + 1
 
   @test isnothing(linkind(MPS(fill(ITensor(), length(psi)), 0, length(psi) + 1), 1))
 
@@ -611,20 +612,20 @@ include(joinpath(@__DIR__, "utils", "util.jl"))
     orthogonalize!(psi, 5)
     phi = psi[5] * psi[6]
     replacebond!(psi, 5, phi; ortho="left")
-    @test ITensors.leftlim(psi) == 5
-    @test ITensors.rightlim(psi) == 7
+    @test ITensorMPS.leftlim(psi) == 5
+    @test ITensorMPS.rightlim(psi) == 7
 
     phi = psi[5] * psi[6]
     replacebond!(psi, 5, phi; ortho="right")
-    @test ITensors.leftlim(psi) == 4
-    @test ITensors.rightlim(psi) == 6
+    @test ITensorMPS.leftlim(psi) == 4
+    @test ITensorMPS.rightlim(psi) == 6
 
-    ITensors.setleftlim!(psi, 3)
-    ITensors.setrightlim!(psi, 7)
+    ITensorMPS.setleftlim!(psi, 3)
+    ITensorMPS.setrightlim!(psi, 7)
     phi = psi[5] * psi[6]
     replacebond!(psi, 5, phi; ortho="left")
-    @test ITensors.leftlim(psi) == 3
-    @test ITensors.rightlim(psi) == 7
+    @test ITensorMPS.leftlim(psi) == 3
+    @test ITensorMPS.rightlim(psi) == 7
 
     # check that replacebond! runs with svd kwargs
     psi = randomMPS(sites)
@@ -648,8 +649,8 @@ end
   init_state = [isodd(n) ? "Up" : "Dn" for n in 1:length(sites)]
   psi0 = MPS(sites, init_state)
   orthogonalize!(psi0, 4)
-  @test ITensors.leftlim(psi0) == 3
-  @test ITensors.rightlim(psi0) == 5
+  @test ITensorMPS.leftlim(psi0) == 3
+  @test ITensorMPS.rightlim(psi0) == 5
 end
 
 # Helper function for making MPS
@@ -691,8 +692,8 @@ end
     M = basicRandomMPS(30)
     orthogonalize!(M, c)
 
-    @test ITensors.leftlim(M) == c - 1
-    @test ITensors.rightlim(M) == c + 1
+    @test ITensorMPS.leftlim(M) == c - 1
+    @test ITensorMPS.rightlim(M) == c + 1
 
     # Test for left-orthogonality
     L = M[1] * prime(M[1], "Link")
@@ -722,7 +723,7 @@ end
     M0 = copy(M)
     truncate!(M; maxdim=5)
 
-    @test ITensors.rightlim(M) == 2
+    @test ITensorMPS.rightlim(M) == 2
 
     # Test for right-orthogonality
     R = M[length(M)] * prime(M[length(M)], "Link")
@@ -779,8 +780,8 @@ end
     sites = siteinds(2, 20)
     M = randomMPS(sites; linkdims=chi)
 
-    @test ITensors.leftlim(M) == 0
-    @test ITensors.rightlim(M) == 2
+    @test ITensorMPS.leftlim(M) == 0
+    @test ITensorMPS.rightlim(M) == 2
 
     @test norm(M[1]) ≈ 1.0
 
@@ -810,8 +811,8 @@ end
     M = randomMPS(sites, state; linkdims=chi)
     @test flux(M) == QN("Sz", 0)
 
-    @test ITensors.leftlim(M) == 0
-    @test ITensors.rightlim(M) == 2
+    @test ITensorMPS.leftlim(M) == 0
+    @test ITensorMPS.rightlim(M) == 2
 
     @test norm(M[1]) ≈ 1.0
     @test inner(M, M) ≈ 1.0
@@ -1100,31 +1101,31 @@ end
 
     # Test map! with limits getting set
     M = orthogonalize(M0, 1)
-    @test ITensors.leftlim(M) == 0
-    @test ITensors.rightlim(M) == 2
+    @test ITensorMPS.leftlim(M) == 0
+    @test ITensorMPS.rightlim(M) == 2
     map!(prime, M)
-    @test ITensors.leftlim(M) == 0
-    @test ITensors.rightlim(M) == length(M0) + 1
+    @test ITensorMPS.leftlim(M) == 0
+    @test ITensorMPS.rightlim(M) == length(M0) + 1
 
     # Test map! without limits getting set
     M = orthogonalize(M0, 1)
     map!(prime, M; set_limits=false)
-    @test ITensors.leftlim(M) == 0
-    @test ITensors.rightlim(M) == 2
+    @test ITensorMPS.leftlim(M) == 0
+    @test ITensorMPS.rightlim(M) == 2
 
     # Test prime! with limits getting set
     M = orthogonalize(M0, 1)
-    @test ITensors.leftlim(M) == 0
-    @test ITensors.rightlim(M) == 2
+    @test ITensorMPS.leftlim(M) == 0
+    @test ITensorMPS.rightlim(M) == 2
     prime!(M; set_limits=true)
-    @test ITensors.leftlim(M) == 0
-    @test ITensors.rightlim(M) == length(M0) + 1
+    @test ITensorMPS.leftlim(M) == 0
+    @test ITensorMPS.rightlim(M) == length(M0) + 1
 
     # Test prime! without limits getting set
     M = orthogonalize(M0, 1)
     prime!(M)
-    @test ITensors.leftlim(M) == 0
-    @test ITensors.rightlim(M) == 2
+    @test ITensorMPS.leftlim(M) == 0
+    @test ITensorMPS.rightlim(M) == 2
   end
 
   @testset "setindex!(::MPS, _, ::Colon)" begin
@@ -1133,15 +1134,15 @@ end
     ϕ = MPS(s, "↑")
     orthogonalize!(ϕ, 1)
     ψ[:] = ϕ
-    @test ITensors.orthocenter(ψ) == 1
+    @test ITensorMPS.orthocenter(ψ) == 1
     @test inner(ψ, ϕ) ≈ 1
 
     ψ = randomMPS(s)
     ϕ = MPS(s, "↑")
     orthogonalize!(ϕ, 1)
-    ψ[:] = ITensors.data(ϕ)
-    @test ITensors.leftlim(ψ) == 0
-    @test ITensors.rightlim(ψ) == length(ψ) + 1
+    ψ[:] = ITensorMPS.data(ϕ)
+    @test ITensorMPS.leftlim(ψ) == 0
+    @test ITensorMPS.rightlim(ψ) == length(ψ) + 1
     @test inner(ψ, ϕ) ≈ 1
   end
 
@@ -1237,21 +1238,21 @@ end
     A = randomITensor(s...)
     ψ = MPS(A, s)
     @test prod(ψ) ≈ A
-    @test ITensors.orthocenter(ψ) == N
+    @test ITensorMPS.orthocenter(ψ) == N
     @test maxlinkdim(ψ) == 4
 
     ψ0 = MPS(s, "↑")
     A = prod(ψ0)
     ψ = MPS(A, s; cutoff=1e-15)
     @test prod(ψ) ≈ A
-    @test ITensors.orthocenter(ψ) == N
+    @test ITensorMPS.orthocenter(ψ) == N
     @test maxlinkdim(ψ) == 1
 
     ψ0 = randomMPS(s; linkdims=2)
     A = prod(ψ0)
     ψ = MPS(A, s; cutoff=1e-15, orthocenter=2)
     @test prod(ψ) ≈ A
-    @test ITensors.orthocenter(ψ) == 2
+    @test ITensorMPS.orthocenter(ψ) == 2
     @test maxlinkdim(ψ) == 2
 
     A = randomITensor(s..., l[1], r[1])
@@ -1260,7 +1261,7 @@ end
     @test hassameinds(ψ[1], (l[1], s[1], ls[1]))
     @test hassameinds(ψ[N], (r[1], s[N], ls[N - 1]))
     @test prod(ψ) ≈ A
-    @test ITensors.orthocenter(ψ) == 3
+    @test ITensorMPS.orthocenter(ψ) == 3
     @test maxlinkdim(ψ) == 12
 
     A = randomITensor(s..., l..., r...)
@@ -1269,7 +1270,7 @@ end
     @test hassameinds(ψ[1], (l..., s[1], ls[1]))
     @test hassameinds(ψ[N], (r..., s[N], ls[N - 1]))
     @test prod(ψ) ≈ A
-    @test ITensors.orthocenter(ψ) == N
+    @test ITensorMPS.orthocenter(ψ) == N
     @test maxlinkdim(ψ) == 36
 
     #
@@ -1290,26 +1291,26 @@ end
     ψ0 = randomMPS(s; linkdims=3)
 
     ψ = orthogonalize(ψ0, 2)
-    A = prod(ITensors.data(ψ)[2:(N - 1)])
+    A = prod(ITensorMPS.data(ψ)[2:(N - 1)])
     randn!(A)
     ϕ = MPS(A, s[2:(N - 1)]; orthocenter=1)
     ψ[2:(N - 1)] = ϕ
     @test prod(ψ) ≈ ψ[1] * A * ψ[N]
     @test maxlinkdim(ψ) == 4
-    @test ITensors.orthocenter(ψ) == 2
+    @test ITensorMPS.orthocenter(ψ) == 2
 
     ψ = orthogonalize(ψ0, 1)
-    A = prod(ITensors.data(ψ)[2:(N - 1)])
+    A = prod(ITensorMPS.data(ψ)[2:(N - 1)])
     randn!(A)
     @test_throws AssertionError ψ[2:(N - 1)] = A
 
     ψ = orthogonalize(ψ0, 2)
-    A = prod(ITensors.data(ψ)[2:(N - 1)])
+    A = prod(ITensorMPS.data(ψ)[2:(N - 1)])
     randn!(A)
     ψ[2:(N - 1), orthocenter=3] = A
     @test prod(ψ) ≈ ψ[1] * A * ψ[N]
     @test maxlinkdim(ψ) == 4
-    @test ITensors.orthocenter(ψ) == 3
+    @test ITensorMPS.orthocenter(ψ) == 3
   end
 
   @testset "movesites reverse sites" begin
@@ -1967,9 +1968,9 @@ end
     ψ = MPS([itensor(randn(ComplexF64, 2), s[n]) for n in 1:N])
     @test all(==(IndexSet()), linkinds(all, ψ))
     ϕ = orthogonalize(ψ, 2)
-    @test ITensors.hasdefaultlinktags(ϕ)
+    @test ITensorMPS.hasdefaultlinktags(ϕ)
     @test ortho_lims(ϕ) == 2:2
-    @test ITensors.dist(ψ, ϕ) ≈ 0 atol = 1e-6
+    @test ITensorMPS.dist(ψ, ϕ) ≈ 0 atol = 1e-6
     # TODO: use this instead?
     # @test lognorm(ψ - ϕ) < -16
     @test norm(ψ - ϕ) ≈ 0 atol = 1e-6
@@ -1980,7 +1981,7 @@ end
     s = siteinds("S=1/2", N)
     ψ = MPS([itensor(randn(ComplexF64, 2), s[n]) for n in 1:N])
     ρ = outer(ψ', ψ)
-    @test !ITensors.hasnolinkinds(ρ)
+    @test !ITensorMPS.hasnolinkinds(ρ)
     @test inner(ρ, ρ) ≈ inner(ψ, ψ)^2
     @test inner(ψ', ρ, ψ) ≈ inner(ψ, ψ)^2
 
@@ -1990,7 +1991,7 @@ end
     @test_deprecated inner(ψ, ρ, ψ)
 
     ρ = @test_deprecated MPO(ψ)
-    @test !ITensors.hasnolinkinds(ρ)
+    @test !ITensorMPS.hasnolinkinds(ρ)
     @test inner(ρ, ρ) ≈ inner(ψ, ψ)^2
     @test inner(ψ', ρ, ψ) ≈ inner(ψ, ψ)^2
   end
@@ -1999,11 +2000,9 @@ end
     N = 4
     s = siteinds("S=1/2", N)
     M = MPO([itensor(randn(ComplexF64, 2, 2), s[n]', dag(s[n])) for n in 1:N])
-    @test ITensors.hasnolinkinds(M)
+    @test ITensorMPS.hasnolinkinds(M)
     Mt = truncate(M; cutoff=1e-15)
-    @test ITensors.hasdefaultlinktags(Mt)
+    @test ITensorMPS.hasdefaultlinktags(Mt)
     @test norm(M - Mt) ≈ 0 atol = 1e-12
   end
 end
-
-nothing
diff --git a/test/ITensorLegacyMPS/base/test_qnmpo.jl b/test/ITensorMPS/base/test_qnmpo.jl
similarity index 99%
rename from test/ITensorLegacyMPS/base/test_qnmpo.jl
rename to test/ITensorMPS/base/test_qnmpo.jl
index 9f3a487f0b..f74710d71c 100644
--- a/test/ITensorLegacyMPS/base/test_qnmpo.jl
+++ b/test/ITensorMPS/base/test_qnmpo.jl
@@ -373,4 +373,3 @@ test_combos = [(make_Heisenberg_AutoMPO, "S=1/2"), (make_Hubbard_AutoMPO, "Elect
     @test norm(dense(Q * dag(prime(Q, q))) - δ(Float64, q, q')) ≈ 0.0 atol = 1e-13
   end
 end
-nothing
diff --git a/test/ITensorLegacyMPS/base/test_readme.jl b/test/ITensorMPS/base/test_readme.jl
similarity index 99%
rename from test/ITensorLegacyMPS/base/test_readme.jl
rename to test/ITensorMPS/base/test_readme.jl
index 95556f1c23..79688626e1 100644
--- a/test/ITensorLegacyMPS/base/test_readme.jl
+++ b/test/ITensorMPS/base/test_readme.jl
@@ -97,5 +97,3 @@ using ITensors, Test
     #println("Final energy = $energy")
   end
 end
-
-nothing
diff --git a/test/ITensorLegacyMPS/base/test_readwrite.jl b/test/ITensorMPS/base/test_readwrite.jl
similarity index 98%
rename from test/ITensorLegacyMPS/base/test_readwrite.jl
rename to test/ITensorMPS/base/test_readwrite.jl
index c1694bcb9f..fe3bb64f61 100644
--- a/test/ITensorLegacyMPS/base/test_readwrite.jl
+++ b/test/ITensorMPS/base/test_readwrite.jl
@@ -36,5 +36,3 @@ include(joinpath(@__DIR__, "utils", "util.jl"))
   #
   rm("data.h5"; force=true)
 end
-
-nothing
diff --git a/test/ITensorLegacyMPS/base/test_symmetrystyle.jl b/test/ITensorMPS/base/test_symmetrystyle.jl
similarity index 100%
rename from test/ITensorLegacyMPS/base/test_symmetrystyle.jl
rename to test/ITensorMPS/base/test_symmetrystyle.jl
diff --git a/test/ITensorLegacyMPS/base/test_threading.jl b/test/ITensorMPS/base/test_threading.jl
similarity index 100%
rename from test/ITensorLegacyMPS/base/test_threading.jl
rename to test/ITensorMPS/base/test_threading.jl
diff --git a/test/ITensorLegacyMPS/base/utils/opsum_hash_bug.jld2 b/test/ITensorMPS/base/utils/opsum_hash_bug.jld2
similarity index 100%
rename from test/ITensorLegacyMPS/base/utils/opsum_hash_bug.jld2
rename to test/ITensorMPS/base/utils/opsum_hash_bug.jld2
diff --git a/test/ITensorLegacyMPS/base/utils/util.jl b/test/ITensorMPS/base/utils/util.jl
similarity index 100%
rename from test/ITensorLegacyMPS/base/utils/util.jl
rename to test/ITensorMPS/base/utils/util.jl
diff --git a/test/ITensorLegacyMPS/runtests.jl b/test/ITensorMPS/runtests.jl
similarity index 100%
rename from test/ITensorLegacyMPS/runtests.jl
rename to test/ITensorMPS/runtests.jl
diff --git a/test/base/test_combiner.jl b/test/base/test_combiner.jl
index c805e29bb1..f668b18be3 100644
--- a/test/base/test_combiner.jl
+++ b/test/base/test_combiner.jl
@@ -261,5 +261,3 @@ using Combinatorics: permutations
     @test B == A
   end
 end
-
-nothing
diff --git a/test/base/test_contract.jl b/test/base/test_contract.jl
index a0f5ec1df9..deeebc0327 100644
--- a/test/base/test_contract.jl
+++ b/test/base/test_contract.jl
@@ -312,5 +312,3 @@ end
     @test array(permute(C, i, j)) ≈ kron(array(A), transpose(array(B)))
   end
 end
-
-nothing
diff --git a/test/base/test_ctmrg.jl b/test/base/test_ctmrg.jl
index ab97498c9e..8d7819a329 100644
--- a/test/base/test_ctmrg.jl
+++ b/test/base/test_ctmrg.jl
@@ -50,5 +50,3 @@ include(joinpath(src_dir, "2d_classical_ising.jl"))
   m = (ACTszₗ * dag(ACₗ))[] / κ
   @test abs(m) ≈ ising_magnetization(β)
 end
-
-nothing
diff --git a/test/base/test_decomp.jl b/test/base/test_decomp.jl
index 3082cc9202..26a21c37f2 100644
--- a/test/base/test_decomp.jl
+++ b/test/base/test_decomp.jl
@@ -500,5 +500,3 @@ end
     @test dim(commonind(U, B)) <= 4
   end
 end
-
-nothing
diff --git a/test/base/test_diagitensor.jl b/test/base/test_diagitensor.jl
index db11eb886a..41e13ccda4 100644
--- a/test/base/test_diagitensor.jl
+++ b/test/base/test_diagitensor.jl
@@ -544,5 +544,3 @@ using Test
     end
   end
 end
-
-nothing
diff --git a/test/base/test_empty.jl b/test/base/test_empty.jl
index abb6b7d823..6ccb0d52d2 100644
--- a/test/base/test_empty.jl
+++ b/test/base/test_empty.jl
@@ -79,5 +79,3 @@ using Test
     @test norm(E - A) < 1E-8
   end
 end
-
-nothing
diff --git a/test/base/test_emptyitensor.jl b/test/base/test_emptyitensor.jl
index 855cda73ae..634b2ee6b6 100644
--- a/test/base/test_emptyitensor.jl
+++ b/test/base/test_emptyitensor.jl
@@ -101,5 +101,3 @@ end
   A = ITensor(i', dag(i))
   @test NDTensors.tensor(zero(A)) isa typeof(NDTensors.tensor(A))
 end
-
-nothing
diff --git a/test/base/test_examples.jl b/test/base/test_examples.jl
index 8d4e44b819..fa381615d3 100644
--- a/test/base/test_examples.jl
+++ b/test/base/test_examples.jl
@@ -10,5 +10,3 @@ using ITensors, Test, Suppressor
     end
   end
 end
-
-nothing
diff --git a/test/base/test_fermions.jl b/test/base/test_fermions.jl
index c3e8edfa74..fae73fbabc 100644
--- a/test/base/test_fermions.jl
+++ b/test/base/test_fermions.jl
@@ -770,5 +770,3 @@ import ITensors: Out, In
 
   ITensors.disable_auto_fermion()
 end
-
-nothing
diff --git a/test/base/test_index.jl b/test/base/test_index.jl
index a72781aa13..03a1264da0 100644
--- a/test/base/test_index.jl
+++ b/test/base/test_index.jl
@@ -183,5 +183,3 @@ import ITensors: In, Out, Neither
     @test hastags(ijk, "test2")
   end
 end
-
-nothing
diff --git a/test/base/test_indexset.jl b/test/base/test_indexset.jl
index c83952e2b1..7fce6ba91b 100644
--- a/test/base/test_indexset.jl
+++ b/test/base/test_indexset.jl
@@ -372,5 +372,3 @@ using Compat
     @test dag(is) == Index[]
   end
 end
-
-nothing
diff --git a/test/base/test_itensor.jl b/test/base/test_itensor.jl
index bcf4a866a7..73a9e47d34 100644
--- a/test/base/test_itensor.jl
+++ b/test/base/test_itensor.jl
@@ -1847,5 +1847,3 @@ end # End Dense ITensor basic functionality
 
 # Disable debug checking once tests are completed
 ITensors.disable_debug_checks()
-
-nothing
diff --git a/test/base/test_itensor_scalar.jl b/test/base/test_itensor_scalar.jl
index 4968559d49..26cf4ebc26 100644
--- a/test/base/test_itensor_scalar.jl
+++ b/test/base/test_itensor_scalar.jl
@@ -65,5 +65,3 @@ using Test
   ITensor(fill(x, (1, 1, 1))) == ITensor(x)
   @test_throws ErrorException ITensor(fill(x, (2, 2)))
 end
-
-nothing
diff --git a/test/base/test_iterativesolvers.jl b/test/base/test_iterativesolvers.jl
deleted file mode 100644
index e704f8cc1d..0000000000
--- a/test/base/test_iterativesolvers.jl
+++ /dev/null
@@ -1,29 +0,0 @@
-using ITensors
-using Test
-
-# Wrap an ITensor with pairs of primed and
-# unprimed indices to pass to davidson
-struct ITensorMap
-  A::ITensor
-end
-Base.eltype(M::ITensorMap) = eltype(M.A)
-Base.size(M::ITensorMap) = dim(IndexSet(filterinds(M.A; plev=0)...))
-(M::ITensorMap)(v::ITensor) = noprime(M.A * v)
-
-@testset "Complex davidson" begin
-  d = 10
-  i = Index(d, "i")
-  A = randomITensor(ComplexF64, i, i')
-  A = mapprime(A * mapprime(dag(A), 0 => 2), 2 => 1)
-  M = ITensorMap(A)
-
-  v = randomITensor(i)
-  λ, v = davidson(M, v; maxiter=10)
-  @test M(v) ≈ λ * v
-
-  v = randomITensor(ComplexF64, i)
-  λ, v = davidson(M, v; maxiter=10)
-  @test M(v) ≈ λ * v
-end
-
-nothing
diff --git a/test/base/test_lattices.jl b/test/base/test_lattices.jl
index ee2a9782f5..0a86ae7a35 100644
--- a/test/base/test_lattices.jl
+++ b/test/base/test_lattices.jl
@@ -12,5 +12,3 @@ end
   tL = triangular_lattice(3, 4; yperiodic=true)
   @test length(tL) == 28 # inc. periodic vertical bonds
 end
-
-nothing
diff --git a/test/base/test_not.jl b/test/base/test_not.jl
index 9278ab186c..f6bdcac81a 100644
--- a/test/base/test_not.jl
+++ b/test/base/test_not.jl
@@ -51,5 +51,3 @@ using ITensors, Test
 
   @test_throws MethodError !1
 end
-
-nothing
diff --git a/test/base/test_phys_site_types.jl b/test/base/test_phys_site_types.jl
index e2be8e9f67..6ff1154893 100644
--- a/test/base/test_phys_site_types.jl
+++ b/test/base/test_phys_site_types.jl
@@ -846,5 +846,3 @@ using ITensors, LinearAlgebra, Test
     @test op(s, "a", 2) ≈ itensor([0 1 0; 0 0 √2; 0 0 0], s[2]', dag(s[2]))
   end
 end
-
-nothing
diff --git a/test/base/test_qn.jl b/test/base/test_qn.jl
index 2b7fb625ec..27063617e0 100644
--- a/test/base/test_qn.jl
+++ b/test/base/test_qn.jl
@@ -148,5 +148,3 @@ import ITensors: nactive
     @test QN("T", -3, 3) == QN("T", 0, 3)
   end
 end
-
-nothing
diff --git a/test/base/test_qncombiner.jl b/test/base/test_qncombiner.jl
index a69b9bb95b..147db15a44 100644
--- a/test/base/test_qncombiner.jl
+++ b/test/base/test_qncombiner.jl
@@ -10,5 +10,3 @@ using ITensors, Test
   Ã = AC * dag(C)
   @test à ≈ A
 end
-
-nothing
diff --git a/test/base/test_qndiagitensor.jl b/test/base/test_qndiagitensor.jl
index fa7b344fa6..178b71175b 100644
--- a/test/base/test_qndiagitensor.jl
+++ b/test/base/test_qndiagitensor.jl
@@ -132,5 +132,3 @@ using ITensors, Test
     @test sU isa String
   end
 end
-
-nothing
diff --git a/test/base/test_qnindex.jl b/test/base/test_qnindex.jl
index b657e0a058..9ba14f5eea 100644
--- a/test/base/test_qnindex.jl
+++ b/test/base/test_qnindex.jl
@@ -68,5 +68,3 @@ import ITensors: In, Out, Neither
     @test_throws ErrorException ITensors.directsum(i, dag(j))
   end
 end
-
-nothing
diff --git a/test/base/test_qnitensor.jl b/test/base/test_qnitensor.jl
index 59d8f5b8fb..1b359afb06 100644
--- a/test/base/test_qnitensor.jl
+++ b/test/base/test_qnitensor.jl
@@ -1895,5 +1895,3 @@ Random.seed!(1234)
     @test !ITensors.have_same_qns([QN(0) => 1, QN(0) => 2, QN(("Sz", 2)) => 1])
   end
 end
-
-nothing
diff --git a/test/base/test_sitetype.jl b/test/base/test_sitetype.jl
index 60cb583f4f..aa2d957bb3 100644
--- a/test/base/test_sitetype.jl
+++ b/test/base/test_sitetype.jl
@@ -568,5 +568,3 @@ end
     @test is_unitary(U; rtol=1e-15)
   end
 end
-
-nothing
diff --git a/test/base/test_smallstring.jl b/test/base/test_smallstring.jl
index 03de75a018..3c40987106 100644
--- a/test/base/test_smallstring.jl
+++ b/test/base/test_smallstring.jl
@@ -71,5 +71,3 @@ import ITensors: SmallString, Tag, isint, isnull, IntChar
     @test sprint(show, t) == "Electron"
   end
 end
-
-nothing
diff --git a/test/base/test_svd.jl b/test/base/test_svd.jl
index bf6a629e00..19b6d2eca0 100644
--- a/test/base/test_svd.jl
+++ b/test/base/test_svd.jl
@@ -230,5 +230,3 @@ include(joinpath(@__DIR__, "utils", "util.jl"))
   ##   # @test F.U * F.S * F.V ≈ T
   ## end
 end
-
-nothing
diff --git a/test/base/test_sweepnext.jl b/test/base/test_sweepnext.jl
index 1cee8cde55..c74008c70b 100644
--- a/test/base/test_sweepnext.jl
+++ b/test/base/test_sweepnext.jl
@@ -49,5 +49,3 @@ using ITensors, Test
     @test count == 2 * (N - 2) + 1
   end
 end
-
-nothing
diff --git a/test/base/test_sweeps.jl b/test/base/test_sweeps.jl
index a5593cf041..49afcc6e3e 100644
--- a/test/base/test_sweeps.jl
+++ b/test/base/test_sweeps.jl
@@ -149,5 +149,3 @@ using Test
     @test maxdim(sw, 5) == typemax(Int)
   end
 end
-
-nothing
diff --git a/test/base/test_tagset.jl b/test/base/test_tagset.jl
index cf692203d1..794a64f258 100644
--- a/test/base/test_tagset.jl
+++ b/test/base/test_tagset.jl
@@ -143,5 +143,3 @@ using ITensors, Test
     ITensors.set_strict_tags!(false)
   end
 end
-
-nothing
diff --git a/test/runtests.jl b/test/runtests.jl
index 958254c98e..c3195191c0 100644
--- a/test/runtests.jl
+++ b/test/runtests.jl
@@ -35,7 +35,7 @@ ITensors.disable_threaded_blocksparse()
     println(
       """\nArguments ARGS = $(test_args) are empty, or contain `"all"` or `"mps"`. Running MPS/MPO ITensors tests.""",
     )
-    dir = "ITensorLegacyMPS"
+    dir = "ITensorMPS"
     println("\nTest $(@__DIR__)/$(dir)")
     @time include(joinpath(@__DIR__, dir, "runtests.jl"))
     if ARGS ≠ test_args