diff --git a/Project.toml b/Project.toml
index 1e7f896..ef37d3b 100644
--- a/Project.toml
+++ b/Project.toml
@@ -19,8 +19,8 @@ SplitApplyCombine = "03a91e81-4c3e-53e1-a0a4-9c0c8f19dd66"
 [compat]
 Compat = "3, 4"
 HDF5 = "0.15, 0.16, 0.17"
-ITensorMPS = "0.1, 0.2"
-ITensors = "0.3.48, 0.4, 0.5, 0.6"
+ITensorMPS = "0.2.2"
+ITensors = "0.6.8"
 Infinities = "0.1"
 IterTools = "1"
 KrylovKit = "0.5, 0.6, 0.7"
diff --git a/examples/development/finite_mps_to_infinite_mps.jl b/examples/development/finite_mps_to_infinite_mps.jl
index 4fbea42..6dd1e37 100644
--- a/examples/development/finite_mps_to_infinite_mps.jl
+++ b/examples/development/finite_mps_to_infinite_mps.jl
@@ -21,7 +21,7 @@ function main(; N, J, h, nsites)
   s = siteinds("S=1/2", N)
 
   H = MPO(ising_opsum_finite(N; J=J, h=h), s)
-  ψ0 = randomMPS(s)
+  ψ0 = random_mps(s)
 
   energy, ψ = dmrg(H, ψ0; nsweeps=10, cutoff=1e-10)
   @show energy / N
diff --git a/examples/development/orthogonalize_infinitemps.jl b/examples/development/orthogonalize_infinitemps.jl
index 101ce85..b3f64ef 100644
--- a/examples/development/orthogonalize_infinitemps.jl
+++ b/examples/development/orthogonalize_infinitemps.jl
@@ -12,7 +12,7 @@ s = siteinds("S=1/2", nsites; conserve_szparity=true)
 space = (("SzParity", 1, 2) => χ ÷ 2) ⊕ (("SzParity", 0, 2) => χ ÷ 2)
 ψ = InfiniteMPS(ComplexF64, s; space=space)
 for n in 1:nsites
-  ψ[n] = randomITensor(inds(ψ[n]))
+  ψ[n] = random_itensor(inds(ψ[n]))
 end
 
 ψ = orthogonalize(ψ, :)
diff --git a/examples/development/vumps/transfer_matrix_spectrum_arpack.jl b/examples/development/vumps/transfer_matrix_spectrum_arpack.jl
index b4694b5..90f1539 100644
--- a/examples/development/vumps/transfer_matrix_spectrum_arpack.jl
+++ b/examples/development/vumps/transfer_matrix_spectrum_arpack.jl
@@ -62,7 +62,7 @@ end
 Nfinite = 100
 sfinite = siteinds("S=1/2", Nfinite; conserve_szparity=true)
 Hfinite = MPO(model, sfinite; model_params...)
-ψfinite = randomMPS(sfinite, initstate)
+ψfinite = random_mps(sfinite, initstate)
 @show flux(ψfinite)
 sweeps = Sweeps(10)
 setmaxdim!(sweeps, maxdim)
@@ -111,8 +111,8 @@ using LinearAlgebra
 
 T = TransferMatrix(ψ.AL)
 Tᵀ = transpose(T)
-vⁱᴿ = randomITensor(dag(input_inds(T)))
-vⁱᴸ = randomITensor(dag(input_inds(Tᵀ)))
+vⁱᴿ = random_itensor(dag(input_inds(T)))
+vⁱᴸ = random_itensor(dag(input_inds(Tᵀ)))
 
 neigs = 10
 tol = 1e-10
diff --git a/examples/development/vumps/vumps_from_dmrg.jl b/examples/development/vumps/vumps_from_dmrg.jl
index cbf42e9..f2e249f 100644
--- a/examples/development/vumps/vumps_from_dmrg.jl
+++ b/examples/development/vumps/vumps_from_dmrg.jl
@@ -11,7 +11,7 @@ J = 1.0
 h = 1.5
 model = Model("ising")
 H = MPO(model, s; J=J, h=h)
-ψ0 = randomMPS(s)
+ψ0 = random_mps(s)
 
 sweeps = Sweeps(10)
 maxdim!(sweeps, 10)
diff --git a/examples/development/vumps/vumps_ising_noncontiguous.jl b/examples/development/vumps/vumps_ising_noncontiguous.jl
index a7787d5..38575a3 100644
--- a/examples/development/vumps/vumps_ising_noncontiguous.jl
+++ b/examples/development/vumps/vumps_ising_noncontiguous.jl
@@ -52,7 +52,7 @@ energy_infinite = expect(ψ_0, H)
 Nfinite = 100
 sfinite = siteinds("S=1/2", Nfinite; conserve_szparity=true)
 Hfinite = MPO(model, sfinite; J=J, h=h)
