[WIP] Start simplifying the necessary OpSum definitions for Models

[mtfishman]

Addresses #61.

[mtfishman]

Here is a minimal example of a VUMPS run with subspace expansion using the changes from this PR:

using ITensors
using ITensorInfiniteMPS

function main()
  initstate(n) = isodd(n) ? "↑" : "↓"
  s = infsiteinds("Electron", 2; initstate, conserve_qns=true)
  ψ = InfMPS(s, initstate)
  H = InfiniteSum{MPO}(Model("hubbard"), s; t=1.0, U=10.0)
  for _ in 1:4
    ψ = subspace_expansion(ψ, H; cutoff=1e-6, maxdim=10)
    ψ = vumps(H, ψ; tol=1e-5, maxiter=20)
  end
  return ψ
end

main()

The main difference here is that the "space shifting" (shifting the QNs of the site indices to make the flux density/flux of the unit cell zero for the state you are interested in targeting) is done automatically by infsiteinds by determining the flux density from the initstate.

In principle, we could add a version infsiteinds("Electron", 2; flux=QN(("Nf",2,-1),("Sz",0)), conserve_qns=true) where you can specify the flux of the unit cell directly, but that isn't included in this PR.

Note that this PR also fixes VUMPS for non-QN conserving states and Hamiltonians.

Another major change in the PR is the interface for defining the Hamiltonian of a model. Here is an example for the new interface:

using ITensors
using ITensorInfiniteMPS

function ITensorInfiniteMPS.unit_cell_terms(::Model"staggered_hubbard"; t₁, U₁, t₂, U₂)
  h₁ = OpSum()
  h₁ += -t₁, "Cdagup", 1, "Cup", 2
  h₁ += -t₁, "Cdagup", 2, "Cup", 1
  h₁ += -t₁, "Cdagdn", 1, "Cdn", 2
  h₁ += -t₁, "Cdagdn", 2, "Cdn", 1
  h₁ += U₁, "Nupdn", 1

  h₂ = OpSum()
  h₂ += -t₂, "Cdagup", 2, "Cup", 3
  h₂ += -t₂, "Cdagup", 3, "Cup", 2
  h₂ += -t₂, "Cdagdn", 2, "Cdn", 3
  h₂ += -t₂, "Cdagdn", 3, "Cdn", 2
  h₂ += U₂, "Nupdn", 2

  return [h₁, h₂]
end

function main()
  initstate(n) = isodd(n) ? "↑" : "↓"
  s = infsiteinds("Electron", 2; initstate, conserve_qns=true)
  ψ = InfMPS(s, initstate)
  H = InfiniteSum{MPO}(Model("staggered_hubbard"), s; t₁=1.0, t₂=1.0, U₁=9.0, U₂=11.0)
  for _ in 1:4
    ψ = subspace_expansion(ψ, H; cutoff=1e-6, maxdim=10)
    ψ = vumps(H, ψ; tol=1e-5, maxiter=20)
  end
  return ψ
end

So instead of the many different overloads required before, you only overload a single function ITensorInfiniteMPS.unit_cell_terms. This should output a single Vector{OpSum} which is the length of the primitive unit cell of the model, where each OpSum corresponds to a term of the Hamiltonian in the primitive unit cell.

This definition is used internally to automatically generate finite MPOs (by repeating the terms of the primitive unit cell and removing any terms that fall outside of the finite system), for example you can directly run DMRG as follows:

function main()
  initstate(n) = isodd(n) ? "↑" : "↓"
  s = siteinds("Electron", 100; conserve_qns=true)
  ψ = randomMPS(s, initstate; linkdims=10)
  H = MPO(Model("staggered_hubbard"), s; t₁=1.0, t₂=1.0, U₁=9.0, U₂=11.0)
  e, ψ = dmrg(H, ψ; nsweeps=10, cutoff=1e-6, maxdim=10)
  return ψ
end

main()

which internally uses the definition ITensorInfiniteMPS.unit_cell_terms above.

@LHerviou I'm interested what you think of the changes. I believe I updated your FQHE code properly but I didn't check carefully.

@hershsingh this should address the issue we were discussed that defining a Model was too complicated before. Happy to take suggestions about the new interface proposed here, i.e. do you think outputting a Vector{OpSum} of the length of the unit cell is an intuitive interface?

[LHerviou]

@mtfishman I will try to have a look this week on the changes on infsiteinds. At worst, it should be possible to just replace the spaces as before, so it should be fine.

I strongly approve of the simplifications for the model creation.

[mtfishman]

@LHerviou thanks for taking a look and glad you approve the new model interface.

In terms of defining custom spaces, I've removed the space keyword argument from infsiteinds for now since it was kind of a bad hack and I'm hoping it won't be necessary in general. The better way to go is to define an overload of ITensors.space for a new SiteType, if the existing site types don't fit your needs. For example, for your FQHE example, I defined a new SiteType:

