diff --git a/ITensorGPU/examples/dmrg.jl b/ITensorGPU/examples/dmrg.jl
index 95d2d43a04..e3448ccace 100644
--- a/ITensorGPU/examples/dmrg.jl
+++ b/ITensorGPU/examples/dmrg.jl
@@ -2,25 +2,23 @@ using ITensors
 using ITensorGPU
 
 # Set to identity to run on CPU
-gpu = cu
+device = cu
 
 N = 50
 sites = siteinds("S=1", N)
 
-ampo = AutoMPO()
+opsum = OpSum()
 for j in 1:(N - 1)
-  ampo .+= 0.5, "S+", j, "S-", j + 1
-  ampo .+= 0.5, "S-", j, "S+", j + 1
-  ampo .+= "Sz", j, "Sz", j + 1
+  opsum .+= 0.5, "S+", j, "S-", j + 1
+  opsum .+= 0.5, "S-", j, "S+", j + 1
+  opsum .+= "Sz", j, "Sz", j + 1
 end
-H = gpu(MPO(ampo, sites))
+H = device(MPO(opsum, sites))
 
-ψ₀ = gpu(randomMPS(sites))
+ψ₀ = device(randomMPS(sites))
 
-sweeps = Sweeps(6)
-maxdim!(sweeps, 10, 20, 40, 100)
-mindim!(sweeps, 1, 10)
-cutoff!(sweeps, 1e-11)
-noise!(sweeps, 1e-10)
-energy, ψ = @time dmrg(H, ψ₀, sweeps)
+dmrg_kwargs = (;
+  nsweeps=6, maxdim=[10, 20, 40, 100], mindim=[1, 10], cutoff=1e-11, noise=1e-10
+)
+energy, ψ = @time dmrg(H, ψ₀; dmrg_kwargs...)
 @show energy
diff --git a/ITensorGPU/test/test_dmrg.jl b/ITensorGPU/test/test_dmrg.jl
index 906bceec75..811343d249 100644
--- a/ITensorGPU/test/test_dmrg.jl
+++ b/ITensorGPU/test/test_dmrg.jl
@@ -34,13 +34,13 @@ end
     N = 10
     sites = siteinds("S=1",N; conserve_qns=true)
 
-    ampo = AutoMPO()
+    opsum = OpSum()
     for j=1:N-1
-      add!(ampo,"Sz",j,"Sz",j+1)
-      add!(ampo,0.5,"S+",j,"S-",j+1)
-      add!(ampo,0.5,"S-",j,"S+",j+1)
+      add!(opsum,"Sz",j,"Sz",j+1)
+      add!(opsum,0.5,"S+",j,"S-",j+1)
+      add!(opsum,0.5,"S-",j,"S+",j+1)
     end
-    H = cuMPO(MPO(ampo,sites))
+    H = cuMPO(MPO(opsum,sites))
 
     state = [isodd(n) ? "Up" : "Dn" for n in 1:N]
     psi = randomCuMPS(sites,state,4)
@@ -63,12 +63,12 @@ end
     Random.seed!(432)
     psi0 = randomCuMPS(sites)
 
-    ampo = AutoMPO()
+    opsum = OpSum()
     for j in 1:N
-      j < N && add!(ampo, -1.0, "Sz", j, "Sz", j + 1)
-      add!(ampo, -0.5, "Sx", j)
+      j < N && add!(opsum, -1.0, "Sz", j, "Sz", j + 1)
+      add!(opsum, -0.5, "Sx", j)
     end
-    H = cuMPO(MPO(ampo, sites))
+    H = cuMPO(MPO(opsum, sites))
 
     sweeps = Sweeps(5)
     maxdim!(sweeps, 10, 20)
@@ -104,14 +104,14 @@ end
     Random.seed!(42)
     psi0 = randomCuMPS(sites)
 
-    ampo = AutoMPO()
+    opsum = OpSum()
     for j = 1:N-1
-      add!(ampo,-1.0,"Sz",j,"Sz",j+1)
+      add!(opsum,-1.0,"Sz",j,"Sz",j+1)
     end
     for j = 1:N
-      add!(ampo,-0.2,"Sx",j)
+      add!(opsum,-0.2,"Sx",j)
     end
-    H = cuMPO(MPO(ampo,sites))
+    H = cuMPO(MPO(opsum,sites))
 
     sweeps = Sweeps(3)
     maxdim!(sweeps,10)
@@ -138,18 +138,18 @@ end
     N = 10
     sites = siteinds("S=1",N)
 
-    ampoZ = AutoMPO()
+    ampoZ = OpSum()
     for j=1:N-1
-      add!(ampoZ,"Sz",j,"Sz",j+1)
+      add!(opsumZ,"Sz",j,"Sz",j+1)
     end
-    HZ = MPO(ampoZ,sites)
+    HZ = MPO(opsumZ,sites)
 
-    ampoXY = AutoMPO()
+    ampoXY = OpSum()
     for j=1:N-1
-      add!(ampoXY,0.5,"S+",j,"S-",j+1)
-      add!(ampoXY,0.5,"S-",j,"S+",j+1)
+      add!(opsumXY,0.5,"S+",j,"S-",j+1)
+      add!(opsumXY,0.5,"S-",j,"S+",j+1)
     end
-    HXY = MPO(ampoXY,sites)
+    HXY = MPO(opsumXY,sites)
 
     psi = randomMPS(sites)
 
@@ -170,13 +170,13 @@ end
     sites[1] = Index(2,"S=1/2,n=1,Site")
     sites[N] = Index(2,"S=1/2,n=$N,Site")
 
-    ampo = AutoMPO()
+    opsum = OpSum()
     for j=1:N-1
-      add!(ampo,"Sz",j,"Sz",j+1)
-      add!(ampo,0.5,"S+",j,"S-",j+1)
-      add!(ampo,0.5,"S-",j,"S+",j+1)
+      add!(opsum,"Sz",j,"Sz",j+1)
+      add!(opsum,0.5,"S+",j,"S-",j+1)
+      add!(opsum,0.5,"S-",j,"S+",j+1)
     end
