diff --git a/.github/workflows/Downstream.yml b/.github/workflows/Downstream.yml
index b6b55bb9bb..9e933c7194 100644
--- a/.github/workflows/Downstream.yml
+++ b/.github/workflows/Downstream.yml
@@ -16,6 +16,7 @@ jobs:
         os: [ubuntu-latest]
         package:
           - {user: ITensor, repo: ITensorGaussianMPS.jl}
+          - {user: ITensor, repo: ITensorMPS.jl}
           - {user: ITensor, repo: ITensorTDVP.jl}
           - {user: ITensor, repo: ITensorUnicodePlots.jl}
           - {user: ITensor, repo: ITensorVisualizationBase.jl}
diff --git a/Project.toml b/Project.toml
index 9fd5c9f574..5420e6128a 100644
--- a/Project.toml
+++ b/Project.toml
@@ -1,7 +1,7 @@
 name = "ITensors"
 uuid = "9136182c-28ba-11e9-034c-db9fb085ebd5"
 authors = ["Matthew Fishman <mfishman@flatironinstitute.org>", "Miles Stoudenmire <mstoudenmire@flatironinstitute.org>"]
-version = "0.6.0"
+version = "0.6.1"
 
 [deps]
 Adapt = "79e6a3ab-5dfb-504d-930d-738a2a938a0e"
diff --git a/docs/settings.jl b/docs/settings.jl
index 8101896cb5..9d9029ceab 100644
--- a/docs/settings.jl
+++ b/docs/settings.jl
@@ -42,7 +42,7 @@ settings = Dict(
         "Observer.md",
         "DMRGObserver.md",
       ],
