Problem using DMRG-X

DMRG and Numerical Methods
1

Hi there,

I’m following the test example of DMRG-X in the TDVP code. As of now the following minimal example does not work:

(Btw by reading the original paper, I’m not sure if this method is completely general as it seems to rely heavily on MBL LIOM properties, but I’d assume your algorithms are not restrictive in terms of physics)

@testset "DMRG-X" begin

  function HB(L)
    ampo = OpSum()
    hop = 1.0
    t = 1.0

    for j=1  :L-1
        
        #println(hop)
        # Hopping
        ampo += -t * hop, "C",j,"Cdag",j+1
        ampo += -t * hop, "C",j+1,"Cdag",j 
    end

    return ampo
  end 

  L = 12
  N = 6

  s = siteinds("Fermion", L, conserve_qns=true)
  H = MPO(HB(L), s)

  initstate = append!([ "Occ" for n=1:N] , ["Emp" for n=1:L -N])
  ψ = randomMPS(s, initstate)

  dmrg_x_kwargs = (
    nsweeps=20, reverse_step=false, normalize=true, maxdim=20, cutoff=1e-10, outputlevel=0
  )

  ϕ = dmrg_x(ProjMPO(H), ψ; nsite=2, dmrg_x_kwargs...)

end

And I get the following errors

