InfiniteMPOMatrix -> InfiniteMPO conversion

src/infinitempo.jl

@@ -11,3 +11,5 @@ mutable struct InfiniteMPO <: AbstractInfiniteMPS
 end
 
 translator(mpo::InfiniteMPO) = mpo.data.translator
+
+InfiniteMPO(data::CelledVector{ITensor}) = InfiniteMPO(data, 0, size(data, 1), false)

src/infinitempomatrix.jl

@@ -7,6 +7,7 @@ mutable struct InfiniteMPOMatrix <: AbstractInfiniteMPS
 end
 
 translator(mpo::InfiniteMPOMatrix) = mpo.data.translator
+data(mpo::InfiniteMPOMatrix) = mpo.data
 
 # TODO better printing?
 function Base.show(io::IO, M::InfiniteMPOMatrix)
@@ -43,4 +44,111 @@ function InfiniteMPOMatrix(data::Vector{Matrix{ITensor}}, translator::Function)
   return InfiniteMPOMatrix(CelledVector(data, translator), 0, size(data)[1], false)
 end
 
-#nrange(H::InfiniteMPOMatrix) = [size(H[j])[1] - 1 for j in 1:nsites(H)]
+#
+#  Hm should have the form below.  Only link indices are shown.
+#
+#    I         0                     0           0      0
+#  M-->l=n     0                     0           0      0
+#    0    l=n-->M-->l=n-1 ...        0           0      0
+#    :         :                     :           :      :
+#    0         0          ...  l=2-->M-->l=1     0      0
+#    0         0          ...        0        l=1-->M   I
+#
+#  We need to capture the corner I's and the sub-diagonal Ms, and paste them together into an ITensor.
+#  This is facilutated by making all elements of Hm into order(4) tensors by adding dummy Dw=1 indices.
+#  We ITensors.directsum() to join all the blocks into the final ITensor.
+#  This code should be 100% dense/blocks-sparse agnostic.
+#
+function cat_to_itensor(Hm::Matrix{ITensor}, inds_array)::Tuple{ITensor,Index,Index}
+  lx, ly = size(Hm)
+  @assert lx == ly
+  #
+  #  Extract the sub diagonal
+  #
+  Ms = map(n -> Hm[lx - n + 1, ly - n], 1:(lx - 1)) #Get the sub diagonal into an array.
+  N = length(Ms)
+  @assert N == lx - 1
+  #
+  # Convert edge Ms to order 4 ITensors using the dummy index.
+  #
+  il0 = inds_array[1][1]
+  T = eltype(Ms[1])
+  Ms[1] *= onehot(T, il0 => 1)
+  Ms[N] *= onehot(T, dag(il0) => 1) #We don't need distinct index IDs for this.  directsum makes all new index anyway.
+  #
+  # Start direct sums to build the single ITensor.  Order is critical to get the desired form.
+  #
+  H = Hm[1, 1] * onehot(T, il0' => 1) * onehot(T, dag(il0) => 1) #Bootstrap with the I op in the top left of Hm
+  H, il = directsum(H => il0', Ms[N] => inds_array[N][1]) #1-D directsum to put M[N] directly below the I op
+  ir = dag(il0) #setup recusion.
+  for i in (N - 1):-1:1 #2-D directsums to place new blocks below and to the right.
+    H, (il, ir) = directsum(
+      H => (il, ir), Ms[i] => inds_array[i]; tags=["Link,left", "Link,right"]
+    )
+  end
+  IN = Hm[N + 1, N + 1] * onehot(T, dag(il0) => 1) * onehot(T, il => dim(il)) #I op in the bottom left of Hm
+  H, ir = directsum(H => ir, IN => il0; tags="Link,right") #1-D directsum to the put I in the bottom row, to the right of M[1]
+
+  return H, il, ir
+end
+
+function find_all_links(Hms::InfiniteMPOMatrix)
+  is = inds(Hms[1][1, 1]) #site inds
+  ir = noncommonind(Hms[1][2, 1], is) #This op should have one link index pointing to the right.
+  #
+  #  Set up return array of 2-tuples
+  #
+  indexT = typeof(ir)
+  TupleT = NTuple{2,indexT}
+  inds_array = Vector{TupleT}[]
+  #
+  #  Make a dummy index
+  #
+  il0 = Index(ITensors.trivial_space(ir); dir=dir(dag(ir)), tags="Link,l=0")
+  #
+  #  Loop over sites and nonzero matrix elements which are linked into the next
+  #  and previous iMPOMatrix.  
+  #
+  for n in 1:nsites(Hms)
+    Hm = Hms[n]
+    inds_n = TupleT[]
+    lx, ly = size(Hm)
+    @assert lx == ly
+    for iy in (ly - 1):-1:1
+      ix = iy + 1
+      il = ix < lx ? commonind(Hm[ix, iy], Hms[n - 1][ix + 1, iy + 1]) : dag(il0)
+      ir = iy > 1 ? commonind(Hm[ix, iy], Hms[n + 1][ix - 1, iy - 1]) : il0
+      push!(inds_n, (il, ir))
+    end
+    push!(inds_array, inds_n)
+  end
+  return inds_array
+end
+
+#
+#  Hm is the InfiniteMPOMatrix
+#  Hlrs is an array of {ITensor,Index,Index}s, one for each site in the unit cell.
+#  Hi is a CelledVector of ITensors.
+#
+function InfiniteMPO(Hm::InfiniteMPOMatrix)
+  inds_array = find_all_links(Hm)
+  Hlrs = cat_to_itensor.(Hm, inds_array) #return an array of {ITensor,Index,Index}
+  #
+  #  Unpack the array of tuples into three arrays.  And also get an array site indices.
+  #
+  Hi = CelledVector([Hlr[1] for Hlr in Hlrs], translator(Hm))
+  ils = CelledVector([Hlr[2] for Hlr in Hlrs], translator(Hm))
+  irs = CelledVector([Hlr[3] for Hlr in Hlrs], translator(Hm))
+  site_inds = [commoninds(Hm[j][1, 1], Hm[j][end, end])[1] for j in 1:nsites(Hm)]
+  #
+  #  Create new tags with proper cell and link numbers.  Also daisy chain
+  #  all the indices so right index at j = dag(left index at j+1)
+  #
+  for j in 1:nsites(Hm)
+    newTag = "Link,c=$(getcell(site_inds[j])),l=$(getsite(site_inds[j]))"
+    ir = replacetags(irs[j], tags(irs[j]), newTag) #new right index
+    Hi[j] = replaceinds(Hi[j], [irs[j]], [ir])
+    Hi[j + 1] = replaceinds(Hi[j + 1], [ils[j + 1]], [dag(ir)])
+  end
+  return InfiniteMPO(Hi)
+end

src/models/models.jl

@@ -82,6 +82,14 @@ end
 import ITensors: op
 op(::OpName"Zero", ::SiteType, s::Index) = ITensor(s', dag(s))
 
+function InfiniteMPO(model::Model, s::CelledVector; kwargs...)
+  return InfiniteMPO(model, s, translator(s); kwargs...)
+end
+
+function InfiniteMPO(model::Model, s::CelledVector, translator::Function; kwargs...)
+  return InfiniteMPO(InfiniteMPOMatrix(model, s, translator; kwargs...))
+end
+
 function InfiniteMPOMatrix(model::Model, s::CelledVector; kwargs...)
   return InfiniteMPOMatrix(model, s, translator(s); kwargs...)
 end

test/test_iMPOConversions.jl

@@ -0,0 +1,141 @@
+using ITensors
+using ITensorInfiniteMPS
+using Test
+#
+# InfiniteMPO has dangling links at the end of the chain.  We contract these on the outside
+#   with l,r terminating vectors, to make a finite lattice MPO.
+#
+function terminate(h::InfiniteMPO)::MPO
+  Ncell = nsites(h)
+  # left termination vector
+  il0 = commonind(h[1], h[0])
+  l = ITensor(0.0, il0)
+  l[il0 => dim(il0)] = 1.0 #assuming lower reg form in h
+  # right termination vector
+  iln = commonind(h[Ncell], h[Ncell + 1])
+  r = ITensor(0.0, iln)
+  r[iln => 1] = 1.0 #assuming lower reg form in h
+  # build up a finite MPO
+  hf = MPO(Ncell)
+  hf[1] = dag(l) * h[1] #left terminate
+  hf[Ncell] = h[Ncell] * dag(r) #right terminate
+  for n in 2:(Ncell - 1)
+    hf[n] = h[n] #fill in the bulk.
+  end
+  return hf
+end
+#
+# Terminate and then call expect
+# for inf ψ and finite h, which is already supported in src/infinitecanonicalmps.jl
+#
+function ITensors.expect(ψ::InfiniteCanonicalMPS, h::InfiniteMPO)
+  return expect(ψ, terminate(h)) #defer to src/infinitecanonicalmps.jl
+end
+
+#H = ΣⱼΣn (½ S⁺ⱼS⁻ⱼ₊n + ½ S⁻ⱼS⁺ⱼ₊n + SᶻⱼSᶻⱼ₊n)
+function ITensorInfiniteMPS.unit_cell_terms(::Model"heisenbergNNN"; NNN::Int64)
+  opsum = OpSum()
+  for n in 1:NNN
+    J = 1.0 / n
+    opsum += J * 0.5, "S+", 1, "S-", 1 + n
+    opsum += J * 0.5, "S-", 1, "S+", 1 + n
+    opsum += J, "Sz", 1, "Sz", 1 + n
+  end
+  return opsum
+end
+
+function ITensorInfiniteMPS.unit_cell_terms(::Model"hubbardNNN"; NNN::Int64)
+  U::Float64 = 0.25
+  t::Float64 = 1.0
+  V::Float64 = 0.5
+  opsum = OpSum()
+  opsum += (U, "Nupdn", 1)
+  for n in 1:NNN
+    tj, Vj = t / n, V / n
+    opsum += -tj, "Cdagup", 1, "Cup", 1 + n
+    opsum += -tj, "Cdagup", 1 + n, "Cup", 1
+    opsum += -tj, "Cdagdn", 1, "Cdn", 1 + n
+    opsum += -tj, "Cdagdn", 1 + n, "Cdn", 1
+    opsum += Vj, "Ntot", 1, "Ntot", 1 + n
+  end
+  return opsum
+end
+
+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
+
+function fermion_momentum_translator(i::Index, n::Integer; N=6)
+  #@show n
+  ts = tags(i)
+  translated_ts = ITensorInfiniteMPS.translatecelltags(ts, n)
+  new_i = replacetags(i, ts => translated_ts)
+  for j in 1:length(new_i.space)
+    ch = new_i.space[j][1][1].val
+    mom = new_i.space[j][1][2].val
+    new_i.space[j] = Pair(QN(("Nf", ch), ("NfMom", mom + n * N * ch)), new_i.space[j][2])
+  end
+  return new_i
+end
+
+function ITensors.op!(Op::ITensor, opname::OpName, ::SiteType"FermionK", s::Index...)
+  return ITensors.op!(Op, opname, SiteType("Fermion"), s...)
+end
+
+@testset verbose = true "InfiniteMPOMatrix -> InfiniteMPO" begin
+  ferro(n) = "↑"
+  antiferro(n) = isodd(n) ? "↑" : "↓"
+
+  models = [(Model"heisenbergNNN"(), "S=1/2"), (Model"hubbardNNN"(), "Electron")]
+  @testset "H=$model, Ncell=$Ncell, NNN=$NNN, Antiferro=$Af, qns=$qns" for (model, site) in
+                                                                           models,
+    qns in [false, true],
+    Ncell in 2:6,
+    NNN in 1:(Ncell - 1),
+    Af in [true, false]
+
+    if isodd(Ncell) && Af #skip test since Af state does fit inside odd cells.
+      continue
+    end
+    initstate(n) = Af ? antiferro(n) : ferro(n)
+    model_kwargs = (NNN=NNN,)
+    s = infsiteinds(site, Ncell; initstate, conserve_qns=qns)
+    ψ = InfMPS(s, initstate)
+    Hi = InfiniteMPO(model, s; model_kwargs...)
+    Hs = InfiniteSum{MPO}(model, s; model_kwargs...)
+    Es = expect(ψ, Hs)
+    Ei = expect(ψ, Hi)
+    #@show Es Ei
+    @test sum(Es[1:(Ncell - NNN)]) ≈ Ei atol = 1e-14
+  end
+
+  @testset "FQHE Hamitonian" begin
+    N = 6
+    model = Model"fqhe_2b_pot"()
+    model_params = (Vs=[1.0, 0.0, 1.0, 0.0, 0.1], Ly=3.0, prec=1e-8)
+    trf = fermion_momentum_translator
+    function initstate(n)
+      mod1(n, 3) == 1 && return 2
+      return 1
+    end
+    p = 1
+    q = 3
+    conserve_momentum = true
+    s = infsiteinds("FermionK", N; translator=trf, initstate, conserve_momentum, p, q)
+    ψ = InfMPS(s, initstate)
+    Hs = InfiniteSum{MPO}(model, s; model_params...)
+    Hi = InfiniteMPO(model, s, trf; model_params...)
+    Es = expect(ψ, Hs)
+    Ei = expect(ψ, Hi)
+    #@show Es Ei
+    @test Es[1] ≈ Ei atol = 1e-14
+  end
+end
+nothing