-ψfinite = randomMPS(sfinite, initstate; linkdims=10)
+ψfinite = random_mps(sfinite, initstate; linkdims=10)
 @show flux(ψfinite)
 sweeps = Sweeps(15)
 setmaxdim!(sweeps, maxdim)
diff --git a/examples/development/vumps/vumps_localham_qns.jl b/examples/development/vumps/vumps_localham_qns.jl
index 5c5c434..5772fa9 100644
--- a/examples/development/vumps/vumps_localham_qns.jl
+++ b/examples/development/vumps/vumps_localham_qns.jl
@@ -23,7 +23,7 @@ initstate(n) = isodd(n) ? "↑" : "↓"
 Nfinite = 100
 sfinite = siteinds("S=1/2", Nfinite; conserve_szparity=true)
 Hfinite = MPO(model, sfinite; model_kwargs...)
-ψfinite = randomMPS(sfinite, initstate)
+ψfinite = random_mps(sfinite, initstate)
 sweeps = Sweeps(20)
 setmaxdim!(sweeps, 1)
 setcutoff!(sweeps, 1E-10)
diff --git a/examples/vumps/transfer_matrix_spectrum.jl b/examples/vumps/transfer_matrix_spectrum.jl
index 881da74..560c443 100644
--- a/examples/vumps/transfer_matrix_spectrum.jl
+++ b/examples/vumps/transfer_matrix_spectrum.jl
@@ -63,7 +63,7 @@ subspace_expansion_kwargs = (cutoff=cutoff, maxdim=maxdim)
 Nfinite = 100
 sfinite = siteinds("S=1/2", Nfinite; conserve_szparity=conserve_qns)
 Hfinite = MPO(model, sfinite; model_params...)
-ψfinite = randomMPS(sfinite, initstate)
+ψfinite = random_mps(sfinite, initstate)
 @show flux(ψfinite)
 sweeps = Sweeps(10)
 setmaxdim!(sweeps, maxdim)
@@ -111,8 +111,8 @@ using LinearAlgebra
 
 T = TransferMatrix(ψ.AL)
 Tᵀ = transpose(T)
-vⁱᴿ = randomITensor(dag(input_inds(T)))
-vⁱᴸ = randomITensor(dag(input_inds(Tᵀ)))
+vⁱᴿ = random_itensor(dag(input_inds(T)))
+vⁱᴸ = random_itensor(dag(input_inds(Tᵀ)))
 
 neigs = 10
 tol = 1e-10
diff --git a/examples/vumps/vumps_heisenberg.jl b/examples/vumps/vumps_heisenberg.jl
index d9da157..c06dc00 100644
--- a/examples/vumps/vumps_heisenberg.jl
+++ b/examples/vumps/vumps_heisenberg.jl
@@ -70,7 +70,7 @@ energy_exact = reference(model, Observable("energy"))
 Nfinite = 100
 sfinite = siteinds("S=1/2", Nfinite; conserve_qns)
 Hfinite = MPO(model, sfinite)
-ψfinite = randomMPS(sfinite, initstate; linkdims=10)
+ψfinite = random_mps(sfinite, initstate; linkdims=10)
 @show flux(ψfinite)
 
 nsweeps = 10
diff --git a/examples/vumps/vumps_hubbard_extended.jl b/examples/vumps/vumps_hubbard_extended.jl
index b28e6bc..1384053 100644
--- a/examples/vumps/vumps_hubbard_extended.jl
+++ b/examples/vumps/vumps_hubbard_extended.jl
@@ -92,7 +92,7 @@ energy_infinite = map(b -> expect_two_site(ψ, H[b], b), bs)
 Nfinite = 100
 sfinite = siteinds("Electron", Nfinite; conserve_qns)
 Hfinite = MPO(model, sfinite; model_params...)
-ψfinite = randomMPS(sfinite, initstate; linkdims=10)
+ψfinite = random_mps(sfinite, initstate; linkdims=10)
 println("\nQN sector of starting finite MPS")
 @show flux(ψfinite)
 
