The ProjMPO-ITensor product error message inside tdvp

Hi,

I am using the tdvp function and it gives me the following error:

ProjMPO-ITensor product P*v is not equal to the order of the ITensor v.

I’ve attached the full error message in the end. Note that I have not tweeked anything inside tdvp. In other words, I am using tdvp straightly.

Here is my julia code. I have a psi loaded from a h5 file. Then I do the following:

n = 9
sites = siteinds(psi)
cdag = op("Cdag", sites, n)
psi[n] = psi[n]*cdag

psi = tdvp(H, t, psi; time_step=dt)

I am a bit stuck on debugging this, and I am wondering if there are any suggestions on what could possibly lead to this error. Maybe it is worth mentioning that if I commented out psi[n] = psi[n]*cdag, then the whole thing runs without error.
The full error message is as follows:

ERROR: LoadError: The order of the ProjMPO-ITensor product P*v is not equal to the order of the ITensor v, this is probably due to an index mismatch.
Common reasons for this error: 
(1) You are trying to multiply the ProjMPO with the 2-site wave-function at the wrong position.
(2) `orthogonalize!` was called, changing the MPS without updating the ProjMPO.

P*v inds: ((dim=2|id=845|"Fermion,Site,n=9"), (dim=2|id=845|"Fermion,Site,n=9")'', (dim=4|id=535|"sum")', (dim=2|id=978|"Fermion,Site,n=2")', (dim=2|id=288|"Fermion,Site,n=1")') 

v inds: ((dim=2|id=288|"Fermion,Site,n=1"), (dim=2|id=978|"Fermion,Site,n=2"), (dim=4|id=535|"sum"))
Stacktrace:
  [1] error(s::String)
    @ Base ./error.jl:35
  [2] product(P::ProjMPO, v::ITensor)
    @ ITensorMPS ~/.julia/packages/ITensorMPS/hghpU/src/abstractprojmpo/abstractprojmpo.jl:73
  [3] AbstractProjMPO
    @ ~/.julia/packages/ITensorMPS/hghpU/src/abstractprojmpo/abstractprojmpo.jl:87 [inlined]
  [4] apply
    @ ~/.julia/packages/KrylovKit/jC5gU/src/apply.jl:2 [inlined]
  [5] expintegrator(A::ProjMPO, t::Float64, u::Tuple{ITensor, ITensor}, alg::KrylovKit.Arnoldi{KrylovKit.ModifiedGramSchmidt2, Float64})
    @ KrylovKit ~/.julia/packages/KrylovKit/jC5gU/src/matrixfun/expintegrator.jl:109
  [6] expintegrator
    @ ~/.julia/packages/KrylovKit/jC5gU/src/matrixfun/expintegrator.jl:102 [inlined]
  [7] expintegrator(::ProjMPO, ::Float64, ::ITensor; kwargs::@Kwargs{})
    @ KrylovKit ~/.julia/packages/KrylovKit/jC5gU/src/matrixfun/expintegrator.jl:98
  [8] expintegrator
    @ ~/.julia/packages/KrylovKit/jC5gU/src/matrixfun/expintegrator.jl:94 [inlined]
  [9] exponentiate
    @ ~/.julia/packages/KrylovKit/jC5gU/src/matrixfun/exponentiate.jl:83 [inlined]
 [10] #exponentiate_updater#765
    @ ~/.julia/packages/ITensorMPS/hghpU/src/solvers/tdvp.jl:5 [inlined]
 [11] exponentiate_updater
    @ ~/.julia/packages/ITensorMPS/hghpU/src/solvers/tdvp.jl:4 [inlined]
 [12] region_update!(nsite_val::Val{2}, reverse_step_val::Val{true}, reduced_operator::ProjMPO, state::MPS, b::Int64; updater::typeof(ITensorMPS.exponentiate_updater), updater_kwargs::@NamedTuple{}, current_time::Float64, time_step::Float64, outputlevel::Int64, normalize::Bool, direction::Base.Order.ForwardOrdering, noise::Bool, which_decomp::Nothing, svd_alg::Nothing, cutoff::Float64, maxdim::Int64, mindim::Int64, maxtruncerr::Float64)
    @ ITensorMPS ~/.julia/packages/ITensorMPS/hghpU/src/solvers/sweep_update.jl:408
 [13] region_update!(reduced_operator::ProjMPO, state::MPS, b::Int64; updater::Function, updater_kwargs::@NamedTuple{}, nsite::Int64, reverse_step::Bool, current_time::Float64, outputlevel::Int64, time_step::Float64, normalize::Bool, direction::Base.Order.ForwardOrdering, noise::Bool, which_decomp::Nothing, svd_alg::Nothing, cutoff::Float64, maxdim::Int64, mindim::Int64, maxtruncerr::Float64)
    @ ITensorMPS ~/.julia/packages/ITensorMPS/hghpU/src/solvers/sweep_update.jl:187
 [14] sub_sweep_update(direction::Base.Order.ForwardOrdering, reduced_operator::ProjMPO, state::MPS; updater::Function, updater_kwargs::@NamedTuple{}, which_decomp::Nothing, svd_alg::Nothing, sweep::Int64, current_time::Float64, time_step::Float64, nsite::Int64, reverse_step::Bool, normalize::Bool, observer!::ITensorMPS.EmptyObserver, outputlevel::Int64, maxdim::Int64, mindim::Int64, cutoff::Float64, noise::Bool)
    @ ITensorMPS ~/.julia/packages/ITensorMPS/hghpU/src/solvers/sweep_update.jl:107
 [15] sweep_update(order::ITensorMPS.TDVPOrder{2, Base.Order.ForwardOrdering()}, reduced_operator::ProjMPO, state::MPS; current_time::Float64, time_step::Float64, kwargs::@Kwargs{updater::typeof(ITensorMPS.exponentiate_updater), updater_kwargs::@NamedTuple{}, nsite::Int64, reverse_step::Bool, sweep::Int64, observer!::ITensorMPS.EmptyObserver, normalize::Bool, outputlevel::Int64, maxdim::Int64, mindim::Int64, cutoff::Float64, noise::Bool})
    @ ITensorMPS ~/.julia/packages/ITensorMPS/hghpU/src/solvers/sweep_update.jl:25
 [16] macro expansion
    @ ~/.julia/packages/ITensorMPS/hghpU/src/solvers/alternating_update.jl:69 [inlined]
 [17] macro expansion
    @ ./timing.jl:421 [inlined]
 [18] alternating_update(operator::MPO, init::MPS; updater::Function, updater_kwargs::@NamedTuple{}, nsweeps::Int64, checkdone::Returns{Bool}, write_when_maxdim_exceeds::Nothing, nsite::Int64, reverse_step::Bool, time_start::Float64, time_step::Float64, order::Int64, observer!::ITensorMPS.EmptyObserver, sweep_observer!::ITensorMPS.EmptyObserver, outputlevel::Int64, normalize::Bool, maxdim::Vector{Int64}, mindim::Int64, cutoff::Vector{Float64}, noise::Bool)
    @ ITensorMPS ~/.julia/packages/ITensorMPS/hghpU/src/solvers/alternating_update.jl:68
 [19] alternating_update
    @ ~/.julia/packages/ITensorMPS/hghpU/src/solvers/alternating_update.jl:23 [inlined]
 [20] #tdvp#767
    @ ~/.julia/packages/ITensorMPS/hghpU/src/solvers/tdvp.jl:86 [inlined]
 [21] top-level scope

Thanks in advance!