-    H = cuMPO(MPO(ampo,sites))
+    H = cuMPO(MPO(opsum,sites))
 
     psi0i = randomCuMPS(sites,10)
 
@@ -211,17 +211,17 @@ end
     state[7] = 2
     psi0 = productMPS(s,state)
 
-    ampo = AutoMPO()
+    opsum = OpSum()
     for j=1:N-1
-      ampo += (-t1, "Cdag", j,   "C", j+1)
-      ampo += (-t1, "Cdag", j+1, "C", j)
-      ampo += (  V, "N",    j,   "N", j+1)
+      opsum += (-t1, "Cdag", j,   "C", j+1)
+      opsum += (-t1, "Cdag", j+1, "C", j)
+      opsum += (  V, "N",    j,   "N", j+1)
     end
     for j=1:N-2
-      ampo += (-t2, "Cdag", j,   "C", j+2)
-      ampo += (-t2, "Cdag", j+2, "C", j)
+      opsum += (-t2, "Cdag", j,   "C", j+2)
+      opsum += (-t2, "Cdag", j+2, "C", j)
     end
-    H = MPO(ampo, s)
+    H = MPO(opsum, s)
 
     sweeps = Sweeps(5)
     maxdim!(sweeps, 10, 20, 100, 100, 200)
diff --git a/ITensorGaussianMPS/README.md b/ITensorGaussianMPS/README.md
index 9b154fe59f..17588e246f 100644
--- a/ITensorGaussianMPS/README.md
+++ b/ITensorGaussianMPS/README.md
@@ -67,16 +67,10 @@ println("\nRandom state starting energy")
 @show inner(ψr, H, ψr)
 
 println("\nRun dmrg with random starting state")
-sweeps = Sweeps(10)
-setmaxdim!(sweeps,10,20,40,60)
-setcutoff!(sweeps,1E-12)
-@time dmrg(H, ψr, sweeps)
+@time dmrg(H, ψr; nsweeps=10, maxdim=[10, 20, 40, 60], cutoff=1e-12)
 
 println("\nRun dmrg with free fermion starting state")
-sweeps = Sweeps(4)
-setmaxdim!(sweeps,60)
-setcutoff!(sweeps,1E-12)
-@time dmrg(H, ψ0, sweeps)
+@time dmrg(H, ψ0; nsweeps=4, maxdim=60, cutoff=1e-12)
 ```
 This will output something like:
 ```julia
diff --git a/ITensorGaussianMPS/examples/Project.toml b/ITensorGaussianMPS/examples/Project.toml
new file mode 100644
index 0000000000..8787f2a31b
--- /dev/null
+++ b/ITensorGaussianMPS/examples/Project.toml
@@ -0,0 +1,3 @@
+[deps]
+ITensorGaussianMPS = "2be41995-7c9f-4653-b682-bfa4e7cebb93"
+ITensors = "9136182c-28ba-11e9-034c-db9fb085ebd5"
diff --git a/ITensorGaussianMPS/examples/hubbard_1d_no_spin_conservation.jl b/ITensorGaussianMPS/examples/broken/hubbard_1d_no_spin_conservation.jl
similarity index 85%
rename from ITensorGaussianMPS/examples/hubbard_1d_no_spin_conservation.jl
rename to ITensorGaussianMPS/examples/broken/hubbard_1d_no_spin_conservation.jl
index 632bbce5ca..c1bb8e3547 100644
--- a/ITensorGaussianMPS/examples/hubbard_1d_no_spin_conservation.jl
+++ b/ITensorGaussianMPS/examples/broken/hubbard_1d_no_spin_conservation.jl
@@ -69,22 +69,16 @@ H = MPO(os, s)
 
 println("Random starting state energy")
 @show flux(ψr)
