Error when using tdvp

oskar-novak · December 5, 2025, 5:33am

I’m getting the following error when using this function:

function linevolve(theta,p)
  # 1. Define L and capture the exact indices L operates on
  os,sites=build_discrete_dephasing_channel_opsum(6,p*theta,theta)
  L_mpo=MPO(os,sites) 
  L_mpo = combine_sites(L_mpo, sites)

  # 2. Prepare the probe state on temporary indices
  s_temp=siteinds("Qubit",3)
  state = fill("0", 3)
  probe_temp = productMPS(s_temp, state) 
  gates::Vector{ITensor} = []
  push!(gates, op("H", s_temp[1])) 
  push!(gates, op("CNOT", s_temp[1], s_temp[2])) 
  push!(gates, op("CNOT", s_temp[2], s_temp[3]))
  probe_temp= apply(gates, probe_temp)
  
  rho0_mpo = outer(probe_temp', probe_temp) 
  
  # 3. Purify, which introduces "doubled" tags/new Index IDs
  rho0_mps = purify_MPO_to_MPS(rho0_mpo)
  
  for i in 1:3
    correct_site=inds(L_mpo[i])
    probe_site=inds(rho0_mps[i])
    for c in correct_site
      if plev(c)==1
        C=delta(c,probe_site[end])
        
        rho0_mps[i] *= C
      end
    end
  end
  
  
  println(siteinds(rho0_mps))
  println(siteinds(L_mpo))
  
 
  nsteps = 100
  t=1.0
  
  # Time evolve the density matrix. Indices are now guaranteed to match.
  probe_evolved = tdvp(L_mpo, t, rho0_mps; nsteps=nsteps, maxdim=50, cutoff=1e-8)
  return probe_evolved
end

ERROR: 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=4|id=356|"CMB,Link"), (dim=4|id=356|"CMB,Link")'', (dim=4|id=971|"Link,l=2")', (dim=4|id=451|"CMB,Link"), (dim=4|id=405|"CMB,Link")) 