diff --git a/examples/vumps/vumps_ising.jl b/examples/vumps/vumps_ising.jl
index 6af5ca9..67d4bc8 100644
--- a/examples/vumps/vumps_ising.jl
+++ b/examples/vumps/vumps_ising.jl
@@ -62,7 +62,7 @@ subspace_expansion_kwargs = (cutoff=cutoff, maxdim=maxdim)
 nsite_finite = 100
 s_finite = siteinds("S=1/2", nsite_finite; conserve_szparity=conserve_qns)
 H_finite = MPO(model, s_finite; model_params...)
-ψ_finite = randomMPS(s_finite, initstate)
+ψ_finite = random_mps(s_finite, initstate)
 @show flux(ψ_finite)
 sweeps = Sweeps(10)
 setmaxdim!(sweeps, maxdim)
diff --git a/examples/vumps/vumps_ising_extended.jl b/examples/vumps/vumps_ising_extended.jl
index cafd12a..1e9e968 100644
--- a/examples/vumps/vumps_ising_extended.jl
+++ b/examples/vumps/vumps_ising_extended.jl
@@ -47,7 +47,7 @@ energy_infinite = expect(ψ, H)
 Nfinite = 100
 sfinite = siteinds("S=1/2", Nfinite; conserve_szparity=conserve_qns)
 Hfinite = MPO(model, sfinite; J=J, J₂=J₂, h=h)
-ψfinite = randomMPS(sfinite, initstate; linkdims=10)
+ψfinite = random_mps(sfinite, initstate; linkdims=10)
 @show flux(ψfinite)
 dmrg_kwargs = (nsweeps=10, maxdim, cutoff)
 energy_finite_total, ψfinite = dmrg(Hfinite, ψfinite; dmrg_kwargs...)
