I see! I attach my code that reproduces the second error. The line produces the error is H = MPO(ampo, sites). And I’ve thought my model more carefully and realized the symmetry I need is particle number and spin SU(2) symmtry (which I think is equivalent to conserve_qns=true?) for electron site, and **No** symmtry for spin1/2 site. I’m not very sure, if I write `sites = siteinds(n->isodd(n) ? "S=1/2" : "Electron", N; conserve_qns=true)`

, does conserve_qns=true only applys to electron or both?

A minimal Code:

```
using ITensors
using HDF5
let
N = 200
Npart = 100
t = 1.0
hz=5
hx=5
error=1E-8
sites = siteinds(n->isodd(n) ? "S=1/2" : "Electron", N; conserve_qns=true)
ampo = OpSum()
# hopping term
for i in 2:2:(N - 2)
ampo += -t, "Cdagup", i, "Cup", i + 2
ampo += -t, "Cdagup", i + 2, "Cup", i
ampo += -t, "Cdagdn", i, "Cdn", i + 2
ampo += -t, "Cdagdn", i + 2, "Cdn", i
end
# pseudospin -hz\sum_i Sz(i) at odd sites
for i in 1:2:N-2
ampo += -hz, "Sz", i
end
# coupling term -hx\sum_i Sx(i)N_up(i+1)
for i in 1:2:N-1
ampo += -hx, "Sx", i, "Nup", i+1
ampo += -hx, "Sx", i, "Ndn", i+1
end
println(ampo)
H = MPO(ampo, sites)
state = ["Emp" for n in 1:N]
p = Npart
for i in N:-2:2
if p > i/2
println("Doubly occupying site $i")
state[i] = "UpDn"
p -= 2
elseif p > 0
println("Singly occupying site $i")
state[i] = (isodd(i÷2) ? "Up" : "Dn")
p -= 1
end
end
for i in N-1:-2:1
state[i]=(isodd(i÷2) ? "Up" : "Dn")
end
# Initialize wavefunction to be bond dimension 10 random MPS with number of particles the same as `state`
psi0 = productMPS(sites, state)
# Start DMRG calculation:
sweeps = Sweeps(50)
setmaxdim!(sweeps,100)
setcutoff!(sweeps, error)
@show sweeps
energy, psi = dmrg(H, psi0, sweeps)
f = h5open("hx$(hx)hz$(hz)_1E-8.h5","w")
write(f,"PGS",psi)
close(f)
println("\nGround State Energy = $energy")
end
```