diff --git a/Project.toml b/Project.toml
index ca542f2944..85ba97eec2 100644
--- a/Project.toml
+++ b/Project.toml
@@ -17,6 +17,7 @@ KrylovKit = "0b1a1467-8014-51b9-945f-bf0ae24f4b77"
 LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
 LinearMaps = "7a12625a-238d-50fd-b39a-03d52299707e"
 NDTensors = "23ae76d9-e61a-49c4-8f12-3f1a16adf9cf"
+Observers = "338f10d5-c7f1-4033-a7d1-f9dec39bcaa0"
 PackageCompiler = "9b87118b-4619-50d2-8e1e-99f35a4d4d9d"
 Pkg = "44cfe95a-1eb2-52ea-b672-e2afdf69b78f"
 Printf = "de0858da-6303-5e67-8744-51eddeeeb8d7"
diff --git a/src/ITensorMPS/ITensorMPS.jl b/src/ITensorMPS/ITensorMPS.jl
index f8fdb9dc07..fd4479e7ac 100644
--- a/src/ITensorMPS/ITensorMPS.jl
+++ b/src/ITensorMPS/ITensorMPS.jl
@@ -1,6 +1,7 @@
 module ITensorMPS
 
-include("usings.jl")
+using ..ITensors
+
 include("imports.jl")
 include("exports.jl")
 include("abstractmps.jl")
@@ -23,4 +24,19 @@ include("autompo/opsum_to_mpo.jl")
 include("autompo/opsum_to_mpo_qn.jl")
 include("deprecated.jl")
 
+include("abstractprojmpo/projmpo_apply.jl")
+include("abstractprojmpo/projmps2.jl")
+include("abstractprojmpo/projmpo_mps2.jl")
+include("defaults.jl")
+include("update_observer.jl")
+include("solver_utils.jl")
+include("tdvporder.jl")
+include("tdvpinfo.jl")
+include("sweep_update.jl")
+include("alternating_update.jl")
+include("tdvp.jl")
+include("dmrg_x.jl")
+include("contract_mpo_mps.jl")
+include("linsolve.jl")
+
 end # module ITensorMPS
diff --git a/src/ITensorMPS/abstractmps.jl b/src/ITensorMPS/abstractmps.jl
index 6e2912085f..c4b98a2716 100644
--- a/src/ITensorMPS/abstractmps.jl
+++ b/src/ITensorMPS/abstractmps.jl
@@ -1,3 +1,6 @@
+using IsApprox: Approx, IsApprox
+using NDTensors: using_auto_fermion, scalartype, tensor
+
 abstract type AbstractMPS end
 
 """
@@ -45,9 +48,9 @@ have type `ComplexF64`, return `ComplexF64`.
 """
 promote_itensor_eltype(m::AbstractMPS) = LinearAlgebra.promote_leaf_eltypes(m)
 
-scalartype(m::AbstractMPS) = LinearAlgebra.promote_leaf_eltypes(m)
-scalartype(m::Array{ITensor}) = LinearAlgebra.promote_leaf_eltypes(m)
-scalartype(m::Array{<:Array{ITensor}}) = LinearAlgebra.promote_leaf_eltypes(m)
+NDTensors.scalartype(m::AbstractMPS) = LinearAlgebra.promote_leaf_eltypes(m)
+NDTensors.scalartype(m::Array{ITensor}) = LinearAlgebra.promote_leaf_eltypes(m)
+NDTensors.scalartype(m::Array{<:Array{ITensor}}) = LinearAlgebra.promote_leaf_eltypes(m)
 
 """
     eltype(m::MPS)
@@ -1141,7 +1144,7 @@ Same as [`inner`](@ref).
 
 See also [`loginner`](@ref), [`logdot`](@ref).
 """
-function dot(M1::MPST, M2::MPST; kwargs...) where {MPST<:AbstractMPS}
+function LinearAlgebra.dot(M1::MPST, M2::MPST; kwargs...) where {MPST<:AbstractMPS}
   return _log_or_not_dot(M1, M2, false; kwargs...)
 end
 
