There is a neat trick for this, which is used internally in SpinDynamica. If the density matrix is represented in the Zeeman ket-bra basis, it can be flattened to a vector corresponding to the coefficients of the Zeeman ket-bra operators. One can then take the dot product with the appropriately-ordered list of basis operators.
To solve your problem in SpinDynamica you can use the ExpressOperator routine, which works more-or-less in the way that Ilya (kuprov) has described:
For more information you can look up the help message by executing