The good news is that this sort of operator tends to be low bond dimension if it is all to all. See the discussion for S^2 here: Total spin operator for Electron sites system
I found the OpSum system quite slow so I built a similar MPO by hand (ala Jan’s post Optimizing dmrg results for expected values of observables - #6 by JanReimers), but you may find better results with e.g. GitHub - ITensor/ITensorMPOConstruction.jl: Alternative backends for constructing MPOs from sums of operators.