VUMPS MPO error

ITensor Julia Questions
1

Hi Everyone,

I am trying to construct an infinite MPO to then construct a VUMPS. I construct a model with custom operators with a bosonic local hilbert space. When I try doing the infinite sum for the MPO, I get the following error, “MethodError: no method matching translatecell(::ITensorInfiniteMPS.var”#_shift_cell#176"{Int64}, ::Scaled{ComplexF64, Prod{Op}}, ::Int64) ". I’ll include my code below. I suspect the problem is arising from my use of custom local operators since I try using a toy model on the same hilbert space that just uses the built in bosonic operators in itensor and I don’t run into any issues.

My model where phi = \frac{1}{\sqrt{2}} (a^\dagger + a) and cpi2 has matrix elements of the square of cpi =\frac{i}{\sqrt{2}} (a^\dagger - a)

function ITensorInfiniteMPS.unit_cell_terms(::Model"Realphi4")
    opsum = OpSum()
    opsum += 0.5, "cpi2",1
    opsum += 0.5*(1+m2), "phi",1,"phi",1
    opsum += 0.5, "phi",2,"phi",2
    opsum += -1, "phi",1,"phi",2
    opsum += ld/24 ,"phi",1,"phi",1,"phi",1,"phi",1
    return opsum
end

My attempt to create an infinite MPO that gives me the mentioned error. This should be runnable if you have the model, which I will include below if the error wants to be reproducible.

using ITensors
using ITensorInfiniteMPS

include(
  joinpath(
    pkgdir(ITensorInfiniteMPS), "examples", "vumps", "src", "vumps_subspace_expansion.jl"
  ),
)

m2= 10
ld =1


maxdim = 100 # Maximum bond dimension
cutoff = 1e-8 # Singular value cutoff when increasing the bond dimension
max_vumps_iters = 100 # Maximum number of iterations of the VUMPS algorithm at each bond dimension
vumps_tol = 1e-6
conserve_qns = false
outer_iters = 8 # Number of times to increase the bond dimension


N = 10 # Number of sites in the unit cell

initstate(n) = isodd(n) ? "1" : "2"
s = infsiteinds("Boson", N; dim=8, conserve_qns, initstate)
ψ = InfMPS(s, initstate)

model = Model("Realphi4")
#model = Model("bosontoy")
# Form the Hamiltonian
H = InfiniteSum{MPO}(model, s)

These are the necessary local operator definitions:

using ITensors

function ITensors.op(::OpName"phi",::SiteType"Boson", d::Int )
    mat = zeros(d,d)
    for i in 0:d-1
        for j in 0:d-1
            mat[i+1 ,j+1]  = (1/ sqrt(2))*( (i-1==j)*sqrt(j+1)+ (i==j-1)*sqrt(i+1) )
        end
    end
    mat[d,d]+= -d/2
    return mat
end

function ITensors.op(::OpName"cpi2",::SiteType"Boson", d::Int )
    mat = zeros(ComplexF64,d,d)
    for i in 0:d-1
        for j in 0:d-1
            mat[i+1 ,j+1] = (1/2)*( -(i-1==j+1)*sqrt((j+1)*(j+2))+ (i==j)*(2*i+1)- (i+1==j-1)*sqrt((i+1)*(i+2)))
        end
    end
    mat[d,d]+= -d/2
    return mat
end

Any help with this would be greatly appreciated. Thank you!

Best,
Jake Montgomery

2

Note that the uncommented bosontoy model is not shown here.
Can you confirm your hamiltonian is correct? Small modifications to it will run (adding a phi^4 term on site 2 for example)

Also the proper way to include m2 and ld is something like the following:

function ITensorInfiniteMPS.unit_cell_terms(::Model"Realphi4"; m2=10, ld=1)
    opsum = OpSum()
    opsum += 0.5, "cpi2",1
    opsum += 0.5*(1+m2), "phi",1,"phi",1
    opsum += 0.5, "phi",2,"phi",2
    opsum += -1, "phi",1,"phi",2
    opsum += ld/24 ,"phi",1,"phi",1,"phi",1,"phi",1
    #opsum += ld/24 ,"phi",2,"phi",2,"phi",2,"phi",2 # ?
    return opsum
end

m2= 10
ld =1
model_params = (m2 = m2, ld = ld)
H = InfiniteSum{MPO}(model, s; model_params...)
3

@ryanlevy Thank you for your reply. That’s my mistake, the “bosontoy” model line should have been commented out. I’ll edit it.

Yes that does seem to work. Do you have an understanding why that would work but not including that line has issues? To me it does not seem like it should matter whether I define a term on one site, or “split it” between two sites. Regardless this shouldn’t affect the model since I can just have the term on both sites 1 and 2 and add a factor of 1/2, but I’d like to understand why the infinite sum has a problem with that. Thank you again.

If you’re particularly curious this is the toy model :

function ITensorInfiniteMPS.unit_cell_terms(::Model"bosontoy")
    opsum = OpSum()
    opsum += 0.5, "n",1
    opsum += 0.5, "adag",1,"adag",1
    opsum += -1, "adag",1,"adag",2
    opsum += 0.5, "a",1,"a",1
    opsum += -1, "a",1,"a",2
    return opsum
end
4

It looks like this actually is an ITensorInfiniteMPS code bug. I have a work around but want to confirm the Hamiltonian you truly want (the bosontoy does what I expect but your model may not)
Essentially because you symmetry broke site 1 and 2 within the unit cell, the OpSum is not going to translate your on site terms for you. If we look at what bosontoy is doing for 4 sites

julia> ITensorInfiniteMPS.opsum_finite(Model("bosontoy"), 4; )

sum(
  0.5 n(1,)
  0.5 adag(1,) adag(1,)
  -1.0 adag(1,) adag(2,)
  0.5 a(1,) a(1,)
  -1.0 a(1,) a(2,)
  0.5 n(2,)
  0.5 adag(2,) adag(2,)
  -1.0 adag(2,) adag(3,)
  0.5 a(2,) a(2,)
  -1.0 a(2,) a(3,)
  0.5 n(3,)
  0.5 adag(3,) adag(3,)
  -1.0 adag(3,) adag(4,)
  0.5 a(3,) a(3,)
  -1.0 a(3,) a(4,)
  0.5 n(4,)
  0.5 adag(4,) adag(4,)
  0.5 a(4,) a(4,)
)

we see the n(1),n(2),n(3), and n(4) terms. But for Realphi4 since you’ve broken the symmetry inside the unit cell, and there are seperate operators on sites 1 and 2, the finite opsum looks like

julia> ITensorInfiniteMPS.opsum_finite(Model("Realphi4"), 4; model_params...)
sum(
  0.5 cpi2(1,)
  5.5 phi(1,) phi(1,)
  -1.0 phi(1,) phi(2,)
  0.041666666666666664 phi(1,) phi(1,) phi(1,) phi(1,)
  0.5 phi(2,) phi(2,)
  0.5 cpi2(3,)
  5.5 phi(3,) phi(3,)
  -1.0 phi(3,) phi(4,)
  0.041666666666666664 phi(3,) phi(3,) phi(3,) phi(3,)
  0.5 phi(4,) phi(4,)
)

Notice only cpi2(1) and cpi2(3) (missing sites 2 and 4).

Edit If this is what you want you’ll need to specify the "phi",2,"phi",3 term to make the unit cell ‘legal’. Otherwise, symmetrizing your unit cell will take care of the bug.

5

I’ve written up what’s going on here: Subtle error message for symmetry broken unitcell · Issue #82 · ITensor/ITensorInfiniteMPS.jl · GitHub

6

Ah yes I see the problem. No this is not what I am after. I’m not quite sure what you mean by “breaking the symmetry” of the unit cell. Do you mean treating site 1 and 2 differently? I figured for constructing an MPO, each site, iteratively, each site i will acquire operators on that are assigned to site 1, and the next site, i+1, will acquire operators placed on site 2. I must be misunderstanding this assignment procedure, this is just how I write the terms when I do normal DMRG.

These terms are arising from the square of a finite difference operator (\phi_{i+1} - \phi_i)^2 = \phi_{i+1}^2 +\phi_i^2 - 2\phi_{i+1} \phi_i , and when specifying the model, I am explicitly keeping the field square term on the next lattice site, and combining the i^{th} site square field with the mass term in my Hamiltonian. Again, in principle, I should be able to lump this i+1 field squared term in with the i^{th} site field square term since when summing the Hamiltonian, that is exactly what happens.

Thank you for looking into this.

7

Simply changing that one term actually does resolve the error :slight_smile: . Thank you for all your help.

Best,
Jake