Density Matrix for Nuclear Spin Evolution

[delta]

I am trying to simulate the evolution of certain nuclear spins in a molecular system at a given temperature. Specifically, I want to obtain the density matrix of these nuclear spins evolving over time. I have attempted to use the evolution function with the 'output' parameter set to 'trajectory', and I obtained state vectors evolving over time. However, I am unsure about the relationship between these state vectors and the density matrix.

For example, when considering the density matrix for a system of 2 spins, a single state vector might be [0,0,0.5,0,0,0,0,0,0.5,0,0,0,0,0,0,0]. If I directly reshape it into a 4x4 matrix, the trace of the density matrix will be all zeros, which does not equal 1.

Is there a correct way to convert state vectors to density matrices? Or is there a better method to handle the problem of obtaining the density matrix of certain nuclear spins evolving over time？

[kuprov]

Set bas.formalism to 'zeeman-hilb', that would switch the kernel to the usual density operator formalism ('zeeman-liouv' and 'sphten-liouv' are more sophisticated options for dissipative dynamics and large spin systems).

As the initial condition, get the thermodynamic equilibrium from equilibrium.m function. Then you'll get the density matrix that has a unit trace; that trace will be preserved by unitary evolution.