diff --git a/Project.toml b/Project.toml
index cfbf7cff..b96d1682 100644
--- a/Project.toml
+++ b/Project.toml
@@ -1,6 +1,6 @@
 name = "ITensorMPS"
 uuid = "0d1a4710-d33b-49a5-8f18-73bdf49b47e2"
-version = "0.3.44"
+version = "0.3.45"
 authors = ["Matthew Fishman <mfishman@flatironinstitute.org>", "Miles Stoudenmire <mstoudenmire@flatironinstitute.org>"]
 
 [workspace]
diff --git a/src/mpo.jl b/src/mpo.jl
index 566a5163..34ec5d2d 100644
--- a/src/mpo.jl
+++ b/src/mpo.jl
@@ -637,6 +637,23 @@ function ITensors.contract(A::MPO, ψ::MPS; alg = nothing, method = alg, kwargs.
     return contract(Algorithm(alg), A, ψ; kwargs...)
 end
 
+function _zipup_docstring(rhs_type::String)
+    rhs_name = rhs_type == "MPO" ? "B" : "ψ"
+    return """    - `"zipup"`: The MPO and $rhs_type tensors are contracted then truncated at each site
+    without enforcing the appropriate orthogonal gauge. Once this sweep is complete a call
+    to `truncate!` occurs. Because the initial truncation is not locally optimal it is
+    recommended to use a loose `cutoff` and `maxdim` and then pass the desired truncation
+    parameters to the locally optimal `truncate!` sweep via the additional keyword argument
+    `truncate_kwargs`. A set of parameters suggested in [^Paeckel2019] is
+    `contract(A, $rhs_name; alg="zipup", cutoff=cutoff / 10, maxdim=2 * maxdim, truncate_kwargs=(; cutoff, maxdim))`."""
+end
+
+function _paeckel2019_citation_docstring()
+    return """
+    [^Paeckel2019]: Time-evolution methods for matrix-product states. Sebastian Paeckel et al. [arXiv:1901.05824](https://arxiv.org/abs/1901.05824)
+    """
+end
+
 contract_mpo_mps_doc = """
     contract(ψ::MPS, A::MPO; kwargs...) -> MPS
     *(::MPS, ::MPO; kwargs...) -> MPS
@@ -659,8 +676,7 @@ convenience function `apply`:
 apply(A, x; kwargs...) = replaceprime(contract(A, x; kwargs...), 2 => 1)`.
 ```
 
-Choose the method with the `method` keyword, for example
-`"densitymatrix"` and `"naive"`.
+Choose the algorithm with the `alg` keyword argument (or the deprecated `method` alias).
 
 # Keywords
 - `cutoff::Float64=1e-13`: the cutoff value for truncating the density matrix
@@ -669,14 +685,18 @@ Choose the method with the `method` keyword, for example
 - `maxdim::Int=maxlinkdim(A) * maxlinkdim(ψ))`: the maximal bond dimension of the results MPS.
 - `mindim::Int=1`: the minimal bond dimension of the resulting MPS.
 - `normalize::Bool=false`: whether or not to normalize the resulting MPS.
-- `method::String="densitymatrix"`: the algorithm to use for the contraction.
-   Currently the options are "densitymatrix", where the network formed by the
-   MPO and MPS is squared and contracted down to a density matrix which is
-   diagonalized iteratively at each site, and "naive", where the MPO and MPS
-   tensor are contracted exactly at each site and then a truncation of the
-   resulting MPS is performed.
+- `alg::String="densitymatrix"`: the algorithm to use for the contraction. Options:
+    - `"densitymatrix"`: The network formed by the MPO and MPS is squared and contracted down
+      to a density matrix which is diagonalized iteratively at each site.
+    - `"naive"`: The MPO and MPS tensors are contracted exactly at each site and then a
+      truncation of the resulting MPS is performed.
+$(_zipup_docstring("MPS"))
+- `truncate_kwargs::NamedTuple`: (for `alg="zipup"` only) truncation parameters passed to
+   the final `truncate!` sweep. Defaults to `(; cutoff, maxdim, mindim)`.
 
 See also [`apply`](@ref).
+
+$(_paeckel2019_citation_docstring())
 """
 
 @doc """
@@ -840,28 +860,34 @@ end
 function ITensors.contract(
         ::Algorithm"zipup",
         A::MPO,
-        B::MPO;
+        B::AbstractMPS;
         cutoff = 1.0e-14,
         maxdim = maxlinkdim(A) * maxlinkdim(B),
         mindim = 1,
+        truncate_kwargs = (; cutoff, maxdim, mindim),
         kwargs...
     )
     if hassameinds(siteinds, A, B)
         error(
-            "In `contract(A::MPO, B::MPO)`, MPOs A and B have the same site indices. The indices of the MPOs in the contraction are taken literally, and therefore they should only share one site index per site so the contraction results in an MPO. You may want to use `replaceprime(contract(A', B), 2 => 1)` or `apply(A, B)` which automatically adjusts the prime levels assuming the input MPOs have pairs of primed and unprimed indices."
+            "In `contract(A::MPO, B::AbstractMPS)`, A and B have the same site indices. " *
+                "The indices in the contraction are taken literally, and therefore they should " *
+                "only share one site index per site so the contraction results in an MPO or MPS. " *
+                "You may want to use `replaceprime(contract(A', B), 2 => 1)` or `apply(A, B)` " *
+                "which automatically adjusts the prime levels assuming the inputs have pairs of " *
+                "primed and unprimed indices."
         )
     end
     N = length(A)
     N != length(B) &&
-        throw(DimensionMismatch("lengths of MPOs A ($N) and B ($(length(B))) do not match"))
+        throw(DimensionMismatch("lengths of A ($N) and B ($(length(B))) do not match"))
     # Special case for a single site
-    N == 1 && return MPO([A[1] * B[1]])
+    N == 1 && return typeof(B)([A[1] * B[1]])
     A = orthogonalize(A, 1)
     B = orthogonalize(B, 1)
     A = sim(linkinds, A)
     sA = siteinds(uniqueinds, A, B)
     sB = siteinds(uniqueinds, B, A)
-    C = MPO(N)
+    C = typeof(B)(N)
     lCᵢ = Index[]
     R = ITensor(true)
     for i in 1:(N - 2)
@@ -892,7 +918,7 @@ function ITensors.contract(
         mindim,
         kwargs...
     )
-    truncate!(C; kwargs...)
+    truncate!(C; truncate_kwargs...)
     return C
 end
 
@@ -963,16 +989,23 @@ C = apply(A, B; alg="naive", truncate=false)
    in general you should set a `cutoff` value.
 - `maxdim::Int=maxlinkdim(A) * maxlinkdim(B))`: the maximal bond dimension of the results MPS.
 - `mindim::Int=1`: the minimal bond dimension of the resulting MPS.
-- `alg="zipup"`: Either `"zipup"` or `"naive"`. `"zipup"` contracts pairs of
-   site tensors and truncates with SVDs in a sweep across the sites, while `"naive"`
-   first contracts pairs of tensor exactly and then truncates at the end if `truncate=true`.
-- `truncate=true`: Enable or disable truncation. If `truncate=false`, ignore
-   other truncation parameters like `cutoff` and `maxdim`. This is most relevant
-   for the `"naive"` version, if you just want to contract the tensors pairwise
-   exactly. This can be useful if you are contracting MPOs that have diverging
-   norms, such as MPOs originating from sums of local operators.
+- `alg::String="zipup"`: the algorithm to use for the contraction. Options:
+    - `"zipup"`: Contracts pairs of site tensors and truncates with SVDs in a sweep across the
+      sites, then performs a final `truncate!` sweep.
+    - `"naive"`: Contracts pairs of tensors exactly and then truncates at the end if
+      `truncate=true`.
+$(_zipup_docstring("MPO"))
+- `truncate_kwargs::NamedTuple`: (for `alg="zipup"` only) truncation parameters passed to
+   the final `truncate!` sweep. Defaults to `(; cutoff, maxdim, mindim)`.
+- `truncate=true`: (for `alg="naive"` only) Enable or disable truncation. If
+   `truncate=false`, ignore other truncation parameters like `cutoff` and `maxdim`. This is
+   most relevant if you just want to contract the tensors pairwise exactly. This can be useful
+   if you are contracting MPOs that have diverging norms, such as MPOs originating from sums
+   of local operators.
 
 See also [`apply`](@ref) for details about the arguments available.
+
+$(_paeckel2019_citation_docstring())
 """
 
 @doc """
diff --git a/test/base/test_mpo.jl b/test/base/test_mpo.jl
index 22a9402a..6e629e59 100644
--- a/test/base/test_mpo.jl
+++ b/test/base/test_mpo.jl
@@ -254,22 +254,23 @@ end
         @test_throws DimensionMismatch error_contract(phi, K, badpsi)
     end
 
-    @testset "contract" begin
+    @testset "contract(::MPO, ::MPS; alg=$alg)" for alg in
+        ["densitymatrix", "naive", "zipup"]
         phi = random_mps(sites)
         K = random_mpo(sites)
         @test maxlinkdim(K) == 1
         psi = random_mps(sites)
-        psi_out = contract(K, psi; maxdim = 1)
+        psi_out = contract(K, psi; alg, maxdim = 1)
         @test inner(phi', psi_out) ≈ inner(phi', K, psi)
-        psi_out = contract(psi, K; maxdim = 1)
+        psi_out = contract(psi, K; alg, maxdim = 1)
         @test inner(phi', psi_out) ≈ inner(phi', K, psi)
         psi_out = psi * K
         @test inner(phi', psi_out) ≈ inner(phi', K, psi)
-        @test_throws MethodError contract(K, psi; method = "fakemethod")
+        @test_throws MethodError contract(K, psi; alg = "fakemethod")
 
         badsites = [Index(2, "Site") for n in 1:(N + 1)]
         badpsi = random_mps(badsites)
-        @test_throws DimensionMismatch contract(K, badpsi)
+        @test_throws DimensionMismatch contract(K, badpsi; alg)
 
         # make bigger random MPO...
         for link_dim in 2:5
@@ -301,8 +302,9 @@ end
             orthogonalize!(psi, 1; maxdim = link_dim)
             orthogonalize!(K, 1; maxdim = link_dim)
             orthogonalize!(phi, 1; normalize = true, maxdim = link_dim)
-            psi_out =
-                contract(deepcopy(K), deepcopy(psi); maxdim = 10 * link_dim, cutoff = 0.0)
+            psi_out = contract(
+                deepcopy(K), deepcopy(psi); alg, maxdim = 10 * link_dim, cutoff = 0.0
+            )
             @test inner(phi', psi_out) ≈ inner(phi', K, psi)
         end
     end
@@ -369,14 +371,14 @@ end
         @test maxlinkdim(H) ≤ maxlinkdim(H₁) + maxlinkdim(H₂)
     end
 
-    @testset "contract(::MPO, ::MPO)" begin
+    @testset "contract(::MPO, ::MPO; alg=$alg)" for alg in ["naive", "zipup"]
         psi = random_mps(sites)
         K = random_mpo(sites)
         L = random_mpo(sites)
         @test maxlinkdim(K) == 1
         @test maxlinkdim(L) == 1
-        KL = contract(prime(K), L; maxdim = 1)
-        psi_kl_out = contract(prime(K), contract(L, psi; maxdim = 1); maxdim = 1)
+        KL = contract(prime(K), L; alg, maxdim = 1)
+        psi_kl_out = contract(prime(K), contract(L, psi; alg, maxdim = 1); alg, maxdim = 1)
         @test inner(psi'', KL, psi) ≈ inner(psi'', psi_kl_out) atol = 5.0e-3
 
         # where both K and L have differently labelled sites
@@ -388,15 +390,15 @@ end
             replaceind!(K[ii], sites[ii]', othersitesk[ii])
             replaceind!(L[ii], sites[ii]', othersitesl[ii])
         end
-        KL = contract(K, L; maxdim = 1)
+        KL = contract(K, L; alg, maxdim = 1)
         psik = random_mps(othersitesk)
         psil = random_mps(othersitesl)
-        psi_kl_out = contract(K, contract(L, psil; maxdim = 1); maxdim = 1)
+        psi_kl_out = contract(K, contract(L, psil; alg, maxdim = 1); alg, maxdim = 1)
         @test inner(psik, KL, psil) ≈ inner(psik, psi_kl_out) atol = 5.0e-3
 
         badsites = [Index(2, "Site") for n in 1:(N + 1)]
         badL = random_mpo(badsites)
-        @test_throws DimensionMismatch contract(K, badL)
+        @test_throws DimensionMismatch contract(K, badL; alg)
     end
 
     @testset "*(::MPO, ::MPO)" begin