The lower bond-dim. of AKLT model

ITensor Julia Questions
4

I’m sorry I didn’t notice before, but are you sure the hamiltonian is correct? The term should be (\vec S_i \vec S_{i+1})^2 = (S^x_iS^x_{i+1}+S^y_iS^y_{i+1}+S^z_iS^z_{i+1})^2 which I don’t think it is equivalent to the one you wrote, the general writing should be this, I’m not sure if the additional terms cancel out.

Also, consider using S^+ and S^- instead of S^x and S^y, in this way the hamiltonian would be real and the computation should be faster.

Furthermore, the S=1 spin with the spin operators is already implemented in ITensors SiteTypes Included with ITensor · ITensorMPS.jl

Edit

I confirm that the error is the implementation of the Hamiltonian, This is the code I have used:

Code Used 
using ITensors, ITensorMPS

os = OpSum()
# OBC
for j = 1:N-1
    global os
    os += 1.0, "Sx", j, "Sx", j + 1
    os -= 1.0, "iSy", j, "iSy", j + 1
    os += 1.0, "Sz", j, "Sz", j + 1
    os += 1 / 3, "Sx", j, "Sx", j + 1, "Sx", j, "Sx", j + 1
    os += 1 / 3, "iSy", j, "iSy", j + 1, "iSy", j, "iSy", j + 1
    os += 1 / 3, "Sz", j, "Sz", j + 1, "Sz", j, "Sz", j + 1
    os -= 1 / 3, "Sx", j, "Sx", j + 1, "iSy", j, "iSy", j + 1
    os += 1 / 3, "Sx", j, "Sx", j + 1, "Sz", j, "Sz", j + 1
    os -= 1 / 3, "iSy", j, "iSy", j + 1, "Sx", j, "Sx", j + 1
    os -= 1 / 3, "iSy", j, "iSy", j + 1, "Sz", j, "Sz", j + 1
    os += 1 / 3, "Sz", j, "Sz", j + 1, "Sx", j, "Sx", j + 1
    os -= 1 / 3, "Sz", j, "Sz", j + 1, "iSy", j, "iSy", j + 1
end

# PBC
# os += 1.0, "Sx", 1, "Sx", N
# os -= 1.0, "iSy", 1, "iSy", N
# os += 1.0, "Sz", 1, "Sz", N
# os += 1 / 3, "Sx", 1, "Sx", N, "Sx", 1, "Sx", N
# os += 1 / 3, "iSy", 1, "iSy", N, "iSy", 1, "iSy", N
# os += 1 / 3, "Sz", 1, "Sz", N, "Sz", 1, "Sz", N
# os -= 1 / 3, "Sx", 1, "Sx", N, "iSy", 1, "iSy", N
# os += 1 / 3, "Sx", 1, "Sx", N, "Sz", 1, "Sz", N
# os -= 1 / 3, "iSy", 1, "iSy", N, "Sx", 1, "Sx", N
# os -= 1 / 3, "iSy", 1, "iSy", N, "Sz", 1, "Sz", N
# os += 1 / 3, "Sz", 1, "Sz", N, "Sx", 1, "Sx", N
# os -= 1 / 3, "Sz", 1, "Sz", N, "iSy", 1, "iSy", N

sites = siteinds("S=1", N)

H = MPO(os, sites)
psi0 = random_mps(sites, linkdims=2)
energy, psi = dmrg(H, psi0; nsweeps=10, maxdim=[10, 10, 20, 40, 80, 100]);

And these are the results obtained with N=128:

julia> include("exp.jl");
After sweep 1 energy=-84.0199504996155  maxlinkdim=10 maxerr=1.01E-04 time=0.383
After sweep 2 energy=-84.65892269585657  maxlinkdim=10 maxerr=5.88E-06 time=0.326
After sweep 3 energy=-84.66601973022019  maxlinkdim=20 maxerr=1.21E-10 time=0.745
After sweep 4 energy=-84.66650798857829  maxlinkdim=40 maxerr=3.71E-13 time=2.431
After sweep 5 energy=-84.66661612521872  maxlinkdim=80 maxerr=1.66E-14 time=6.331
After sweep 6 energy=-84.66665106329236  maxlinkdim=100 maxerr=3.63E-15 time=11.978
After sweep 7 energy=-84.6666626385208  maxlinkdim=100 maxerr=1.23E-15 time=14.806
After sweep 8 energy=-84.66666560968093  maxlinkdim=100 maxerr=1.66E-15 time=13.052
After sweep 9 energy=-84.66666639004092  maxlinkdim=100 maxerr=1.00E-15 time=9.586
After sweep 10 energy=-84.6666666048294  maxlinkdim=69 maxerr=2.22E-16 time=5.896

julia> energy, psi = dmrg(H, psi; nsweeps = 10, maxdim=100);
After sweep 1 energy=-84.66666665655427  maxlinkdim=51 maxerr=2.22E-16 time=2.880
After sweep 2 energy=-84.66666666534802  maxlinkdim=37 maxerr=2.21E-16 time=1.291
After sweep 3 energy=-84.66666666653006  maxlinkdim=25 maxerr=2.22E-16 time=0.592
After sweep 4 energy=-84.66666666665955  maxlinkdim=14 maxerr=2.20E-16 time=0.495
After sweep 5 energy=-84.66666666666681  maxlinkdim=6 maxerr=2.15E-16 time=0.200
After sweep 6 energy=-84.66666666666667  maxlinkdim=2 maxerr=1.50E-17 time=0.203
After sweep 7 energy=-84.66666666666609  maxlinkdim=2 maxerr=3.48E-22 time=0.149
After sweep 8 energy=-84.66666666666649  maxlinkdim=2 maxerr=1.27E-26 time=0.158
After sweep 9 energy=-84.66666666666688  maxlinkdim=2 maxerr=8.51E-29 time=0.154
After sweep 10 energy=-84.6666666666664  maxlinkdim=2 maxerr=9.59E-29 time=0.135

As you can see in less than 20 sweeps it has reached the convergence with a maxlinkdim of 2 (which corresponds to the VBS state)

For the future, I suggest you to use exact diagonalization for small systems and compare the results and the hamiltonian, so that you can spot these mistakes previously.