diff --git a/Project.toml b/Project.toml
index 9207780..1e1b99e 100644
--- a/Project.toml
+++ b/Project.toml
@@ -1,7 +1,7 @@
 name = "ITensorTDVP"
 uuid = "25707e16-a4db-4a07-99d9-4d67b7af0342"
 authors = ["Matthew Fishman <mfishman@flatironinstitute.org> and contributors"]
-version = "0.3.1"
+version = "0.4.0"
 
 [deps]
 ITensors = "9136182c-28ba-11e9-034c-db9fb085ebd5"
diff --git a/examples/01_tdvp.jl b/examples/01_tdvp.jl
index 9297d47..1c48327 100644
--- a/examples/01_tdvp.jl
+++ b/examples/01_tdvp.jl
@@ -21,9 +21,9 @@ function main()
 
   ϕ = tdvp(
     H,
-    -1.0,
+    -20.0,
     ψ;
-    nsweeps=20,
+    time_step=-1.0,
     reverse_step=false,
     normalize=true,
     maxdim=30,
diff --git a/examples/02_dmrg-x.jl b/examples/02_dmrg-x.jl
index 61b2aa4..cdd207e 100644
--- a/examples/02_dmrg-x.jl
+++ b/examples/02_dmrg-x.jl
@@ -31,9 +31,7 @@ function main()
   initstate = rand(["↑", "↓"], n)
   ψ = MPS(s, initstate)
 
-  dmrg_x_kwargs = (
-    nsweeps=10, reverse_step=false, normalize=true, maxdim=20, cutoff=1e-10, outputlevel=1
-  )
+  dmrg_x_kwargs = (nsweeps=10, normalize=true, maxdim=20, cutoff=1e-10, outputlevel=1)
 
   e, ϕ = dmrg_x(H, ψ; dmrg_x_kwargs...)
 
diff --git a/src/alternating_update.jl b/src/alternating_update.jl
index 12701a0..351b1b3 100644
--- a/src/alternating_update.jl
+++ b/src/alternating_update.jl
@@ -1,30 +1,36 @@
-using ITensors: permute
+using ITensors: ITensors, permute
 using ITensors.ITensorMPS:
-  AbstractObserver,
+  ## AbstractObserver,
   MPO,
   MPS,
   ProjMPO,
   ProjMPOSum,
   check_hascommoninds,
-  checkdone!,
+  ## checkdone!,
   disk,
   linkind,
   maxlinkdim,
   siteinds
 
-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 _compute_nsweeps(t; time_step=default_time_step(t), nsweeps=default_nsweeps())
+## 
+##   @show t, time_step, 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
+## 
+##   @show nsweeps
+## 
+##   return nsweeps
+## end
 
 function _extend_sweeps_param(param, nsweeps)
   if param isa Number
@@ -48,16 +54,15 @@ end
 
 function alternating_update(
   solver,
-  PH,
-  t::Number,
-  psi0::MPS;
+  reduced_operator,
+  init::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),
+  time_step=default_time_step(),
   order=default_order(),
   (observer!)=default_observer!(),
   (step_observer!)=default_step_observer!(),
@@ -65,13 +70,13 @@ function alternating_update(
   normalize=default_normalize(),
   maxdim=default_maxdim(),
   mindim=default_mindim(),
-  cutoff=default_cutoff(Float64),
+  cutoff=default_cutoff(ITensors.scalartype(init)),
   noise=default_noise(),
 )
-  nsweeps = _compute_nsweeps(t; time_step, nsweeps)
+  ## 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)
+  state = copy(init)
   # Keep track of the start of the current time step.
   # Helpful for tracking the total time, for example
   # when using time-dependent solvers.
@@ -86,17 +91,17 @@ function alternating_update(
           "write_when_maxdim_exceeds = $write_when_maxdim_exceeds and maxdim(sweeps, sw) = $(maxdim(sweeps, sweep)), writing environment tensors to disk",
         )
       end
-      PH = disk(PH)
+      reduced_operator = disk(reduced_operator)
     end
-    sweep_time = @elapsed begin
-      psi, PH, info = sweep_update(
+    sweep_elapsed_time = @elapsed begin
+      state, reduced_operator, info = sweep_update(
         forward_order,
         solver,
-        PH,
-        time_step,
-        psi;
+        reduced_operator,
+        state;
         nsite,
         current_time,
+        time_step,
         reverse_step,
         sweep,
         observer!,
@@ -107,51 +112,54 @@ function alternating_update(
         noise=noise[sweep],
       )
     end
-    current_time += time_step
-    update_observer!(step_observer!; psi, sweep, outputlevel, current_time)
+    if !isnothing(time_step)
+      current_time += time_step
+    end
+    update_observer!(step_observer!; state, sweep, outputlevel, current_time)
     if outputlevel >= 1
       print("After sweep ", sweep, ":")
-      print(" maxlinkdim=", maxlinkdim(psi))
+      print(" maxlinkdim=", maxlinkdim(state))
       @printf(" maxerr=%.2E", info.maxtruncerr)
       print(" current_time=", round(current_time; digits=3))
-      print(" time=", round(sweep_time; digits=3))
+      print(" time=", round(sweep_elapsed_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 = checkdone(; state, sweep, outputlevel)
+    ## isdone = false
+    ## if !isnothing(checkdone)
+    ##  isdone = checkdone(; state, sweep, outputlevel)
+    ## elseif observer! isa AbstractObserver
+    ##   isdone = checkdone!(observer!; state, sweep, outputlevel)
+    ## end
     isdone && break
   end
-  return psi
+  return state
 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, ITensors.scalartype(psi0)(Inf), psi0; kwargs...)
-end
+## function alternating_update(solver, operator, init::MPS; kwargs...)
+##   return alternating_update(solver, operator, ITensors.scalartype(init)(Inf), init; kwargs...)
+## end
 
-function alternating_update(solver, H::MPO, t::Number, psi0::MPS; kwargs...)
-  check_hascommoninds(siteinds, H, psi0)
-  check_hascommoninds(siteinds, H, psi0')
+function alternating_update(solver, operator::MPO, init::MPS; kwargs...)
+  check_hascommoninds(siteinds, operator, init)
+  check_hascommoninds(siteinds, operator, init')
   # 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...)
+  operator = permute(operator, (linkind, siteinds, linkind))
+  reduced_operator = ProjMPO(operator)
+  return alternating_update(solver, reduced_operator, init; 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')
+function alternating_update(solver, operators::Vector{MPO}, init::MPS; kwargs...)
+  for operator in operators
+    check_hascommoninds(siteinds, operator, init)
+    check_hascommoninds(siteinds, operator, init')
   end
-  Hs .= ITensors.permute.(Hs, Ref((linkind, siteinds, linkind)))
-  PHs = ProjMPOSum(Hs)
-  return alternating_update(solver, PHs, t, psi0; kwargs...)
+  operators .= ITensors.permute.(operators, Ref((linkind, siteinds, linkind)))
+  reduced_operator = ProjMPOSum(operators)
+  return alternating_update(solver, reduced_operator, init; kwargs...)
 end
diff --git a/src/contract_mpo_mps.jl b/src/contract_mpo_mps.jl
index 8f2dac4..c40ac27 100644
--- a/src/contract_mpo_mps.jl
+++ b/src/contract_mpo_mps.jl
@@ -12,48 +12,55 @@ using ITensors:
   siteinds
 
 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]
+  function solver(reduced_operator, psi; kws...)
+    reduced_state = ITensor(true)
+    for j in (reduced_operator.lpos + 1):(reduced_operator.rpos - 1)
+      reduced_state *= reduced_operator.input_state[j]
     end
-    Hpsi0 = contract(PH, v)
-    return Hpsi0, nothing
+    reduced_state = contract(reduced_operator, reduced_state)
+    return reduced_state, nothing
   end
   return solver
 end
 
+# `init_mps` is for backwards compatibility.
 function ITensors.contract(
-  ::Algorithm"fit", A::MPO, psi0::MPS; init_mps=psi0, nsweeps=1, kwargs...
+  ::Algorithm"fit",
+  operator::MPO,
+  input_state::MPS;
+  init=input_state,
+  init_mps=init,
+  kwargs...,
 )::MPS
-  n = length(A)
-  n != length(psi0) &&
-    throw(DimensionMismatch("lengths of MPO ($n) and MPS ($(length(psi0))) do not match"))
+  n = length(operator)
+  n != length(input_state) && throw(
+    DimensionMismatch("lengths of MPO ($n) and MPS ($(length(input_state))) do not match")
+  )
   if n == 1
-    return MPS([A[1] * psi0[1]])
+    return MPS([operator[1] * input_state[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)
+  any(i -> isempty(i), siteinds(commoninds, operator, input_state)) && error(
+    "In `contract(operator::MPO, x::MPS)`, `operator` and `x` must share a set of site indices",
+  )
+  # In case operator and input_state have the same link indices
+  operator = sim(linkinds, operator)
+  # Fix site and link inds of init
+  init = deepcopy(init)
+  init = sim(linkinds, init)
+  siteinds_operator = siteinds(operator)
   ti = Vector{Index}(undef, n)
   for j in 1:n
-    for i in Ai[j]
-      if !hasind(psi0[j], i)
+    for i in siteinds_operator[j]
+      if !hasind(input_state[j], i)
         ti[j] = i
         break
       end
     end
   end
-  replace_siteinds!(init_mps, ti)
-  reverse_step = false
-  PH = ProjMPOApply(psi0, A)
+  replace_siteinds!(init, ti)
+  reduced_operator = ProjMPOApply(input_state, operator)
   psi = alternating_update(
-    contractmpo_solver(; kwargs...), PH, init_mps; nsweeps, reverse_step, kwargs...
+    contractmpo_solver(; kwargs...), reduced_operator, init; kwargs...
   )
   return psi
 end
diff --git a/src/defaults.jl b/src/defaults.jl
index dd745a2..df037a4 100644
--- a/src/defaults.jl
+++ b/src/defaults.jl
@@ -1,16 +1,16 @@
-using ITensors: NoObserver
+## using ITensors: NoObserver
 using KrylovKit: eigsolve, exponentiate
 
 default_nsweeps() = nothing
-default_checkdone() = nothing
+default_checkdone() = Returns(false)
 default_write_when_maxdim_exceeds() = nothing
 default_nsite() = 2
-default_reverse_step() = true
-default_time_start() = 0
-default_time_step(t) = t
+default_reverse_step() = false
+default_time_start() = nothing
+default_time_step() = nothing
 default_order() = 2
-default_observer!() = NoObserver()
-default_step_observer!() = NoObserver()
+default_observer!() = EmptyObserver()
+default_step_observer!() = EmptyObserver()
 default_outputlevel() = 0
 default_normalize() = false
 default_sweep() = 1
diff --git a/src/dmrg.jl b/src/dmrg.jl
index ef6f412..4fd2ec9 100644
--- a/src/dmrg.jl
+++ b/src/dmrg.jl
@@ -7,12 +7,12 @@ function dmrg_solver(
   solver_maxiter,
   solver_verbosity,
 )
-  function solver(H, t, psi0; current_time, outputlevel)
+  function solver(operator, init; current_time, time_step, outputlevel)
     howmany = 1
     which = solver_which_eigenvalue
     vals, vecs, info = f(
-      H,
-      psi0,
+      operator,
+      init,
       howmany,
       which;
       ishermitian=default_ishermitian(),
@@ -27,8 +27,8 @@ function dmrg_solver(
 end
 
 function dmrg(
-  H,
-  psi0::MPS;
+  operator,
+  init::MPS;
   ishermitian=default_ishermitian(),
   solver_which_eigenvalue=default_solver_which_eigenvalue(eigsolve),
   solver_tol=default_solver_tol(eigsolve),
@@ -38,11 +38,10 @@ function dmrg(
   (observer!)=default_observer!(),
   kwargs...,
 )
-  reverse_step = false
   info_ref! = Ref{Any}()
   info_observer! = values_observer(; info=info_ref!)
   observer! = compose_observers(observer!, info_observer!)
-  psi = alternating_update(
+  state = alternating_update(
     dmrg_solver(
       eigsolve;
       solver_which_eigenvalue,
@@ -52,11 +51,10 @@ function dmrg(
       solver_maxiter,
       solver_verbosity,
     ),
-    H,
-    psi0;
-    reverse_step,
+    operator,
+    init;
     observer!,
     kwargs...,
   )
-  return info_ref![].eigval, psi
+  return info_ref![].eigval, state
 end
diff --git a/src/dmrg_x.jl b/src/dmrg_x.jl
index d640237..37fce33 100644
--- a/src/dmrg_x.jl
+++ b/src/dmrg_x.jl
@@ -12,12 +12,10 @@ function dmrg_x_solver(PH, t, psi0; current_time, outputlevel)
   return U_max, (; eigval=D_max)
 end
 
-function dmrg_x(
-  PH, psi0::MPS; reverse_step=false, (observer!)=default_observer!(), kwargs...
-)
+function dmrg_x(PH, psi0::MPS; (observer!)=default_observer!(), kwargs...)
   info_ref = Ref{Any}()
   info_observer! = values_observer(; info=info_ref)
   observer! = compose_observers(observer!, info_observer!)
-  psi = alternating_update(dmrg_x_solver, PH, psi0; reverse_step, observer!, kwargs...)
+  psi = alternating_update(dmrg_x_solver, PH, psi0; observer!, kwargs...)
   return info_ref[].eigval, psi
 end
diff --git a/src/linsolve.jl b/src/linsolve.jl
index 88dd403..f943963 100644
--- a/src/linsolve.jl
+++ b/src/linsolve.jl
@@ -24,13 +24,7 @@ Keyword arguments:
     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...,
+  A::MPO, b::MPS, x₀::MPS, a₀::Number=false, a₁::Number=true; solver_kwargs=(;), kwargs...
 )
   function linsolve_solver(P::ProjMPO_MPS2, t, x₀; current_time, outputlevel)
     b = dag(only(proj_mps(P)))
@@ -38,5 +32,5 @@ function KrylovKit.linsolve(
     return x, nothing
   end
   P = ProjMPO_MPS2(A, b)
-  return alternating_update(linsolve_solver, P, x₀; reverse_step=false, tdvp_kwargs...)
+  return alternating_update(linsolve_solver, P, x₀; kwargs...)
 end
diff --git a/src/projmpo_apply.jl b/src/projmpo_apply.jl
index 0493e48..637f0be 100644
--- a/src/projmpo_apply.jl
+++ b/src/projmpo_apply.jl
@@ -3,88 +3,125 @@ using ITensors.ITensorMPS: ITensorMPS, AbstractProjMPO, MPO, MPS
 
 """
 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>).
+MPO `operator` onto an MPS `input_state` but "projected" by
+the basis of a different MPS `state` (which
+could be an approximation to operator|state>).
 
 As an implementation of the AbstractProjMPO
 type, it supports multiple `nsite` values for
 one- and two-site algorithms.
 
 ```
-     *--*--*-      -*--*--*--*--*--* <psi|
+     *--*--*-      -*--*--*--*--*--* <state|
      |  |  |  |  |  |  |  |  |  |  |
-     h--h--h--h--h--h--h--h--h--h--h H  
+     h--h--h--h--h--h--h--h--h--h--h operator  
      |  |  |  |  |  |  |  |  |  |  |
-     o--o--o-      -o--o--o--o--o--o |psi0>
+     o--o--o-      -o--o--o--o--o--o |input_state>
 ```
 """
 mutable struct ProjMPOApply <: AbstractProjMPO
   lpos::Int
   rpos::Int
   nsite::Int
-  psi0::MPS
-  H::MPO
+  input_state::MPS
+  operator::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)))
+function ProjMPOApply(input_state::MPS, operator::MPO)
+  lpos = 0
+  rpos = length(operator) + 1
+  nsite = 2
+  LR = Vector{ITensor}(undef, length(operator))
+  return ProjMPOApply(
+    lpos,
+    rpos,
+    nsite,
+    input_state,
+    operator,
+    LR,
+  )
 end
 
-function Base.copy(P::ProjMPOApply)
-  return ProjMPOApply(P.lpos, P.rpos, P.nsite, copy(P.psi0), copy(P.H), copy(P.LR))
+function Base.getproperty(reduced_operator::ProjMPOApply, sym::Symbol)
+  if sym === :H
+    # This is for compatibility with `AbstractProjMPO`.
+    # TODO: Don't use `reduced_operator.H` in `AbstractProjMPO`.
+    return getfield(reduced_operator, :operator)
+  end
+  return getfield(reduced_operator, sym)
+end
+
+function Base.copy(reduced_operator::ProjMPOApply)
+  return ProjMPOApply(
+    reduced_operator.lpos,
+    reduced_operator.rpos,
+    reduced_operator.nsite,
+    copy(reduced_operator.input_state),
+    copy(reduced_operator.operator),
+    copy(reduced_operator.LR),
+  )
 end
 
-function ITensorMPS.set_nsite!(P::ProjMPOApply, nsite)
-  P.nsite = nsite
-  return P
+Base.length(reduced_operator::ProjMPOApply) = length(reduced_operator.operator)
+
+function ITensorMPS.set_nsite!(reduced_operator::ProjMPOApply, nsite)
+  reduced_operator.nsite = nsite
+  return reduced_operator
 end
 
-function ITensorMPS.makeL!(P::ProjMPOApply, psi::MPS, k::Int)
+function ITensorMPS.makeL!(reduced_operator::ProjMPOApply, state::MPS, k::Int)
   # Save the last `L` that is made to help with caching
   # for DiskProjMPO
-  ll = P.lpos
+  ll = reduced_operator.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
+    reduced_operator.lpos = k
     return nothing
   end
   # Make sure ll is at least 0 for the generic logic below
   ll = max(ll, 0)
-  L = lproj(P)
+  L = lproj(reduced_operator)
   while ll < k
-    L = L * P.psi0[ll + 1] * P.H[ll + 1] * dag(psi[ll + 1])
-    P.LR[ll + 1] = L
+    L =
+      L *
+      reduced_operator.input_state[ll + 1] *
+      reduced_operator.operator[ll + 1] *
+      dag(state[ll + 1])
+    reduced_operator.LR[ll + 1] = L
     ll += 1
   end
   # Needed when moving lproj backward.
-  P.lpos = k
-  return P
+  reduced_operator.lpos = k
+  return reduced_operator
 end
 
-function ITensorMPS.makeR!(P::ProjMPOApply, psi::MPS, k::Int)
+function ITensorMPS.makeR!(reduced_operator::ProjMPOApply, state::MPS, k::Int)
   # Save the last `R` that is made to help with caching
   # for DiskProjMPO
-  rl = P.rpos
+  rl = reduced_operator.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
+    reduced_operator.rpos = k
     return nothing
   end
-  N = length(P.H)
+  N = length(reduced_operator.operator)
   # Make sure rl is no bigger than `N + 1` for the generic logic below
   rl = min(rl, N + 1)
-  R = rproj(P)
+  R = rproj(reduced_operator)
   while rl > k
-    R = R * P.psi0[rl - 1] * P.H[rl - 1] * dag(psi[rl - 1])
-    P.LR[rl - 1] = R
+    R =
+      R *
+      reduced_operator.input_state[rl - 1] *
+      reduced_operator.operator[rl - 1] *
+      dag(state[rl - 1])
+    reduced_operator.LR[rl - 1] = R
     rl -= 1
   end
-  P.rpos = k
-  return P
+  reduced_operator.rpos = k
+  return reduced_operator
 end
diff --git a/src/sweep_update.jl b/src/sweep_update.jl
index e950911..cbfb947 100644
--- a/src/sweep_update.jl
+++ b/src/sweep_update.jl
@@ -1,29 +1,44 @@
-using ITensors: uniqueinds
+using ITensors: ITensors, uniqueinds
 using ITensors.ITensorMPS:
   ITensorMPS, MPS, isortho, orthocenter, orthogonalize!, position!, replacebond!, set_nsite!
 using LinearAlgebra: norm, normalize!, svd
 using Printf: @printf
 
 function sweep_update(
-  order::TDVPOrder, solver, PH, time_step::Number, psi::MPS; current_time=0, kwargs...
+  order::TDVPOrder,
+  solver,
+  reduce_operator,
+  state::MPS;
+  current_time=nothing,
+  time_step=nothing,
+  kwargs...,
 )
   order_orderings = orderings(order)
-  order_sub_time_steps = eltype(time_step).(sub_time_steps(order))
-  order_sub_time_steps *= time_step
+  order_sub_time_steps = sub_time_steps(order)
+  if !isnothing(time_step)
+    order_sub_time_steps = eltype(time_step).(order_sub_time_steps)
+    order_sub_time_steps *= time_step
+  end
   info = nothing
+  sub_time_step = nothing
   for substep in 1:length(order_sub_time_steps)
-    psi, PH, info = sub_sweep_update(
+    if !isnothing(time_step)
+      sub_time_step = order_sub_time_steps[substep]
+    end
+    state, reduce_operator, info = sub_sweep_update(
       order_orderings[substep],
       solver,
-      PH,
-      order_sub_time_steps[substep],
-      psi;
+      reduce_operator,
+      state;
       current_time,
+      time_step=sub_time_step,
       kwargs...,
     )
-    current_time += order_sub_time_steps[substep]
+    if !isnothing(time_step)
+      current_time += sub_time_step
+    end
   end
-  return psi, PH, info
+  return state, reduce_operator, info
 end
 
 isforward(direction::Base.ForwardOrdering) = true
@@ -50,13 +65,13 @@ end
 function sub_sweep_update(
   direction::Base.Ordering,
   solver,
-  PH,
-  time_step::Number,
-  psi::MPS;
+  reduce_operator,
+  state::MPS;
   which_decomp=nothing,
   svd_alg=nothing,
   sweep=default_sweep(),
-  current_time=default_current_time(),
+  current_time=nothing,
+  time_step=nothing,
   nsite=default_nsite(),
   reverse_step=default_reverse_step(),
   normalize=default_normalize(),
@@ -64,38 +79,38 @@ function sub_sweep_update(
   outputlevel=default_outputlevel(),
   maxdim=default_maxdim(),
   mindim=default_mindim(),
-  cutoff=default_cutoff(time_step),
+  cutoff=default_cutoff(ITensors.scalartype(state)),
   noise=default_noise(),
 )
-  PH = copy(PH)
-  psi = copy(psi)
-  if length(psi) == 1
+  reduce_operator = copy(reduce_operator)
+  state = copy(state)
+  if length(state) == 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)
+  N = length(state)
+  set_nsite!(reduce_operator, nsite)
   if isforward(direction)
-    if !isortho(psi) || orthocenter(psi) != 1
-      orthogonalize!(psi, 1)
+    if !isortho(state) || orthocenter(state) != 1
+      orthogonalize!(state, 1)
     end
-    @assert isortho(psi) && orthocenter(psi) == 1
-    position!(PH, psi, 1)
+    @assert isortho(state) && orthocenter(state) == 1
+    position!(reduce_operator, state, 1)
   elseif isreverse(direction)
-    if !isortho(psi) || orthocenter(psi) != N - nsite + 1
-      orthogonalize!(psi, N - nsite + 1)
+    if !isortho(state) || orthocenter(state) != N - nsite + 1
+      orthogonalize!(state, N - nsite + 1)
     end
-    @assert(isortho(psi) && (orthocenter(psi) == N - nsite + 1))
-    position!(PH, psi, N - nsite + 1)
+    @assert(isortho(state) && (orthocenter(state) == N - nsite + 1))
+    position!(reduce_operator, state, 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,
+      reduce_operator,
+      state,
       b;
       nsite,
       reverse_step,
@@ -126,14 +141,14 @@ function sub_sweep_update(
       println()
       if spec != nothing
         @printf(
-          "  Trunc. err=%.2E, bond dimension %d\n", spec.truncerr, dim(linkind(psi, b))
+          "  Trunc. err=%.2E, bond dimension %d\n", spec.truncerr, dim(linkind(state, b))
         )
       end
       flush(stdout)
     end
     update_observer!(
       observer!;
-      psi,
+      state,
       bond=b,
       sweep,
       half_sweep=isforward(direction) ? 1 : 2,
@@ -145,14 +160,14 @@ function sub_sweep_update(
     )
   end
   # Just to be sure:
-  normalize && normalize!(psi)
-  return psi, PH, TDVPInfo(maxtruncerr)
+  normalize && normalize!(state)
+  return state, reduce_operator, TDVPInfo(maxtruncerr)
 end
 
 function region_update!(
   solver,
-  PH,
-  psi,
+  reduce_operator,
+  state,
   b;
   nsite,
   reverse_step,
@@ -173,8 +188,8 @@ function region_update!(
     Val(nsite),
     Val(reverse_step),
     solver,
-    PH,
-    psi,
+    reduce_operator,
+    state,
     b;
     current_time,
     outputlevel,
@@ -195,8 +210,8 @@ function region_update!(
   nsite_val::Val{1},
   reverse_step_val::Val{false},
   solver,
-  PH,
-  psi,
+  reduce_operator,
+  state,
   b;
   current_time,
   outputlevel,
@@ -211,21 +226,25 @@ function region_update!(
   mindim,
   maxtruncerr,
 )
-  N = length(psi)
+  N = length(state)
   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))
+  set_nsite!(reduce_operator, nsite)
+  position!(reduce_operator, state, b)
+  reduced_state = state[b]
+  reduced_state, info = solver(
+    reduce_operator, reduced_state; current_time, time_step, outputlevel
+  )
+  if !isnothing(time_step)
+    current_time += time_step
+  end
+  normalize && (reduced_state /= norm(reduced_state))
   spec = nothing
-  psi[b] = phi1
+  state[b] = reduced_state
   if !is_half_sweep_done(direction, b, N; ncenter=nsite)
     # Move ortho center
     Δ = (isforward(direction) ? +1 : -1)
-    orthogonalize!(psi, b + Δ)
+    orthogonalize!(state, b + Δ)
   end
   return current_time, maxtruncerr, spec, info
 end
@@ -234,8 +253,8 @@ function region_update!(
   nsite_val::Val{1},
   reverse_step_val::Val{true},
   solver,
-  PH,
-  psi,
+  reduce_operator,
+  state,
   b;
   current_time,
   outputlevel,
@@ -250,42 +269,46 @@ function region_update!(
   mindim,
   maxtruncerr,
 )
-  N = length(psi)
+  N = length(state)
   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)
+  set_nsite!(reduce_operator, nsite)
+  position!(reduce_operator, state, b)
+  reduced_state = state[b]
+  reduced_state, info = solver(
+    reduce_operator, time_step, reduced_state; current_time, outputlevel
+  )
   current_time += time_step
-  normalize && (phi1 /= norm(phi1))
+  normalize && (reduced_state /= norm(reduced_state))
   spec = nothing
-  psi[b] = phi1
+  state[b] = reduced_state
   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
+    uinds = uniqueinds(reduced_state, state[b + Δ])
+    U, S, V = svd(reduced_state, uinds)
+    state[b] = U
+    bond_reduced_state = S * V
     if isforward(direction)
-      ITensorMPS.setleftlim!(psi, b)
+      ITensorMPS.setleftlim!(state, b)
     elseif isreverse(direction)
-      ITensorMPS.setrightlim!(psi, b)
+      ITensorMPS.setrightlim!(state, b)
     end
-    set_nsite!(PH, nsite - 1)
-    position!(PH, psi, b1)
-    phi0, info = solver(PH, -time_step, phi0; current_time, outputlevel)
+    set_nsite!(reduce_operator, nsite - 1)
+    position!(reduce_operator, state, b1)
+    bond_reduced_state, info = solver(
+      reduce_operator, -time_step, bond_reduced_state; current_time, outputlevel
+    )
     current_time -= time_step
-    normalize && (phi0 ./= norm(phi0))
-    psi[b + Δ] = phi0 * psi[b + Δ]
+    normalize && (bond_reduced_state ./= norm(bond_reduced_state))
+    state[b + Δ] = bond_reduced_state * state[b + Δ]
     if isforward(direction)
-      ITensorMPS.setrightlim!(psi, b + Δ + 1)
+      ITensorMPS.setrightlim!(state, b + Δ + 1)
     elseif isreverse(direction)
-      ITensorMPS.setleftlim!(psi, b + Δ - 1)
+      ITensorMPS.setleftlim!(state, b + Δ - 1)
     end
-    set_nsite!(PH, nsite)
+    set_nsite!(reduce_operator, nsite)
   end
   return current_time, maxtruncerr, spec, info
 end
@@ -294,12 +317,12 @@ function region_update!(
   nsite_val::Val{2},
   reverse_step_val::Val{false},
   solver,
-  PH,
-  psi,
+  reduce_operator,
+  state,
   b;
   current_time,
-  outputlevel,
   time_step,
+  outputlevel,
   normalize,
   direction,
   noise,
@@ -310,25 +333,29 @@ function region_update!(
   mindim,
   maxtruncerr,
 )
-  N = length(psi)
+  N = length(state)
   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))
+  set_nsite!(reduce_operator, nsite)
+  position!(reduce_operator, state, b)
+  reduced_state = state[b] * state[b + 1]
+  reduced_state, info = solver(
+    reduce_operator, reduced_state; current_time, time_step, outputlevel
+  )
+  if !isnothing(time_step)
+    current_time += time_step
+  end
+  normalize && (reduced_state /= norm(reduced_state))
   spec = nothing
   ortho = isforward(direction) ? "left" : "right"
   drho = nothing
   if noise > 0.0 && isforward(direction)
-    drho = noise * noiseterm(PH, phi, ortho)
+    drho = noise * noiseterm(reduce_operator, reduced_state, ortho)
   end
   spec = replacebond!(
-    psi,
+    state,
     b,
-    phi1;
+    reduced_state;
     maxdim,
     mindim,
     cutoff,
@@ -346,12 +373,12 @@ function region_update!(
   nsite_val::Val{2},
   reverse_step_val::Val{true},
   solver,
-  PH,
-  psi,
+  reduce_operator,
+  state,
   b;
   current_time,
-  outputlevel,
   time_step,
+  outputlevel,
   normalize,
   direction,
   noise,
@@ -362,25 +389,27 @@ function region_update!(
   mindim,
   maxtruncerr,
 )
-  N = length(psi)
+  N = length(state)
   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)
+  set_nsite!(reduce_operator, nsite)
+  position!(reduce_operator, state, b)
+  reduced_state = state[b] * state[b + 1]
+  reduced_state, info = solver(
+    reduce_operator, reduced_state; current_time, time_step, outputlevel
+  )
   current_time += time_step
-  normalize && (phi1 /= norm(phi1))
+  normalize && (reduced_state /= norm(reduced_state))
   spec = nothing
   ortho = isforward(direction) ? "left" : "right"
   drho = nothing
   if noise > 0.0 && isforward(direction)
-    drho = noise * noiseterm(PH, phi, ortho)
+    drho = noise * noiseterm(reduce_operator, phi, ortho)
   end
   spec = replacebond!(
-    psi,
+    state,
     b,
-    phi1;
+    reduced_state;
     maxdim,
     mindim,
     cutoff,
@@ -395,14 +424,16 @@ function region_update!(
     # 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)
+    bond_reduced_state = state[b1]
+    set_nsite!(reduce_operator, nsite - 1)
+    position!(reduce_operator, state, b1)
+    bond_reduced_state, info = solver(
+      reduce_operator, bond_reduced_state; current_time, time_step=-time_step, outputlevel
+    )
     current_time -= time_step
-    normalize && (phi0 ./= norm(phi0))
-    psi[b1] = phi0
-    set_nsite!(PH, nsite)
+    normalize && (bond_reduced_state /= norm(bond_reduced_state))
+    state[b1] = bond_reduced_state
+    set_nsite!(reduce_operator, nsite)
   end
   return current_time, maxtruncerr, spec, info
 end
@@ -411,8 +442,8 @@ function region_update!(
   ::Val{nsite},
   ::Val{reverse_step},
   solver,
-  PH,
-  psi,
+  reduce_operator,
+  state,
   b;
   current_time,
   outputlevel,
diff --git a/src/tdvp.jl b/src/tdvp.jl
index fedfffd..b5759ef 100644
--- a/src/tdvp.jl
+++ b/src/tdvp.jl
@@ -20,11 +20,11 @@ function tdvp_solver(
   solver_maxiter,
   solver_outputlevel,
 )
-  function solver(H, t, psi0; current_time, outputlevel)
-    psi, info = f(
-      H,
+  function solver(operator, t, init; current_time, outputlevel)
+    return f(
+      operator,
       t,
-      psi0;
+      init;
       ishermitian,
       issymmetric,
       tol=solver_tol,
@@ -33,15 +33,59 @@ function tdvp_solver(
       verbosity=solver_outputlevel,
       eager=true,
     )
-    return psi, info
   end
   return solver
 end
 
+function time_step_and_nsteps(t, time_step::Nothing, nsteps::Nothing)
+  return error("Must specify either `time_step`, `nsteps`, or both.")
+end
+
+function time_step_and_nsteps(t, time_step::Nothing, nsteps)
+  return t / nsteps, nsteps
+end
+
+function time_step_and_nsteps(t, time_step, nsteps::Nothing)
+  nsteps, rem = divrem(t, time_step)
+  if rem ≉ 0
+    return error("`t / time_step = $t / $time_step = $(t / time_step)` must be an integer.")
+  end
+  return time_step, Int(nsteps)
+end
+
+function time_step_and_nsteps(t, time_step, nsteps)
+  if time_step * nsteps ≠ t
+    return error(
+      "`t = $t`, `time_step = $time_step`, and `nsteps = $nsteps` must satisfy `time_steps * nsteps == t`, while `time_steps * nsteps = $time_steps * $nsteps`.",
+    )
+  end
+  return time_step, nsteps
+end
+
+"""
+    tdvp(operator, t::Number, init::MPS; time_step, nsteps, kwargs...)
+
+Use the time dependent variational principle (TDVP) algorithm
+to compute `exp(t * operator) * init` using an efficient algorithm based
+on alternating optimization of the MPS tensors and local Krylov
+exponentiation of `operator`.
+
+Specify one of `time_step` or `nsteps`. If they are both specified, they
+must satisfy `time_step * nsteps == t`.
+
+Returns:
+* `state::MPS` - time-evolved MPS
+
+"""
 function tdvp(
-  H,
+  operator,
   t::Number,
-  psi0::MPS;
+  init::MPS;
+  reverse_step=true,
+  time_step=nothing,
+  time_start=0,
+  nsweeps=nothing,
+  nsteps=nsweeps,
   ishermitian=default_ishermitian(),
   issymmetric=default_issymmetric(),
   solver_backend=default_tdvp_solver_backend(),
@@ -52,7 +96,10 @@ function tdvp(
   solver_outputlevel=default_solver_outputlevel(solver_function),
   kwargs...,
 )
-  return tdvp(
+  @show t, time_step, nsteps
+  time_step, nsteps = time_step_and_nsteps(t, time_step, nsteps)
+  @show time_step, nsweeps
+  return alternating_update(
     tdvp_solver(
       solver_function;
       ishermitian,
@@ -62,69 +109,71 @@ function tdvp(
       solver_maxiter,
       solver_outputlevel,
     ),
-    H,
-    t,
-    psi0;
+    operator,
+    init;
+    reverse_step,
+    nsweeps=nsteps,
+    time_step,
     kwargs...,
   )
 end
 
-function tdvp(t::Number, H, psi0::MPS; kwargs...)
-  return tdvp(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
 
-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
 
-"""
-    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.
+## function tdvp(solver, t::Number, H, psi0::MPS; kwargs...)
+##   return tdvp(solver, H, t, psi0; kwargs...)
+## end
 
-Returns:
-* `psi::MPS` - time-evolved MPS
+## function tdvp(solver, H, psi0::MPS, t::Number; kwargs...)
+##   return tdvp(solver, H, t, psi0; kwargs...)
+## end
 
-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
+## """
+##     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/tdvp_sweeps.jl b/src/tdvp_sweeps.jl
index d39ea92..0c9b23c 100644
--- a/src/tdvp_sweeps.jl
+++ b/src/tdvp_sweeps.jl
@@ -1,17 +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
+## 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/update_observer.jl b/src/update_observer.jl
index beed55d..18358f3 100644
--- a/src/update_observer.jl
+++ b/src/update_observer.jl
@@ -2,6 +2,9 @@ function update_observer!(observer; kwargs...)
   return error("Not implemented")
 end
 
+struct EmptyObserver end
+update_observer!(observer::EmptyObserver; kwargs...) = observer
+
 using ITensors.ITensorMPS: ITensorMPS
 function update_observer!(observer::ITensorMPS.AbstractObserver; kwargs...)
   return ITensorMPS.measure!(observer; kwargs...)
diff --git a/test/test_contract_mpo.jl b/test/test_contract_mpo.jl
index 5e085fb..5e974a1 100644
--- a/test/test_contract_mpo.jl
+++ b/test/test_contract_mpo.jl
@@ -50,12 +50,12 @@ using Test: @test, @testset
   # Test with less good initial guess MPS not equal to psi
   psi_guess = copy(psi)
   truncate!(psi_guess; maxdim=2)
-  Hpsi = apply(H, psi; alg="fit", nsweeps=4, init_mps=psi_guess)
+  Hpsi = apply(H, psi; alg="fit", nsweeps=4, init=psi_guess)
   @test ITensors.scalartype(Hpsi) == elt
   @test inner(psit, Hpsi) ≈ inner(psit, H, psi) rtol = 20 * √eps(real(elt))
   # Test with nsite=1
   Hpsi_guess = apply(H, psi; alg="naive", cutoff=1e-4)
-  Hpsi = apply(H, psi; alg="fit", init_mps=Hpsi_guess, nsite=1, nsweeps=2)
+  Hpsi = apply(H, psi; alg="fit", init=Hpsi_guess, nsite=1, nsweeps=2)
   @test ITensors.scalartype(Hpsi) == elt
   scale(::Type{Float32}) = 10^2
   scale(::Type{Float64}) = 10^6
diff --git a/test/test_dmrg_x.jl b/test/test_dmrg_x.jl
index 34ce26b..d8467d7 100644
--- a/test/test_dmrg_x.jl
+++ b/test/test_dmrg_x.jl
@@ -31,9 +31,7 @@ using Test: @test, @testset
   H = MPO(elt, heisenberg(n; h), s)
   initstate = rand(["↑", "↓"], n)
   ψ = MPS(elt, s, initstate)
-  dmrg_x_kwargs = (
-    nsweeps=20, reverse_step=false, normalize=true, maxdim=20, cutoff=1e-10, outputlevel=0
-  )
+  dmrg_x_kwargs = (nsweeps=20, normalize=true, maxdim=20, cutoff=1e-10, outputlevel=0)
   e, ϕ = dmrg_x(ProjMPO(H), ψ; nsite=2, dmrg_x_kwargs...)
   @test ITensors.scalartype(ϕ) == elt
   @test inner(ϕ', H, ϕ) / inner(ϕ, ϕ) ≈ e