v inds: ((dim=4|id=405|"CMB,Link")', (dim=4|id=451|"CMB,Link")', (dim=4|id=971|"Link,l=2"))
Stacktrace:
  [1] error(s::String)
    @ Base ./error.jl:44
  [2] product(P::ProjMPO, v::ITensor)
    @ ITensorMPS ~/.julia/packages/ITensorMPS/kjso2/src/abstractprojmpo/abstractprojmpo.jl:73
  [3] AbstractProjMPO
    @ ~/.julia/packages/ITensorMPS/kjso2/src/abstractprojmpo/abstractprojmpo.jl:87 [inlined]
  [4] apply
    @ ~/.julia/packages/KrylovKit/ZcdRg/src/apply.jl:2 [inlined]
  [5] expintegrator(A::ProjMPO, t::Float64, u::Tuple{…}, alg::KrylovKit.Arnoldi{…})
    @ KrylovKit ~/.julia/packages/KrylovKit/ZcdRg/src/matrixfun/expintegrator.jl:109
  [6] expintegrator
    @ ~/.julia/packages/KrylovKit/ZcdRg/src/matrixfun/expintegrator.jl:102 [inlined]
  [7] expintegrator(::ProjMPO, ::Float64, ::ITensor; kwargs::@Kwargs{})
    @ KrylovKit ~/.julia/packages/KrylovKit/ZcdRg/src/matrixfun/expintegrator.jl:98
  [8] expintegrator
    @ ~/.julia/packages/KrylovKit/ZcdRg/src/matrixfun/expintegrator.jl:94 [inlined]
  [9] exponentiate
    @ ~/.julia/packages/KrylovKit/ZcdRg/src/matrixfun/exponentiate.jl:83 [inlined]
 [10] #exponentiate_updater#922
    @ ~/.julia/packages/ITensorMPS/kjso2/src/solvers/tdvp.jl:5 [inlined]
 [11] exponentiate_updater
    @ ~/.julia/packages/ITensorMPS/kjso2/src/solvers/tdvp.jl:4 [inlined]
 [12] region_update!(nsite_val::Val{…}, reverse_step_val::Val{…}, 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/kjso2/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/kjso2/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/kjso2/src/solvers/sweep_update.jl:107
 [15] sweep_update(order::ITensorMPS.TDVPOrder{…}, reduced_operator::ProjMPO, state::MPS; current_time::Float64, time_step::Float64, kwargs::@Kwargs{…})
    @ ITensorMPS ~/.julia/packages/ITensorMPS/kjso2/src/solvers/sweep_update.jl:25
 [16] macro expansion
    @ ~/.julia/packages/ITensorMPS/kjso2/src/solvers/alternating_update.jl:69 [inlined]
 [17] macro expansion
    @ ./timing.jl:461 [inlined]
 [18] alternating_update(operator::MPO, init::MPS; updater::Function, updater_kwargs::@NamedTuple{}, nsweeps::Int64, checkdone::Returns{…}, 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::Int64, mindim::Int64, cutoff::Float64, noise::Bool)
    @ ITensorMPS ~/.julia/packages/ITensorMPS/kjso2/src/solvers/alternating_update.jl:68
 [19] alternating_update
    @ ~/.julia/packages/ITensorMPS/kjso2/src/solvers/alternating_update.jl:23 [inlined]
 [20] #tdvp#924
    @ ~/.julia/packages/ITensorMPS/kjso2/src/solvers/tdvp.jl:86 [inlined]
 [21] linevolve(theta::Float64, p::Float64)
    @ Main ~/Documents/UMASS summer 2025/ITensorHandsOn/Quantum Metrology Project/ghzsim.jl:213
 [22] top-level scope
    @ ~/Documents/UMASS summer 2025/ITensorHandsOn/Quantum Metrology Project/ghzsim.jl:223
Some type information was truncated. Use `show(err)` to see complete types.

I’m curious what I’m doing wrong, since I thought the delta tensor would help me match indices between the MPS and MPO.
Thanks!

oskar-novak · December 5, 2025, 5:36am

Here are the helper functions if needed:

function build_discrete_dephasing_channel_opsum(N, P_flip,theta)
    sites = siteinds("Qubit", N)
    os = OpSum()
    range=Int(N/2)
    I_minus_P = 1.0 - P_flip
    for i in 0:(2^range - 1)
        coefficient = 1.0
        ops_tuple = ()
        current_combination_id = i 
        for j in 1:2:N
            is_flip = (current_combination_id >> (j-1)) & 1 == 1
            if is_flip
                ops_tuple = (ops_tuple..., "Z", j, "Z", j+1)
                coefficient *= P_flip
            else
                coefficient *= I_minus_P
            end
        end
        if isempty(ops_tuple)
           for j in 1:2:N
            os += -1im*theta, "Z", j
            os += +1im*theta, "Z", j+1
           end
           
        else
             os += coefficient, ops_tuple... 
        end
    end
    os += -1, "Id", 1
    return os, sites
end

function combine_sites(L_mpo::MPO, sites::Vector{Index{Int64}})
    N = length(sites)
    
  
    L_tensor = contract(L_mpo) 
    
  
    combiner_tensors = ITensor[]
    for j in 1:2:N
        s1 = sites[j]
        s2 = sites[j+1]
        
        C = combiner(s1, s2)
        push!(combiner_tensors, C)
        
       
        L_tensor = L_tensor * C
      
        L_tensor = L_tensor * prime(C)
    end
    
  
    combined_site_indices = [combinedind(C) for C in combiner_tensors]

  
   
    L_combined_mpo = MPO(L_tensor, combined_site_indices)
    
    return L_combined_mpo
end

function purify_MPO_to_MPS(M::MPO)::MPS
    N = length(M)
    s_unprimed_vec = [siteind(M, j, plev=0) for j in 1:N]
    s_primed_vec = [siteind(M, j, plev=1) for j in 1:N]

    
    doubled_sites = [addtags(s_unprimed_vec[j], "doubled") for j in 1:N]

    doubled_mps_tensors = ITensor[]
    
    for j in 1:N
        mpo_tensor = M[j]
        s_unprimed = s_unprimed_vec[j]
        s_primed = s_primed_vec[j]
        s_doubled_new = doubled_sites[j]

        
        bell_state = ITensor(s_doubled_new, s_primed)
        for k in 1:dim(s_doubled_new)
            bell_state[s_doubled_new(k), s_primed(k)] = 1.0
        end

        
        combined_tensor = mpo_tensor * bell_state
        
        
        C = combiner(s_unprimed, s_doubled_new)
      
        
        output_tensor = combined_tensor * C
         
        
        push!(doubled_mps_tensors, output_tensor)
        
    end
    
   
    purified_psi = MPS(doubled_mps_tensors)
    
    return normalize(purified_psi)
end

Topic		Replies	Views
The ProjMPO-ITensor product error message inside tdvp ITensor Julia Questions mpo , tdvp	6	79	May 4, 2025
Generalization of the TDVP function DMRG and Numerical Methods julia , mpo	5	665	May 24, 2023
tdvp time developing ITensor Julia Questions julia , mps	4	779	December 28, 2022
Combined re-orthogonalization of an MPS and repositioning of a ProjMPO ITensor Julia Questions julia , mps , mpo	2	370	April 11, 2023
BoundsError: attempt to access 0-element Vector{Int64} at index [1] ITensor Julia Questions julia , mps , mpo	3	84	March 25, 2025

Error when using tdvp

Related topics