-      "OpSum (AutoMPO)" => "OpSum.md",
+      "OpSum" => "OpSum.md",
     ],
     "Frequently Asked Questions" => [
       "Programming Language (Julia, C++, ...) FAQs" => "faq/JuliaAndCpp.md",
diff --git a/docs/src/AdvancedUsageGuide.md b/docs/src/AdvancedUsageGuide.md
index 3cfbb04e25..7441c1952a 100644
--- a/docs/src/AdvancedUsageGuide.md
+++ b/docs/src/AdvancedUsageGuide.md
@@ -141,6 +141,8 @@ NDTensors.Dense{Float64,Array{Float64,1}}
 A common place you might accidentally come across this is when
 you are creating a Hamiltonian with `OpSum`:
 ```julia
+julia> using ITensors, ITensorMPS
+
 julia> N = 4;
 
 julia> sites = siteinds("S=1/2",N);
@@ -180,7 +182,7 @@ Julia provides many tools for searching for documentation interactively at the R
 julia> using ITensors
 
 julia> ?ITensor
-search: ITensor ITensors itensor emptyITensor randomITensor
+search: ITensor ITensors itensor randomITensor
 
   An ITensor is a tensor whose interface is independent of its
   memory layout. Therefore it is not necessary to know the ordering
@@ -220,27 +222,7 @@ julia> fieldnames(ITensor)
 ```
 which shows the fields of a type. Note that in general the specific names of the fields and structures of types may change (we consider those to be internal details), however we often make functions to access the fields of a type that have the same name as the field, so it is a good place to get started. For example, you can access the storage and indices of an ITensor `A` with the functions `store(A)` and `inds(A)`.
 
-Another helpful function is `apropos`, which search through all documentation for a string (ignoring the case) and prints a list of all types and methods with documentation that contain the string. For example:
-```julia
-julia> apropos("IndexSet")
-ITensors.IndexSet
-ITensors.push
-ITensors.insertat
-ITensors.getfirst
-ITensors.commoninds
-ITensors.pushfirst
-NDTensors.mindim
-[...]
-```
-This can often return too much information. A helpful way to narrow down the search is with regular expressions, for example:
-```julia
-julia> apropos(r"ITensor.*IndexSet")
-ITensors.block
-ITensors.hasinds
-ITensors.ITensor
-NDTensors.inds
-```
-where the notation `r"..."` is Julia notation for making a string that will be interpreted as a [regular expression](https://docs.julialang.org/en/v1/manual/strings/#Regular-Expressions). Here, we are searching for any documentation that contains the string "ITensor" followed at some point by "IndexSet". The notation `.*` is regular expression notation for matching any number of any type of character.
+Another helpful function is `apropos`, which search through all documentation for a string (ignoring the case) and prints a list of all types and methods with documentation that contain the string.
 
 Based on the `apropos` function, we can make some helper functions that may be useful. For example:
 ```julia
@@ -949,7 +931,7 @@ This allows us to have code like:
 julia> i = Index(2, "i")
 (dim=2|id=811|"i")
 
-julia> A = emptyITensor(i', i);
+julia> A = ITensor(i', i);
 
 julia> @show A;
 A = ITensor ord=2
diff --git a/docs/src/ContractionSequenceOptimization.md b/docs/src/ContractionSequenceOptimization.md
index a441fdf9a1..5991946e80 100644
--- a/docs/src/ContractionSequenceOptimization.md
+++ b/docs/src/ContractionSequenceOptimization.md
@@ -38,11 +38,11 @@ h₃ = Index(k, "h₃")
 s₁ = Index(d, "s₁")
 s₂ = Index(d, "s₂")
 
-H₁ = emptyITensor(dag(s₁), s₁', dag(h₁), h₂)
-H₂ = emptyITensor(dag(s₂), s₂', dag(h₂), h₃)
-L = emptyITensor(dag(l), l', h₁)
-R = emptyITensor(dag(r), r', h₃)
-ψ = emptyITensor(l, s₁, s₂, r)
+H₁ = ITensor(dag(s₁), s₁', dag(h₁), h₂)
+H₂ = ITensor(dag(s₂), s₂', dag(h₂), h₃)
+L = ITensor(dag(l), l', h₁)
+R = ITensor(dag(r), r', h₃)
+ψ = ITensor(l, s₁, s₂, r)
 
 TN = [ψ, L, H₁, H₂, R]
 sequence1 = Any[2, Any[3, Any[4, Any[1, 5]]]]
@@ -80,11 +80,11 @@ function tensor_network(; m, k, d)
   s₁ = Index(d, "s₁")
   s₂ = Index(d, "s₂")
 
-  ψ = emptyITensor(l, s₁, s₂, r)
-  L = emptyITensor(dag(l), l', h₁)
-  H₁ = emptyITensor(dag(s₁), s₁', dag(h₁), h₂)
-  H₂ = emptyITensor(dag(s₂), s₂', dag(h₂), h₃)
-  R = emptyITensor(dag(r), r', h₃)
+  ψ = ITensor(l, s₁, s₂, r)
+  L = ITensor(dag(l), l', h₁)
+  H₁ = ITensor(dag(s₁), s₁', dag(h₁), h₂)
+  H₂ = ITensor(dag(s₂), s₂', dag(h₂), h₃)
+  R = ITensor(dag(r), r', h₃)
   return [ψ, L, H₁, H₂, R]
 end
 
diff --git a/docs/src/Multithreading.md b/docs/src/Multithreading.md
index c7c31613b3..982da86a5c 100644
--- a/docs/src/Multithreading.md
+++ b/docs/src/Multithreading.md
@@ -34,8 +34,6 @@ julia> using LinearAlgebra
 julia> BLAS.vendor()  # Check which BLAS you are using
 :mkl
 
-julia> using ITensors
-
 julia> BLAS.get_num_threads()
 6
 
@@ -60,17 +58,19 @@ $ julia -t 4
 julia> Threads.nthreads()
 4
 
-julia> ITensors.Strided.get_num_threads()
+julia> using Strided
+
+julia> Strided.get_num_threads()
 4
 ```
 We find that this threading competes with BLAS threading as well as ITensors.jl's own block sparse
 multithreading, so if you are using Julia with multiple threads you may want to disable
 Strided.jl's threading with:
 ```julia
-julia> ITensors.Strided.disable_threads()
+julia> Strided.disable_threads()
 1
 
-julia> ITensors.Strided.get_num_threads()
+julia> Strided.get_num_threads()
 1
 ```
 in favor of either BLAS threading or ITensors.jl's block sparse threading.
@@ -88,12 +88,13 @@ Here is a simple example of using block sparse multithreading to speed up a spar
 tensor contraction:
 ```julia
 using BenchmarkTools
-using ITensors
+using ITensors, ITensorMPS
 using LinearAlgebra
+using Strided
 
 function main(; d = 20, order = 4)
   BLAS.set_num_threads(1)
-  ITensors.Strided.set_num_threads(1)
+  Strided.set_num_threads(1)
 
   println("#################################################")
   println("# order = ", order)
@@ -139,7 +140,7 @@ julia> main(d = 20, order = 4)
 Threads.nthreads() = 5
 Sys.CPU_THREADS = 6
 BLAS.get_num_threads() = 1
-ITensors.Strided.get_num_threads() = 1
+Strided.get_num_threads() = 1
 
 Serial contract:
   21.558 ms (131 allocations: 7.34 MiB)
diff --git a/docs/src/Observer.md b/docs/src/Observer.md
index 37cb985062..57bce78631 100644
--- a/docs/src/Observer.md
+++ b/docs/src/Observer.md
@@ -30,6 +30,8 @@ For example, let's make a type called `DemoObserver`
 as:
 
 ```julia
+using ITensors, ITensorMPS
+
 mutable struct DemoObserver <: AbstractObserver
    energy_tol::Float64
    last_energy::Float64
@@ -137,7 +139,7 @@ energy, psi = dmrg(H,psi0,sweeps; observer=obs, outputlevel=1)
 ## Complete Sample Code
 
 ```julia
-using ITensors
+using ITensors, ITensorMPS
 
 mutable struct DemoObserver <: AbstractObserver
    energy_tol::Float64
diff --git a/docs/src/examples/DMRG.md b/docs/src/examples/DMRG.md
index e5f48ad0fb..44a6227705 100644
--- a/docs/src/examples/DMRG.md
+++ b/docs/src/examples/DMRG.md
@@ -7,6 +7,7 @@ or a wavefunction we need to construct a "site set" which will hold the site ind
 the physical Hilbert space:
 
 ```julia
+using ITensors, ITensorMPS
 N = 100
 sites = siteinds("S=1",N)
 ```
@@ -66,7 +67,7 @@ approximation to the ground state as the variable `psi`.
 Below you can find a complete working code that includes all of these steps:
 
 ```julia
-using ITensors
+using ITensors, ITensorMPS
 
 let
   N = 100
@@ -118,7 +119,7 @@ These tags tell the OpSum system which local operators to use for these
 sites when building the Hamiltonian MPO.
 
 ```julia
-using ITensors
+using ITensors, ITensorMPS
 
 let
   N = 100
@@ -208,7 +209,7 @@ the geometry a cylinder.
 **Full example code:**
 
 ```julia
-using ITensors
+using ITensors, ITensorMPS
 
 let
   Ny = 6
@@ -285,7 +286,7 @@ within the same quantum number sector.
 **Full Example code:**
 
 ```julia
-using ITensors
+using ITensors, ITensorMPS
 
 let
   N = 20
@@ -377,6 +378,8 @@ First we define our custom observer type, `EntanglementObserver`, and overload t
 for it:
 
 ```julia
+using ITensors, ITensorMPS
+
 mutable struct EntanglementObserver <: AbstractObserver
 end
 
@@ -398,7 +401,7 @@ has just finished optimizing.
 Here is a complete sample code including constructing the observer and passing it to DMRG:
 
 ```julia
-using ITensors
+using ITensors, ITensorMPS
 
 mutable struct EntanglementObserver <: AbstractObserver
 end
@@ -473,6 +476,8 @@ First we define our custom observer type, `SizeObserver`, and overload the `meas
 for it:
 
 ```julia
+using ITensors, ITensorMPS
+
 mutable struct SizeObserver <: AbstractObserver
 end
 
@@ -493,7 +498,7 @@ When it runs, it calls `Base.summarysize` on the wavefunction `psi` object and t
 Here is a complete sample code including constructing the observer and passing it to DMRG:
 
 ```julia
-using ITensors
+using ITensors, ITensorMPS
 
 mutable struct SizeObserver <: AbstractObserver
 end
diff --git a/docs/src/examples/MPSandMPO.md b/docs/src/examples/MPSandMPO.md
index a9677dc25f..9f1891dd2e 100644
--- a/docs/src/examples/MPSandMPO.md
+++ b/docs/src/examples/MPSandMPO.md
@@ -83,7 +83,7 @@ In the example code below we will obtain the element ``T^{1,2,1,1,2,1,2,2,2,1}``
 which is (implicitly) defined by the MPS psi:
 
 ```@example mps_element
-using ITensors # hide
+using ITensors, ITensorMPS
 let # hide
 N = 10
 s = siteinds(2,N)
@@ -139,7 +139,7 @@ following example we will use a random MPS of bond dimension ``\chi=4`` but the
 MPS could be obtained other ways such as through a DMRG calculation.
 
 ```@example expect
-using ITensors # hide
+using ITensors, ITensorMPS
 N = 10
 chi = 4
 sites = siteinds("S=1/2",N)
@@ -306,6 +306,8 @@ as well as limits on the maximum bond dimension (`maxdim` keyword argument).
 **Complete code example**
 
 ```julia
+using ITensors, ITensorMPS
+
 psi = orthogonalize(psi, 3)
 
 wf = (psi[3] * psi[4]) * G
@@ -355,7 +357,8 @@ the algorithm for sampling MPS is a `perfect' sampling algorithm with no autocor
 In more detail, say we have a set of `N` site indices `s` and define a random MPS
 with these sites:
 ```@example sample_mps; continued=true
-using ITensors # hide
+using ITensors, ITensorMPS
+
 N = 10 # number of sites
 d = 3  # dimension of each site
 chi = 16 # bond dimension of the MPS
@@ -384,7 +387,7 @@ can be brought into orthogonal form by calling `psi = orthogonalize(psi, 1)`.
 
 ## Write and Read an MPS or MPO to Disk with HDF5
 
-!!! info 
+!!! info
     Make sure to install the HDF5 package to use this feature. (Run `julia> ] add HDF5` in the Julia REPL console.)
 
 **Writing an MPS to an HDF5 File**
@@ -427,6 +430,8 @@ So for example, to create an MPO from an OpSum which has the same site indices
 as your MPS `psi`, do the following:
 
 ```julia
+using ITensors, ITensorMPS
+
 os = OpSum()
 # Then put operators into os...
 
diff --git a/docs/src/examples/Physics.md b/docs/src/examples/Physics.md
index 99855bd294..751e76e4f4 100644
--- a/docs/src/examples/Physics.md
+++ b/docs/src/examples/Physics.md
@@ -9,15 +9,21 @@ can be found in the section [SiteTypes Included with ITensor](@ref).
 If you have an Index `s` carrying a "S=1/2" tag, for example, you can obtain the "Sz"
 operator like this:
 ```julia
+using ITensors, ITensorMPS
+
 op("Sz",s)
 ```
 
 Usually indices with physical tags come from an array of indices returned from the `siteinds` function
 ```julia
+using ITensors, ITensorMPS
+
+N = 10
 sites = siteinds("S=1/2",N)
 ```
 in which case one might want the "Sz" operator on site 4
 ```julia
+using ITensors, ITensorMPS
 Sz4 = op("Sz",sites[4])
 ```
 
@@ -34,6 +40,7 @@ M = [1/2 0 ; 0 -1/2]
 ```
 by calling
 ```julia
+using ITensors, ITensorMPS
 Sz = op(M,s)
 ```
 
@@ -51,6 +58,7 @@ Perhaps the most common part of the site type system one wishes to extend
 are the various `op` or `op!` function overloads which allow code like
 
 ```julia
+using ITensors, ITensorMPS
 s = siteind("S=1/2")
 Sz = op("Sz",s)
 ```
@@ -84,8 +92,7 @@ To add this operator to the ITensor `op` system, we just need to introduce the f
 code
 
 ```julia
-using ITensors
-
+using ITensors, ITensorMPS
 ITensors.op(::OpName"Pup",::SiteType"S=1/2") =
  [1 0
   0 0]
@@ -100,6 +107,7 @@ other functions of the name `op` inside the ITensors module.
 Having defined the above code, we can now do things like
 
 ```julia
+using ITensors, ITensorMPS
 s = siteind("S=1/2")
 Pup = op("Pup",s)
 ```
@@ -108,6 +116,7 @@ to obtain the `"Pup"` operator for our `"S=1/2"` Index `s`. Or we can do a simil
 thing for an array of site indices:
 
 ```julia
+using ITensors, ITensorMPS
 N = 40
 s = siteinds("S=1/2",N)
 Pup1 = op("Pup",s[1])
@@ -117,8 +126,7 @@ Pup3 = op("Pup",s[3])
 Note that for the `"Qudit"`/`"Boson"` site types, you have to define your overload
 of `op` with the dimension of the local Hilbert space, for example:
 ```julia
-using ITensors
-
+using ITensors, ITensorMPS
 function ITensors.op(::OpName"P1", ::SiteType"Boson", d::Int)
   o = zeros(d, d)
   o[1, 1] = 1
@@ -127,6 +135,7 @@ end
 ```
 Alternatively you could use Julia's [array comprehension](https://docs.julialang.org/en/v1/manual/arrays/#man-comprehensions) syntax:
 ```julia
+using ITensors, ITensorMPS
 ITensors.op(::OpName"P1", ::SiteType"Boson", d::Int) =
   [(i == j == 1) ? 1.0 : 0.0 for i in 1:d, j in 1:d]
 ```
@@ -142,6 +151,7 @@ indices.
 For example, we could now make an OpSum involving our custom operator such as:
 
 ```julia
+using ITensors, ITensorMPS
 N = 100
 sites = siteinds("S=1/2",N)
 os = OpSum()
@@ -166,6 +176,7 @@ Say we want to define a new state for the "Electron" site type called "+", which
 the meaning of one electron with its spin in the +X direction. First let's review
 the existing state definitions:
 ```julia
+using ITensors, ITensorMPS
 ITensors.state(::StateName"Emp", ::SiteType"Electron") = [1.0, 0, 0, 0]
 ITensors.state(::StateName"Up", ::SiteType"Electron") = [0.0, 1, 0, 0]
 ITensors.state(::StateName"Dn", ::SiteType"Electron") = [0.0, 0, 1, 0]
@@ -186,6 +197,7 @@ which makes the state
 Having defined this overload of `state`, if we have an Index of type "Electron"
 we can obtain our new state for it by doing
 ```julia
+using ITensors, ITensorMPS
 s = siteind("Electron")
 plus = state("+",s)
 ```
@@ -255,7 +267,7 @@ First let's see the minimal code needed to define and use this new
 the code is doing.
 
 ```julia
-using ITensors
+using ITensors, ITensorMPS
 
 ITensors.space(::SiteType"S=3/2") = 4
 
@@ -308,6 +320,7 @@ describe the vector space corresponding to that site type. For our
 carrying the `"S=3/2"` tag, the definition is
 
 ```julia
+using ITensors, ITensorMPS
 ITensors.space(::SiteType"S=3/2") = 4
 ```
 
@@ -325,12 +338,14 @@ numbers.
 After defining this `space` function, you can just write code like:
 
 ```julia
+using ITensors, ITensorMPS
 s = siteind("S=3/2")
 ```
 
 to obtain a single `"S=3/2"` Index, or write code like
 
 ```julia
+using ITensors, ITensorMPS
 N = 100
 sites = siteinds("S=3/2",N)
 ```
@@ -352,8 +367,7 @@ In our example above, we defined this function for the case of the `"Sz"`
 operator as:
 
 ```@example S32
-using ITensors # hide
-
+using ITensors, ITensorMPS
 ITensors.op(::OpName"Sz",::SiteType"S=3/2") =
   [+3/2   0    0    0
      0  +1/2   0    0
@@ -379,6 +393,7 @@ then, for example, you can get the `"Sz"` operator for this Index
 and print it out by doing:
 
 ```@example S32
+using ITensors, ITensorMPS
 Sz = op("Sz",s)
 println(Sz)
 ```
@@ -392,6 +407,7 @@ You can use the `op` function yourself with a set of site indices created from
 the `siteinds` function like this:
 
 ```julia
+using ITensors, ITensorMPS
 N = 100
 sites = siteinds("S=3/2",N)
 Sz1 = op("Sz",sites[1])
@@ -439,7 +455,7 @@ quantum numbers of our ``S=3/2`` site type, then we will discuss what each part
 the code is doing.
 
 ```julia
-using ITensors
+using ITensors, ITensorMPS
 
 function ITensors.space(::SiteType"S=3/2";
                         conserve_qns=false)
@@ -491,6 +507,7 @@ it indicates the dimension of the subspace.
 After defining the `space` function this way, you can write code like:
 
 ```julia
+using ITensors, ITensorMPS
 s = siteind("S=3/2"; conserve_qns=true)
 ```
 
@@ -502,6 +519,7 @@ which are conventional for all site indices.
 You can now also call code like:
 
 ```julia
+using ITensors, ITensorMPS
 N = 100
 sites = siteinds("S=3/2",N; conserve_qns=true)
 ```
diff --git a/docs/src/index.md b/docs/src/index.md
index 550d5ea4ec..3e1e5adddb 100644
--- a/docs/src/index.md
+++ b/docs/src/index.md
@@ -29,6 +29,8 @@ Development of ITensor is supported by the Flatiron Institute, a division of the
 
 ## News
 
+- May 9, 2024: A new package ITensorMPS.jl has been released. We plan to move all of the MPS/MPO code in ITensors.jl to ITensorMPS.jl. For now, ITensorMPS.jl just re-exports the MPS/MPO functionality of ITensors.jl (as well as of [ITensorTDVP.jl](https://github.com/ITensor/ITensorTDVP.jl)), such as `dmrg`, `siteinds`, `MPS`, `MPO`, etc. To prepare for the change over to ITensorMPS.jl, please change `using ITensors` to `using ITensors, ITensorMPS` in any code that makes use of MPS/MPO functionality, and if you are using ITensorTDVP.jl change `using ITensorTDVP` to `using ITensorMPS` in your code.
+
 - May 8, 2024: ITensors.jl v0.6 has been released. This version deletes the experimental "combine-contract" contraction backend, which was enabled by `ITensors.enable_combine_contract()`. This feature enabled performing ITensor contractions by first combining indices and then performing contractions as matrix multiplications, which potentially could lead to speedups for certain contractions involving higher-order QN-conserving tensors. However, the speedups weren't consistent with the current implementation, and this feature will be incorporated into the library in a more systematic way when we release our new non-abelian symmetric tensor backend.
 
 - May 2, 2024: ITensors.jl v0.5 has been released. This version removes PackageCompiler.jl as a dependency and moves the package compilation functionality into a package extension. In order to use the `ITensors.compile()` function going forward, you need to install the PackageCompiler.jl package with `using Pkg: Pkg; Pkg.add("PackageCompiler")` and put `using PackageCompiler` together with `using ITensors` in your code.
@@ -273,7 +275,7 @@ eigenvectors of Hamiltonian (energy) matrices which model
 the behavior of quantum systems.
 
 ```julia
-using ITensors
+using ITensors, ITensorMPS
 let
   # Create 100 spin-one indices
   N = 100
diff --git a/docs/src/tutorials/DMRG.md b/docs/src/tutorials/DMRG.md
index 69d89cd594..428274ff1f 100644
--- a/docs/src/tutorials/DMRG.md
+++ b/docs/src/tutorials/DMRG.md
@@ -25,7 +25,7 @@ discuss the main steps. If you need help running the code below, see the getting
 started page on [Running ITensor and Julia Codes](@ref).
 
 ```julia
-using ITensors
+using ITensors, ITensorMPS
 let
   N = 100
   sites = siteinds("S=1",N)
@@ -56,7 +56,7 @@ end
 The first two lines
 
 ```@example siteinds; continued=true
-using ITensors # hide
+using ITensors, ITensorMPS
 N = 100
 sites = siteinds("S=1",N)
 ```
diff --git a/docs/src/tutorials/MPSTimeEvolution.md b/docs/src/tutorials/MPSTimeEvolution.md
index 31bbd8622f..58886a6c6b 100644
--- a/docs/src/tutorials/MPSTimeEvolution.md
+++ b/docs/src/tutorials/MPSTimeEvolution.md
@@ -63,7 +63,7 @@ discuss the main steps. (If you need help running the code below, see the gettin
 started page on running ITensor codes.)
 
 ```julia
-using ITensors
+using ITensors, ITensorMPS
 
 let
   N = 100
diff --git a/docs/src/tutorials/QN_DMRG.md b/docs/src/tutorials/QN_DMRG.md
index bb3b988058..390e3a04ba 100644
--- a/docs/src/tutorials/QN_DMRG.md
+++ b/docs/src/tutorials/QN_DMRG.md
@@ -153,7 +153,7 @@ to turn it into a code which conserves the total ``S^z`` quantum
 number throughout the DMRG calculation. The resulting code is:
 
 ```julia
-using ITensors
+using ITensors, ITensorMPS
 let
   N = 100
   sites = siteinds("S=1",N;conserve_qns=true)
diff --git a/docs/src/tutorials/tebd.jl b/docs/src/tutorials/tebd.jl
index e811633001..7dfb452b7d 100644
--- a/docs/src/tutorials/tebd.jl
+++ b/docs/src/tutorials/tebd.jl
@@ -1,4 +1,4 @@
-using ITensors
+using ITensors, ITensorMPS
 
 let
   N = 100
diff --git a/src/lib/ITensorMPS/examples/autodiff/circuit_optimization/op.jl b/src/lib/ITensorMPS/examples/autodiff/circuit_optimization/op.jl
index 10a3cde7b9..bdb40e6334 100644
--- a/src/lib/ITensorMPS/examples/autodiff/circuit_optimization/op.jl
+++ b/src/lib/ITensorMPS/examples/autodiff/circuit_optimization/op.jl
@@ -1,4 +1,4 @@
-using ITensors
+using ITensors, ITensorMPS
 using Zygote
 
 s = siteind("Qubit")
diff --git a/src/lib/ITensorMPS/examples/autodiff/circuit_optimization/state_preparation.jl b/src/lib/ITensorMPS/examples/autodiff/circuit_optimization/state_preparation.jl
index d8050d108f..bb09325f83 100644
--- a/src/lib/ITensorMPS/examples/autodiff/circuit_optimization/state_preparation.jl
+++ b/src/lib/ITensorMPS/examples/autodiff/circuit_optimization/state_preparation.jl
@@ -1,4 +1,4 @@
-using ITensors
+using ITensors, ITensorMPS
 using OptimKit
 using Random
 using Zygote
diff --git a/src/lib/ITensorMPS/examples/autodiff/circuit_optimization/vqe.jl b/src/lib/ITensorMPS/examples/autodiff/circuit_optimization/vqe.jl
index 35b1751ed5..5ac27f66a5 100644
--- a/src/lib/ITensorMPS/examples/autodiff/circuit_optimization/vqe.jl
+++ b/src/lib/ITensorMPS/examples/autodiff/circuit_optimization/vqe.jl
@@ -1,4 +1,4 @@
-using ITensors
+using ITensors, ITensorMPS
 using OptimKit
 using Random
 using Zygote
diff --git a/src/lib/ITensorMPS/examples/autodiff/mps_autodiff.jl b/src/lib/ITensorMPS/examples/autodiff/mps_autodiff.jl
index 5130a8a641..41a16f05ec 100644
--- a/src/lib/ITensorMPS/examples/autodiff/mps_autodiff.jl
+++ b/src/lib/ITensorMPS/examples/autodiff/mps_autodiff.jl
@@ -1,4 +1,4 @@
-using ITensors
+using ITensors, ITensorMPS
 using OptimKit
 using Zygote
 
diff --git a/src/lib/ITensorMPS/examples/autodiff/ops/ops_ad.jl b/src/lib/ITensorMPS/examples/autodiff/ops/ops_ad.jl
index 024ea37dca..47bb051ebd 100644
--- a/src/lib/ITensorMPS/examples/autodiff/ops/ops_ad.jl
+++ b/src/lib/ITensorMPS/examples/autodiff/ops/ops_ad.jl
@@ -1,4 +1,4 @@
-using ITensors
+using ITensors, ITensorMPS
 using ITensors.Ops
 using Zygote
 
diff --git a/src/lib/ITensorMPS/examples/autodiff/ops/trotter_ad_1.jl b/src/lib/ITensorMPS/examples/autodiff/ops/trotter_ad_1.jl
index eabb9db64f..c6c8fc2eb5 100644
--- a/src/lib/ITensorMPS/examples/autodiff/ops/trotter_ad_1.jl
+++ b/src/lib/ITensorMPS/examples/autodiff/ops/trotter_ad_1.jl
@@ -1,4 +1,4 @@
-using ITensors
+using ITensors, ITensorMPS
 using Zygote
 using OptimKit
 
diff --git a/src/lib/ITensorMPS/examples/autodiff/ops/trotter_ad_2.jl b/src/lib/ITensorMPS/examples/autodiff/ops/trotter_ad_2.jl
index b880f0adb1..163dcc7d67 100644
--- a/src/lib/ITensorMPS/examples/autodiff/ops/trotter_ad_2.jl
+++ b/src/lib/ITensorMPS/examples/autodiff/ops/trotter_ad_2.jl
@@ -1,4 +1,4 @@
-using ITensors
+using ITensors, ITensorMPS
 using Zygote
 using OptimKit
 
diff --git a/src/lib/ITensorMPS/examples/autodiff/ops/vqe.jl b/src/lib/ITensorMPS/examples/autodiff/ops/vqe.jl
index e1ad7979aa..6757a8d3ce 100644
--- a/src/lib/ITensorMPS/examples/autodiff/ops/vqe.jl
+++ b/src/lib/ITensorMPS/examples/autodiff/ops/vqe.jl
@@ -1,4 +1,4 @@
-using ITensors
+using ITensors, ITensorMPS
 using OptimKit
 using Random
 using Zygote
diff --git a/src/lib/ITensorMPS/examples/dmrg/1d_heisenberg.jl b/src/lib/ITensorMPS/examples/dmrg/1d_heisenberg.jl
index 57bca58dc6..1b9111b2f3 100644
--- a/src/lib/ITensorMPS/examples/dmrg/1d_heisenberg.jl
+++ b/src/lib/ITensorMPS/examples/dmrg/1d_heisenberg.jl
@@ -1,4 +1,4 @@
-using ITensors
+using ITensors, ITensorMPS
 using Printf
 using Random
 
diff --git a/src/lib/ITensorMPS/examples/dmrg/1d_heisenberg_conserve_spin.jl b/src/lib/ITensorMPS/examples/dmrg/1d_heisenberg_conserve_spin.jl
index 33c60383f4..e9b44b2be4 100644
--- a/src/lib/ITensorMPS/examples/dmrg/1d_heisenberg_conserve_spin.jl
+++ b/src/lib/ITensorMPS/examples/dmrg/1d_heisenberg_conserve_spin.jl
@@ -1,8 +1,8 @@
-using ITensors
-using ITensors.Strided
+using ITensors, ITensorMPS
 using LinearAlgebra
 using Printf
 using Random
+using Strided
 
 Random.seed!(1234)
 BLAS.set_num_threads(1)
diff --git a/src/lib/ITensorMPS/examples/dmrg/1d_hubbard_extended.jl b/src/lib/ITensorMPS/examples/dmrg/1d_hubbard_extended.jl
index 063e59c3e0..9c2deaa730 100644
--- a/src/lib/ITensorMPS/examples/dmrg/1d_hubbard_extended.jl
+++ b/src/lib/ITensorMPS/examples/dmrg/1d_hubbard_extended.jl
@@ -1,4 +1,4 @@
-using ITensors
+using ITensors, ITensorMPS
 
 #
 # DMRG calculation of the extended Hubbard model
diff --git a/src/lib/ITensorMPS/examples/dmrg/1d_ising_with_observer.jl b/src/lib/ITensorMPS/examples/dmrg/1d_ising_with_observer.jl
index 7166f6423f..f4b9847e2c 100644
--- a/src/lib/ITensorMPS/examples/dmrg/1d_ising_with_observer.jl
+++ b/src/lib/ITensorMPS/examples/dmrg/1d_ising_with_observer.jl
@@ -1,7 +1,7 @@
 # In this example we show how to pass a DMRGObserver to
 # the dmrg function which allows tracking energy convergence and
 # convergence of local operators.
-using ITensors
+using ITensors, ITensorMPS
 
 """
   Get MPO of transverse field Ising model Hamiltonian with field strength h
diff --git a/src/lib/ITensorMPS/examples/dmrg/2d_heisenberg_conserve_spin.jl b/src/lib/ITensorMPS/examples/dmrg/2d_heisenberg_conserve_spin.jl
index acdab1d0e9..6278658ffe 100644
--- a/src/lib/ITensorMPS/examples/dmrg/2d_heisenberg_conserve_spin.jl
+++ b/src/lib/ITensorMPS/examples/dmrg/2d_heisenberg_conserve_spin.jl
@@ -1,4 +1,4 @@
-using ITensors
+using ITensors, ITensorMPS
 
 let
   Ny = 6
diff --git a/src/lib/ITensorMPS/examples/dmrg/2d_hubbard_conserve_momentum.jl b/src/lib/ITensorMPS/examples/dmrg/2d_hubbard_conserve_momentum.jl
index d0afe6f25f..2ed348c7f9 100644
--- a/src/lib/ITensorMPS/examples/dmrg/2d_hubbard_conserve_momentum.jl
+++ b/src/lib/ITensorMPS/examples/dmrg/2d_hubbard_conserve_momentum.jl
@@ -1,6 +1,7 @@
-using ITensors
+using ITensors, ITensorMPS
 using LinearAlgebra
 using Random
+using Strided
 
 include(joinpath(@__DIR__, "..", "src", "electronk.jl"))
 include(joinpath(@__DIR__, "..", "src", "hubbard.jl"))
@@ -13,7 +14,7 @@ energy, H, psi = main(; Nx=8, Ny=4, U=4.0, t=1.0, nsweeps=10, maxdim=3000, threa
 energy, H, psi = main(; Nx=8, Ny=4, U=4.0, t=1.0, nsweeps=10, maxdim=3000, random_init=false, threaded_blocksparse=false);
 energy, H, psi = main(; Nx=8, Ny=4, U=4.0, t=1.0, nsweeps=10, maxdim=3000, random_init=false, threaded_blocksparse=true);
 
-using ITensors.HDF5
+using HDF5
 h5open("2d_hubbard_conserve_momentum.h5", "w") do fid
   fid["energy"] = energy
   fid["H"] = H
@@ -49,7 +50,7 @@ function main(;
 
   # Disable other threading
   BLAS.set_num_threads(1)
-  ITensors.Strided.set_num_threads(1)
+  Strided.set_num_threads(1)
 
   ITensors.enable_threaded_blocksparse(threaded_blocksparse)
   @show ITensors.using_threaded_blocksparse()
diff --git a/src/lib/ITensorMPS/examples/dmrg/2d_hubbard_conserve_particles.jl b/src/lib/ITensorMPS/examples/dmrg/2d_hubbard_conserve_particles.jl
index 8c6129a622..928d21f7f0 100644
--- a/src/lib/ITensorMPS/examples/dmrg/2d_hubbard_conserve_particles.jl
+++ b/src/lib/ITensorMPS/examples/dmrg/2d_hubbard_conserve_particles.jl
@@ -1,4 +1,4 @@
-using ITensors
+using ITensors, ITensorMPS
 
 function main(; Nx=6, Ny=3, U=4.0, t=1.0)
   N = Nx * Ny
diff --git a/src/lib/ITensorMPS/examples/dmrg/Project.toml b/src/lib/ITensorMPS/examples/dmrg/Project.toml
index d163dd05ed..e6b91775d5 100644
--- a/src/lib/ITensorMPS/examples/dmrg/Project.toml
+++ b/src/lib/ITensorMPS/examples/dmrg/Project.toml
@@ -1,4 +1,6 @@
 [deps]
+ITensorMPS = "0d1a4710-d33b-49a5-8f18-73bdf49b47e2"
 ITensors = "9136182c-28ba-11e9-034c-db9fb085ebd5"
 LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
 Random = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c"
+Strided = "5e0ebb24-38b0-5f93-81fe-25c709ecae67"
diff --git a/src/lib/ITensorMPS/examples/dmrg/input_files/ArgParse/1d_hubbard_extended.jl b/src/lib/ITensorMPS/examples/dmrg/input_files/ArgParse/1d_hubbard_extended.jl
index 8034b26869..ebe26b9180 100644
--- a/src/lib/ITensorMPS/examples/dmrg/input_files/ArgParse/1d_hubbard_extended.jl
+++ b/src/lib/ITensorMPS/examples/dmrg/input_files/ArgParse/1d_hubbard_extended.jl
@@ -1,4 +1,4 @@
-using ITensors
+using ITensors, ITensorMPS
 
 #
 # DMRG calculation of the extended Hubbard model
diff --git a/src/lib/ITensorMPS/examples/dmrg/input_files/argsdict/1d_hubbard_extended.jl b/src/lib/ITensorMPS/examples/dmrg/input_files/argsdict/1d_hubbard_extended.jl
index fed0a41db8..26681cf346 100644
--- a/src/lib/ITensorMPS/examples/dmrg/input_files/argsdict/1d_hubbard_extended.jl
+++ b/src/lib/ITensorMPS/examples/dmrg/input_files/argsdict/1d_hubbard_extended.jl
@@ -1,4 +1,4 @@
-using ITensors
+using ITensors, ITensorMPS
 
 #
 # DMRG calculation of the extended Hubbard model
diff --git a/src/lib/ITensorMPS/examples/dmrg/input_files/script/1d_hubbard_extended.jl b/src/lib/ITensorMPS/examples/dmrg/input_files/script/1d_hubbard_extended.jl
index 82c5e55a86..51ff4f99c5 100644
--- a/src/lib/ITensorMPS/examples/dmrg/input_files/script/1d_hubbard_extended.jl
+++ b/src/lib/ITensorMPS/examples/dmrg/input_files/script/1d_hubbard_extended.jl
@@ -1,4 +1,4 @@
-using ITensors
+using ITensors, ITensorMPS
 
 #
 # DMRG calculation of the extended Hubbard model
diff --git a/src/lib/ITensorMPS/examples/dmrg/threaded_blocksparse/2d_hubbard_conserve_momentum.jl b/src/lib/ITensorMPS/examples/dmrg/threaded_blocksparse/2d_hubbard_conserve_momentum.jl
index b9b7a01ee9..2df178da49 100644
--- a/src/lib/ITensorMPS/examples/dmrg/threaded_blocksparse/2d_hubbard_conserve_momentum.jl
+++ b/src/lib/ITensorMPS/examples/dmrg/threaded_blocksparse/2d_hubbard_conserve_momentum.jl
@@ -1,4 +1,4 @@
-using ITensors
+using ITensors, ITensorMPS
 using LinearAlgebra
 using Random
 
diff --git a/src/lib/ITensorMPS/examples/dmrg/write_to_disk/1d_heisenberg.jl b/src/lib/ITensorMPS/examples/dmrg/write_to_disk/1d_heisenberg.jl
index 2b9d831f67..ff0373cc9d 100644
--- a/src/lib/ITensorMPS/examples/dmrg/write_to_disk/1d_heisenberg.jl
+++ b/src/lib/ITensorMPS/examples/dmrg/write_to_disk/1d_heisenberg.jl
@@ -1,8 +1,8 @@
-using ITensors
-using ITensors.Strided
+using ITensors, ITensorMPS
 using LinearAlgebra
 using Printf
 using Random
+using Strided
 
 Random.seed!(1234)
 Strided.set_num_threads(1)
@@ -29,7 +29,6 @@ let
   maxdim = [10, 20, 100, 100, 200]
   # Set maximum truncation error allowed when adapting bond dimensions
   cutoff = 1E-10
-  @show sweeps
 
   # Run the DMRG algorithm, returning energy and optimized MPS
   energy, psi = dmrg(H, psi0; nsweeps, cutoff, maxdim, write_when_maxdim_exceeds=25)
diff --git a/src/lib/ITensorMPS/examples/dmrg/write_to_disk/Project.toml b/src/lib/ITensorMPS/examples/dmrg/write_to_disk/Project.toml
new file mode 100644
index 0000000000..abb011c624
--- /dev/null
+++ b/src/lib/ITensorMPS/examples/dmrg/write_to_disk/Project.toml
@@ -0,0 +1,7 @@
+[deps]
+ITensorMPS = "0d1a4710-d33b-49a5-8f18-73bdf49b47e2"
+ITensors = "9136182c-28ba-11e9-034c-db9fb085ebd5"
+LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
+Printf = "de0858da-6303-5e67-8744-51eddeeeb8d7"
+Random = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c"
+Strided = "5e0ebb24-38b0-5f93-81fe-25c709ecae67"
diff --git a/src/lib/ITensorMPS/examples/exact_diagonalization/Project.toml b/src/lib/ITensorMPS/examples/exact_diagonalization/Project.toml
new file mode 100644
index 0000000000..226addea8c
--- /dev/null
+++ b/src/lib/ITensorMPS/examples/exact_diagonalization/Project.toml
@@ -0,0 +1,7 @@
+[deps]
+ITensorMPS = "0d1a4710-d33b-49a5-8f18-73bdf49b47e2"
+ITensors = "9136182c-28ba-11e9-034c-db9fb085ebd5"
+KrylovKit = "0b1a1467-8014-51b9-945f-bf0ae24f4b77"
+LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
+MKL = "33e6dc65-8f57-5167-99aa-e5a354878fb2"
+Strided = "5e0ebb24-38b0-5f93-81fe-25c709ecae67"
diff --git a/src/lib/ITensorMPS/examples/exact_diagonalization/exact_diagonalization.jl b/src/lib/ITensorMPS/examples/exact_diagonalization/exact_diagonalization.jl
index 1b883086c9..389809e4c9 100644
--- a/src/lib/ITensorMPS/examples/exact_diagonalization/exact_diagonalization.jl
+++ b/src/lib/ITensorMPS/examples/exact_diagonalization/exact_diagonalization.jl
@@ -1,11 +1,12 @@
-using ITensors
+using ITensors, ITensorMPS
 using KrylovKit
 using LinearAlgebra
 using MKL
+using Strided
 
 include("fuse_inds.jl")
 
-ITensors.Strided.disable_threads()
+Strided.disable_threads()
 ITensors.disable_threaded_blocksparse()
 
 function heisenberg(n)
diff --git a/src/lib/ITensorMPS/examples/exact_diagonalization/fuse_inds.jl b/src/lib/ITensorMPS/examples/exact_diagonalization/fuse_inds.jl
index 4813559de2..82f90e91da 100644
--- a/src/lib/ITensorMPS/examples/exact_diagonalization/fuse_inds.jl
+++ b/src/lib/ITensorMPS/examples/exact_diagonalization/fuse_inds.jl
@@ -1,4 +1,4 @@
-using ITensors
+using ITensors, ITensorMPS
 
 function fusion_tree(s::Vector{<:Index})
   n = length(s)
diff --git a/src/lib/ITensorMPS/examples/finite_temperature/metts.jl b/src/lib/ITensorMPS/examples/finite_temperature/metts.jl
index cfeead5e0a..ec53c5da5d 100644
--- a/src/lib/ITensorMPS/examples/finite_temperature/metts.jl
+++ b/src/lib/ITensorMPS/examples/finite_temperature/metts.jl
@@ -1,4 +1,4 @@
-using ITensors
+using ITensors, ITensorMPS
 using Printf
 
 #=
diff --git a/src/lib/ITensorMPS/examples/finite_temperature/purification.jl b/src/lib/ITensorMPS/examples/finite_temperature/purification.jl
index 4bba52b12b..fd374e47dd 100644
--- a/src/lib/ITensorMPS/examples/finite_temperature/purification.jl
+++ b/src/lib/ITensorMPS/examples/finite_temperature/purification.jl
@@ -1,4 +1,4 @@
-using ITensors
+using ITensors, ITensorMPS
 using Printf
 
 #=
diff --git a/src/lib/ITensorMPS/examples/gate_evolution/mpo_gate_evolution.jl b/src/lib/ITensorMPS/examples/gate_evolution/mpo_gate_evolution.jl
index 30ffb9af0a..04c16c7e14 100644
--- a/src/lib/ITensorMPS/examples/gate_evolution/mpo_gate_evolution.jl
+++ b/src/lib/ITensorMPS/examples/gate_evolution/mpo_gate_evolution.jl
@@ -1,8 +1,6 @@
-using ITensors
+using ITensors, ITensorMPS
 
-import ITensors: op
-
-function op(::OpName"expτSS", ::SiteType"S=1/2", s1::Index, s2::Index; τ)
+function ITensors.op(::OpName"expτSS", ::SiteType"S=1/2", s1::Index, s2::Index; τ)
   h =
     1 / 2 * op("S+", s1) * op("S-", s2) +
     1 / 2 * op("S-", s1) * op("S+", s2) +
diff --git a/src/lib/ITensorMPS/examples/gate_evolution/quantum_simulator.jl b/src/lib/ITensorMPS/examples/gate_evolution/quantum_simulator.jl
index 18aec97ee3..91e2d57c00 100644
--- a/src/lib/ITensorMPS/examples/gate_evolution/quantum_simulator.jl
+++ b/src/lib/ITensorMPS/examples/gate_evolution/quantum_simulator.jl
@@ -1,4 +1,4 @@
-using ITensors
+using ITensors, ITensorMPS
 
 N = 10
 s = siteinds("Qubit", N)
diff --git a/src/lib/ITensorMPS/examples/gate_evolution/trotter_suzuki_decomposition.jl b/src/lib/ITensorMPS/examples/gate_evolution/trotter_suzuki_decomposition.jl
index cd0c3ed4b1..6222470a91 100644
--- a/src/lib/ITensorMPS/examples/gate_evolution/trotter_suzuki_decomposition.jl
+++ b/src/lib/ITensorMPS/examples/gate_evolution/trotter_suzuki_decomposition.jl
@@ -1,4 +1,4 @@
-using ITensors
+using ITensors, ITensorMPS
 
 function heisenberg(N)
   os = Sum{Op}()
diff --git a/src/lib/ITensorMPS/examples/mps_mpo_algebra/mps_density_matrix.jl b/src/lib/ITensorMPS/examples/mps_mpo_algebra/mps_density_matrix.jl
index e1738ed23c..bea370842d 100644
--- a/src/lib/ITensorMPS/examples/mps_mpo_algebra/mps_density_matrix.jl
+++ b/src/lib/ITensorMPS/examples/mps_mpo_algebra/mps_density_matrix.jl
@@ -1,4 +1,4 @@
-using ITensors
+using ITensors, ITensorMPS
 
 N = 4
 nmps = 3
diff --git a/src/lib/ITensorMPS/examples/src/electronk.jl b/src/lib/ITensorMPS/examples/src/electronk.jl
index 0326aec23d..22a73f1e38 100644
--- a/src/lib/ITensorMPS/examples/src/electronk.jl
+++ b/src/lib/ITensorMPS/examples/src/electronk.jl
@@ -1,7 +1,6 @@
+using ITensors, ITensorMPS
 
-import ITensors: space, state, op!, has_fermion_string
-
-function space(
+function ITensors.space(
   ::SiteType"ElecK",
   n::Int;
   conserve_qns=false,
@@ -61,111 +60,115 @@ function space(
   return 4
 end
 
-state(::StateName"Emp", ::SiteType"ElecK") = [1.0 0.0 0.0 0.0]
-state(::StateName"Up", ::SiteType"ElecK") = [0.0 1.0 0.0 0.0]
-state(::StateName"Dn", ::SiteType"ElecK") = [0.0 0.0 1.0 0.0]
-state(::StateName"UpDn", ::SiteType"ElecK") = [0.0 0.0 0.0 1.0]
-state(::StateName"0", st::SiteType"ElecK") = state(StateName("Emp"), st)
-state(::StateName"↑", st::SiteType"ElecK") = state(StateName("Up"), st)
-state(::StateName"↓", st::SiteType"ElecK") = state(StateName("Dn"), st)
-state(::StateName"↑↓", st::SiteType"ElecK") = state(StateName("UpDn"), st)
+ITensors.state(::StateName"Emp", ::SiteType"ElecK") = [1.0 0.0 0.0 0.0]
+ITensors.state(::StateName"Up", ::SiteType"ElecK") = [0.0 1.0 0.0 0.0]
+ITensors.state(::StateName"Dn", ::SiteType"ElecK") = [0.0 0.0 1.0 0.0]
+ITensors.state(::StateName"UpDn", ::SiteType"ElecK") = [0.0 0.0 0.0 1.0]
+ITensors.state(::StateName"0", st::SiteType"ElecK") = state(StateName("Emp"), st)
+ITensors.state(::StateName"↑", st::SiteType"ElecK") = state(StateName("Up"), st)
+ITensors.state(::StateName"↓", st::SiteType"ElecK") = state(StateName("Dn"), st)
+ITensors.state(::StateName"↑↓", st::SiteType"ElecK") = state(StateName("UpDn"), st)
 
-function op!(Op::ITensor, ::OpName"Nup", ::SiteType"ElecK", s::Index)
+function ITensors.op!(Op::ITensor, ::OpName"Nup", ::SiteType"ElecK", s::Index)
   Op[s' => 2, s => 2] = 1.0
   return Op[s' => 4, s => 4] = 1.0
 end
 
-function op!(Op::ITensor, ::OpName"Ndn", ::SiteType"ElecK", s::Index)
+function ITensors.op!(Op::ITensor, ::OpName"Ndn", ::SiteType"ElecK", s::Index)
   Op[s' => 3, s => 3] = 1.0
   return Op[s' => 4, s => 4] = 1.0
 end
 
-function op!(Op::ITensor, ::OpName"Nupdn", ::SiteType"ElecK", s::Index)
+function ITensors.op!(Op::ITensor, ::OpName"Nupdn", ::SiteType"ElecK", s::Index)
   return Op[s' => 4, s => 4] = 1.0
 end
 
-function op!(Op::ITensor, ::OpName"Ntot", ::SiteType"ElecK", s::Index)
+function ITensors.op!(Op::ITensor, ::OpName"Ntot", ::SiteType"ElecK", s::Index)
   Op[s' => 2, s => 2] = 1.0
   Op[s' => 3, s => 3] = 1.0
   return Op[s' => 4, s => 4] = 2.0
 end
 
-function op!(Op::ITensor, ::OpName"Cup", ::SiteType"ElecK", s::Index)
+function ITensors.op!(Op::ITensor, ::OpName"Cup", ::SiteType"ElecK", s::Index)
   Op[s' => 1, s => 2] = 1.0
   return Op[s' => 3, s => 4] = 1.0
 end
 
-function op!(Op::ITensor, ::OpName"Cdagup", ::SiteType"ElecK", s::Index)
+function ITensors.op!(Op::ITensor, ::OpName"Cdagup", ::SiteType"ElecK", s::Index)
   Op[s' => 2, s => 1] = 1.0
   return Op[s' => 4, s => 3] = 1.0
 end
 
-function op!(Op::ITensor, ::OpName"Cdn", ::SiteType"ElecK", s::Index)
+function ITensors.op!(Op::ITensor, ::OpName"Cdn", ::SiteType"ElecK", s::Index)
   Op[s' => 1, s => 3] = 1.0
   return Op[s' => 2, s => 4] = -1.0
 end
 
-function op!(Op::ITensor, ::OpName"Cdagdn", ::SiteType"ElecK", s::Index)
+function ITensors.op!(Op::ITensor, ::OpName"Cdagdn", ::SiteType"ElecK", s::Index)
   Op[s' => 3, s => 1] = 1.0
   return Op[s' => 4, s => 2] = -1.0
 end
 
-function op!(Op::ITensor, ::OpName"Aup", ::SiteType"ElecK", s::Index)
+function ITensors.op!(Op::ITensor, ::OpName"Aup", ::SiteType"ElecK", s::Index)
   Op[s' => 1, s => 2] = 1.0
   return Op[s' => 3, s => 4] = 1.0
 end
 
-function op!(Op::ITensor, ::OpName"Adagup", ::SiteType"ElecK", s::Index)
+function ITensors.op!(Op::ITensor, ::OpName"Adagup", ::SiteType"ElecK", s::Index)
   Op[s' => 2, s => 1] = 1.0
   return Op[s' => 4, s => 3] = 1.0
 end
 
-function op!(Op::ITensor, ::OpName"Adn", ::SiteType"ElecK", s::Index)
+function ITensors.op!(Op::ITensor, ::OpName"Adn", ::SiteType"ElecK", s::Index)
   Op[s' => 1, s => 3] = 1.0
   return Op[s' => 2, s => 4] = 1.0
 end
 
-function op!(Op::ITensor, ::OpName"Adagdn", ::SiteType"ElecK", s::Index)
+function ITensors.op!(Op::ITensor, ::OpName"Adagdn", ::SiteType"ElecK", s::Index)
   Op[s' => 3, s => 1] = 1.0
   return Op[s' => 2, s => 4] = 1.0
 end
 
-function op!(Op::ITensor, ::OpName"F", ::SiteType"ElecK", s::Index)
+function ITensors.op!(Op::ITensor, ::OpName"F", ::SiteType"ElecK", s::Index)
   Op[s' => 1, s => 1] = +1.0
   Op[s' => 2, s => 2] = -1.0
   Op[s' => 3, s => 3] = -1.0
   return Op[s' => 4, s => 4] = +1.0
 end
 
-function op!(Op::ITensor, ::OpName"Fup", ::SiteType"ElecK", s::Index)
+function ITensors.op!(Op::ITensor, ::OpName"Fup", ::SiteType"ElecK", s::Index)
   Op[s' => 1, s => 1] = +1.0
   Op[s' => 2, s => 2] = -1.0
   Op[s' => 3, s => 3] = +1.0
   return Op[s' => 4, s => 4] = -1.0
 end
 
-function op!(Op::ITensor, ::OpName"Fdn", ::SiteType"ElecK", s::Index)
+function ITensors.op!(Op::ITensor, ::OpName"Fdn", ::SiteType"ElecK", s::Index)
   Op[s' => 1, s => 1] = +1.0
   Op[s' => 2, s => 2] = +1.0
   Op[s' => 3, s => 3] = -1.0
   return Op[s' => 4, s => 4] = -1.0
 end
 
-function op!(Op::ITensor, ::OpName"Sz", ::SiteType"ElecK", s::Index)
+function ITensors.op!(Op::ITensor, ::OpName"Sz", ::SiteType"ElecK", s::Index)
   Op[s' => 2, s => 2] = +0.5
   return Op[s' => 3, s => 3] = -0.5
 end
 
-op!(Op::ITensor, ::OpName"Sᶻ", st::SiteType"ElecK", s::Index) = op!(Op, OpName("Sz"), st, s)
+function ITensors.op!(Op::ITensor, ::OpName"Sᶻ", st::SiteType"ElecK", s::Index)
+  return op!(Op, OpName("Sz"), st, s)
+end
 
-function op!(Op::ITensor, ::OpName"Sx", ::SiteType"ElecK", s::Index)
+function ITensors.op!(Op::ITensor, ::OpName"Sx", ::SiteType"ElecK", s::Index)
   Op[s' => 2, s => 3] = 0.5
   return Op[s' => 3, s => 2] = 0.5
 end
 
-op!(Op::ITensor, ::OpName"Sˣ", st::SiteType"ElecK", s::Index) = op!(Op, OpName("Sx"), st, s)
+function ITensors.op!(Op::ITensor, ::OpName"Sˣ", st::SiteType"ElecK", s::Index)
+  return op!(Op, OpName("Sx"), st, s)
+end
 
-function op!(Op::ITensor, ::OpName"S+", ::SiteType"ElecK", s::Index)
+function ITensors.op!(Op::ITensor, ::OpName"S+", ::SiteType"ElecK", s::Index)
   return Op[s' => 2, s => 3] = 1.0
 end
 
@@ -179,13 +182,17 @@ function op!(Op::ITensor, ::OpName"S-", ::SiteType"ElecK", s::Index)
   return Op[s' => 3, s => 2] = 1.0
 end
 
-op!(Op::ITensor, ::OpName"S⁻", st::SiteType"ElecK", s::Index) = op!(Op, OpName("S-"), st, s)
-op!(Op::ITensor, ::OpName"Sm", st::SiteType"ElecK", s::Index) = op!(Op, OpName("S-"), st, s)
-function op!(Op::ITensor, ::OpName"Sminus", st::SiteType"ElecK", s::Index)
+function ITensors.op!(Op::ITensor, ::OpName"S⁻", st::SiteType"ElecK", s::Index)
+  return op!(Op, OpName("S-"), st, s)
+end
+function ITensors.op!(Op::ITensor, ::OpName"Sm", st::SiteType"ElecK", s::Index)
+  return op!(Op, OpName("S-"), st, s)
+end
+function ITensors.op!(Op::ITensor, ::OpName"Sminus", st::SiteType"ElecK", s::Index)
   return op!(Op, OpName("S-"), st, s)
 end
 
-has_fermion_string(::OpName"Cup", ::SiteType"ElecK") = true
-has_fermion_string(::OpName"Cdagup", ::SiteType"ElecK") = true
-has_fermion_string(::OpName"Cdn", ::SiteType"ElecK") = true
-has_fermion_string(::OpName"Cdagdn", ::SiteType"ElecK") = true
+ITensors.has_fermion_string(::OpName"Cup", ::SiteType"ElecK") = true
+ITensors.has_fermion_string(::OpName"Cdagup", ::SiteType"ElecK") = true
+ITensors.has_fermion_string(::OpName"Cdn", ::SiteType"ElecK") = true
+ITensors.has_fermion_string(::OpName"Cdagdn", ::SiteType"ElecK") = true
diff --git a/src/lib/ITensorMPS/src/abstractmps.jl b/src/lib/ITensorMPS/src/abstractmps.jl
index 140573b1bf..1e614dccd0 100644
--- a/src/lib/ITensorMPS/src/abstractmps.jl
+++ b/src/lib/ITensorMPS/src/abstractmps.jl
@@ -252,7 +252,7 @@ Use [`deepcopy`](@ref) for an alternative that copies the ITensors as well.
 
 # Examples
 ```julia
-julia> using ITensors
+julia> using ITensors, ITensorMPS
 
 julia> s = siteinds("S=1/2", 3);
 
@@ -300,7 +300,7 @@ copying the ITensor data.
 
 # Examples
 ```julia
-julia> using ITensors
+julia> using ITensors, ITensorMPS
 
 julia> s = siteinds("S=1/2", 3);
 
diff --git a/src/lib/ITensorMPS/src/mpo.jl b/src/lib/ITensorMPS/src/mpo.jl
index 654d3105bd..a829577457 100644
--- a/src/lib/ITensorMPS/src/mpo.jl
+++ b/src/lib/ITensorMPS/src/mpo.jl
@@ -36,7 +36,7 @@ function MPO(::Type{ElT}, sites::Vector{<:Index}) where {ElT<:Number}
   if N == 0
     return MPO()
   elseif N == 1
-    v[1] = emptyITensor(ElT, dag(sites[1]), sites[1]')
+    v[1] = ITensor(ElT, dag(sites[1]), sites[1]')
     return MPO(v)
   end
   space_ii = all(hasqns, sites) ? [QN() => 1] : 1
@@ -44,11 +44,11 @@ function MPO(::Type{ElT}, sites::Vector{<:Index}) where {ElT<:Number}
   for ii in eachindex(sites)
     s = sites[ii]
     if ii == 1
-      v[ii] = emptyITensor(ElT, dag(s), s', l[ii])
+      v[ii] = ITensor(ElT, dag(s), s', l[ii])
     elseif ii == N
-      v[ii] = emptyITensor(ElT, dag(l[ii - 1]), dag(s), s')
+      v[ii] = ITensor(ElT, dag(l[ii - 1]), dag(s), s')
     else
-      v[ii] = emptyITensor(ElT, dag(l[ii - 1]), dag(s), s', l[ii])
+      v[ii] = ITensor(ElT, dag(l[ii - 1]), dag(s), s', l[ii])
     end
   end
   return MPO(v)
diff --git a/src/lib/ITensorMPS/src/opsum_to_mpo/opsum_to_mpo_generic.jl b/src/lib/ITensorMPS/src/opsum_to_mpo/opsum_to_mpo_generic.jl
index aea3a3e47c..e21806b4aa 100644
--- a/src/lib/ITensorMPS/src/opsum_to_mpo/opsum_to_mpo_generic.jl
+++ b/src/lib/ITensorMPS/src/opsum_to_mpo/opsum_to_mpo_generic.jl
@@ -204,7 +204,8 @@ function sorteachterm(os::OpSum, sites)
       end
     end
 
-    (t_parity == -1) && error("Parity-odd fermionic terms not yet supported by AutoMPO")
+    (t_parity == -1) &&
+      error("Parity-odd fermionic terms not yet supported by OpSum to MPO conversion")
 
     # Keep only fermionic op positions (non-zero entries)
     filter!(!iszero, perm)
diff --git a/src/lib/ITensorMPS/test/base/test_mps.jl b/src/lib/ITensorMPS/test/base/test_mps.jl
index 110da45908..bfe39da013 100644
--- a/src/lib/ITensorMPS/test/base/test_mps.jl
+++ b/src/lib/ITensorMPS/test/base/test_mps.jl
@@ -973,7 +973,7 @@ end
         ("S+", "S-"),
       ],
     )
-    # can't test ,("Cdn","Cdn") yet, because AutoMPO thinks this is antisymmetric
+    # can't test ,("Cdn","Cdn") yet, because OpSum to MPO thinks this is antisymmetric
 
     #trigger unsupported error
     let err = nothing
diff --git a/src/lib/ITensorMPS/test/base/test_qnmpo.jl b/src/lib/ITensorMPS/test/base/test_qnmpo.jl
index 8f72e6577e..1f417d5db8 100644
--- a/src/lib/ITensorMPS/test/base/test_qnmpo.jl
+++ b/src/lib/ITensorMPS/test/base/test_qnmpo.jl
@@ -296,7 +296,7 @@ end
 #
 #  Build up Hamiltonians with non trival QN spaces in the link indices and further neighbour interactions.
 #
-function make_Heisenberg_AutoMPO(sites, NNN::Int64; J::Float64=1.0, kwargs...)::MPO
+function make_heisenberg_opsum(sites, NNN::Int64; J::Float64=1.0, kwargs...)::MPO
   N = length(sites)
   @assert N >= NNN
   opsum = OpSum()
@@ -311,7 +311,7 @@ function make_Heisenberg_AutoMPO(sites, NNN::Int64; J::Float64=1.0, kwargs...)::
   return MPO(opsum, sites; kwargs...)
 end
 
-function make_Hubbard_AutoMPO(
+function make_hubbard_opsum(
   sites, NNN::Int64; U::Float64=1.0, t::Float64=1.0, V::Float64=0.5, kwargs...
 )::MPO
   N = length(sites)
@@ -333,7 +333,7 @@ function make_Hubbard_AutoMPO(
   return MPO(os, sites; kwargs...)
 end
 
-test_combos = [(make_Heisenberg_AutoMPO, "S=1/2"), (make_Hubbard_AutoMPO, "Electron")]
+test_combos = [(make_heisenberg_opsum, "S=1/2"), (make_hubbard_opsum, "Electron")]
 
 @testset "QR/QL MPO tensors with complex block structures, H=$(test_combo[1])" for test_combo in
                                                                                    test_combos
diff --git a/src/lib/LazyApply/src/LazyApply.jl b/src/lib/LazyApply/src/LazyApply.jl
index 77d4aca34b..f0e69fa71f 100644
--- a/src/lib/LazyApply/src/LazyApply.jl
+++ b/src/lib/LazyApply/src/LazyApply.jl
@@ -306,7 +306,7 @@ lastindex(a::Scaled{C,<:Sum}) where {C} = lastindex(argument(a))
 lastindex(a::Scaled{C,<:Prod}) where {C} = lastindex(argument(a))
 
 #
-# Functions convenient for AutoMPO code
+# Functions convenient for OpSum code
 #
 
 terms(a::Union{Sum,Prod}) = only(a.args)
diff --git a/test/Project.toml b/test/Project.toml
index 3d7da6330d..34fe930acb 100644
--- a/test/Project.toml
+++ b/test/Project.toml
@@ -5,6 +5,7 @@ Combinatorics = "861a8166-3701-5b0c-9a16-15d98fcdc6aa"
 Compat = "34da2185-b29b-5c13-b0c7-acf172513d20"
 FiniteDifferences = "26cc04aa-876d-5657-8c51-4c34ba976000"
 HDF5 = "f67ccb44-e63f-5c2f-98bd-6dc0ccc4ba2f"
+ITensorMPS = "0d1a4710-d33b-49a5-8f18-73bdf49b47e2"
 ITensors = "9136182c-28ba-11e9-034c-db9fb085ebd5"
 JLD2 = "033835bb-8acc-5ee8-8aae-3f567f8a3819"
 KrylovKit = "0b1a1467-8014-51b9-945f-bf0ae24f4b77"
diff --git a/test/base/test_indices.jl b/test/base/test_indices.jl
index a15eca9f56..ebde10a6b0 100644
--- a/test/base/test_indices.jl
+++ b/test/base/test_indices.jl
@@ -37,12 +37,12 @@ using ITensors.NDTensors
   @test hassameinds(ITensor(undef, is1), is)
   @test hassameinds(ITensor(undef, is2), is)
   @test hassameinds(ITensor(undef, is1...), is)
-  @test hassameinds(emptyITensor(is1), is)
-  @test hassameinds(emptyITensor(is2), is)
-  @test hassameinds(emptyITensor(is1...), is)
-  @test hassameinds(emptyITensor(Float64, is1), is)
-  @test hassameinds(emptyITensor(Float64, is2), is)
-  @test hassameinds(emptyITensor(Float64, is1...), is)
+  @test hassameinds(ITensor(is1), is)
+  @test hassameinds(ITensor(is2), is)
+  @test hassameinds(ITensor(is1...), is)
+  @test hassameinds(ITensor(Float64, is1), is)
+  @test hassameinds(ITensor(Float64, is2), is)
+  @test hassameinds(ITensor(Float64, is1...), is)
   @test hassameinds(ITensor(A, is1), is)
   @test hassameinds(ITensor(A, is2), is)
   @test hassameinds(ITensor(A, is1...), is)
@@ -125,12 +125,12 @@ end
   @test hassameinds(ITensor(undef, is1), is)
   @test hassameinds(ITensor(undef, is2), is)
   @test hassameinds(ITensor(undef, is1...), is)
-  @test hassameinds(emptyITensor(is1), is)
-  @test hassameinds(emptyITensor(is2), is)
-  @test hassameinds(emptyITensor(is1...), is)
-  @test hassameinds(emptyITensor(Float64, is1), is)
-  @test hassameinds(emptyITensor(Float64, is2), is)
-  @test hassameinds(emptyITensor(Float64, is1...), is)
+  @test hassameinds(ITensor(is1), is)
+  @test hassameinds(ITensor(is2), is)
+  @test hassameinds(ITensor(is1...), is)
+  @test hassameinds(ITensor(Float64, is1), is)
+  @test hassameinds(ITensor(Float64, is2), is)
+  @test hassameinds(ITensor(Float64, is1...), is)
   @test hassameinds(ITensor(A, is1), is)
   @test hassameinds(ITensor(A, is2), is)
   @test hassameinds(ITensor(A, is1...), is)
diff --git a/test/base/test_itensor_scalar_contract.jl b/test/base/test_itensor_scalar_contract.jl
index c322139bba..c91e07dc1d 100644
--- a/test/base/test_itensor_scalar_contract.jl
+++ b/test/base/test_itensor_scalar_contract.jl
@@ -18,11 +18,11 @@ Random.seed!(1234)
   C = A * B
   @test C ≈ B[1, 1, 1] * A * ITensor(1, inds(B))
 
-  C = emptyITensor(is..., α', α'')
+  C = ITensor(is..., α', α'')
   C .= A .* B
   @test C ≈ B[1, 1, 1] * A * ITensor(1, inds(B))
 
-  C = emptyITensor(shuffle([(is..., α', α'')...])...)
+  C = ITensor(shuffle([(is..., α', α'')...])...)
   C .= A .* B
   @test C ≈ B[1, 1, 1] * A * ITensor(1, inds(B))
 end
diff --git a/test/base/test_qndiagitensor.jl b/test/base/test_qndiagitensor.jl
index 178b71175b..c2ec0c2075 100644
--- a/test/base/test_qndiagitensor.jl
+++ b/test/base/test_qndiagitensor.jl
@@ -83,7 +83,7 @@ using ITensors, Test
     r = addtags(s, "right")
     u = addtags(s, "up")
     d = dag(addtags(s, "down"))
-    A = emptyITensor(l, r, u, d)
+    A = ITensor(l, r, u, d)
     A[1, 1, 1, 1] = 1.0
     A[1, 1, 2, 2] = 1.0
     A[2, 2, 1, 1] = 1.0
diff --git a/test/base/test_qnitensor.jl b/test/base/test_qnitensor.jl
index 1b359afb06..61ebcf5c1d 100644
--- a/test/base/test_qnitensor.jl
+++ b/test/base/test_qnitensor.jl
@@ -292,8 +292,8 @@ Random.seed!(1234)
     @test T[] == 0
 
     s = Index(QN(-1) => 1, QN(1) => 1)
-    A = emptyITensor(s, dag(s'))
-    B = emptyITensor(s', dag(s))
+    A = ITensor(s, dag(s'))
+    B = ITensor(s', dag(s))
     A[1, 1] = 1
     B[2, 2] = 1
     C = A * B
@@ -306,7 +306,7 @@ Random.seed!(1234)
   @testset "Empty constructor" begin
     i = Index([QN(0) => 1, QN(1) => 2], "i")
 
-    A = emptyITensor(i, dag(i'))
+    A = ITensor(i, dag(i'))
 
     @test nnzblocks(A) == 0
     @test nnz(A) == 0
@@ -426,7 +426,7 @@ Random.seed!(1234)
     @testset "Test 1" begin
       s1 = Index([QN("N", 0, -1) => 1, QN("N", 1, -1) => 1], "s1")
       s2 = Index([QN("N", 0, -1) => 1, QN("N", 1, -1) => 1], "s2")
-      A = emptyITensor(s1, s2)
+      A = ITensor(s1, s2)
 
       @test nnzblocks(A) == 0
       @test nnz(A) == 0
@@ -452,7 +452,7 @@ Random.seed!(1234)
     @testset "Test 2" begin
       s1 = Index([QN("N", 0, -1) => 1, QN("N", 1, -1) => 1], "s1")
       s2 = Index([QN("N", 0, -1) => 1, QN("N", 1, -1) => 1], "s2")
-      A = emptyITensor(s1, s2)
+      A = ITensor(s1, s2)
 
       @test nnzblocks(A) == 0
       @test nnz(A) == 0
@@ -872,7 +872,7 @@ Random.seed!(1234)
   @testset "Combiner for block deficient ITensor" begin
     i = Index(QN(0, 2) => 2, QN(1, 2) => 2; tags="i")
     j = settags(i, "j")
-    A = emptyITensor(i, j, dag(i'))
+    A = ITensor(i, j, dag(i'))
     A[1, 1, 1] = 1.0
     C = combiner(i, j; tags="c")
     AC = A * C
@@ -1529,7 +1529,7 @@ Random.seed!(1234)
         QN("Nf", 1, -1) => 2, QN("Nf", 0, -1) => 1, QN("Nf", 1, -1) => 2; tags="Link,u"
       )
 
-      A = emptyITensor(ElT, l, s, dag(r))
+      A = ITensor(ElT, l, s, dag(r))
 
       insertblock!(A, Block(2, 1, 2))
       insertblock!(A, Block(1, 2, 2))
@@ -1561,7 +1561,7 @@ Random.seed!(1234)
         QN("Sz", 2) => 4,
         QN("Sz", 4) => 1,
       )
-      A = emptyITensor(ElT, s, s')
+      A = ITensor(ElT, s, s')
       insertblock!(A, Block(5, 2))
       insertblock!(A, Block(4, 3))
       insertblock!(A, Block(3, 4))
@@ -1579,7 +1579,7 @@ Random.seed!(1234)
         QN("Sz", 2) => 4,
         QN("Sz", 4) => 1,
       )
-      A = emptyITensor(ElT, s, s')
+      A = ITensor(ElT, s, s')
       insertblock!(A, Block(5, 1))
       insertblock!(A, Block(4, 2))
       insertblock!(A, Block(3, 3))
@@ -1598,7 +1598,7 @@ Random.seed!(1234)
         QN("Sz", 2) => 4,
         QN("Sz", 4) => 1,
       )
-      A = emptyITensor(ElT, s, s')
+      A = ITensor(ElT, s, s')
       insertblock!(A, Block(5, 1))
       insertblock!(A, Block(4, 2))
       insertblock!(A, Block(3, 3))
diff --git a/test/base/test_sitetype.jl b/test/base/test_sitetype.jl
index aa2d957bb3..9dfd794c72 100644
--- a/test/base/test_sitetype.jl
+++ b/test/base/test_sitetype.jl
@@ -95,7 +95,7 @@ end
     # is for S=3/2, because we might define the
     # "S=3/2" TagType inside ITensors.jl later
     function ITensors.op(::OpName"Sz", ::SiteType"_Custom_", s::Index)
-      Op = emptyITensor(s', dag(s))
+      Op = ITensor(s', dag(s))
       Op[s' => 1, s => 1] = +3 / 2
       Op[s' => 2, s => 2] = +1 / 2
       Op[s' => 3, s => 3] = -1 / 2
@@ -104,7 +104,7 @@ end
     end
 
     function ITensors.op(::OpName"α", ::SiteType"_Custom_", s1::Index, s2::Index)
-      Op = emptyITensor(s1', s2', dag(s1), dag(s2))
+      Op = ITensor(s1', s2', dag(s1), dag(s2))
       Op[s1' => 1, s2' => 2, s1 => 1, s2 => 2] = +3 / 2
       Op[s1' => 2, s2' => 1, s1 => 2, s2 => 2] = +1 / 2
       Op[s1' => 3, s2' => 3, s1 => 3, s2 => 4] = -1 / 2
@@ -115,7 +115,7 @@ end
     function ITensors.op(
       ::OpName"β", ::SiteType"_Custom1", ::SiteType"_Custom2", s1::Index, s2::Index
     )
-      Op = emptyITensor(s1', s2', dag(s1), dag(s2))
+      Op = ITensor(s1', s2', dag(s1), dag(s2))
       Op[s1' => 1, s2' => 2, s1 => 1, s2 => 2] = +5 / 2
       Op[s1' => 2, s2' => 1, s1 => 2, s2 => 2] = +3 / 2
       Op[s1' => 3, s2' => 3, s1 => 3, s2 => 4] = -3 / 2
@@ -151,7 +151,7 @@ end
 
   @testset "Custom OpName with long name" begin
     function ITensors.op(::OpName"my_favorite_operator", ::SiteType"S=1/2", s::Index)
-      Op = emptyITensor(s', dag(s))
+      Op = ITensor(s', dag(s))
       Op[s' => 1, s => 1] = 0.11
       Op[s' => 1, s => 2] = 0.12
       Op[s' => 2, s => 1] = 0.21
@@ -291,7 +291,7 @@ end
     # is for S=3/2, because we might define the
     # "S=3/2" TagType inside ITensors.jl later
     function ITensors.op(::SiteType"_Custom_", s::Index, opname::AbstractString)
-      Op = emptyITensor(s', dag(s))
+      Op = ITensor(s', dag(s))
       if opname == "S+"
         Op[s' => 1, s => 2] = sqrt(3)
         Op[s' => 2, s => 3] = 2
diff --git a/test/threading/test_threading.jl b/test/threading/test_threading.jl
index 3add3d3e75..d63f24c2cd 100644
--- a/test/threading/test_threading.jl
+++ b/test/threading/test_threading.jl
@@ -60,8 +60,8 @@ end
 
   @testset "Contraction resulting in no blocks with threading bug" begin
     i = Index([QN(0) => 1, QN(1) => 1])
-    A = emptyITensor(i', dag(i))
-    B = emptyITensor(i', dag(i))
+    A = ITensor(i', dag(i))
+    B = ITensor(i', dag(i))
     A[i' => 1, i => 1] = 11.0
     B[i' => 2, i => 2] = 22.0