Hello Miles,

I have been trying to calculate correlations from the ground state MPS saved as mps.h5 file. Here is a minimal example for free fermions:

Step 1: Obtaining and saving the ground state MPS:

```
using ITensors, HDF5, JLD, DelimitedFiles
let
N = 4
numb_sweeps = 5
sweeps = Sweeps(numb_sweeps)
maxdim!(sweeps,10,20,100)
cutoff!(sweeps,1e-10)
sites = siteinds("Fermion",N)
states = ["Occ","Occ","Emp","Emp"]
psi0 = productMPS(sites,states)
ampo = AutoMPO()
for j = 1:N-1
ampo += "C",j,"Cdag",j+1
ampo += "C",j+1,"Cdag",j
end
ampo += "C",N,"Cdag",1
ampo += "C",1,"Cdag",N
H = MPO(ampo,sites)
energy,psi = dmrg(H,psi0,sweeps)
f = h5open("mps.h5","w")
write(f,"mps",psi)
close(f)
return
end
```

Step 2: use the save mps.h5 file to calculate correlation

```
using ITensors, HDF5, JLD, DelimitedFiles
let
N = 4
sites = siteinds("Fermion",N)
ampo = AutoMPO()
ampo += "C",1,"Cdag",2
corr = MPO(ampo,sites)
f = h5open("mps.h5","r")
psi=read(f,"mps",MPS)
close(f)
println("correlation = ",inner(psi,corr,psi))
return
end
```

Everything is pretty standard, but I got the error: "ERROR: LoadError: MPO A and MPS B must share site indices. On site 1, A has site indices IndexSet{2} (dim=2|id=347|"Fermion,Site,n=1")' (dim=2|id=347|"Fermion,Site,n=1") while B has site indices IndexSet{1} (dim=2|id=806|"Fermion,Site,n=1")."

It seems that when I define the correlation operator corr in the 2nd code, the site indices are regenerated and no longer match with that of the MPS obtained in the 1st code. If this is indeed the case, then is there a way around this issue?

Thank you for your time,

-Mason