I never thought about it. Do you mean that under certain
circumstances, the quadrupolar term has to be accounted for in the
computation of the equilibrium density matrix (so that the matrix is
longer proportional to Iz)?
At this point, there is no way to automatically generate such a
matrix. But you still can do this manually. You just need to figure
out the correct expression for your initial density matrix that
includes the quadrupolar term, and specify it in the input file as
usual (at the rho0 line). You would have to use the bra-ket notation
for this though, because I did not implement any matrix algebra in the
interpreter of the rho0 expression (so, for example, I1z^2 will not be
understood by the program).