Error with conserved quantum numbers

ITensor Julia Questions
,
1

I am running into an error where if I create an MPO, say one with the SiteSet “Fermion”, and I set conserved_nf=true, I get the following error if the MPO gets too large in number of terms:

ERROR: LoadError: DimensionMismatch: array could not be broadcast to match destination
Stacktrace:
  [1] check_broadcast_shape(shp::Tuple{Base.OneTo{Int64}}, Ashp::Tuple{Base.OneTo{Int64}})
    @ Base.Broadcast ./broadcast.jl:540
  [2] check_broadcast_shape(shp::Tuple{Base.OneTo{Int64}, Base.OneTo{Int64}}, Ashp::Tuple{Base.OneTo{Int64}, Base.OneTo{Int64}})
    @ Base.Broadcast ./broadcast.jl:541
  [3] check_broadcast_axes(shp::Tuple{Base.OneTo{Int64}, Base.OneTo{Int64}}, A::Matrix{Float64})
    @ Base.Broadcast ./broadcast.jl:543
  [4] instantiate(bc::Base.Broadcast.Broadcasted{Base.Broadcast.DefaultArrayStyle{2}, Tuple{Base.OneTo{Int64}, Base.OneTo{Int64}}, typeof(identity), Tuple{Matrix{Float64}}})
    @ Base.Broadcast ./broadcast.jl:284
  [5] materialize!(#unused#::Base.Broadcast.ArrayStyle{NDTensors.DenseTensor{Float64, 2, Tuple{Int64, Int64}, NDTensors.Dense{Float64, SubArray{Float64, 1, Vector{Float64}, Tuple{UnitRange{Int64}}, true}}}}, dest::NDTensors.DenseTensor{Float64, 2, Tuple{Int64, Int64}, NDTensors.Dense{Float64, SubArray{Float64, 1, Vector{Float64}, Tuple{UnitRange{Int64}}, true}}}, bc::Base.Broadcast.Broadcasted{Base.Broadcast.DefaultArrayStyle{2}, Nothing, typeof(identity), Tuple{Matrix{Float64}}})
    @ Base.Broadcast ./broadcast.jl:871
  [6] materialize!(dest::NDTensors.DenseTensor{Float64, 2, Tuple{Int64, Int64}, NDTensors.Dense{Float64, SubArray{Float64, 1, Vector{Float64}, Tuple{UnitRange{Int64}}, true}}}, bc::Base.Broadcast.Broadcasted{Base.Broadcast.DefaultArrayStyle{2}, Nothing, typeof(identity), Tuple{Matrix{Float64}}})
    @ Base.Broadcast ./broadcast.jl:868
  [7] qn_svdMPO(os::Sum{Scaled{ComplexF64, Prod{Op}}}, sites::Vector{Index{Vector{Pair{QN, Int64}}}}; kwargs::Base.Pairs{Symbol, Int64, Tuple{Symbol}, NamedTuple{(:maxdim,), Tuple{Int64}}})
    @ ITensors ~/.julia/packages/ITensors/5CAqA/src/physics/autompo/opsum_to_mpo_qn.jl:223
  [8] MPO(os::Sum{Scaled{ComplexF64, Prod{Op}}}, sites::Vector{Index{Vector{Pair{QN, Int64}}}}; splitblocks::Bool, kwargs::Base.Pairs{Symbol, Int64, Tuple{Symbol}, NamedTuple{(:maxdim,), Tuple{Int64}}})
    @ ITensors ~/.julia/packages/ITensors/5CAqA/src/physics/autompo/opsum_to_mpo_generic.jl:246

Does anyone have any idea what might be causing this error? Keep in mind that for a smaller Hamiltonian, I am able to create the MPO. I should also note that I am using AutoMPO to create my Hamiltonian and that I am basically transcribing a molecular Hamiltonian from OpenFermion.

2

Hm, I’m really not sure. It’s puzzling that it works for fewer terms then gives this message for larger numbers of terms, because the message doesn’t seem to say anything about a large size or use of memory or anything like that. Sometimes the real error can be different from what error messages you get.

Two next steps:

  • if it might be possible to provide some steps to reproduce and/or file a bug report on our Issues page, we could try to look into it
  • a workaround could be to split the terms up into two or more groups (especially if there is some natural grouping you can think of) and then make a few separate MPOs. Then you can pass all of these MPOs into our DMRG code like this: dmrg([H1,H2,H3], ...) where I’m showing the case of three MPOs, and it will work the same or even faster than if you used a single MPO. One natural way to split up terms can be hopping versus interaction, and also up electron processes versus down electron processes.
3

Hi Miles,

Thanks for responding to my problem so promptly. Also, congratulations to you and the people at Flatiron on the tensor network Eagle simulation paper! It’s been making some noise in the community.

Anyway, to answer your first point, the Hamiltonian I turning into an MPO is shown below:

ampo += (0.3379505065631, "Cdag", 2, "Cdag", 2, "C", 2, "C", 2)
ampo += (0.0904377249491051, "Cdag", 1, "Cdag", 1, "C", 3, "C", 3)
ampo += (0.3493792455889514, "Cdag", 4, "Cdag", 4, "C", 4, "C", 4)
ampo += (0.0904377249491051, "Cdag", 4, "Cdag", 2, "C", 4, "C", 2)
ampo += (0.3323780270343749, "Cdag", 4, "Cdag", 2, "C", 2, "C", 4)
ampo += (0.0904377249491051, "Cdag", 1, "Cdag", 4, "C", 2, "C", 3)
ampo += (0.0904377249491051, "Cdag", 2, "Cdag", 4, "C", 2, "C", 4)
ampo += (0.0904377249491051, "Cdag", 3, "Cdag", 3, "C", 1, "C", 1)
ampo += (0.3323780270343749, "Cdag", 2, "Cdag", 4, "C", 4, "C", 2)
ampo += (0.3323780270343749, "Cdag", 4, "Cdag", 1, "C", 1, "C", 4)
ampo += (0.0904377249491051, "Cdag", 2, "Cdag", 2, "C", 4, "C", 4)
ampo += (0.3493792455889514, "Cdag", 3, "Cdag", 4, "C", 4, "C", 3)
ampo += (0.0904377249491051, "Cdag", 4, "Cdag", 3, "C", 1, "C", 2)
ampo += (-0.47125568651010585, "Cdag", 4, "C", 4)
ampo += (0.0904377249491051, "Cdag", 1, "Cdag", 3, "C", 1, "C", 3)
ampo += (0.0904377249491051, "Cdag", 3, "Cdag", 4, "C", 2, "C", 1)
ampo += (0.0904377249491051, "Cdag", 4, "Cdag", 1, "C", 3, "C", 2)
ampo += (-0.47125568651010585, "Cdag", 3, "C", 3)
ampo += (0.3323780270343749, "Cdag", 1, "Cdag", 4, "C", 4, "C", 1)
ampo += (0.3379505065631, "Cdag", 1, "Cdag", 1, "C", 1, "C", 1)
ampo += (0.0904377249491051, "Cdag", 2, "Cdag", 1, "C", 3, "C", 4)
ampo += (0.3493792455889514, "Cdag", 4, "Cdag", 3, "C", 3, "C", 4)
ampo += (0.3379505065631, "Cdag", 1, "Cdag", 2, "C", 2, "C", 1)
ampo += (0.3379505065631, "Cdag", 2, "Cdag", 1, "C", 1, "C", 2)
ampo += (-1.2569462599509786, "Cdag", 1, "C", 1)
ampo += (0.0904377249491051, "Cdag", 2, "Cdag", 3, "C", 1, "C", 4)
ampo += (0.3323780270343749, "Cdag", 2, "Cdag", 3, "C", 3, "C", 2)
ampo += (0.3323780270343749, "Cdag", 3, "Cdag", 2, "C", 2, "C", 3)
ampo += (0.0904377249491051, "Cdag", 3, "Cdag", 2, "C", 4, "C", 1)
ampo += (0.0904377249491051, "Cdag", 4, "Cdag", 4, "C", 2, "C", 2)
ampo += (-1.2569462599509786, "Cdag", 2, "C", 2)
ampo += (0.3323780270343749, "Cdag", 1, "Cdag", 3, "C", 3, "C", 1)
ampo += (0.3493792455889514, "Cdag", 3, "Cdag", 3, "C", 3, "C", 3)
ampo += (0.0904377249491051, "Cdag", 3, "Cdag", 1, "C", 3, "C", 1)
ampo += (0.0904377249491051, "Cdag", 1, "Cdag", 2, "C", 4, "C", 3)
ampo += (0.3323780270343749, "Cdag", 3, "Cdag", 1, "C", 1, "C", 3)

It’s basically the Hamiltonian for molecular Hydrogen using 4 orbitals, or sites. To construct this MPO, I am using the FermionSites site type with the following parameters:

sites = siteinds(
            "Fermion",
            n_sites;
            conserve_qns=true,
          )

and then I turn the autoMPO into an MPO using the MPO(ampo, sites) command. What’s even weirder is that the MPO successfully builds if I set conserve_qns=false.

With regards to your second point, I will see if splitting the MPO will fix the issue.

Thanks again,
Jacob

4

Hi Jacob,
Thanks! Joey Tindall and Matt Fishman here at CCQ were already working on that belief propagation method and had it all coded up (in ITensor of course – actually a new overlay library called ITensorNetworks.jl) so Joey just inputted the heavy-hex Ising system and got great results very quickly.

About the MPO, it’s not too surprising to me that the QN versus non-QN case might behave differently since the backend code for those is basically separate. I’m more surprised by how it works up to a point (I guess for smaller numbers of terms) then fails for a larger number, because I don’t see where in the code it would be too sensitive to the number of terms, apart from needing more memory, say.

By the way, are you getting the error even for the explicit case you wrote above? Or only when you have even more terms? Thanks –

Miles

5

Hi Miles,

I am getting the error for the explicit case above. However, I just discovered that I do not get an error if I split the MPO up in the following way:

ampo = AutoMPO()
ampo += (0.3379505065631, "Cdag", 2, "Cdag", 2, "C", 2, "C", 2)
ampo += (0.0904377249491051, "Cdag", 1, "Cdag", 1, "C", 3, "C", 3)
ampo += (0.3493792455889514, "Cdag", 4, "Cdag", 4, "C", 4, "C", 4)
ampo += (0.0904377249491051, "Cdag", 4, "Cdag", 2, "C", 4, "C", 2)
ampo += (0.3323780270343749, "Cdag", 4, "Cdag", 2, "C", 2, "C", 4)
# ampo += (0.0904377249491051, "Cdag", 1, "Cdag", 4, "C", 2, "C", 3)
ampo += (0.0904377249491051, "Cdag", 2, "Cdag", 4, "C", 2, "C", 4)
ampo += (0.0904377249491051, "Cdag", 3, "Cdag", 3, "C", 1, "C", 1)
ampo += (0.3323780270343749, "Cdag", 2, "Cdag", 4, "C", 4, "C", 2)
ampo += (0.3323780270343749, "Cdag", 4, "Cdag", 1, "C", 1, "C", 4)
ampo += (0.0904377249491051, "Cdag", 2, "Cdag", 2, "C", 4, "C", 4)
ampo += (0.3493792455889514, "Cdag", 3, "Cdag", 4, "C", 4, "C", 3)
# ampo += (0.0904377249491051, "Cdag", 4, "Cdag", 3, "C", 1, "C", 2)
ampo += (-0.47125568651010585, "Cdag", 4, "C", 4)
ampo += (0.0904377249491051, "Cdag", 1, "Cdag", 3, "C", 1, "C", 3)
# ampo += (0.0904377249491051, "Cdag", 3, "Cdag", 4, "C", 2, "C", 1)
# ampo += (0.0904377249491051, "Cdag", 4, "Cdag", 1, "C", 3, "C", 2)
ampo += (-0.47125568651010585, "Cdag", 3, "C", 3)
ampo += (0.3323780270343749, "Cdag", 1, "Cdag", 4, "C", 4, "C", 1)
ampo += (0.3379505065631, "Cdag", 1, "Cdag", 1, "C", 1, "C", 1)
# ampo += (0.0904377249491051, "Cdag", 2, "Cdag", 1, "C", 3, "C", 4)
ampo += (0.3493792455889514, "Cdag", 4, "Cdag", 3, "C", 3, "C", 4)
ampo += (0.3379505065631, "Cdag", 1, "Cdag", 2, "C", 2, "C", 1)
ampo += (0.3379505065631, "Cdag", 2, "Cdag", 1, "C", 1, "C", 2)
ampo += (-1.2569462599509786, "Cdag", 1, "C", 1)
# ampo += (0.0904377249491051, "Cdag", 2, "Cdag", 3, "C", 1, "C", 4)
ampo += (0.3323780270343749, "Cdag", 2, "Cdag", 3, "C", 3, "C", 2)
ampo += (0.3323780270343749, "Cdag", 3, "Cdag", 2, "C", 2, "C", 3)
# ampo += (0.0904377249491051, "Cdag", 3, "Cdag", 2, "C", 4, "C", 1)
ampo += (0.0904377249491051, "Cdag", 4, "Cdag", 4, "C", 2, "C", 2)
ampo += (-1.2569462599509786, "Cdag", 2, "C", 2)
ampo += (0.3323780270343749, "Cdag", 1, "Cdag", 3, "C", 3, "C", 1)
ampo += (0.3493792455889514, "Cdag", 3, "Cdag", 3, "C", 3, "C", 3)
ampo += (0.0904377249491051, "Cdag", 3, "Cdag", 1, "C", 3, "C", 1)
# ampo += (0.0904377249491051, "Cdag", 1, "Cdag", 2, "C", 4, "C", 3)
ampo += (0.3323780270343749, "Cdag", 3, "Cdag", 1, "C", 1, "C", 3)

ampo2 = AutoMPO()
ampo2 += (0.0904377249491051, "Cdag", 1, "Cdag", 4, "C", 2, "C", 3)
ampo2 += (0.0904377249491051, "Cdag", 4, "Cdag", 3, "C", 1, "C", 2)
ampo2 += (0.0904377249491051, "Cdag", 3, "Cdag", 4, "C", 2, "C", 1)
ampo2 += (0.0904377249491051, "Cdag", 4, "Cdag", 1, "C", 3, "C", 2)
ampo2 += (0.0904377249491051, "Cdag", 2, "Cdag", 1, "C", 3, "C", 4)
ampo2 += (0.0904377249491051, "Cdag", 2, "Cdag", 3, "C", 1, "C", 4)
ampo2 += (0.0904377249491051, "Cdag", 3, "Cdag", 2, "C", 4, "C", 1)
ampo2 += (0.0904377249491051, "Cdag", 1, "Cdag", 2, "C", 4, "C", 3)

Here, I the terms that are greyed out in ampo are added to a separate MPO called ampo2. I am not sure why splitting the MPO this way works, but what I can say is that those terms that I consolidated into ampo2 are the only terms that contain more than 2 sites, e.g. the fermions scatter from 1 and 4 to 2 and 3 instead of either undergoing electrostatic interaction (for example <1,2|V|2,1>) or exchange between sites (for example <1,2|V|1,2>). Maybe this might give more insight into the problem.

Thank you again Miles for your help. It seems that splitting the MPO as you suggested fixed my problem and I was able to compute the ground state.

6

Glad it’s working! But splitting it is of course a workaround so I’d like to get your whole case working.

I was able to reproduce the same error as you once I input all of the terms. In fact, I could get the error just inputting five specific terms, so I will look further into this minimal case to identify the bug. It’s definitely a bug on our side – you aren’t doing anything wrong. Thanks for pointing it out to us.

Miles

7

Hi @manalo , it turns out this case above was a subtle bug with OpSum. In your input, it had to do with the fact that some of the terms put two \hat{c}_j operators on the same site, which gives zero, but that should not result in a bug, only just that such a term won’t contribute to the Hamiltonian.

We’ve now fixed this issue in a recent patch to ITensors.jl and it will be fixed in the next tagged version.

(Fixed in this pull request).