NF = 16
IndexDn_2p = {0,2,4}
IndexUp_2p = {1,3,5}
IndexDn_3d = {6,8,10,12,14}
IndexUp_3d = {7,9,11,13,15}
Shell_2p = {0,1,2,3,4,5}
Shell_3d = {6,7,8,9,10,11,12,13,14,15}

H = NewOperator("U", NF, IndexUp_3d, IndexDn_3d, {5.0, 8.0, 6.0})
  + NewOperator("U", NF, IndexUp_2p, IndexDn_2p, IndexUp_3d, IndexDn_3d, {6.0, 7.0}, {5.0, 3.0})
  + NewOperator("CF", NF, IndexUp_3d, IndexDn_3d, PotentialExpandedOnClm("Oh", 2, {0.6, -0.4}))
  + 0.4 * NewOperator("ldots", NF, IndexUp_2p, IndexDn_2p)

Hav = ShellAverage(ShellAverage(H, Shell_2p), Shell_3d)