@@ -1319,7 +1322,7 @@ lognorm_ψ[1] == -Inf # There was an infinite norm
 
 See also [`normalize`](@ref), [`norm`](@ref), [`lognorm`](@ref).
 """
-function normalize!(M::AbstractMPS; (lognorm!)=[])
+function LinearAlgebra.normalize!(M::AbstractMPS; (lognorm!)=[])
   c = ortho_lims(M)
   lognorm_M = lognorm(M)
   push!(lognorm!, lognorm_M)
diff --git a/src/ITensorMPS/abstractprojmpo/projmpo_apply.jl b/src/ITensorMPS/abstractprojmpo/projmpo_apply.jl
new file mode 100644
index 0000000000..9134aad563
--- /dev/null
+++ b/src/ITensorMPS/abstractprojmpo/projmpo_apply.jl
@@ -0,0 +1,89 @@
+using ITensors: ITensor
+
+"""
+A ProjMPOApply represents the application of an
+MPO `H` onto an MPS `psi0` but "projected" by
+the basis of a different MPS `psi` (which
+could be an approximation to H|psi>).
+
+As an implementation of the AbstractProjMPO
+type, it supports multiple `nsite` values for
+one- and two-site algorithms.
+
+```
+     *--*--*-      -*--*--*--*--*--* <psi|
+     |  |  |  |  |  |  |  |  |  |  |
+     h--h--h--h--h--h--h--h--h--h--h H  
+     |  |  |  |  |  |  |  |  |  |  |
+     o--o--o-      -o--o--o--o--o--o |psi0>
+```
+"""
+mutable struct ProjMPOApply <: AbstractProjMPO
+  lpos::Int
+  rpos::Int
+  nsite::Int
+  psi0::MPS
+  H::MPO
+  LR::Vector{ITensor}
+end
+
+function ProjMPOApply(psi0::MPS, H::MPO)
+  return ProjMPOApply(0, length(H) + 1, 2, psi0, H, Vector{ITensor}(undef, length(H)))
+end
+
+function Base.copy(P::ProjMPOApply)
+  return ProjMPOApply(P.lpos, P.rpos, P.nsite, copy(P.psi0), copy(P.H), copy(P.LR))
+end
+
+function set_nsite!(P::ProjMPOApply, nsite)
+  P.nsite = nsite
+  return P
+end
+
+function makeL!(P::ProjMPOApply, psi::MPS, k::Int)
+  # Save the last `L` that is made to help with caching
+  # for DiskProjMPO
+  ll = P.lpos
+  if ll ≥ k
+    # Special case when nothing has to be done.
+    # Still need to change the position if lproj is
+    # being moved backward.
+    P.lpos = k
+    return nothing
+  end
+  # Make sure ll is at least 0 for the generic logic below
+  ll = max(ll, 0)
+  L = lproj(P)
+  while ll < k
+    L = L * P.psi0[ll + 1] * P.H[ll + 1] * dag(psi[ll + 1])
+    P.LR[ll + 1] = L
+    ll += 1
+  end
+  # Needed when moving lproj backward.
+  P.lpos = k
+  return P
+end
+
+function makeR!(P::ProjMPOApply, psi::MPS, k::Int)
+  # Save the last `R` that is made to help with caching
+  # for DiskProjMPO
+  rl = P.rpos
+  if rl ≤ k
+    # Special case when nothing has to be done.
+    # Still need to change the position if rproj is
+    # being moved backward.
+    P.rpos = k
+    return nothing
+  end
+  N = length(P.H)
+  # Make sure rl is no bigger than `N + 1` for the generic logic below
+  rl = min(rl, N + 1)
+  R = rproj(P)
+  while rl > k
+    R = R * P.psi0[rl - 1] * P.H[rl - 1] * dag(psi[rl - 1])
+    P.LR[rl - 1] = R
+    rl -= 1
+  end
+  P.rpos = k
+  return P
+end
diff --git a/src/ITensorMPS/abstractprojmpo/projmpo_mps2.jl b/src/ITensorMPS/abstractprojmpo/projmpo_mps2.jl
new file mode 100644
index 0000000000..c051819304
--- /dev/null
+++ b/src/ITensorMPS/abstractprojmpo/projmpo_mps2.jl
@@ -0,0 +1,46 @@
+using ITensors: ITensors, contract
+
+mutable struct ProjMPO_MPS2 <: AbstractProjMPO
+  PH::ProjMPO
+  Ms::Vector{ProjMPS2}
+end
+
+function ProjMPO_MPS2(H::MPO, M::MPS)
+  return ProjMPO_MPS2(ProjMPO(H), [ProjMPS2(M)])
+end
+
+function ProjMPO_MPS2(H::MPO, Mv::Vector{MPS})
+  return ProjMPO_MPS2(ProjMPO(H), [ProjMPS2(m) for m in Mv])
+end
+
+Base.copy(P::ProjMPO_MPS2) = ProjMPO_MPS2(copy(P.PH), copy(P.Ms))
+
+nsite(P::ProjMPO_MPS2) = nsite(P.PH)
+
+function set_nsite!(P::ProjMPO_MPS2, nsite)
+  set_nsite!(P.PH, nsite)
+  for m in P.Ms
+    set_nsite!(m, nsite)
+  end
+  return P
+end
+
+function makeL!(P::ProjMPO_MPS2, psi::MPS, k::Int)
+  makeL!(P.PH, psi, k)
+  for m in P.Ms
+    makeL!(m, psi, k)
+  end
+  return P
+end
+
+function makeR!(P::ProjMPO_MPS2, psi::MPS, k::Int)
+  makeR!(P.PH, psi, k)
+  for m in P.Ms
+    makeR!(m, psi, k)
+  end
+  return P
+end
+
+ITensors.contract(P::ProjMPO_MPS2, v::ITensor) = contract(P.PH, v)
+
+proj_mps(P::ProjMPO_MPS2) = [proj_mps(m) for m in P.Ms]
diff --git a/src/ITensorMPS/abstractprojmpo/projmps2.jl b/src/ITensorMPS/abstractprojmpo/projmps2.jl
new file mode 100644
index 0000000000..12a25f6cae
--- /dev/null
+++ b/src/ITensorMPS/abstractprojmpo/projmps2.jl
@@ -0,0 +1,122 @@
+using ITensors: ITensors, ITensor, dag, dim, prime
+
+"""
+Holds the following data where psi
+is the MPS being optimized and M is the 
+MPS held constant by the ProjMPS.
+```
+     o--o--o--o--o--o--o--o--o--o--o <M|
+     |  |  |  |  |  |  |  |  |  |  |
+     *--*--*-      -*--*--*--*--*--* |psi>
+```
+"""
+mutable struct ProjMPS2 <: AbstractProjMPO
+  lpos::Int
+  rpos::Int
+  nsite::Int
+  M::MPS
+  LR::Vector{ITensor}
+end
+
+function ProjMPS2(M::MPS)
+  return ProjMPS2(0, length(M) + 1, 2, M, Vector{ITensor}(undef, length(M)))
+end
+
+Base.length(P::ProjMPS2) = length(P.M)
+
+function Base.copy(P::ProjMPS2)
+  return ProjMPS2(P.lpos, P.rpos, P.nsite, copy(P.M), copy(P.LR))
+end
+
+function set_nsite!(P::ProjMPS2, nsite)
+  P.nsite = nsite
+  return P
+end
+
+function makeL!(P::ProjMPS2, psi::MPS, k::Int)
+  # Save the last `L` that is made to help with caching
+  # for DiskProjMPO
+  ll = P.lpos
+  if ll ≥ k
+    # Special case when nothing has to be done.
+    # Still need to change the position if lproj is
+    # being moved backward.
+    P.lpos = k
+    return nothing
+  end
+  # Make sure ll is at least 0 for the generic logic below
+  ll = max(ll, 0)
+  L = lproj(P)
+  while ll < k
+    L = L * psi[ll + 1] * dag(prime(P.M[ll + 1], "Link"))
+    P.LR[ll + 1] = L
+    ll += 1
+  end
+  # Needed when moving lproj backward.
+  P.lpos = k
+  return P
+end
+
+function makeR!(P::ProjMPS2, psi::MPS, k::Int)
+  # Save the last `R` that is made to help with caching
+  # for DiskProjMPO
+  rl = P.rpos
+  if rl ≤ k
+    # Special case when nothing has to be done.
+    # Still need to change the position if rproj is
+    # being moved backward.
+    P.rpos = k
+    return nothing
+  end
+  N = length(P.M)
+  # Make sure rl is no bigger than `N + 1` for the generic logic below
+  rl = min(rl, N + 1)
+  R = rproj(P)
+  while rl > k
+    R = R * psi[rl - 1] * dag(prime(P.M[rl - 1], "Link"))
+    P.LR[rl - 1] = R
+    rl -= 1
+  end
+  P.rpos = k
+  return P
+end
+
+function ITensors.contract(P::ProjMPS2, v::ITensor)
+  itensor_map = Union{ITensor,OneITensor}[lproj(P)]
+  append!(itensor_map, [prime(t, "Link") for t in P.M[site_range(P)]])
+  push!(itensor_map, rproj(P))
+
+  # Reverse the contraction order of the map if
+  # the first tensor is a scalar (for example we
+  # are at the left edge of the system)
+  if dim(first(itensor_map)) == 1
+    reverse!(itensor_map)
+  end
+
+  # Apply the map
+  Mv = v
+  for it in itensor_map
+    Mv *= it
+  end
+  return Mv
+end
+
+function proj_mps(P::ProjMPS2)
+  itensor_map = Union{ITensor,OneITensor}[lproj(P)]
+  append!(itensor_map, [dag(prime(t, "Link")) for t in P.M[site_range(P)]])
+  push!(itensor_map, rproj(P))
+
+  # Reverse the contraction order of the map if
+  # the first tensor is a scalar (for example we
+  # are at the left edge of the system)
+  if dim(first(itensor_map)) == 1
+    reverse!(itensor_map)
+  end
+
+  # Apply the map
+  m = ITensor(true)
+  for it in itensor_map
+    m *= it
+  end
+  return m
+end
diff --git a/src/ITensorMPS/adapt.jl b/src/ITensorMPS/adapt.jl
index 9c0192fac8..b73ab27bd5 100644
--- a/src/ITensorMPS/adapt.jl
+++ b/src/ITensorMPS/adapt.jl
@@ -1 +1,2 @@
-adapt_structure(to, x::Union{MPS,MPO}) = map(xᵢ -> adapt(to, xᵢ), x)
+using Adapt: Adapt
+Adapt.adapt_structure(to, x::Union{MPS,MPO}) = map(xᵢ -> adapt(to, xᵢ), x)
diff --git a/src/ITensorMPS/alternating_update.jl b/src/ITensorMPS/alternating_update.jl
new file mode 100644
index 0000000000..3770f8b86e
--- /dev/null
+++ b/src/ITensorMPS/alternating_update.jl
@@ -0,0 +1,147 @@
+using ITensors: permute
+using NDTensors: scalartype
+using Printf: @printf
+
+function _compute_nsweeps(t; time_step=default_time_step(t), nsweeps=default_nsweeps())
+  if isinf(t) && isnothing(nsweeps)
+    nsweeps = 1
+  elseif !isnothing(nsweeps) && time_step != t
+    error("Cannot specify both time_step and nsweeps in alternating_update")
+  elseif isfinite(time_step) && abs(time_step) > 0 && isnothing(nsweeps)
+    nsweeps = convert(Int, ceil(abs(t / time_step)))
+    if !(nsweeps * time_step ≈ t)
+      error("Time step $time_step not commensurate with total time t=$t")
+    end
+  end
+  return nsweeps
+end
+
+function _extend_sweeps_param(param, nsweeps)
+  if param isa Number
+    eparam = fill(param, nsweeps)
+  else
+    length(param) == nsweeps && return param
+    eparam = Vector(undef, nsweeps)
+    eparam[1:length(param)] = param
+    eparam[(length(param) + 1):end] .= param[end]
+  end
+  return eparam
+end
+
+function process_sweeps(; nsweeps, maxdim, mindim, cutoff, noise)
+  maxdim = _extend_sweeps_param(maxdim, nsweeps)
+  mindim = _extend_sweeps_param(mindim, nsweeps)
+  cutoff = _extend_sweeps_param(cutoff, nsweeps)
+  noise = _extend_sweeps_param(noise, nsweeps)
+  return (; maxdim, mindim, cutoff, noise)
+end
+
+function alternating_update(
+  solver,
+  PH,
+  t::Number,
+  psi0::MPS;
+  nsweeps=default_nsweeps(),
+  checkdone=default_checkdone(),
+  write_when_maxdim_exceeds=default_write_when_maxdim_exceeds(),
+  nsite=default_nsite(),
+  reverse_step=default_reverse_step(),
+  time_start=default_time_start(),
+  time_step=default_time_step(t),
+  order=default_order(),
+  (observer!)=default_observer!(),
+  (step_observer!)=default_step_observer!(),
+  outputlevel=default_outputlevel(),
+  normalize=default_normalize(),
+  maxdim=default_maxdim(),
+  mindim=default_mindim(),
+  cutoff=default_cutoff(Float64),
+  noise=default_noise(),
+)
+  nsweeps = _compute_nsweeps(t; time_step, nsweeps)
+  maxdim, mindim, cutoff, noise = process_sweeps(; nsweeps, maxdim, mindim, cutoff, noise)
+  forward_order = TDVPOrder(order, Base.Forward)
+  psi = copy(psi0)
+  # Keep track of the start of the current time step.
+  # Helpful for tracking the total time, for example
+  # when using time-dependent solvers.
+  # This will be passed as a keyword argument to the
+  # `solver`.
+  current_time = time_start
+  info = nothing
+  for sweep in 1:nsweeps
+    if !isnothing(write_when_maxdim_exceeds) && maxdim[sweep] > write_when_maxdim_exceeds
+      if outputlevel >= 2
+        println(
+          "write_when_maxdim_exceeds = $write_when_maxdim_exceeds and maxdim(sweeps, sw) = $(maxdim(sweeps, sweep)), writing environment tensors to disk",
+        )
+      end
+      PH = disk(PH)
+    end
+    sweep_time = @elapsed begin
+      psi, PH, info = sweep_update(
+        forward_order,
+        solver,
+        PH,
+        time_step,
+        psi;
+        nsite,
+        current_time,
+        reverse_step,
+        sweep,
+        observer!,
+        normalize,
+        maxdim=maxdim[sweep],
+        mindim=mindim[sweep],
+        cutoff=cutoff[sweep],
+        noise=noise[sweep],
+      )
+    end
+    current_time += time_step
+    update!(step_observer!; psi, sweep, outputlevel, current_time)
+    if outputlevel >= 1
+      print("After sweep ", sweep, ":")
+      print(" maxlinkdim=", maxlinkdim(psi))
+      @printf(" maxerr=%.2E", info.maxtruncerr)
+      print(" current_time=", round(current_time; digits=3))
+      print(" time=", round(sweep_time; digits=3))
+      println()
+      flush(stdout)
+    end
+    isdone = false
+    if !isnothing(checkdone)
+      isdone = checkdone(; psi, sweep, outputlevel)
+    elseif observer! isa AbstractObserver
+      isdone = checkdone!(observer!; psi, sweep, outputlevel)
+    end
+    isdone && break
+  end
+  return psi
+end
+
+# Convenience wrapper to not have to specify time step.
+# Use a time step of `Inf` as a convention, since TDVP
+# with an infinite time step corresponds to DMRG.
+function alternating_update(solver, H, psi0::MPS; kwargs...)
+  return alternating_update(solver, H, scalartype(psi0)(Inf), psi0; kwargs...)
+end
+
+function alternating_update(solver, H::MPO, t::Number, psi0::MPS; kwargs...)
+  check_hascommoninds(siteinds, H, psi0)
+  check_hascommoninds(siteinds, H, psi0')
+  # Permute the indices to have a better memory layout
+  # and minimize permutations
+  H = permute(H, (linkind, siteinds, linkind))
+  PH = ProjMPO(H)
+  return alternating_update(solver, PH, t, psi0; kwargs...)
+end
+
+function alternating_update(solver, Hs::Vector{MPO}, t::Number, psi0::MPS; kwargs...)
+  for H in Hs
+    check_hascommoninds(siteinds, H, psi0)
+    check_hascommoninds(siteinds, H, psi0')
+  end
+  Hs .= ITensors.permute.(Hs, Ref((linkind, siteinds, linkind)))
+  PHs = ProjMPOSum(Hs)
+  return alternating_update(solver, PHs, t, psi0; kwargs...)
+end
diff --git a/src/ITensorMPS/autompo/opsum_to_mpo.jl b/src/ITensorMPS/autompo/opsum_to_mpo.jl
index ad63aa2151..114eec3bb8 100644
--- a/src/ITensorMPS/autompo/opsum_to_mpo.jl
+++ b/src/ITensorMPS/autompo/opsum_to_mpo.jl
@@ -1,3 +1,5 @@
+using NDTensors: using_auto_fermion
+
 # `ValType::Type{<:Number}` is used instead of `ValType::Type` for efficiency, possibly due to increased method specialization.
 # See https://github.com/ITensor/ITensors.jl/pull/1183.
 function svdMPO(
diff --git a/src/ITensorMPS/autompo/opsum_to_mpo_generic.jl b/src/ITensorMPS/autompo/opsum_to_mpo_generic.jl
index 0cd310993c..7a855ea213 100644
--- a/src/ITensorMPS/autompo/opsum_to_mpo_generic.jl
+++ b/src/ITensorMPS/autompo/opsum_to_mpo_generic.jl
@@ -1,3 +1,5 @@
+using NDTensors: using_auto_fermion
+
 # TODO: Deprecate.
 const AutoMPO = OpSum
 
diff --git a/src/ITensorMPS/autompo/opsum_to_mpo_qn.jl b/src/ITensorMPS/autompo/opsum_to_mpo_qn.jl
index 0149ea7634..25f75b30ba 100644
--- a/src/ITensorMPS/autompo/opsum_to_mpo_qn.jl
+++ b/src/ITensorMPS/autompo/opsum_to_mpo_qn.jl
@@ -1,3 +1,5 @@
+using NDTensors: using_auto_fermion
+
 # `ValType::Type{<:Number}` is used instead of `ValType::Type` for efficiency, possibly due to increased method specialization.
 # See https://github.com/ITensor/ITensors.jl/pull/1183.
 function qn_svdMPO(
diff --git a/src/ITensorMPS/contract_mpo_mps.jl b/src/ITensorMPS/contract_mpo_mps.jl
new file mode 100644
index 0000000000..4156328fc8
--- /dev/null
+++ b/src/ITensorMPS/contract_mpo_mps.jl
@@ -0,0 +1,48 @@
+using ITensors: ITensors, Index, ITensor, @Algorithm_str, commoninds, contract, hasind, sim
+
+function contractmpo_solver(; kwargs...)
+  function solver(PH, t, psi; kws...)
+    v = ITensor(true)
+    for j in (PH.lpos + 1):(PH.rpos - 1)
+      v *= PH.psi0[j]
+    end
+    Hpsi0 = contract(PH, v)
+    return Hpsi0, nothing
+  end
+  return solver
+end
+
+function ITensors.contract(
+  ::Algorithm"fit", A::MPO, psi0::MPS; init_mps=psi0, nsweeps=1, kwargs...
+)::MPS
+  n = length(A)
+  n != length(psi0) &&
+    throw(DimensionMismatch("lengths of MPO ($n) and MPS ($(length(psi0))) do not match"))
+  if n == 1
+    return MPS([A[1] * psi0[1]])
+  end
+  any(i -> isempty(i), siteinds(commoninds, A, psi0)) &&
+    error("In `contract(A::MPO, x::MPS)`, `A` and `x` must share a set of site indices")
+  # In case A and psi0 have the same link indices
+  A = sim(linkinds, A)
+  # Fix site and link inds of init_mps
+  init_mps = deepcopy(init_mps)
+  init_mps = sim(linkinds, init_mps)
+  Ai = siteinds(A)
+  ti = Vector{Index}(undef, n)
+  for j in 1:n
+    for i in Ai[j]
+      if !hasind(psi0[j], i)
+        ti[j] = i
+        break
+      end
+    end
+  end
+  replace_siteinds!(init_mps, ti)
+  reverse_step = false
+  PH = ProjMPOApply(psi0, A)
+  psi = alternating_update(
+    contractmpo_solver(; kwargs...), PH, init_mps; nsweeps, reverse_step, kwargs...
+  )
+  return psi
+end
diff --git a/src/ITensorMPS/defaults.jl b/src/ITensorMPS/defaults.jl
new file mode 100644
index 0000000000..72ffa408f9
--- /dev/null
+++ b/src/ITensorMPS/defaults.jl
@@ -0,0 +1,49 @@
+using KrylovKit: eigsolve, exponentiate
+
+default_nsweeps() = nothing
+default_checkdone() = nothing
+default_write_when_maxdim_exceeds() = nothing
+default_nsite() = 2
+default_reverse_step() = true
+default_time_start() = 0
+default_time_step(t) = t
+default_order() = 2
+default_observer!() = NoObserver()
+default_step_observer!() = NoObserver()
+default_outputlevel() = 0
+default_normalize() = false
+default_sweep() = 1
+default_current_time() = 0
+
+# Truncation
+default_maxdim() = typemax(Int)
+default_mindim() = 1
+default_cutoff(type::Type{<:Number}) = eps(real(type))
+default_noise() = 0
+
+# Solvers
+default_tdvp_solver_backend() = "exponentiate"
+default_ishermitian() = true
+default_issymmetric() = true
+default_solver_verbosity() = 0
+
+# Customizable based on solver function
+default_solver_outputlevel(::Function) = 0
+
+default_solver_tol(::Function) = error("Not implemented")
+default_solver_which_eigenvalue(::Function) = error("Not implemented")
+default_solver_krylovdim(::Function) = error("Not implemented")
+default_solver_verbosity(::Function) = error("Not implemented")
+
+## Solver-specific keyword argument defaults
+
+# dmrg/eigsolve
+default_solver_tol(::typeof(eigsolve)) = 1e-14
+default_solver_krylovdim(::typeof(eigsolve)) = 3
+default_solver_maxiter(::typeof(eigsolve)) = 1
+default_solver_which_eigenvalue(::typeof(eigsolve)) = :SR
+
+# tdvp/exponentiate
+default_solver_tol(::typeof(exponentiate)) = 1e-12
+default_solver_krylovdim(::typeof(exponentiate)) = 30
+default_solver_maxiter(::typeof(exponentiate)) = 100
diff --git a/src/ITensorMPS/dmrg.jl b/src/ITensorMPS/dmrg.jl
index 37885fcbb4..78d2ff79f8 100644
--- a/src/ITensorMPS/dmrg.jl
+++ b/src/ITensorMPS/dmrg.jl
@@ -1,4 +1,8 @@
+using Adapt: adapt
 using KrylovKit: eigsolve
+using NDTensors: scalartype, timer
+using Printf: @printf
+using TupleTools: TupleTools
 
 function permute(
   M::AbstractMPS, ::Tuple{typeof(linkind),typeof(siteinds),typeof(linkind)}
@@ -8,7 +12,7 @@ function permute(
     lₙ₋₁ = linkind(M, n - 1)
     lₙ = linkind(M, n)
     s⃗ₙ = TupleTools.sort(Tuple(siteinds(M, n)); by=plev)
-    M̃[n] = permute(M[n], filter(!isnothing, (lₙ₋₁, s⃗ₙ..., lₙ)))
+    M̃[n] = ITensors.permute(M[n], filter(!isnothing, (lₙ₋₁, s⃗ₙ..., lₙ)))
   end
   set_ortho_lims!(M̃, ortho_lims(M))
   return M̃
@@ -340,16 +344,11 @@ function dmrg(
   return (energy, psi)
 end
 
-default_maxdim() = typemax(Int)
-default_mindim() = 1
-default_cutoff() = 1e-8
-default_noise() = false
-
 function _dmrg_sweeps(;
   nsweeps,
   maxdim=default_maxdim(),
   mindim=default_mindim(),
-  cutoff=default_cutoff(),
+  cutoff=default_cutoff(Float64),
   noise=default_noise(),
 )
   sweeps = Sweeps(nsweeps)
@@ -367,7 +366,7 @@ function dmrg(
   nsweeps,
   maxdim=default_maxdim(),
   mindim=default_mindim(),
-  cutoff=default_cutoff(),
+  cutoff=default_cutoff(Float64),
   noise=default_noise(),
   kwargs...,
 )
@@ -382,9 +381,70 @@ function dmrg(
   nsweeps,
   maxdim=default_maxdim(),
   mindim=default_mindim(),
-  cutoff=default_cutoff(),
+  cutoff=default_cutoff(Float64),
   noise=default_noise(),
   kwargs...,
 )
   return dmrg(x1, psi0, _dmrg_sweeps(; nsweeps, maxdim, mindim, cutoff, noise); kwargs...)
 end
+
+# Implementation of DMRG originally from ITensorTDVP package
+
+function dmrg_solver(
+  f::typeof(eigsolve);
+  solver_which_eigenvalue,
+  ishermitian,
+  solver_tol,
+  solver_krylovdim,
+  solver_maxiter,
+  solver_verbosity,
+)
+  function solver(H, t, psi0; current_time, outputlevel)
+    howmany = 1
+    which = solver_which_eigenvalue
+    vals, vecs, info = f(
+      H,
+      psi0,
+      howmany,
+      which;
+      ishermitian=default_ishermitian(),
+      tol=solver_tol,
+      krylovdim=solver_krylovdim,
+      maxiter=solver_maxiter,
+      verbosity=solver_verbosity,
+    )
+    psi = vecs[1]
+    return psi, info
+  end
+  return solver
+end
+
+function alternating_update_dmrg(
+  H,
+  psi0::MPS;
+  ishermitian=default_ishermitian(),
+  solver_which_eigenvalue=default_solver_which_eigenvalue(eigsolve),
+  solver_tol=default_solver_tol(eigsolve),
+  solver_krylovdim=default_solver_krylovdim(eigsolve),
+  solver_maxiter=default_solver_maxiter(eigsolve),
+  solver_verbosity=default_solver_verbosity(),
+  kwargs...,
+)
+  reverse_step = false
+  psi = alternating_update(
+    dmrg_solver(
+      eigsolve;
+      solver_which_eigenvalue,
+      ishermitian,
+      solver_tol,
+      solver_krylovdim,
+      solver_maxiter,
+      solver_verbosity,
+    ),
+    H,
+    psi0;
+    reverse_step,
+    kwargs...,
+  )
+  return psi
+end
diff --git a/src/ITensorMPS/dmrg_x.jl b/src/ITensorMPS/dmrg_x.jl
new file mode 100644
index 0000000000..50911599a4
--- /dev/null
+++ b/src/ITensorMPS/dmrg_x.jl
@@ -0,0 +1,16 @@
+using ITensors: array, contract, dag, uniqueind, onehot
+using LinearAlgebra: eigen
+
+function dmrg_x_solver(PH, t, psi0; current_time, outputlevel)
+  H = contract(PH, ITensor(true))
+  D, U = eigen(H; ishermitian=true)
+  u = uniqueind(U, H)
+  max_overlap, max_ind = findmax(abs, array(psi0 * dag(U)))
+  U_max = U * dag(onehot(eltype(U), u => max_ind))
+  return U_max, nothing
+end
+
+function dmrg_x(PH, psi0::MPS; reverse_step=false, kwargs...)
+  psi = alternating_update(dmrg_x_solver, PH, psi0; reverse_step, kwargs...)
+  return psi
+end
diff --git a/src/ITensorMPS/imports.jl b/src/ITensorMPS/imports.jl
index 937f431ba5..8bb438cfe2 100644
--- a/src/ITensorMPS/imports.jl
+++ b/src/ITensorMPS/imports.jl
@@ -80,10 +80,6 @@ import Base.Broadcast:
   broadcastable,
   instantiate
 
-import Adapt: adapt_structure, adapt_storage
-
-import LinearAlgebra: dot, normalize!, tr
-
 import ..ITensors.NDTensors:
   Algorithm,
   @Algorithm_str,
@@ -102,12 +98,7 @@ import ..ITensors.NDTensors:
   fill!!,
   randn!!,
   permutedims,
-  permutedims!,
-  scalartype,
-  single_precision,
-  tensor,
-  timer,
-  using_auto_fermion
+  permutedims!
 
 import ..ITensors:
   AbstractRNG,
diff --git a/src/ITensorMPS/linsolve.jl b/src/ITensorMPS/linsolve.jl
new file mode 100644
index 0000000000..2fa5d8869a
--- /dev/null
+++ b/src/ITensorMPS/linsolve.jl
@@ -0,0 +1,41 @@
+using KrylovKit: KrylovKit, linsolve
+
+"""
+Compute a solution x to the linear system:
+
+(a₀ + a₁ * A)*x = b
+
+using starting guess x₀. Leaving a₀, a₁
+set to their default values solves the 
+system A*x = b.
+
+To adjust the balance between accuracy of solution
+and speed of the algorithm, it is recommed to first try
+adjusting the solver keyword arguments as descibed below.
+
+Keyword arguments:
+  - `nsweeps`, `cutoff`, `maxdim`, etc. (like for other MPO/MPS solvers).
+  - `solver_kwargs=(;)` - a `NamedTuple` containing keyword arguments that will get forwarded to the local solver,
+    in this case `KrylovKit.linsolve` which is a GMRES linear solver. For example:
+    ```juli
+    linsolve(A, b, x; maxdim=100, cutoff=1e-8, nsweeps=10, solver_kwargs=(; ishermitian=true, tol=1e-6, maxiter=20, krylovdim=30))
+    ```
+    See `KrylovKit.jl` documentation for more details on available keyword arguments.
+"""
+function KrylovKit.linsolve(
+  A::MPO,
+  b::MPS,
+  x₀::MPS,
+  a₀::Number=false,
+  a₁::Number=true;
+  solver_kwargs=(;),
+  tdvp_kwargs...,
+)
+  function linsolve_solver(P::ProjMPO_MPS2, t, x₀; current_time, outputlevel)
+    b = dag(only(proj_mps(P)))
+    x, info = KrylovKit.linsolve(P, b, x₀, a₀, a₁; solver_kwargs...)
+    return x, nothing
+  end
+  P = ProjMPO_MPS2(A, b)
+  return alternating_update(linsolve_solver, P, x₀; reverse_step=false, tdvp_kwargs...)
+end
diff --git a/src/ITensorMPS/mpo.jl b/src/ITensorMPS/mpo.jl
index 9be5974cfe..215d882363 100644
--- a/src/ITensorMPS/mpo.jl
+++ b/src/ITensorMPS/mpo.jl
@@ -1,3 +1,6 @@
+using Adapt: adapt
+using HDF5: attributes, create_group, open_group
+using Random: Random
 
 """
     MPO
