How to avoid index collisions in ITensors.jl?

ITensor Julia Questions
1

I have a fixed lattice (300-site spin-1/2 chain), but different Hamiltonians. I want to be able to do the following.

  • Store states computed using DMRG on different machines.
  • Load them later and, e.g., compute inner products between them.
  • Avoid spurious issues due to unintended index collisions (see example below).
  • Ideally, have reproducible computations (which return the same result if reran).
    What is the recommended way to handle the indices to achieve this?

I currently use a method that works more than 90% of the time, but sometimes fails spectacularly, as illustrated by the following example (requires a file with a saved state to reproduce):

# Reproducible 100% of the time due to seeding. Happens rarely otherwise.
import HDF5
using ITensors
import Random
HDF5.h5open("example.hdf5") do f
    global loaded_state = read(f, "elben_xxz_dmrg_res/states/2", MPS)
end;

# Ensure same site indices:
Random.seed!(index_id_rng(), 300)
sites = siteinds("S=1/2", 300)

sz_id = MPO(ITensor[op("Id", s) for s in sites])
z_prod = MPO(ITensor[2 * op("Sz", s) for s in sites])
zprod_proj1 = (sz_id - z_prod) / 2;

# Some unrelated computation (not used in the end result,
# but the issue does not reproduce if it is removed):
ops = OpSum()
for i in 1:(length(sites)-1)
    J = isodd(i) ? 1.0 : 0.0
    ops .+= J, "S+", i, "S-", i+1
    ops .+= J, "S-", i, "S+", i+1
end
ham = MPO(ops, sites)
_ = sz_id + z_prod;
ops = OpSum()
ops .+= 1.0, "Id", 1
obs1 = MPO(ops, sites);
_ = apply(obs1, obs1)
psi1 = productMPS(sites, [isodd(i) ? "Dn" : "Up" for i in 1:length(sites)])
sweeps = Sweeps(10)
noise!(sweeps, 1e-6, 1e-7, 1e-8, 1e-9, 0)
psi1 = dmrg(ham, psi1, sweeps)[2]
_ = inner(psi1', ham, psi1)
_ = inner(psi1, psi1)
_ = apply(zprod_proj1, psi1)
_ = dmrg(ham, psi1; nsweeps=1, cutoff=1e-12)[2]
_ = inner(psi1', obs1, psi1)
_ = apply(zprod_proj1, psi1)

# The issue:
psi2 = loaded_state
psi2_a = apply(zprod_proj1, psi2) # norm-1 MPS, correct
psi2_b = apply(zprod_proj1, psi2) # 0 MPS, incorrect, even though produced by the same formula
psi2_a_norm = norm(psi2_a)
psi2_b_norm = norm(psi2_b)
@assert (psi2_a_norm == psi2_b_norm) "$(psi2_a_norm) != $(psi2_b_norm)"
# AssertionError: 1.0000000000000469 != 0.0
2

I have to admit, that I don’t fully understand your example, but normally if you have states from different runs defined over the same Hilbertspace, you can use the replace_siteinds function to match the indices for the different states as in the example:

using ITensors

fix_hilbertspace(psi1, psi2) = replace_siteinds(psi1, siteinds(psi2))
N = 10
s1 = siteinds("SpinHalf", N)
s2 = siteinds("SpinHalf", N)

states = [isodd(j) ? "Up" : "Dn" for j in 1:N]
psi1 = MPS(s1, states)
psi2 = MPS(s2, states)


inner(psi1, fix_hilbertspace(psi2, psi1)) # -> returns 1