DMRG-X: Error During Test at /path/behavior_test.jl:8
  Got exception outside of a @test
  Attempting to contract IndexSet:
  
  ((dim=2|id=346|"Fermion,Site,n=1") <Out>
   1: QN("Nf",0,-1) => 1
   2: QN("Nf",1,-1) => 1, (dim=2|id=823|"Fermion,Site,n=2") <Out>
   1: QN("Nf",0,-1) => 1
   2: QN("Nf",1,-1) => 1, (dim=1|id=859|"Link,l=2") <Out>
   1: QN("Nf",4,-1) => 1)
  
  with IndexSet:
  
  ((dim=1|id=859|"Link,l=2") <Out>
   1: QN("Nf",4,-1) => 1, (dim=2|id=823|"Fermion,Site,n=2") <Out>
   1: QN("Nf",0,-1) => 1
   2: QN("Nf",1,-1) => 1, (dim=2|id=346|"Fermion,Site,n=1") <Out>
   1: QN("Nf",0,-1) => 1
   2: QN("Nf",1,-1) => 1, (dim=2|id=780|"Link,eigen") <In>
   1: QN("Nf",5,-1) => 2)
  
  QN indices must have opposite direction to contract, but indices:
  
  (dim=2|id=346|"Fermion,Site,n=1") <Out>
   1: QN("Nf",0,-1) => 1
   2: QN("Nf",1,-1) => 1
  
  and:
  
  (dim=2|id=346|"Fermion,Site,n=1") <Out>
   1: QN("Nf",0,-1) => 1
   2: QN("Nf",1,-1) => 1
  
  do not have opposite directions.
  Stacktrace:
    [1] error(s::String)
      @ Base ./error.jl:33
    [2] compute_contraction_labels(Ais::Tuple{Index{Vector{Pair{QN, Int64}}}, Index{Vector{Pair{QN, Int64}}}, Index{Vector{Pair{QN, Int64}}}}, Bis::NTuple{4, Index{Vector{Pair{QN, Int64}}}})
      @ ITensors ~/.julia/packages/ITensors/OjQuG/src/indexset.jl:648
    [3] _contract(A::NDTensors.BlockSparseTensor{Float64, 3, Tuple{Index{Vector{Pair{QN, Int64}}}, Index{Vector{Pair{QN, Int64}}}, Index{Vector{Pair{QN, Int64}}}}, NDTensors.BlockSparse{Float64, Vector{Float64}, 3}}, B::NDTensors.BlockSparseTensor{Float64, 4, NTuple{4, Index{Vector{Pair{QN, Int64}}}}, NDTensors.BlockSparse{Float64, Vector{Float64}, 4}})
      @ ITensors ~/.julia/packages/ITensors/OjQuG/src/itensor.jl:1863
    [4] _contract(A::ITensor, B::ITensor)
      @ ITensors ~/.julia/packages/ITensors/OjQuG/src/itensor.jl:1870
    [5] contract(A::ITensor, B::ITensor)
      @ ITensors ~/.julia/packages/ITensors/OjQuG/src/itensor.jl:1974
    [6] *
      @ ~/.julia/packages/ITensors/OjQuG/src/itensor.jl:1961 [inlined]
    [7] dmrg_x_solver(PH::ProjMPO, t::Float64, psi0::ITensor; kwargs::Base.Pairs{Symbol, Real, Tuple{Symbol, Symbol}, NamedTuple{(:current_time, :outputlevel), Tuple{Float64, Int64}}})
      @ ITensorTDVP ~/.julia/packages/ITensorTDVP/p6CK4/src/dmrg_x.jl:5
    [8] tdvp_site_update!(nsite_val::Val{2}, reverse_step_val::Val{false}, solver::typeof(ITensorTDVP.dmrg_x_solver), PH::ProjMPO, psi::MPS, b::Int64; current_time::Float64, outputlevel::Int64, time_step::Float64, normalize::Bool, direction::Base.Order.ForwardOrdering, noise::Float64, which_decomp::Nothing, svd_alg::String, cutoff::Float64, maxdim::Int64, mindim::Int64, maxtruncerr::Float64)
      @ ITensorTDVP ~/.julia/packages/ITensorTDVP/p6CK4/src/tdvp_step.jl:302
    [9] tdvp_site_update!(solver::typeof(ITensorTDVP.dmrg_x_solver), PH::ProjMPO, psi::MPS, b::Int64; nsite::Int64, reverse_step::Bool, current_time::Float64, outputlevel::Int64, time_step::Float64, normalize::Bool, direction::Base.Order.ForwardOrdering, noise::Float64, which_decomp::Nothing, svd_alg::String, cutoff::Float64, maxdim::Int64, mindim::Int64, maxtruncerr::Float64)
      @ ITensorTDVP ~/.julia/packages/ITensorTDVP/p6CK4/src/tdvp_step.jl:155
   [10] tdvp_sweep(direction::Base.Order.ForwardOrdering, solver::Function, PH::ProjMPO, time_step::Float64, psi::MPS; kwargs::Base.Pairs{Symbol, Real, NTuple{11, Symbol}, NamedTuple{(:current_time, :reverse_step, :nsite, :nsweeps, :normalize, :maxdim, :cutoff, :outputlevel, :sweep, :mindim, :noise), Tuple{Float64, Bool, Int64, Int64, Bool, Int64, Float64, Int64, Int64, Int64, Float64}}})
      @ ITensorTDVP ~/.julia/packages/ITensorTDVP/p6CK4/src/tdvp_step.jl:79
   [11] tdvp_step(order::ITensorTDVP.TDVPOrder{2, Base.Order.ForwardOrdering()}, solver::Function, PH::ProjMPO, time_step::Float64, psi::MPS; current_time::Float64, kwargs::Base.Pairs{Symbol, Real, NTuple{10, Symbol}, NamedTuple{(:reverse_step, :nsite, :nsweeps, :normalize, :maxdim, :cutoff, :outputlevel, :sweep, :mindim, :noise), Tuple{Bool, Int64, Int64, Bool, Int64, Float64, Int64, Int64, Int64, Float64}}})
      @ ITensorTDVP ~/.julia/packages/ITensorTDVP/p6CK4/src/tdvp_step.jl:9
   [12] macro expansion
      @ ~/.julia/packages/ITensorTDVP/p6CK4/src/tdvp_generic.jl:84 [inlined]
   [13] macro expansion
      @ ./timing.jl:299 [inlined]
   [14] tdvp(solver::Function, PH::ProjMPO, t::Float64, psi0::MPS; kwargs::Base.Pairs{Symbol, Real, NTuple{7, Symbol}, NamedTuple{(:reverse_step, :nsite, :nsweeps, :normalize, :maxdim, :cutoff, :outputlevel), Tuple{Bool, Int64, Int64, Bool, Int64, Float64, Int64}}})
      @ ITensorTDVP ~/.julia/packages/ITensorTDVP/p6CK4/src/tdvp_generic.jl:83
   [15] #dmrg_x#50
      @ ~/.julia/packages/ITensorTDVP/p6CK4/src/dmrg_x.jl:12 [inlined]
   [16] macro expansion
      @ ~/Desktop/Code/iTensor/mysrc/behavior_test.jl:39 [inlined]
   [17] macro expansion
      @ /Applications/Julia-1.7.app/Contents/Resources/julia/share/julia/stdlib/v1.7/Test/src/Test.jl:1283 [inlined]
   [18] top-level scope
      @ ~/Desktop/Code/iTensor/mysrc/behavior_test.jl:10
   [19] eval
      @ ./boot.jl:373 [inlined]
   [20] include_string(mapexpr::typeof(identity), mod::Module, code::String, filename::String)
      @ Base ./loading.jl:1196
   [21] include_string(m::Module, txt::String, fname::String)
      @ Base ./loading.jl:1206
   [22] invokelatest(::Any, ::Any, ::Vararg{Any}; kwargs::Base.Pairs{Symbol, Union{}, Tuple{}, NamedTuple{(), Tuple{}}})
      @ Base ./essentials.jl:716
   [23] invokelatest(::Any, ::Any, ::Vararg{Any})
      @ Base ./essentials.jl:714
  + further VSCode env errors
Test Summary: | Error  Total
DMRG-X        |     1      1
ERROR: Some tests did not pass: 0 passed, 0 failed, 1 errored, 0 broken.
2

It looks like we just hadn’t tested out the DMRG-X code for QN conserving tensors, could you raise an issue about it at Issues · ITensor/ITensorTDVP.jl · GitHub?

The part that is particular to the physics in some more abstract way is that currently the solver works by diagonalizing the local Hamiltonian and picking the eigenstate that has the largest overlap with the previous state. This is inspired by the l-bit view of MBL, where you try to find an eigenstate that is close to a certain product state.

The particular details of the solver may depend on the application (for example, it could be modified to target a particular eigenvalue, which I have done for other applications).

The implementation is very simple: ITensorTDVP.jl/dmrg_x.jl at ae85e5e854559f59482b94ff5e762e634ca4a797 · ITensor/ITensorTDVP.jl · GitHub so you could easily define your own solver for your application (though of course we should fix the bug that you’ve come across with QN conserving tensors)!

3

Thank you! I have raised an issue and will look into specific applications.