@@ -398,7 +401,7 @@ end
 
 Same as [`inner`](@ref).
 """
-function dot(y::MPS, A::MPO, x::MPS; make_inds_match::Bool=true, kwargs...)
+function LinearAlgebra.dot(y::MPS, A::MPO, x::MPS; make_inds_match::Bool=true, kwargs...)
   return _log_or_not_dot(y, A, x, false; make_inds_match=make_inds_match, kwargs...)
 end
 
@@ -457,7 +460,9 @@ loginner(y::MPS, A::MPO, x::MPS; kwargs...) = logdot(y, A, x; kwargs...)
 
 Same as [`inner`](@ref).
 """
-function dot(B::MPO, y::MPS, A::MPO, x::MPS; make_inds_match::Bool=true, kwargs...)::Number
+function LinearAlgebra.dot(
+  B::MPO, y::MPS, A::MPO, x::MPS; make_inds_match::Bool=true, kwargs...
+)::Number
   !make_inds_match && error(
     "make_inds_match = false not currently supported in dot(::MPO, ::MPS, ::MPO, ::MPS)"
   )
@@ -515,7 +520,7 @@ Same as [`dot`](@ref).
 """
 inner(B::MPO, y::MPS, A::MPO, x::MPS) = dot(B, y, A, x)
 
-function dot(M1::MPO, M2::MPO; make_inds_match::Bool=false, kwargs...)
+function LinearAlgebra.dot(M1::MPO, M2::MPO; make_inds_match::Bool=false, kwargs...)
   if make_inds_match
     error("In dot(::MPO, ::MPO), make_inds_match is not currently supported")
   end
@@ -531,7 +536,7 @@ function logdot(M1::MPO, M2::MPO; make_inds_match::Bool=false, kwargs...)
   return _log_or_not_dot(M1, M2, true; make_inds_match=make_inds_match)
 end
 
-function tr(M::MPO; plev::Pair{Int,Int}=0 => 1, tags::Pair=ts"" => ts"")
+function LinearAlgebra.tr(M::MPO; plev::Pair{Int,Int}=0 => 1, tags::Pair=ts"" => ts"")
   N = length(M)
   #
   # TODO: choose whether to contract or trace
diff --git a/src/ITensorMPS/mps.jl b/src/ITensorMPS/mps.jl
index 612f531797..d5bb6a1c06 100644
--- a/src/ITensorMPS/mps.jl
+++ b/src/ITensorMPS/mps.jl
@@ -1,3 +1,7 @@
+using Adapt: adapt
+using HDF5: attributes, create_group, open_group
+using NDTensors: using_auto_fermion
+using Random: Random
 
 """
     MPS
