I am seeking an algorithm to generate a random wavefunction = $\sum {c_i |\varphi _i\rangle }$ from a thermal ensemble, whose density matrix $\rho \sim e^{-\beta H}$, without the need to diagonalize the Hamiltonian. One possible method is to generate a random number for each basis and compare it with $e^{-\beta \langle\varphi _i|H|\varphi _i \rangle }$ and if it is lower we assign the corresponding $c_i$ to complex number with norm = 1 and totally random phase and otherwise make it just zero. Another possible algorithm to assign all $c_i$ to random complex numbers of norm one and then evaluate $e^{- \frac{\beta}{2} H}|\psi \rangle$. Are there other algorithms?
Edit: I now tried the second method with evaluating $e^{- \frac{\beta}{2} H}|\psi \rangle$ by propagating the wavefunction in imaginary time and it works fine! So I recommend this one.
