Problems in benchmarking DMRG + measurements - XX chain

Hey,

I am trying to benchmark my exact free fermion numerics for a setup (that I shall now explain) on DMRG. My setup involves taking the ground state of the XX chain and making random 0/1 projections/measurements on a subset of sites of the system and then probing the resultant leftover state. I pick the 0/1 projection for a given site with 50-50 probability. Now, in the exact free fermion numerics, I observe that my state is sometimes killed (which makes sense since there would be outcomes chosen such that the overlap with the ground state is 0). However, that does not seem to happen in the DMRG code. Here is my measurement code:

function ITensors.op(::OpName"projupz", ::SiteType"S=1/2", s::Index)
    mat = [1 0
           0 0]
    return itensor(mat,  s', s)

end

function ITensors.op(::OpName"projdnz", ::SiteType"S=1/2", s::Index)
    mat = [0 0
           0 1]
    return itensor(mat,  s', s)
end

outcomes = rand(0:1, length(meas_loc))

function measure_avg(psi, loc, outcomes, sites)
    for (i, site_index) in enumerate(loc)
        op_tuple = (outcomes[i] == 0 ? ("projdnz", site_index) : ("projupz", site_index))
        psi = apply(ITensors.ops([op_tuple], sites), psi)
        psi /= norm(psi)
    end
    return psi
end

  sites = siteinds("S=1/2",L)
  nsweeps = 40
  maxdims = [200,300, 350, 400]
  cutoff = [1E-15]

  os = OpSum()
  
  for j=1:L-1
      os -= "Sy",j,"Sy",j+1
      os -= "Sx",j,"Sx",j+1
  end

  os -= "Sy",L,"Sy",1
  os -= "Sx",L,"Sx",1
  H = MPO(os,sites)

psi0 = randomMPS(sites;linkdims=4)
observe = DMRGObserver(; energy_tol=1e-12, minsweeps=10, energy_type=Float64)
energy,psi_gs = dmrg(H,psi0;nsweeps,maxdim,cutoff, observer=observe)
result = measure_avg(psi_gs, meas_loc, outcomes, sites)

The above is just a skeleton of my code. I tried debugging by taking the input state to the measurement function to be the product states of S_x and S_z eigenstates and I get the expected results. Not sure where the issue could be? Any help would be appreciated.

You’ve just been lucky, but it actually happens. I just tried to force it using the bell state \varphi^+=(\ket{00}+\ket{11})/\sqrt2 and projecting on opposite spins. Here the modified code:

using ITensors, ITensorMPS

function measure_avg(psi, loc, outcomes, sites)
   for (i, site_index) in enumerate(loc)
   	op_tuple = (outcomes[i] == 0 ? ("ProjUp", site_index) : ("ProjDn", site_index))
   	psi = apply(ITensors.op(op_tuple, sites), psi)
   	normψ = norm(psi)
   	if isapprox(normψ, 0.0; atol=1E-12)
   		throw("The norm of the state is zero")
   	end
   	psi /= normψ
   	truncate!(psi;cutoff=1E-12)
   end
   return psi
end

# Construct the bell state |ϕ⁺⟩ = 1/√2(|00⟩ + |11⟩)
sites = siteinds("S=1/2",2)
ϕ⁺ = MPS([1,0,0,1]/sqrt(2.), sites)

result = measure_avg(ϕ⁺, [1,2], [0,1], sites)	# Force to project on opposite spins

I also made some changes like using the buil-ins “ProjUp” and “ProjDn” and calling directly ITesnors.op since you are applying only 1 operator. I also added a truncation of the state, since projections are going to reduce the entanglement.