-@show inner(ψr, H, ψr)
+@show inner(ψr', H, ψr)
 println()
 println("Free fermion starting state energy")
 @show flux(ψ0)
-@show inner(ψ0, H, ψ0)
+@show inner(ψ0', H, ψ0)
 
 println("\nStart from product state")
-sweeps = Sweeps(10)
-setmaxdim!(sweeps, 10, 20, _maxlinkdim)
-setcutoff!(sweeps, _cutoff)
-@time dmrg(H, ψr, sweeps)
+@time dmrg(H, ψr; nsweeps=10, maxdim=[10, 20, _maxlinkdim], cutoff=_cutoff)
 
 println("\nStart from free fermion state")
-sweeps = Sweeps(5)
-setmaxdim!(sweeps, _maxlinkdim)
-setcutoff!(sweeps, _cutoff)
-@time dmrg(H, ψ0, sweeps)
+@time dmrg(H, ψ0; nsweeps=5, maxdim=_maxlinkdim, cutoff=_cutoff)
 
 nothing
diff --git a/ITensorGaussianMPS/examples/hubbard_2d_no_spin_conservation.jl b/ITensorGaussianMPS/examples/broken/hubbard_2d_no_spin_conservation.jl
similarity index 86%
rename from ITensorGaussianMPS/examples/hubbard_2d_no_spin_conservation.jl
rename to ITensorGaussianMPS/examples/broken/hubbard_2d_no_spin_conservation.jl
index 6da63b41b3..f69a551c5e 100644
--- a/ITensorGaussianMPS/examples/hubbard_2d_no_spin_conservation.jl
+++ b/ITensorGaussianMPS/examples/broken/hubbard_2d_no_spin_conservation.jl
@@ -82,16 +82,15 @@ println("\nFree fermion MPS starting state energy")
 @show inner(ψ0, H, ψ0)
 
 println("\nStart from random product state")
-sweeps = Sweeps(10)
-setmaxdim!(sweeps, 10, 20, 100, 200, _maxlinkdim)
-setcutoff!(sweeps, _cutoff)
-setnoise!(sweeps, 1e-7, 1e-8, 1e-10, 0.0)
-@time dmrg(H, ψr, sweeps)
+dmrg_kwargs = (;
+  nsweeps=10,
+  maxdim=[10, 20, 100, 200, _maxlinkdim],
+  cutoff=_cutoff,
+  noise=[1e-7, 1e-8, 1e-10, 0.0],
+)
+@time dmrg(H, ψr; dmrg_kwargs...)
 
 println("\nStart from free fermion state")
-sweeps = Sweeps(10)
-setmaxdim!(sweeps, _maxlinkdim)
-setcutoff!(sweeps, _cutoff)
-@time dmrg(H, ψ0, sweeps)
+@time dmrg(H, ψ0; nsweeps=10, maxdim=_maxlinkdim, cutoff=_cutoff)
 
 nothing
diff --git a/ITensorGaussianMPS/examples/hubbard_1d_spin_conservation.jl b/ITensorGaussianMPS/examples/hubbard_1d_spin_conservation.jl
index 3ff3c33868..460bdd4df9 100644
--- a/ITensorGaussianMPS/examples/hubbard_1d_spin_conservation.jl
+++ b/ITensorGaussianMPS/examples/hubbard_1d_spin_conservation.jl
@@ -68,7 +68,7 @@ for n in 1:(N - 1)
 end
 
 H_noninteracting = MPO(os_noninteracting, s)
-@show inner(ψ0, H_noninteracting, ψ0)
+@show inner(ψ0', H_noninteracting, ψ0)
 @show sum(diag(Φ_up' * h_up * Φ_up)) + sum(diag(Φ_dn' * h_dn * Φ_dn))
 
 # The total interacting Hamiltonian
@@ -90,25 +90,19 @@ H = MPO(os_interacting, s)
 
 println("Random starting state energy")
 @show flux(ψr)
-@show inner(ψr, H, ψr)
+@show inner(ψr', H, ψr)
 println()
 println("Free fermion starting state energy")
 @show flux(ψ0)
-@show inner(ψ0, H, ψ0)
+@show inner(ψ0', H, ψ0)
 
 println("\nStart from random product state")
-sweeps = Sweeps(10)
-setmaxdim!(sweeps, 10, 20, _maxlinkdim)
-setcutoff!(sweeps, _cutoff)
-er, ψ̃r = @time dmrg(H, ψr, sweeps)
+er, ψ̃r = @time dmrg(H, ψr; nsweeps=10, maxdim=[10, 20, _maxlinkdim], cutoff=_cutoff)
 @show er
 @show flux(ψ̃r)
 
 println("\nStart from free fermion state")
-sweeps = Sweeps(5)
-setmaxdim!(sweeps, _maxlinkdim)
-setcutoff!(sweeps, _cutoff)
-e0, ψ̃0 = @time dmrg(H, ψ0, sweeps)
+e0, ψ̃0 = @time dmrg(H, ψ0; nsweeps=5, maxdim=_maxlinkdim, cutoff=_cutoff)
 @show e0
 @show flux(ψ̃0)
 
diff --git a/ITensorGaussianMPS/examples/hubbard_2d_spin_conservation.jl b/ITensorGaussianMPS/examples/hubbard_2d_spin_conservation.jl
index 2984f8dc23..41274b5a43 100644
--- a/ITensorGaussianMPS/examples/hubbard_2d_spin_conservation.jl
+++ b/ITensorGaussianMPS/examples/hubbard_2d_spin_conservation.jl
@@ -85,23 +85,22 @@ H = MPO(os, s)
 
 println("\nRandom starting state energy")
 @show flux(ψr)
-@show inner(ψr, H, ψr)
+@show inner(ψr', H, ψr)
 
 println("\nFree fermion MPS starting state energy")
 @show flux(ψ0)
-@show inner(ψ0, H, ψ0)
+@show inner(ψ0', H, ψ0)
 
 println("\nStart from random product state")
-sweeps = Sweeps(10)
-setmaxdim!(sweeps, 10, 20, 100, 200, _maxlinkdim)
-setcutoff!(sweeps, _cutoff)
-setnoise!(sweeps, 1e-7, 1e-8, 1e-10, 0.0)
-@time dmrg(H, ψr, sweeps)
+dmrg_kwargs = (;
+  nsweeps=10,
+  maxdim=[10, 20, 100, 200, _maxlinkdim],
+  cutoff=_cutoff,
+  noise=[1e-7, 1e-8, 1e-10, 0.0],
+)
+@time dmrg(H, ψr; dmrg_kwargs...)
 
 println("\nStart from free fermion state")
-sweeps = Sweeps(10)
-setmaxdim!(sweeps, _maxlinkdim)
-setcutoff!(sweeps, _cutoff)
-@time dmrg(H, ψ0, sweeps)
+@time dmrg(H, ψ0; nsweeps=10, maxdim=_maxlinkdim, cutoff=_cutoff)
 
 nothing
diff --git a/ITensorGaussianMPS/examples/mps_to_determinants.jl b/ITensorGaussianMPS/examples/mps_to_determinants.jl
index 335ac536f8..c3590b4bc7 100644
--- a/ITensorGaussianMPS/examples/mps_to_determinants.jl
+++ b/ITensorGaussianMPS/examples/mps_to_determinants.jl
@@ -65,13 +65,10 @@ H = MPO(os, s)
 
 println("Free fermion starting state energy")
 @show flux(ψ0)
-@show inner(ψ0, H, ψ0)
+@show inner(ψ0', H, ψ0)
 
 println("\nStart from free fermion state")
-sweeps = Sweeps(5)
-setmaxdim!(sweeps, _maxlinkdim)
-setcutoff!(sweeps, _cutoff)
-e, ψ = @time dmrg(H, ψ0, sweeps)
+e, ψ = @time dmrg(H, ψ0; nsweeps=5, maxdim=_maxlinkdim, cutoff=_cutoff)
 @show e
 @show flux(ψ)
 
@@ -81,4 +78,4 @@ using ITensorGaussianMPS: correlation_matrix_to_gmps, correlation_matrix_to_mps,
 Λ_dn = correlation_matrix(ψ, "Cdagdn", "Cdn")
 ψ̃0 = correlation_matrix_to_mps(s, Λ_up, Λ_dn; eigval_cutoff=1e-2, maxblocksize=4)
 @show inner(ψ̃0, ψ)
-@show inner(ψ̃0, H, ψ̃0)
+@show inner(ψ̃0', H, ψ̃0)
diff --git a/ITensorGaussianMPS/examples/spinless_fermion.jl b/ITensorGaussianMPS/examples/spinless_fermion.jl
index 63db2e3dbe..7d2f74b9af 100644
--- a/ITensorGaussianMPS/examples/spinless_fermion.jl
+++ b/ITensorGaussianMPS/examples/spinless_fermion.jl
@@ -56,22 +56,16 @@ H = MPO(os, s)
 
 println("\nRandom state starting energy")
 @show flux(ψr)
-@show inner(ψr, H, ψr)
+@show inner(ψr', H, ψr)
 
 println("\nFree fermion starting energy")
 @show flux(ψ0)
-@show inner(ψ0, H, ψ0)
+@show inner(ψ0', H, ψ0)
 
 println("\nRun dmrg with random starting state")
-sweeps = Sweeps(20)
-setmaxdim!(sweeps, 10, 20, 40, _maxlinkdim)
-setcutoff!(sweeps, _cutoff)
-@time dmrg(H, ψr, sweeps)
+@time dmrg(H, ψr; nsweeps=20, maxdim=[10, 20, 40, _maxlinkdim], cutoff=_cutoff)
 
 println("\nRun dmrg with free fermion starting state")
-sweeps = Sweeps(4)
-setmaxdim!(sweeps, _maxlinkdim)
-setcutoff!(sweeps, _cutoff)
-@time dmrg(H, ψ0, sweeps)
+@time dmrg(H, ψ0; nsweeps=4, maxdim=_maxlinkdim, cutoff=_cutoff)
 
 nothing
diff --git a/ITensorGaussianMPS/examples/spinless_fermion_pairing.jl b/ITensorGaussianMPS/examples/spinless_fermion_pairing.jl
index b1d477ce2e..145b042e41 100644
--- a/ITensorGaussianMPS/examples/spinless_fermion_pairing.jl
+++ b/ITensorGaussianMPS/examples/spinless_fermion_pairing.jl
@@ -67,19 +67,16 @@ let
 
   println("\nFree fermion starting energy")
   @show flux(psi)
-  @show inner(psi, H, psi)
+  @show inner(psi', H, psi)
   println("\nRun dmrg with GMPS starting state")
-  sweeps = Sweeps(12)
-  setmaxdim!(sweeps, 10, 20, 40, _maxlinkdim)
-  setcutoff!(sweeps, _cutoff)
-  _, psidmrg = dmrg(H, psi, sweeps)
+  _, psidmrg = dmrg(H, psi; nsweeps=12, maxdim=[10, 20, 40, _maxlinkdim], cutoff=_cutoff)
   cdagc_dmrg = correlation_matrix(psidmrg, "C", "Cdag")
   cc_dmrg = correlation_matrix(psidmrg, "C", "C")
 
   @show norm(cdagc_dmrg - cdagc)
   @show norm(cc_dmrg - cc)
 
-  @show inner(psidmrg, H, psidmrg)
+  @show inner(psidmrg', H, psidmrg)
   @show(abs(inner(psidmrg, psi)))
 
   #return
diff --git a/ITensorUnicodePlots/examples/ex_qn_mps.jl b/ITensorUnicodePlots/examples/ex_qn_mps.jl
index d6c4edf370..11a3b59b92 100644
--- a/ITensorUnicodePlots/examples/ex_qn_mps.jl
+++ b/ITensorUnicodePlots/examples/ex_qn_mps.jl
@@ -3,7 +3,7 @@ using ITensorUnicodePlots
 
 s = siteinds("S=1/2", 5; conserve_qns=true)
 ψ = randomMPS(s, n -> isodd(n) ? "↑" : "↓"; linkdims=2)
-orthogonalize!(ψ, 2)
+ψ = orthogonalize(ψ, 2)
 ψdag = prime(linkinds, dag(ψ))
 tn = [ψ..., ψdag...]
 
diff --git a/docs/src/examples/MPSandMPO.md b/docs/src/examples/MPSandMPO.md
index 54ca7f226e..8feb12a3ca 100644
--- a/docs/src/examples/MPSandMPO.md
+++ b/docs/src/examples/MPSandMPO.md
@@ -223,7 +223,7 @@ removing the prime from the ``s'_3`` index afterward:
 
 ```julia
 newA = G * psi[3]
-noprime!(newA)
+newA = noprime(newA)
 ```
 
 Finally, put the new tensor back into MPS `psi` to update its third MPS tensor:
@@ -238,11 +238,11 @@ Afterward, we can visualize the modified MPS as:
 
 As a technical note, if you are working in a context where gauge or orthogonality
 properties of the MPS are important, such as in time evolution using two-site gates,
-then you may want to call `orthogonalize!(psi,3)`
+then you may want to call `psi = orthogonalize(psi, 3)`
 before modifying the tensor at site 3, which will ensure that the MPS remains in a
 well-defined orthogonal gauge centered on site 3. Modifying a tensor which is left- or right-orthogonal
 (i.e. not the "center" tensor of the gauge) will destroy the gauge condition and
-require extra operations to restore it. (Calling `orthogonalize!` method will automatically
+require extra operations to restore it. (Calling `orthogonalize` method will automatically
 fix this but will have to do extra work to do so.)
 
 
@@ -264,7 +264,7 @@ that either site 3 or 4 is the *orthogonality center*. Here we make site 3
 the center:
 
 ```julia
-orthogonalize!(psi,3)
+psi = orthogonalize(psi, 3)
 ```
 
 ![](twosite_figures/gate_gauge.png)
@@ -278,7 +278,7 @@ Next, contract the gate tensor G with the MPS tensors for sites 3 and 4
 
 ```julia
 wf = (psi[3] * psi[4]) * G
-noprime!(wf)
+wf = noprime(wf)
 ```
 
 Finally, use the singular value decomposition (SVD) to factorize the
@@ -306,15 +306,15 @@ as well as limits on the maximum bond dimension (`maxdim` keyword argument).
 **Complete code example**
 
 ```julia
-orthogonalize!(psi,3)
+psi = orthogonalize(psi, 3)
 
 wf = (psi[3] * psi[4]) * G
-noprime!(wf)
+wf = noprime(wf)
 
-inds3 = uniqueinds(psi[3],psi[4])
-U,S,V = svd(wf,inds3,cutoff=1E-8)
+inds3 = uniqueinds(psi[3], psi[4])
+U, S, V = svd(wf, inds3; cutoff=1E-8)
 psi[3] = U
-psi[4] = S*V
+psi[4] = S * V
 ```
 
 ## Computing the Entanglement Entropy of an MPS
@@ -326,8 +326,8 @@ Say that we have obtained an MPS `psi` of length N and we wish to compute the en
 Then the following code formula can be used to accomplish this task:
 
 ```julia
-orthogonalize!(psi, b)
-U,S,V = svd(psi[b], (linkind(psi, b-1), siteind(psi,b)))
+psi = orthogonalize(psi, b)
+U,S,V = svd(psi[b], (linkinds(psi, b-1)..., siteinds(psi, b)...))
 SvN = 0.0
 for n=1:dim(S, 1)
   p = S[n,n]^2
@@ -335,7 +335,7 @@ for n=1:dim(S, 1)
 end
 ```
 
-As a brief explanation of the code above, the call to `orthogonalize!(psi,b)`
+As a brief explanation of the code above, the call to `psi = orthogonalize(psi, b)`
 shifts the orthogonality center to site `b` of the MPS.
 
 The call to the `svd` routine says to treat the link (virtual or bond) Index connecting the b'th MPS tensor `psi[b]` and the b'th physical Index as "row" indices for the purposes of the SVD (these indices will end up on `U`, along with the Index connecting `U` to `S`).
@@ -379,7 +379,7 @@ in the "computational basis".
 
 For reasons of efficiency, the `sample` function requires the MPS to be in orthogonal
 form, orthogonalized to the first site. If it is not already in this form, it
-can be brought into orthogonal form by calling `orthogonalize!(psi,1)`.
+can be brought into orthogonal form by calling `psi = orthogonalize(psi, 1)`.
 
 
 ## Write and Read an MPS or MPO to Disk with HDF5
diff --git a/examples/dmrg/input_files/argsdict/1d_hubbard_extended.jl b/examples/dmrg/input_files/argsdict/1d_hubbard_extended.jl
index 917f6d290f..fed0a41db8 100644
--- a/examples/dmrg/input_files/argsdict/1d_hubbard_extended.jl
+++ b/examples/dmrg/input_files/argsdict/1d_hubbard_extended.jl
@@ -128,7 +128,7 @@ energy, psi = dmrg(H, psi0, sweeps)
 upd = fill(0.0, N)
 dnd = fill(0.0, N)
 for j in 1:N
-  orthogonalize!(psi, j)
+  psi = orthogonalize(psi, j)
   psidag_j = dag(prime(psi[j], "Site"))
   upd[j] = scalar(psidag_j * op(sites, "Nup", j) * psi[j])
   dnd[j] = scalar(psidag_j * op(sites, "Ndn", j) * psi[j])
diff --git a/examples/src/hubbard.jl b/examples/src/hubbard.jl
index a8ad9a00b0..2810b4e37e 100644
--- a/examples/src/hubbard.jl
+++ b/examples/src/hubbard.jl
@@ -1,15 +1,15 @@
 
 function hubbard_1d(; N::Int, t=1.0, U=0.0)
-  ampo = OpSum()
+  opsum = OpSum()
   for b in 1:(N - 1)
-    ampo .+= -t, "Cdagup", b, "Cup", b + 1
-    ampo .+= -t, "Cdagup", b + 1, "Cup", b
-    ampo .+= -t, "Cdagdn", b, "Cdn", b + 1
-    ampo .+= -t, "Cdagdn", b + 1, "Cdn", b
+    opsum .+= -t, "Cdagup", b, "Cup", b + 1
+    opsum .+= -t, "Cdagup", b + 1, "Cup", b
+    opsum .+= -t, "Cdagdn", b, "Cdn", b + 1
+    opsum .+= -t, "Cdagdn", b + 1, "Cdn", b
   end
   if U ≠ 0
     for n in 1:N
-      ampo .+= U, "Nupdn", n
+      opsum .+= U, "Nupdn", n
     end
   end
   return ampo
@@ -18,30 +18,30 @@ end
 function hubbard_2d(; Nx::Int, Ny::Int, t=1.0, U=0.0, yperiodic::Bool=true)
   N = Nx * Ny
   lattice = square_lattice(Nx, Ny; yperiodic=yperiodic)
-  ampo = OpSum()
+  opsum = OpSum()
   for b in lattice
-    ampo .+= -t, "Cdagup", b.s1, "Cup", b.s2
-    ampo .+= -t, "Cdagup", b.s2, "Cup", b.s1
-    ampo .+= -t, "Cdagdn", b.s1, "Cdn", b.s2
-    ampo .+= -t, "Cdagdn", b.s2, "Cdn", b.s1
+    opsum .+= -t, "Cdagup", b.s1, "Cup", b.s2
+    opsum .+= -t, "Cdagup", b.s2, "Cup", b.s1
+    opsum .+= -t, "Cdagdn", b.s1, "Cdn", b.s2
+    opsum .+= -t, "Cdagdn", b.s2, "Cdn", b.s1
   end
   if U ≠ 0
     for n in 1:N
-      ampo .+= U, "Nupdn", n
+      opsum .+= U, "Nupdn", n
     end
   end
   return ampo
 end
 
 function hubbard_2d_ky(; Nx::Int, Ny::Int, t=1.0, U=0.0)
-  ampo = OpSum()
+  opsum = OpSum()
   for x in 0:(Nx - 1)
     for ky in 0:(Ny - 1)
       s = x * Ny + ky + 1
       disp = -2 * t * cos((2 * π / Ny) * ky)
       if abs(disp) > 1e-12
-        ampo .+= disp, "Nup", s
-        ampo .+= disp, "Ndn", s
+        opsum .+= disp, "Nup", s
+        opsum .+= disp, "Ndn", s
       end
     end
   end
@@ -49,10 +49,10 @@ function hubbard_2d_ky(; Nx::Int, Ny::Int, t=1.0, U=0.0)
     for ky in 0:(Ny - 1)
       s1 = x * Ny + ky + 1
       s2 = (x + 1) * Ny + ky + 1
-      ampo .+= -t, "Cdagup", s1, "Cup", s2
-      ampo .+= -t, "Cdagup", s2, "Cup", s1
-      ampo .+= -t, "Cdagdn", s1, "Cdn", s2
-      ampo .+= -t, "Cdagdn", s2, "Cdn", s1
+      opsum .+= -t, "Cdagup", s1, "Cup", s2
+      opsum .+= -t, "Cdagup", s2, "Cup", s1
+      opsum .+= -t, "Cdagdn", s1, "Cdn", s2
+      opsum .+= -t, "Cdagdn", s2, "Cdn", s1
     end
   end
   if U ≠ 0
@@ -64,7 +64,7 @@ function hubbard_2d_ky(; Nx::Int, Ny::Int, t=1.0, U=0.0)
             s2 = x * Ny + (py - qy + Ny) % Ny + 1
             s3 = x * Ny + py + 1
             s4 = x * Ny + ky + 1
-            ampo .+= (U / Ny), "Cdagdn", s1, "Cdagup", s2, "Cup", s3, "Cdn", s4
+            opsum .+= (U / Ny), "Cdagdn", s1, "Cdagup", s2, "Cup", s3, "Cdn", s4
           end
         end
       end
@@ -75,11 +75,11 @@ end
 
 function hubbard(; Nx::Int, Ny::Int=1, t=1.0, U=0.0, yperiodic::Bool=true, ky::Bool=false)
   if Ny == 1
-    ampo = hubbard_1d(; N=Nx, t=t, U=U)
+    opsum = hubbard_1d(; N=Nx, t=t, U=U)
   elseif ky
-    ampo = hubbard_2d_ky(; Nx=Nx, Ny=Ny, t=t, U=U)
+    opsum = hubbard_2d_ky(; Nx=Nx, Ny=Ny, t=t, U=U)
   else
-    ampo = hubbard_2d(; Nx=Nx, Ny=Ny, yperiodic=yperiodic, t=t, U=U)
+    opsum = hubbard_2d(; Nx=Nx, Ny=Ny, yperiodic=yperiodic, t=t, U=U)
   end
   return ampo
 end
diff --git a/precompile/snoop/dmrg.jl b/precompile/snoop/dmrg.jl
index 974510cb36..e6cab70284 100644
--- a/precompile/snoop/dmrg.jl
+++ b/precompile/snoop/dmrg.jl
@@ -5,13 +5,13 @@ let
   # to precompile.jl by hand.
   N = 4
   sites = siteinds("S=1", N)
-  ampo = AutoMPO()
+  opsum = OpSum()
   for j in 1:(N - 1)
-    ampo .+= ("Sz", j, "Sz", j + 1)
-    ampo .+= (0.5, "S+", j, "S-", j + 1)
-    ampo .+= (0.5, "S-", j, "S+", j + 1)
+    opsum .+= ("Sz", j, "Sz", j + 1)
+    opsum .+= (0.5, "S+", j, "S-", j + 1)
+    opsum .+= (0.5, "S-", j, "S+", j + 1)
   end
-  H = MPO(ampo, sites)
+  H = MPO(opsum, sites)
   psi0 = randomMPS(sites, 10)
   sweeps = Sweeps(1)
   maxdim!(sweeps, 10)
diff --git a/src/ITensorMPS/abstractmps.jl b/src/ITensorMPS/abstractmps.jl
index d5d38531d0..6e2912085f 100644
--- a/src/ITensorMPS/abstractmps.jl
+++ b/src/ITensorMPS/abstractmps.jl
@@ -103,7 +103,7 @@ Returns the range of sites of the orthogonality center of the MPS/MPO.
 ```julia
 s = siteinds("S=½", 5)
 ψ = randomMPS(s)
-orthogonalize!(ψ, 3)
+ψ = orthogonalize(ψ, 3)
 
 # ortho_lims(ψ) = 3:3
 @show ortho_lims(ψ)
@@ -174,8 +174,8 @@ s = siteinds("S=1/2", 4)
 # Make random MPS with bond dimension 2
 ψ₁ = randomMPS(s, "↑", 2)
 ψ₂ = randomMPS(s, "↑", 2)
-orthogonalize!(ψ₁, 1)
-orthogonalize!(ψ₂, 1)
+ψ₁ = orthogonalize(ψ₁, 1)
+ψ₂ = orthogonalize(ψ₂, 1)
 
 # ortho_lims(ψ₁) = 1:1
 @show ortho_lims(ψ₁)
@@ -1948,7 +1948,7 @@ function (::Type{MPST})(
   M = MPST(ψ)
   setleftlim!(M, N - 1)
   setrightlim!(M, N + 1)
-  orthogonalize!(M, orthocenter)
+  M = orthogonalize(M, orthocenter)
   return M
 end
 
@@ -1969,9 +1969,9 @@ function swapbondsites(ψ::AbstractMPS, b::Integer; ortho="right", kwargs...)
     orthocenter = b
   end
   if leftlim(ψ) < b - 1
-    orthogonalize!(ψ, b)
+    ψ = orthogonalize(ψ, b)
   elseif rightlim(ψ) > b + 2
-    orthogonalize!(ψ, b + 1)
+    ψ = orthogonalize(ψ, b + 1)
   end
   ψ[b:(b + 1), orthocenter=orthocenter, perm=[2, 1], kwargs...] = ψ[b] * ψ[b + 1]
   return ψ
diff --git a/src/ITensorMPS/mps.jl b/src/ITensorMPS/mps.jl
index a6bec60f96..612f531797 100644
--- a/src/ITensorMPS/mps.jl
+++ b/src/ITensorMPS/mps.jl
@@ -784,8 +784,7 @@ function correlation_matrix(
     end
   end
 
-  psi = copy(psi)
-  orthogonalize!(psi, start_site)
+  psi = orthogonalize(psi, start_site)
   norm2_psi = norm(psi[start_site])^2
 
   # Nb = size of block of correlation matrix
@@ -994,13 +993,13 @@ function expect(psi::MPS, ops; sites=1:length(psi), site_range=nothing)
 
   el_types = map(o -> ishermitian(op(o, s[start_site])) ? real(ElT) : ElT, ops)
 
-  orthogonalize!(psi, start_site)
+  psi = orthogonalize(psi, start_site)
   norm2_psi = norm(psi)^2
   iszero(norm2_psi) && error("MPS has zero norm in function `expect`")
 
   ex = map((o, el_t) -> zeros(el_t, Ns), ops, el_types)
   for (entry, j) in enumerate(site_range)
-    orthogonalize!(psi, j)
+    psi = orthogonalize(psi, j)
     for (n, opname) in enumerate(ops)
       oⱼ = adapt(datatype(psi[j]), op(opname, s[j]))
       val = inner(psi[j], apply(oⱼ, psi[j])) / norm2_psi
diff --git a/src/packagecompile/precompile_itensors.jl b/src/packagecompile/precompile_itensors.jl
index dfb6c120db..d99242873d 100644
--- a/src/packagecompile/precompile_itensors.jl
+++ b/src/packagecompile/precompile_itensors.jl
@@ -10,26 +10,24 @@ using ITensors
 #                  "runtests.jl"))
 
 N = 6
-sweeps = Sweeps(3)
-maxdim!(sweeps, 10)
-cutoff!(sweeps, 1E-13)
+dmrg_kwargs = (; nsweeps=3, maxdim=10, cutoff=1e-13)
 
-ampo = OpSum()
+opsum = OpSum()
 for j in 1:(N - 1)
-  ampo .+= "Sz", j, "Sz", j + 1
-  ampo .+= 0.5, "S+", j, "S-", j + 1
-  ampo .+= 0.5, "S-", j, "S+", j + 1
+  opsum .+= "Sz", j, "Sz", j + 1
+  opsum .+= 0.5, "S+", j, "S-", j + 1
+  opsum .+= 0.5, "S-", j, "S+", j + 1
 end
 
 sites = siteinds("S=1", N)
-H = MPO(ampo, sites)
+H = MPO(opsum, sites)
 psi0 = randomMPS(sites; linkdims=2)
-dmrg(H, psi0, sweeps)
+dmrg(H, psi0; dmrg_kwargs...)
 
 sites_qn = siteinds("S=1", N; conserve_qns=true)
 if !hasqns(sites_qn[1])
   throw(ErrorException("Index does not have QNs in part of precompile script"))
 end
-H_qn = MPO(ampo, sites_qn)
+H_qn = MPO(opsum, sites_qn)
 psi0_qn = randomMPS(sites_qn, [isodd(n) ? "Up" : "Dn" for n in 1:N]; linkdims=2)
-dmrg(H_qn, psi0_qn, sweeps)
+dmrg(H_qn, psi0_qn; dmrg_kwargs...)
diff --git a/test/ITensorMPS/base/test_mps.jl b/test/ITensorMPS/base/test_mps.jl
index d44e0edf1e..319c213309 100644
--- a/test/ITensorMPS/base/test_mps.jl
+++ b/test/ITensorMPS/base/test_mps.jl
@@ -1797,13 +1797,13 @@ end
 
       t = 1.0
       U = 1.0
-      ampo = OpSum()
+      opsum = OpSum()
       for b in 1:(N - 1)
-        ampo .+= -t, "Cdag", b, "C", b + 1
-        ampo .+= -t, "Cdag", b + 1, "C", b
-        ampo .+= U, "N", b, "N", b + 1
+        opsum .+= -t, "Cdag", b, "C", b + 1
+        opsum .+= -t, "Cdag", b + 1, "C", b
+        opsum .+= U, "N", b, "N", b + 1
       end
-      H = MPO(ampo, s)
+      H = MPO(opsum, s)
 
       sweeps = Sweeps(6)
       maxdim!(sweeps, 10, 20, 40)
@@ -1831,9 +1831,9 @@ end
           G2 *= op("C", s, j)
         end
 
-        ampo = OpSum()
-        ampo += "Cdag", i, "C", j
-        G3 = MPO(ampo, s)
+        opsum = OpSum()
+        opsum += "Cdag", i, "C", j
+        G3 = MPO(opsum, s)
 
         A_OP = prod(product(G1, ψ0; cutoff=1e-6))
 
@@ -1862,9 +1862,9 @@ end
           end
           G2 *= op("C", s, l)
 
-          ampo = OpSum()
-          ampo += "Cdag", i, "Cdag", j, "C", k, "C", l
-          G3 = MPO(ampo, s)
+          opsum = OpSum()
+          opsum += "Cdag", i, "Cdag", j, "C", k, "C", l
+          G3 = MPO(opsum, s)
 
           A_OP = prod(product(G1, ψ0; cutoff=1e-16))
 
@@ -1883,17 +1883,17 @@ end
       ψ0 = randomMPS(s, n -> isodd(n) ? "↑" : "↓")
       t = 1.0
       U = 1.0
-      ampo = OpSum()
+      opsum = OpSum()
       for b in 1:(N - 1)
-        ampo .+= -t, "Cdagup", b, "Cup", b + 1
-        ampo .+= -t, "Cdagup", b + 1, "Cup", b
-        ampo .+= -t, "Cdagdn", b, "Cdn", b + 1
-        ampo .+= -t, "Cdagdn", b + 1, "Cdn", b
+        opsum .+= -t, "Cdagup", b, "Cup", b + 1
+        opsum .+= -t, "Cdagup", b + 1, "Cup", b
+        opsum .+= -t, "Cdagdn", b, "Cdn", b + 1
+        opsum .+= -t, "Cdagdn", b + 1, "Cdn", b
       end
       for n in 1:N
-        ampo .+= U, "Nupdn", n
+        opsum .+= U, "Nupdn", n
       end
-      H = MPO(ampo, s)
+      H = MPO(opsum, s)
       sweeps = Sweeps(6)
       maxdim!(sweeps, 10, 20, 40)
       cutoff!(sweeps, 1E-12)
@@ -1904,9 +1904,9 @@ end
       end
 
       for i in 1:(N - 1), j in (i + 1):N
-        ampo = OpSum()
-        ampo += "Cdagup", i, "Cup", j
-        G1 = MPO(ampo, s)
+        opsum = OpSum()
+        opsum += "Cdagup", i, "Cup", j
+        G1 = MPO(opsum, s)
         G2 = op("CCup", s, i, j)
         A_MPO = prod(noprime(contract(G1, ψ; cutoff=1e-8)))
         A_OP = prod(product(G2, ψ; cutoff=1e-8))
diff --git a/test/ITensorMPS/base/test_qnmpo.jl b/test/ITensorMPS/base/test_qnmpo.jl
index f74710d71c..8f72e6577e 100644
--- a/test/ITensorMPS/base/test_qnmpo.jl
+++ b/test/ITensorMPS/base/test_qnmpo.jl
@@ -197,22 +197,22 @@ end
 @testset "splitblocks" begin
   N = 4
   sites = siteinds("S=1", N; conserve_qns=true)
-  ampo = OpSum()
+  opsum = OpSum()
   for j in 1:(N - 1)
-    ampo .+= 0.5, "S+", j, "S-", j + 1
-    ampo .+= 0.5, "S-", j, "S+", j + 1
-    ampo .+= "Sz", j, "Sz", j + 1
+    opsum .+= 0.5, "S+", j, "S-", j + 1
+    opsum .+= 0.5, "S-", j, "S+", j + 1
+    opsum .+= "Sz", j, "Sz", j + 1
   end
-  H = MPO(ampo, sites; splitblocks=false)
+  H = MPO(opsum, sites; splitblocks=false)
 
   # Split the tensors to make them more sparse
   # Drops zero blocks by default
   H̃ = splitblocks(linkinds, H)
 
-  H̃2 = MPO(ampo, sites; splitblocks=true)
+  H̃2 = MPO(opsum, sites; splitblocks=true)
 
   # Defaults to true
-  H̃3 = MPO(ampo, sites)
+  H̃3 = MPO(opsum, sites)
 
   @test prod(H) ≈ prod(H̃)
   @test prod(H) ≈ prod(H̃2)
@@ -299,16 +299,16 @@ end
 function make_Heisenberg_AutoMPO(sites, NNN::Int64; J::Float64=1.0, kwargs...)::MPO
   N = length(sites)
   @assert N >= NNN
-  ampo = OpSum()
+  opsum = OpSum()
   for dj in 1:NNN
     f = J / dj
     for j in 1:(N - dj)
-      add!(ampo, f, "Sz", j, "Sz", j + dj)
-      add!(ampo, f * 0.5, "S+", j, "S-", j + dj)
-      add!(ampo, f * 0.5, "S-", j, "S+", j + dj)
+      add!(opsum, f, "Sz", j, "Sz", j + dj)
+      add!(opsum, f * 0.5, "S+", j, "S-", j + dj)
+      add!(opsum, f * 0.5, "S-", j, "S+", j + dj)
     end
   end
-  return MPO(ampo, sites; kwargs...)
+  return MPO(opsum, sites; kwargs...)
 end
 
 function make_Hubbard_AutoMPO(
diff --git a/test/ITensorMPS/base/test_readme.jl b/test/ITensorMPS/base/test_readme.jl
index 79688626e1..414a8d4872 100644
--- a/test/ITensorMPS/base/test_readme.jl
+++ b/test/ITensorMPS/base/test_readme.jl
@@ -76,13 +76,13 @@ using ITensors, Test
     # Input operator terms which define
     # a Hamiltonian matrix, and convert
     # these terms to an MPO tensor network
-    ampo = OpSum()
+    opsum = OpSum()
     for j in 1:(N - 1)
-      add!(ampo, "Sz", j, "Sz", j + 1)
-      add!(ampo, 0.5, "S+", j, "S-", j + 1)
-      add!(ampo, 0.5, "S-", j, "S+", j + 1)
+      add!(opsum, "Sz", j, "Sz", j + 1)
+      add!(opsum, 0.5, "S+", j, "S-", j + 1)
+      add!(opsum, 0.5, "S-", j, "S+", j + 1)
     end
-    H = MPO(ampo, sites)
+    H = MPO(opsum, sites)
 
     # Create an initial random matrix product state
     psi0 = randomMPS(sites)
diff --git a/test/ITensorMPS/base/test_threading.jl b/test/ITensorMPS/base/test_threading.jl
index 4af99a3a3d..885923f40e 100644
--- a/test/ITensorMPS/base/test_threading.jl
+++ b/test/ITensorMPS/base/test_threading.jl
@@ -31,17 +31,17 @@ end
     noise!(sweeps, 1e-6, 1e-7, 1e-8, 0.0)
     sites = siteinds("Electron", N; conserve_qns=true)
     lattice = square_lattice(Nx, Ny; yperiodic=true)
-    ampo = OpSum()
+    opsum = OpSum()
     for b in lattice
-      ampo .+= -t, "Cdagup", b.s1, "Cup", b.s2
-      ampo .+= -t, "Cdagup", b.s2, "Cup", b.s1
-      ampo .+= -t, "Cdagdn", b.s1, "Cdn", b.s2
-      ampo .+= -t, "Cdagdn", b.s2, "Cdn", b.s1
+      opsum .+= -t, "Cdagup", b.s1, "Cup", b.s2
+      opsum .+= -t, "Cdagup", b.s2, "Cup", b.s1
+      opsum .+= -t, "Cdagdn", b.s1, "Cdn", b.s2
+      opsum .+= -t, "Cdagdn", b.s2, "Cdn", b.s1
     end
     for n in 1:N
-      ampo .+= U, "Nupdn", n
+      opsum .+= U, "Nupdn", n
     end
-    H = MPO(ampo, sites)
+    H = MPO(opsum, sites)
     Hsplit = splitblocks(linkinds, H)
     state = [isodd(n) ? "↑" : "↓" for n in 1:N]
     ψ0 = productMPS(sites, state)