function ITensors.space(::SiteType"FermionK", pos::Int; p=1, q=1, conserve_momentum=true)
  if !conserve_momentum
    return [QN("Nf", -p) => 1, QN("Nf", q - p) => 1]
  else
    return [
      QN(("Nf", -p), ("NfMom", -p * pos)) => 1,
      QN(("Nf", q - p), ("NfMom", (q - p) * pos)) => 1,
    ]
  end
end

and then shifting to zero flux density is handled by inputting an initstate to infsiteinds, though I don't know if I got the details correct.

[mtfishman]

@LHerviou I just checked that infsiteinds on main outputs this for the FQHE test:

s = Index{Vector{Pair{QN, Int64}}}[(dim=2|id=261|"Fermion,Site,c=1,n=1") <Out>
 1: QN(("Nf",-1),("NfMom",0)) => 1
 2: QN(("Nf",2),("NfMom",3)) => 1, (dim=2|id=112|"Fermion,Site,c=1,n=2") <Out>
 1: QN(("Nf",-1),("NfMom",-1)) => 1
 2: QN(("Nf",2),("NfMom",5)) => 1, (dim=2|id=77|"Fermion,Site,c=1,n=3") <Out>
 1: QN(("Nf",-1),("NfMom",-2)) => 1
 2: QN(("Nf",2),("NfMom",7)) => 1, (dim=2|id=816|"Fermion,Site,c=1,n=4") <Out>
 1: QN(("Nf",-1),("NfMom",-3)) => 1
 2: QN(("Nf",2),("NfMom",9)) => 1, (dim=2|id=867|"Fermion,Site,c=1,n=5") <Out>
 1: QN(("Nf",-1),("NfMom",-4)) => 1
 2: QN(("Nf",2),("NfMom",11)) => 1, (dim=2|id=115|"Fermion,Site,c=1,n=6") <Out>
 1: QN(("Nf",-1),("NfMom",-5)) => 1
 2: QN(("Nf",2),("NfMom",13)) => 1]

while the PR outputs:

s = Index{Vector{Pair{QN, Int64}}}[(dim=2|id=717|"FermionK,Site,c=1,n=1") <Out>
 1: QN(("Nf",-1),("NfMom",0)) => 1
 2: QN(("Nf",2),("NfMom",3)) => 1, (dim=2|id=394|"FermionK,Site,c=1,n=2") <Out>
 1: QN(("Nf",-1),("NfMom",-1)) => 1
 2: QN(("Nf",2),("NfMom",5)) => 1, (dim=2|id=22|"FermionK,Site,c=1,n=3") <Out>
 1: QN(("Nf",-1),("NfMom",-2)) => 1
 2: QN(("Nf",2),("NfMom",7)) => 1, (dim=2|id=306|"FermionK,Site,c=1,n=4") <Out>
 1: QN(("Nf",-1),("NfMom",-3)) => 1
 2: QN(("Nf",2),("NfMom",9)) => 1, (dim=2|id=434|"FermionK,Site,c=1,n=5") <Out>
 1: QN(("Nf",-1),("NfMom",-4)) => 1
 2: QN(("Nf",2),("NfMom",11)) => 1, (dim=2|id=559|"FermionK,Site,c=1,n=6") <Out>
 1: QN(("Nf",-1),("NfMom",-5)) => 1
 2: QN(("Nf",2),("NfMom",13)) => 1]

so it looks like they are the same up to the new SiteType.

[mtfishman]

In the latest commit, I simplified the overload requirements for new models even more. The new interface is that ITensorInfiniteMPS.unit_cell_terms should output a single OpSum containing all of the terms in the unit cell, so the staggered Hubbard model example above becomes:

function ITensorInfiniteMPS.unit_cell_terms(::Model"staggered_hubbard"; t₁, U₁, t₂, U₂)
  h = OpSum()
  h += -t₁, "Cdagup", 1, "Cup", 2
  h += -t₁, "Cdagup", 2, "Cup", 1
  h += -t₁, "Cdagdn", 1, "Cdn", 2
  h += -t₁, "Cdagdn", 2, "Cdn", 1
  h += -t₂, "Cdagup", 2, "Cup", 3
  h += -t₂, "Cdagup", 3, "Cup", 2
  h += -t₂, "Cdagdn", 2, "Cdn", 3
  h += -t₂, "Cdagdn", 3, "Cdn", 2
  h += U₁, "Nupdn", 1
  h += U₂, "Nupdn", 2
  return h
end

The number of sites in the unit cell is deduced from the OpSum automatically, by grouping the terms based on the location of the first site of each term. A corner case that isn't handled right now would be if there were sites in the unit cell that didn't have any terms of the Hamiltonian that had support on those sites, but I think that is pretty rare and could be handled by a future PR, for example with another function you can optionally overload to specify the number of sites in the unit cell explicitly.