NF = 14
IndexDn = {0,2,4,6,8,10,12}
IndexUp = {1,3,5,7,9,11,13}

f = AngularMomentumOperators(NF, IndexUp, IndexDn)

-- the last two terms are tiny, they are there only to lift the accidental
-- degeneracy between terms that the Coulomb interaction puts at the same energy
H = 1.0*f.SlaterF[2] + 0.6*f.SlaterF[4] + 0.4*f.SlaterF[6]
  + 1E-7*f.Casimir.G2 + 1E-9*f.Casimir.SO7

psiList = Eigensystem(H, {NF, 0, {"11111111111111", 3, 3}}, 364, {{"DenseBorder", 400}})

-- keep one state per multiplet, so that each term shows up once
psiTerm = {}
Eprev = nil
for i = 1, #psiList do
  local E = psiList[i] * H * psiList[i]
  if (Eprev == nil) or (math.abs(E - Eprev) > 1E-8) then
    psiTerm[#psiTerm + 1] = psiList[i]
    Eprev = E
  end
end

PrintExpectationValues(psiTerm, {f.Ssqr, f.Lsqr, f.Casimir.SO3, f.Casimir.G2, f.Casimir.SO7}, H)