diff --git a/src/infinitemps_approx.jl b/src/infinitemps_approx.jl
index 8009d3d..be692fb 100644
--- a/src/infinitemps_approx.jl
+++ b/src/infinitemps_approx.jl
@@ -39,7 +39,7 @@ function infinitemps_approx(
   # XXX: Why doesn't this work?
   # randn!.(A)
   for n in 1:N
-    A[n] = randomITensor(inds(A[n]))
+    A[n] = random_itensor(inds(A[n]))
   end
   ψ∞ = orthogonalize(A, :)
 
@@ -62,14 +62,14 @@ function infinitemps_approx(
       TR = TR * ψ[n] * ψ∞.AR[inf_range[n]]
     end
 
-    L0 = randomITensor(linkind(ψ∞.AL, first(inf_range) - 1 => first(inf_range)))
+    L0 = random_itensor(linkind(ψ∞.AL, first(inf_range) - 1 => first(inf_range)))
     vals, vecs, _ = eigsolve(transpose(ITensorMap(ψ∞.AL[inf_range], ψ)), L0)
     L0 = vecs[1]
     λL0 = vals[1]
     LTL = L0 * TL
     LTL = replaceinds(LTL, inds(LTL) => inds(L0))
 
-    RN = randomITensor(linkind(ψ∞.AR, last(inf_range) => last(inf_range) + 1))
+    RN = random_itensor(linkind(ψ∞.AR, last(inf_range) => last(inf_range) + 1))
     vals, vecs, _ = eigsolve(ITensorMap(ψ∞.AR[inf_range], ψ), RN)
     RN = vecs[1]
     λRN = vals[1]
diff --git a/src/orthogonalize.jl b/src/orthogonalize.jl
index 5d1ba22..4ed000a 100644
--- a/src/orthogonalize.jl
+++ b/src/orthogonalize.jl
@@ -7,7 +7,7 @@ function right_orthogonalize(
   T = TransferMatrix(ψ)
 
   # TODO: make an optional initial state
-  v₁ᴿᴺ = randomITensor(dag(input_inds(T)))
+  v₁ᴿᴺ = random_itensor(dag(input_inds(T)))
 
   # Start by getting the right eivenvector/eigenvalue of T
   # TODO: make a function `right_environments(::InfiniteMPS)` that computes
diff --git a/test/test_itensormap.jl b/test/test_itensormap.jl
index ee3b0d3..c0308ed 100644
--- a/test/test_itensormap.jl
+++ b/test/test_itensormap.jl
@@ -5,16 +5,16 @@ using Random
 
 @testset "ITensorMap" begin
   i = Index(2)
-  A1 = randomITensor(i'', dag(i))
-  A2 = randomITensor(i', dag(i''))
-  B1 = randomITensor(i'', dag(i))
-  B2 = randomITensor(i', dag(i''))
+  A1 = random_itensor(i'', dag(i))
+  A2 = random_itensor(i', dag(i''))
+  B1 = random_itensor(i'', dag(i))
+  B2 = random_itensor(i', dag(i''))
   c1 = 2.3
   c2 = 3.4
 
   MA = ITensorMap([A1, A2])
   MB = ITensorMap([B1, B2])
-  v = randomITensor(i)
+  v = random_itensor(i)
 
   @test MA(v) ≈ noprime(A2 * A1 * v)
   @test MB(v) ≈ noprime(B2 * B1 * v)
diff --git a/test/test_vumps.jl b/test/test_vumps.jl
index ad16e49..9921697 100644
--- a/test/test_vumps.jl
+++ b/test/test_vumps.jl
@@ -33,7 +33,7 @@ end
   Nfinite = 100
   sfinite = siteinds("S=1/2", Nfinite; conserve_szparity=true)
   Hfinite = MPO(model, sfinite; model_kwargs...)
-  ψfinite = randomMPS(sfinite, initstate)
+  ψfinite = random_mps(sfinite, initstate)
   nsweeps = 20
   energy_finite_total, ψfinite = dmrg(
     Hfinite, ψfinite; nsweeps, maxdim=10, cutoff=1e-10, outputlevel=0
@@ -139,7 +139,7 @@ end
   Nfinite = 100
   sfinite = siteinds("S=1/2", Nfinite; conserve_szparity=true)
   Hfinite = MPO(model, sfinite; model_kwargs...)
-  ψfinite = randomMPS(sfinite, initstate)
+  ψfinite = random_mps(sfinite, initstate)
   nsweeps = 20
   energy_finite_total, ψfinite = dmrg(
     Hfinite, ψfinite; nsweeps, maxdim=10, cutoff=1e-10, outputlevel=0
diff --git a/test/test_vumps_extendedising.jl b/test/test_vumps_extendedising.jl
index a5ac2b1..e6a8de1 100644
--- a/test/test_vumps_extendedising.jl
+++ b/test/test_vumps_extendedising.jl
@@ -19,7 +19,7 @@ end
   Nfinite = 100
   sfinite = siteinds("S=1/2", Nfinite; conserve_szparity=true)
   Hfinite = MPO(model, sfinite; model_kwargs...)
-  ψfinite = randomMPS(sfinite, initstate)
+  ψfinite = random_mps(sfinite, initstate)
   energy_finite_total, ψfinite = dmrg(
     Hfinite, ψfinite; outputlevel=0, nsweeps=30, maxdim=30, cutoff=1e-10
   )
@@ -89,7 +89,7 @@ end
   space = (("SzParity", 1, 2) => χ ÷ 2) ⊕ (("SzParity", 0, 2) => χ ÷ 2)
   ψ = InfiniteMPS(ComplexF64, s; space=space)
   for n in 1:N
-    ψ[n] = randomITensor(inds(ψ[n]))
+    ψ[n] = random_itensor(inds(ψ[n]))
   end
   ψ = orthogonalize(ψ, :)
   @test contract(ψ.AL[1:N]) * ψ.C[N] ≈ ψ.C[0] * contract(ψ.AR[1:N])
diff --git a/test/test_vumpsmpo.jl b/test/test_vumpsmpo.jl
index 83583a9..817103a 100644
--- a/test/test_vumpsmpo.jl
+++ b/test/test_vumpsmpo.jl
@@ -28,7 +28,7 @@ end
   Nfinite = 100
   sfinite = siteinds("S=1/2", Nfinite; conserve_szparity=true)
   Hfinite = MPO(model, sfinite; model_kwargs...)
-  ψfinite = randomMPS(sfinite, initstate)
+  ψfinite = random_mps(sfinite, initstate)
 
   sweeps = Sweeps(20)
   setmaxdim!(sweeps, 10)
@@ -102,7 +102,7 @@ end
   Nfinite = 100
   sfinite = siteinds("S=1/2", Nfinite; conserve_szparity=true)
   Hfinite = MPO(model, sfinite; model_kwargs...)
-  ψfinite = randomMPS(sfinite, initstate)
+  ψfinite = random_mps(sfinite, initstate)
   nsweeps = 20
   energy_finite_total, ψfinite = dmrg(
     Hfinite, ψfinite; nsweeps, maxdim=10, cutoff=1e-10, outputlevel=0
@@ -173,7 +173,7 @@ end
   Nfinite = 100
   sfinite = siteinds("S=1/2", Nfinite; conserve_szparity=true)
   Hfinite = MPO(model, sfinite; model_kwargs...)
-  ψfinite = randomMPS(sfinite, initstate)
+  ψfinite = random_mps(sfinite, initstate)
   sweeps = Sweeps(20)
   setmaxdim!(sweeps, 10)
   setcutoff!(sweeps, 1E-10)