Hello Miles,

After fixing the single-flavor fermion chain following your advice at http://itensor.org/support/3237/strange-inconsistency-in-free-fermion-benchmark, I went on to add in a second flavor of free fermion. I would expect the GS energy would simply double if the number of fermions is the same for the two flavors. However, this is not what I got. Is this because of the fermionic statistics between the two flavors? For single-flavor, I know that it has been taken care of by "has*fermion*string", but for two flavors I'm not sure whether it has been included properly or not. If not, I would like to know how to incorporate the statistics.

Here is a minimal version of the code for your reference. Thanks a lot. - Mason

```
using ITensors
function ITensors.space(
::SiteType"CFermion";
conserve_qns=false,
conserve_Q1=conserve_qns,
conserve_nfparity=conserve_qns,
qnname_Q1="Q1",
qnname_nfparity="NfParity",
)
if conserve_nfparity && conserve_Q1
return [
QN((qnname_nfparity, 0, -2), (qnname_Q1, 0)) => 1,
QN((qnname_nfparity, 1, -2), (qnname_Q1, +1)) => 1,
QN((qnname_nfparity, 1, -2), (qnname_Q1, +2)) => 1,
QN((qnname_nfparity, 0, -2), (qnname_Q1, +3)) => 1
]
end
return 4
end
ITensors.state(::SiteType"CFermion",::StateName"Emp" ) = 1
ITensors.state(::SiteType"CFermion",::StateName"f1" ) = 2
ITensors.state(::SiteType"CFermion",::StateName"f2" ) = 3
ITensors.state(::SiteType"CFermion",::StateName"f1f2" ) = 4
function ITensors.op!(Op::ITensor, ::OpName"C1", ::SiteType"CFermion", s::Index)
Op[s' => 1, s => 2] = 1.0
return Op[s' => 3, s => 4] = 1.0
end
function ITensors.op!(Op::ITensor, ::OpName"C1d", ::SiteType"CFermion", s::Index)
Op[s' => 2, s => 1] = 1.0
return Op[s' => 4, s => 3] = 1.0
end
function ITensors.op!(Op::ITensor, ::OpName"C2", ::SiteType"CFermion", s::Index)
Op[s' => 1, s => 3] = 1.0
return Op[s' => 2, s => 4] = 1.0
end
function ITensors.op!(Op::ITensor, ::OpName"C2d", ::SiteType"CFermion", s::Index)
Op[s' => 3, s => 1] = 1.0
return Op[s' => 4, s => 2] = 1.0
end
function ITensors.op!(Op::ITensor, ::OpName"F", ::SiteType"CFermion", s::Index)
Op[s'=>1,s=>1] = +1.0
Op[s'=>2,s=>2] = -1.0
Op[s'=>3,s=>3] = -1.0
Op[s'=>4,s=>4] = +1.0
end
ITensors.has_fermion_string(::OpName"C1", ::SiteType"CFermion") = true
ITensors.has_fermion_string(::OpName"C1d", ::SiteType"CFermion") = true
ITensors.has_fermion_string(::OpName"C2", ::SiteType"CFermion") = true
ITensors.has_fermion_string(::OpName"C2d", ::SiteType"CFermion") = true
let
N = 4
sweeps = Sweeps(10)
maxdim!(sweeps,10,20,50)
cutoff!(sweeps,1e-10)
sites = siteinds("CFermion",N; conserve_qns = true)
states = ["f1f2","Emp","Emp","Emp"]
ampo = AutoMPO()
for i = 1:(N-1)
ampo += -1, "C1d", i, "C1", i+1
ampo += -1, "C2d", i, "C2", i+1
ampo += -1, "C1d", i+1, "C1", i
ampo += -1, "C2d", i+1, "C2", i
end
# PBC
ampo += -1, "C1d", N, "C1", 1
ampo += -1, "C2d", N, "C2", 1
ampo += -1, "C1d", 1, "C1", N
ampo += -1, "C2d", 1, "C2", N
H = MPO(ampo,sites)
psi0 = productMPS(sites,states)
energy0,psi = dmrg(H,psi0,sweeps)
return
end
```