Using the packages developed for fermions, I wrote a very simple code to calculate the GS energy of free fermions.
N = 4
numbsweeps = 5
sweeps = Sweeps(numbsweeps)
sites = siteinds("Fermion",N; conserve_qns = true)
# initialize a state at half-filling
states = vcat(["Occ" for n = 1:N/2],["Emp" for n = 1:N/2])
psi0 = productMPS(sites,states)
ampo = AutoMPO()
for j = 1:N-1
ampo += "C",j,"Cdag",j+1
ampo += "Cdag",j,"C",j+1
H = MPO(ampo,sites)
energy,psi = dmrg(H,psi0,sweeps)
However, the GS energy calculated from inner(psi,*(H,psi)) is quite different from the result from DMRG. I also tried this for Ising spin chain, and it checks out there. For the fermion case, I don't know where went wrong. Any help is greatly appreciated.