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

F0, F2, F4 = 5.0, 10.0, 6.0
zeta, Bz   = 0.1, 0.01

d = AngularMomentumOperators(NF, IndexUp, IndexDn)

H = F0*d.SlaterF[0] + F2*d.SlaterF[2] + F4*d.SlaterF[4] + zeta*d.ldots + Bz*(2*d.Sz + d.Lz)