diff --git a/src/ITensorMPS/solver_utils.jl b/src/ITensorMPS/solver_utils.jl
new file mode 100644
index 0000000000..6372e56473
--- /dev/null
+++ b/src/ITensorMPS/solver_utils.jl
@@ -0,0 +1,69 @@
+using ITensors: ITensor, apply, array, inds, itensor, permute
+
+# Utilities for making it easier
+# to define solvers (like ODE solvers)
+# for TDVP
+
+"""
+    to_vec(x)
+Transform `x` into a `Vector`, and return the vector, and a closure which inverts the
+transformation.
+
+Modeled after FiniteDifferences.to_vec:
+
+https://github.com/JuliaDiff/FiniteDifferences.jl/blob/main/src/to_vec.jl
+"""
+to_vec(x) = error("Not implemented")
+
+function to_vec(x::ITensor)
+  function to_itensor(x_vec)
+    return itensor(x_vec, inds(x))
+  end
+  return vec(array(x)), to_itensor
+end
+
+# Represents a time-dependent sum of terms:
+#
+# H(t) = f[1](t) * H0[1] + f[2](t) * H0[2] + …
+#
+struct TimeDependentSum{S,T}
+  f::Vector{S}
+  H0::T
+end
+TimeDependentSum(f::Vector, H0::ProjMPOSum) = TimeDependentSum(f, ITensors.terms(H0))
+Base.length(H::TimeDependentSum) = length(H.f)
+
+function Base.:*(c::Number, H::TimeDependentSum)
+  return TimeDependentSum([t -> c * fₙ(t) for fₙ in H.f], H.H0)
+end
+Base.:*(H::TimeDependentSum, c::Number) = c * H
+
+# Calling a `TimeDependentOpSum` at a certain time like:
+#
+# H(t)
+#
+# Returns a `ScaledSum` at that time.
+(H::TimeDependentSum)(t::Number) = ScaledSum([fₙ(t) for fₙ in H.f], H.H0)
+
+# Represents the sum of scaled terms:
+#
+# H = coefficient[1] * H[1] + coefficient * H[2] + …
+#
+struct ScaledSum{S,T}
+  coefficients::Vector{S}
+  H::T
+end
+Base.length(H::ScaledSum) = length(H.coefficients)
+
+# Apply the scaled sum of terms:
+#
+# H(ψ₀) = coefficient[1] * H[1](ψ₀) + coefficient[2] * H[2](ψ₀) + …
+#
+# onto ψ₀.
+function (H::ScaledSum)(ψ₀)
+  ψ = ITensor(inds(ψ₀))
+  for n in 1:length(H)
+    ψ += H.coefficients[n] * apply(H.H[n], ψ₀)
+  end
+  return permute(ψ, inds(ψ₀))
+end
diff --git a/src/ITensorMPS/sweep_update.jl b/src/ITensorMPS/sweep_update.jl
new file mode 100644
index 0000000000..3679e7a748
--- /dev/null
+++ b/src/ITensorMPS/sweep_update.jl
@@ -0,0 +1,434 @@
+using ITensors: uniqueinds
+using ITensors.ITensorMPS:
+  ITensorMPS, MPS, isortho, orthocenter, orthogonalize!, position!, replacebond!, set_nsite!
+using LinearAlgebra: norm, normalize!, svd
+using Observers: update!
+using Printf: @printf
+
+function sweep_update(
+  order::TDVPOrder, solver, PH, time_step::Number, psi::MPS; current_time=0, kwargs...
+)
+  order_orderings = orderings(order)
+  order_sub_time_steps = eltype(time_step).(sub_time_steps(order))
+  order_sub_time_steps *= time_step
+  info = nothing
+  for substep in 1:length(order_sub_time_steps)
+    psi, PH, info = sub_sweep_update(
+      order_orderings[substep],
+      solver,
+      PH,
+      order_sub_time_steps[substep],
+      psi;
+      current_time,
+      kwargs...,
+    )
+    current_time += order_sub_time_steps[substep]
+  end
+  return psi, PH, info
+end
+
+isforward(direction::Base.ForwardOrdering) = true
+isforward(direction::Base.ReverseOrdering) = false
+isreverse(direction) = !isforward(direction)
+
+function sweep_bonds(direction::Base.ForwardOrdering, n::Int; ncenter::Int)
+  return 1:(n - ncenter + 1)
+end
+
+function sweep_bonds(direction::Base.ReverseOrdering, n::Int; ncenter::Int)
+  return reverse(sweep_bonds(Base.Forward, n; ncenter))
+end
+
+is_forward_done(direction::Base.ForwardOrdering, b, n; ncenter) = (b + ncenter - 1 == n)
+is_forward_done(direction::Base.ReverseOrdering, b, n; ncenter) = false
+is_reverse_done(direction::Base.ForwardOrdering, b, n; ncenter) = false
+is_reverse_done(direction::Base.ReverseOrdering, b, n; ncenter) = (b == 1)
+function is_half_sweep_done(direction, b, n; ncenter)
+  return is_forward_done(direction, b, n; ncenter) ||
+         is_reverse_done(direction, b, n; ncenter)
+end
+
+function sub_sweep_update(
+  direction::Base.Ordering,
+  solver,
+  PH,
+  time_step::Number,
+  psi::MPS;
+  which_decomp=nothing,
+  svd_alg=nothing,
+  sweep=default_sweep(),
+  current_time=default_current_time(),
+  nsite=default_nsite(),
+  reverse_step=default_reverse_step(),
+  normalize=default_normalize(),
+  (observer!)=default_observer!(),
+  outputlevel=default_outputlevel(),
+  maxdim=default_maxdim(),
+  mindim=default_mindim(),
+  cutoff=default_cutoff(time_step),
+  noise=default_noise(),
+)
+  PH = copy(PH)
+  psi = copy(psi)
+  if length(psi) == 1
+    error(
+      "`tdvp`, `dmrg`, `linsolve`, etc. currently does not support system sizes of 1. You can diagonalize the MPO tensor directly with tools like `LinearAlgebra.eigen`, `KrylovKit.exponentiate`, etc.",
+    )
+  end
+  N = length(psi)
+  set_nsite!(PH, nsite)
+  if isforward(direction)
+    if !isortho(psi) || orthocenter(psi) != 1
+      orthogonalize!(psi, 1)
+    end
+    @assert isortho(psi) && orthocenter(psi) == 1
+    position!(PH, psi, 1)
+  elseif isreverse(direction)
+    if !isortho(psi) || orthocenter(psi) != N - nsite + 1
+      orthogonalize!(psi, N - nsite + 1)
+    end
+    @assert(isortho(psi) && (orthocenter(psi) == N - nsite + 1))
+    position!(PH, psi, N - nsite + 1)
+  end
+  maxtruncerr = 0.0
+  info = nothing
+  for b in sweep_bonds(direction, N; ncenter=nsite)
+    current_time, maxtruncerr, spec, info = region_update!(
+      solver,
+      PH,
+      psi,
+      b;
+      nsite,
+      reverse_step,
+      current_time,
+      outputlevel,
+      time_step,
+      normalize,
+      direction,
+      noise,
+      which_decomp,
+      svd_alg,
+      cutoff,
+      maxdim,
+      mindim,
+      maxtruncerr,
+    )
+    if outputlevel >= 2
+      if nsite == 1
+        @printf("Sweep %d, direction %s, bond (%d,) \n", sweep, direction, b)
+      elseif nsite == 2
+        @printf("Sweep %d, direction %s, bond (%d,%d) \n", sweep, direction, b, b + 1)
+      end
+      print("  Truncated using")
+      @printf(" cutoff=%.1E", cutoff)
+      @printf(" maxdim=%.1E", maxdim)
+      print(" mindim=", mindim)
+      print(" current_time=", round(current_time; digits=3))
+      println()
+      if spec != nothing
+        @printf(
+          "  Trunc. err=%.2E, bond dimension %d\n", spec.truncerr, dim(linkind(psi, b))
+        )
+      end
+      flush(stdout)
+    end
+    update!(
+      observer!;
+      psi,
+      bond=b,
+      sweep,
+      half_sweep=isforward(direction) ? 1 : 2,
+      spec,
+      outputlevel,
+      half_sweep_is_done=is_half_sweep_done(direction, b, N; ncenter=nsite),
+      current_time,
+      info,
+    )
+  end
+  # Just to be sure:
+  normalize && normalize!(psi)
+  return psi, PH, TDVPInfo(maxtruncerr)
+end
+
+function region_update!(
+  solver,
+  PH,
+  psi,
+  b;
+  nsite,
+  reverse_step,
+  current_time,
+  outputlevel,
+  time_step,
+  normalize,
+  direction,
+  noise,
+  which_decomp,
+  svd_alg,
+  cutoff,
+  maxdim,
+  mindim,
+  maxtruncerr,
+)
+  return region_update!(
+    Val(nsite),
+    Val(reverse_step),
+    solver,
+    PH,
+    psi,
+    b;
+    current_time,
+    outputlevel,
+    time_step,
+    normalize,
+    direction,
+    noise,
+    which_decomp,
+    svd_alg,
+    cutoff,
+    maxdim,
+    mindim,
+    maxtruncerr,
+  )
+end
+
+function region_update!(
+  nsite_val::Val{1},
+  reverse_step_val::Val{false},
+  solver,
+  PH,
+  psi,
+  b;
+  current_time,
+  outputlevel,
+  time_step,
+  normalize,
+  direction,
+  noise,
+  which_decomp,
+  svd_alg,
+  cutoff,
+  maxdim,
+  mindim,
+  maxtruncerr,
+)
+  N = length(psi)
+  nsite = 1
+  # Do 'forwards' evolution step
+  set_nsite!(PH, nsite)
+  position!(PH, psi, b)
+  phi1 = psi[b]
+  phi1, info = solver(PH, time_step, phi1; current_time, outputlevel)
+  current_time += time_step
+  normalize && (phi1 /= norm(phi1))
+  spec = nothing
+  psi[b] = phi1
+  if !is_half_sweep_done(direction, b, N; ncenter=nsite)
+    # Move ortho center
+    Δ = (isforward(direction) ? +1 : -1)
+    orthogonalize!(psi, b + Δ)
+  end
+  return current_time, maxtruncerr, spec, info
+end
+
+function region_update!(
+  nsite_val::Val{1},
+  reverse_step_val::Val{true},
+  solver,
+  PH,
+  psi,
+  b;
+  current_time,
+  outputlevel,
+  time_step,
+  normalize,
+  direction,
+  noise,
+  which_decomp,
+  svd_alg,
+  cutoff,
+  maxdim,
+  mindim,
+  maxtruncerr,
+)
+  N = length(psi)
+  nsite = 1
+  # Do 'forwards' evolution step
+  set_nsite!(PH, nsite)
+  position!(PH, psi, b)
+  phi1 = psi[b]
+  phi1, info = solver(PH, time_step, phi1; current_time, outputlevel)
+  current_time += time_step
+  normalize && (phi1 /= norm(phi1))
+  spec = nothing
+  psi[b] = phi1
+  if !is_half_sweep_done(direction, b, N; ncenter=nsite)
+    # Do backwards evolution step
+    b1 = (isforward(direction) ? b + 1 : b)
+    Δ = (isforward(direction) ? +1 : -1)
+    uinds = uniqueinds(phi1, psi[b + Δ])
+    U, S, V = svd(phi1, uinds)
+    psi[b] = U
+    phi0 = S * V
+    if isforward(direction)
+      ITensorMPS.setleftlim!(psi, b)
+    elseif isreverse(direction)
+      ITensorMPS.setrightlim!(psi, b)
+    end
+    set_nsite!(PH, nsite - 1)
+    position!(PH, psi, b1)
+    phi0, info = solver(PH, -time_step, phi0; current_time, outputlevel)
+    current_time -= time_step
+    normalize && (phi0 ./= norm(phi0))
+    psi[b + Δ] = phi0 * psi[b + Δ]
+    if isforward(direction)
+      ITensorMPS.setrightlim!(psi, b + Δ + 1)
+    elseif isreverse(direction)
+      ITensorMPS.setleftlim!(psi, b + Δ - 1)
+    end
+    set_nsite!(PH, nsite)
+  end
+  return current_time, maxtruncerr, spec, info
+end
+
+function region_update!(
+  nsite_val::Val{2},
+  reverse_step_val::Val{false},
+  solver,
+  PH,
+  psi,
+  b;
+  current_time,
+  outputlevel,
+  time_step,
+  normalize,
+  direction,
+  noise,
+  which_decomp,
+  svd_alg,
+  cutoff,
+  maxdim,
+  mindim,
+  maxtruncerr,
+)
+  N = length(psi)
+  nsite = 2
+  # Do 'forwards' evolution step
+  set_nsite!(PH, nsite)
+  position!(PH, psi, b)
+  phi1 = psi[b] * psi[b + 1]
+  phi1, info = solver(PH, time_step, phi1; current_time, outputlevel)
+  current_time += time_step
+  normalize && (phi1 /= norm(phi1))
+  spec = nothing
+  ortho = isforward(direction) ? "left" : "right"
+  drho = nothing
+  if noise > 0.0 && isforward(direction)
+    drho = noise * noiseterm(PH, phi, ortho)
+  end
+  spec = replacebond!(
+    psi,
+    b,
+    phi1;
+    maxdim,
+    mindim,
+    cutoff,
+    eigen_perturbation=drho,
+    ortho=ortho,
+    normalize,
+    which_decomp,
+    svd_alg,
+  )
+  maxtruncerr = max(maxtruncerr, spec.truncerr)
+  return current_time, maxtruncerr, spec, info
+end
+
+function region_update!(
+  nsite_val::Val{2},
+  reverse_step_val::Val{true},
+  solver,
+  PH,
+  psi,
+  b;
+  current_time,
+  outputlevel,
+  time_step,
+  normalize,
+  direction,
+  noise,
+  which_decomp,
+  svd_alg,
+  cutoff,
+  maxdim,
+  mindim,
+  maxtruncerr,
+)
+  N = length(psi)
+  nsite = 2
+  # Do 'forwards' evolution step
+  set_nsite!(PH, nsite)
+  position!(PH, psi, b)
+  phi1 = psi[b] * psi[b + 1]
+  phi1, info = solver(PH, time_step, phi1; current_time, outputlevel)
+  current_time += time_step
+  normalize && (phi1 /= norm(phi1))
+  spec = nothing
+  ortho = isforward(direction) ? "left" : "right"
+  drho = nothing
+  if noise > 0.0 && isforward(direction)
+    drho = noise * noiseterm(PH, phi, ortho)
+  end
+  spec = replacebond!(
+    psi,
+    b,
+    phi1;
+    maxdim,
+    mindim,
+    cutoff,
+    eigen_perturbation=drho,
+    ortho=ortho,
+    normalize,
+    which_decomp,
+    svd_alg,
+  )
+  maxtruncerr = max(maxtruncerr, spec.truncerr)
+  if !is_half_sweep_done(direction, b, N; ncenter=nsite)
+    # Do backwards evolution step
+    b1 = (isforward(direction) ? b + 1 : b)
+    Δ = (isforward(direction) ? +1 : -1)
+    phi0 = psi[b1]
+    set_nsite!(PH, nsite - 1)
+    position!(PH, psi, b1)
+    phi0, info = solver(PH, -time_step, phi0; current_time, outputlevel)
+    current_time -= time_step
+    normalize && (phi0 ./= norm(phi0))
+    psi[b1] = phi0
+    set_nsite!(PH, nsite)
+  end
+  return current_time, maxtruncerr, spec, info
+end
+
+function region_update!(
+  ::Val{nsite},
+  ::Val{reverse_step},
+  solver,
+  PH,
+  psi,
+  b;
+  current_time,
+  outputlevel,
+  time_step,
+  normalize,
+  direction,
+  noise,
+  which_decomp,
+  svd_alg,
+  cutoff,
+  maxdim,
+  mindim,
+  maxtruncerr,
+) where {nsite,reverse_step}
+  return error(
+    "`tdvp`, `dmrg`, `linsolve`, etc. with `nsite=$nsite` and `reverse_step=$reverse_step` not implemented.",
+  )
+end
diff --git a/src/ITensorMPS/sweeps.jl b/src/ITensorMPS/sweeps.jl
index a1db337dea..1f56633dab 100644
--- a/src/ITensorMPS/sweeps.jl
+++ b/src/ITensorMPS/sweeps.jl
@@ -1,3 +1,4 @@
+using Printf: @printf
 
 """
 A Sweeps objects holds information
diff --git a/src/ITensorMPS/tdvp.jl b/src/ITensorMPS/tdvp.jl
new file mode 100644
index 0000000000..331e9acefe
--- /dev/null
+++ b/src/ITensorMPS/tdvp.jl
@@ -0,0 +1,136 @@
+using ITensors: Algorithm, @Algorithm_str
+
+# Select solver function
+solver_function(solver_backend::String) = solver_function(Algorithm(solver_backend))
+solver_function(::Algorithm"exponentiate") = exponentiate
+function solver_function(solver_backend::Algorithm)
+  return error(
+    "solver_backend=$(String(solver_backend)) not recognized (only \"exponentiate\" is supported)",
+  )
+end
+
+# Kept for backwards compatibility
+function solver_function(::Algorithm"applyexp")
+  println(
+    "Warning: the `solver_backend` option `\"applyexp\"` in `tdvp` has been removed. `\"exponentiate\"` will be used instead. To remove this warning, don't specify the `solver_backend` keyword argument.",
+  )
+  return solver_function(Algorithm"exponentiate"())
+end
+
+function tdvp_solver(
+  f::typeof(exponentiate);
+  ishermitian,
+  issymmetric,
+  solver_tol,
+  solver_krylovdim,
+  solver_maxiter,
+  solver_outputlevel,
+)
+  function solver(H, t, psi0; current_time, outputlevel)
+    psi, info = f(
+      H,
+      t,
+      psi0;
+      ishermitian,
+      issymmetric,
+      tol=solver_tol,
+      krylovdim=solver_krylovdim,
+      maxiter=solver_maxiter,
+      verbosity=solver_outputlevel,
+      eager=true,
+    )
+    return psi, info
+  end
+  return solver
+end
+
+function tdvp(
+  H,
+  t::Number,
+  psi0::MPS;
+  ishermitian=default_ishermitian(),
+  issymmetric=default_issymmetric(),
+  solver_backend=default_tdvp_solver_backend(),
+  solver_function=solver_function(solver_backend),
+  solver_tol=default_solver_tol(solver_function),
+  solver_krylovdim=default_solver_krylovdim(solver_function),
+  solver_maxiter=default_solver_maxiter(solver_function),
+  solver_outputlevel=default_solver_outputlevel(solver_function),
+  kwargs...,
+)
+  return tdvp(
+    tdvp_solver(
+      solver_function;
+      ishermitian,
+      issymmetric,
+      solver_tol,
+      solver_krylovdim,
+      solver_maxiter,
+      solver_outputlevel,
+    ),
+    H,
+    t,
+    psi0;
+    kwargs...,
+  )
+end
+
+function tdvp(t::Number, H, psi0::MPS; kwargs...)
+  return tdvp(H, t, psi0; kwargs...)
+end
+
+function tdvp(H, psi0::MPS, t::Number; kwargs...)
+  return tdvp(H, t, psi0; kwargs...)
+end
+
+"""
+    tdvp(H::MPO,psi0::MPS,t::Number; kwargs...)
+    tdvp(H::MPO,psi0::MPS,t::Number; kwargs...)
+
+Use the time dependent variational principle (TDVP) algorithm
+to compute `exp(t*H)*psi0` using an efficient algorithm based
+on alternating optimization of the MPS tensors and local Krylov
+exponentiation of H.
+
+Returns:
+* `psi::MPS` - time-evolved MPS
+
+Optional keyword arguments:
+* `outputlevel::Int = 1` - larger outputlevel values resulting in printing more information and 0 means no output
+* `observer` - object implementing the [Observer](@ref observer) interface which can perform measurements and stop early
+* `write_when_maxdim_exceeds::Int` - when the allowed maxdim exceeds this value, begin saving tensors to disk to free memory in large calculations
+"""
+function tdvp(solver, H::MPO, t::Number, psi0::MPS; kwargs...)
+  return alternating_update(solver, H, t, psi0; kwargs...)
+end
+
+function tdvp(solver, t::Number, H, psi0::MPS; kwargs...)
+  return tdvp(solver, H, t, psi0; kwargs...)
+end
+
+function tdvp(solver, H, psi0::MPS, t::Number; kwargs...)
+  return tdvp(solver, H, t, psi0; kwargs...)
+end
+
+"""
+    tdvp(Hs::Vector{MPO},psi0::MPS,t::Number; kwargs...)
+    tdvp(Hs::Vector{MPO},psi0::MPS,t::Number, sweeps::Sweeps; kwargs...)
+
+Use the time dependent variational principle (TDVP) algorithm
+to compute `exp(t*H)*psi0` using an efficient algorithm based
+on alternating optimization of the MPS tensors and local Krylov
+exponentiation of H.
+
+This version of `tdvp` accepts a representation of H as a
+Vector of MPOs, Hs = [H1,H2,H3,...] such that H is defined
+as H = H1+H2+H3+...
+Note that this sum of MPOs is not actually computed; rather
+the set of MPOs [H1,H2,H3,..] is efficiently looped over at
+each step of the algorithm when optimizing the MPS.
+
+Returns:
+* `psi::MPS` - time-evolved MPS
+"""
+function tdvp(solver, Hs::Vector{MPO}, t::Number, psi0::MPS; kwargs...)
+  return alternating_update(solver, Hs, t, psi0; kwargs...)
+end
diff --git a/src/ITensorMPS/tdvp_sweeps.jl b/src/ITensorMPS/tdvp_sweeps.jl
new file mode 100644
index 0000000000..d39ea92ac7
--- /dev/null
+++ b/src/ITensorMPS/tdvp_sweeps.jl
@@ -0,0 +1,17 @@
+function process_sweeps(s::Sweeps)
+  return (;
+    nsweeps=s.nsweep, maxdim=s.maxdim, mindim=s.mindim, cutoff=s.cutoff, noise=s.noise
+  )
+end
+
+function tdvp(H, t::Number, psi0::MPS, sweeps::Sweeps; kwargs...)
+  return tdvp(H, t, psi0; process_sweeps(sweeps)..., kwargs...)
+end
+
+function tdvp(solver, H, t::Number, psi0::MPS, sweeps::Sweeps; kwargs...)
+  return tdvp(solver, H, t, psi0; process_sweeps(sweeps)..., kwargs...)
+end
+
+function dmrg(H, psi0::MPS, sweeps::Sweeps; kwargs...)
+  return dmrg(H, psi0; process_sweeps(sweeps)..., kwargs...)
+end
diff --git a/src/ITensorMPS/tdvpinfo.jl b/src/ITensorMPS/tdvpinfo.jl
new file mode 100644
index 0000000000..fb4c128527
--- /dev/null
+++ b/src/ITensorMPS/tdvpinfo.jl
@@ -0,0 +1,7 @@
+"""
+#fields
+- `maxtruncerr::Float64`: the maximum tuncation error
+"""
+struct TDVPInfo
+  maxtruncerr::Float64
+end
diff --git a/src/ITensorMPS/tdvporder.jl b/src/ITensorMPS/tdvporder.jl
new file mode 100644
index 0000000000..e37e00e251
--- /dev/null
+++ b/src/ITensorMPS/tdvporder.jl
@@ -0,0 +1,24 @@
+struct TDVPOrder{order,direction} end
+
+TDVPOrder(order::Int, direction::Base.Ordering) = TDVPOrder{order,direction}()
+
+orderings(::TDVPOrder) = error("Not implemented")
+sub_time_steps(::TDVPOrder) = error("Not implemented")
+
+function orderings(::TDVPOrder{1,direction}) where {direction}
+  return [direction, Base.ReverseOrdering(direction)]
+end
+sub_time_steps(::TDVPOrder{1}) = [1, 0]
+
+function orderings(::TDVPOrder{2,direction}) where {direction}
+  return [direction, Base.ReverseOrdering(direction)]
+end
+sub_time_steps(::TDVPOrder{2}) = [1 / 2, 1 / 2]
+
+function orderings(::TDVPOrder{4,direction}) where {direction}
+  return [direction, Base.ReverseOrdering(direction)]
+end
+function sub_time_steps(::TDVPOrder{4})
+  s = 1 / (2 - 2^(1 / 3))
+  return [s / 2, s / 2, (1 - 2 * s) / 2, (1 - 2 * s) / 2, s / 2, s / 2]
+end
diff --git a/src/ITensorMPS/update_observer.jl b/src/ITensorMPS/update_observer.jl
new file mode 100644
index 0000000000..149f3d485d
--- /dev/null
+++ b/src/ITensorMPS/update_observer.jl
@@ -0,0 +1,5 @@
+using Observers: Observers
+
+function Observers.update!(observer::AbstractObserver; kwargs...)
+  return measure!(observer; kwargs...)
+end
diff --git a/src/ITensorMPS/usings.jl b/src/ITensorMPS/usings.jl
deleted file mode 100644
index 697e8d6870..0000000000
--- a/src/ITensorMPS/usings.jl
+++ /dev/null
@@ -1,7 +0,0 @@
-using Adapt
-using ..ITensors
-using HDF5: attributes, create_group, open_group
-using IsApprox
-using Printf
-using Random
-using TupleTools
diff --git a/src/imports.jl b/src/imports.jl
index d4ce90d82b..87f089ea1b 100644
--- a/src/imports.jl
+++ b/src/imports.jl
@@ -171,7 +171,6 @@ import ITensors.NDTensors:
   permuteblocks,
   polar,
   ql,
-  scalartype,
   scale!,
   setblock!,
   setblockdim!,
diff --git a/src/itensor.jl b/src/itensor.jl
index cd3b12d3e5..5b13376371 100644
--- a/src/itensor.jl
+++ b/src/itensor.jl
@@ -795,7 +795,7 @@ symmetrystyle(T::Tensor) = symmetrystyle(inds(T))
 symmetrystyle(T::ITensor)::SymmetryStyle = symmetrystyle(tensor(T))
 
 eltype(T::ITensor) = eltype(tensor(T))
-scalartype(x::ITensor) = eltype(x)
+NDTensors.scalartype(x::ITensor) = eltype(x)
 
 """
     order(A::ITensor)
