using ITensors
let
N = 10
s = siteinds("S=1/2",N)
hterms = OpSum()
for j=1:2:N
hterms += "Sz",j,"Sz",j+1
hterms += 1/2,"S+",j,"S-",j+1
hterms += 1/2,"S-",j,"S+",j+1
end
H = MPO(hterms,s)
energy, psi = dmrg(H,randomMPS(s); nsweeps=2, maxdim=2)
@show inner(psi',H,psi)
@show (3/4)*(N÷2)
return
end

Another way would be to explicitly make each singlet state as a two-site tensor (for each bond), then use singular value decompositions (SVD) to split the tensors into an MPS.