diff --git a/test/ITensorMPS/base/test_autompo.jl b/test/ITensorMPS/base/test_autompo.jl
index b44ef903ef..bc1e8f0d1b 100644
--- a/test/ITensorMPS/base/test_autompo.jl
+++ b/test/ITensorMPS/base/test_autompo.jl
@@ -1,4 +1,5 @@
 using ITensors, Test, Random, JLD2
+using NDTensors: scalartype
 
 include(joinpath(@__DIR__, "utils", "util.jl"))
 
@@ -294,26 +295,26 @@ end
     H1 = MPO(elt, O1, sites)
     H2 = MPO(elt, O2, sites)
     H = H1 + 2 * H2
-    @test ITensors.scalartype(H1) == elt
-    @test ITensors.scalartype(H2) == elt
-    @test ITensors.scalartype(H) == elt
+    @test scalartype(H1) == elt
+    @test scalartype(H2) == elt
+    @test scalartype(H) == elt
     @test prod(MPO(O, sites)) ≈ prod(H)
 
-    @test ITensors.scalartype(MPO(elt, Op("Sz", 1), sites)) == elt
-    @test ITensors.scalartype(MPO(elt, Op("Sz", 1) + Op("Sz", 2), sites)) == elt
-    @test ITensors.scalartype(MPO(elt, 2 * Op("Sz", 1) + 3 * Op("Sz", 2), sites)) == elt
-    @test ITensors.scalartype(MPO(elt, 2 * Op("Sz", 1), sites)) == elt
-    @test ITensors.scalartype(MPO(elt, Op("Sz", 1) * Op("Sz", 2), sites)) == elt
-    @test ITensors.scalartype(MPO(elt, 2 * Op("Sz", 1) * Op("Sz", 2), sites)) == elt
+    @test scalartype(MPO(elt, Op("Sz", 1), sites)) == elt
+    @test scalartype(MPO(elt, Op("Sz", 1) + Op("Sz", 2), sites)) == elt
+    @test scalartype(MPO(elt, 2 * Op("Sz", 1) + 3 * Op("Sz", 2), sites)) == elt
+    @test scalartype(MPO(elt, 2 * Op("Sz", 1), sites)) == elt
+    @test scalartype(MPO(elt, Op("Sz", 1) * Op("Sz", 2), sites)) == elt
+    @test scalartype(MPO(elt, 2 * Op("Sz", 1) * Op("Sz", 2), sites)) == elt
 
     O = O1 - 2 * O2
     @test length(O) == 2 * n - 1
     H1 = MPO(elt, O1, sites)
     H2 = MPO(elt, O2, sites)
     H = H1 - 2 * H2
-    @test ITensors.scalartype(H1) == elt
-    @test ITensors.scalartype(H2) == elt
-    @test ITensors.scalartype(H) == elt
+    @test scalartype(H1) == elt
+    @test scalartype(H2) == elt
+    @test scalartype(H) == elt
     @test prod(MPO(O, sites)) ≈ prod(H)
 
     O = O1 - O2 / 2
@@ -321,9 +322,9 @@ end
     H1 = MPO(elt, O1, sites)
     H2 = MPO(elt, O2, sites)
     H = H1 - H2 / 2
-    @test ITensors.scalartype(H1) == elt
-    @test ITensors.scalartype(H2) == elt
-    @test ITensors.scalartype(H) == elt
+    @test scalartype(H1) == elt
+    @test scalartype(H2) == elt
+    @test scalartype(H) == elt
     @test prod(MPO(O, sites)) ≈ prod(H)
   end
 
@@ -345,7 +346,7 @@ end
     end
     sites = siteinds("S=1/2", N)
     Ha = MPO(os, sites)
-    @test ITensors.scalartype(Ha) <: Float64
+    @test scalartype(Ha) <: Float64
     He = isingMPO(sites)
     psi = makeRandomMPS(sites)
     Oa = inner(psi', Ha, psi)
@@ -353,7 +354,7 @@ end
     @test Oa ≈ Oe
 
     H_complex = MPO(ComplexF64, os, sites)
-    @test ITensors.scalartype(H_complex) <: ComplexF64
+    @test scalartype(H_complex) <: ComplexF64
     @test H_complex ≈ Ha
   end
 
diff --git a/test/ITensorMPS/base/test_mpo.jl b/test/ITensorMPS/base/test_mpo.jl
index cbe96e205b..9813a4697f 100644
--- a/test/ITensorMPS/base/test_mpo.jl
+++ b/test/ITensorMPS/base/test_mpo.jl
@@ -1,5 +1,6 @@
 using Combinatorics
 using ITensors
+using NDTensors: scalartype
 using Test
 
 include(joinpath(@__DIR__, "utils", "util.jl"))
@@ -800,7 +801,7 @@ end
     sites = siteinds("S=1/2", 4)
     A = randn(T) * convert_leaf_eltype(T, randomMPO(sites))
     B = randn(T) * convert_leaf_eltype(T, randomMPO(sites))
-    @test ITensors.scalartype(apply(A, B)) == T
+    @test scalartype(apply(A, B)) == T
   end
   @testset "sample" begin
     N = 6