Right, so the reason you are getting complex energies in that case is because apply(gates,rho) is defined to only apply the gates to one side of rho. So it is doing e^{-i H \tau} \rho.
In contrast, the correct way to evolve a density matrix in real time is e^{-i H \tau} \rho e^{i H \tau} for a time step \tau. For this purpose, there is a keyword argument called apply_dag you can use like this:
rho = apply(gates,rho; apply_dag=true)
The documentation for apply_dag can currently be hard to find, so I can see why you did not know about it!
For the imaginary time case of what you are doing you don’t need apply_dag however, because computing just e^{- H \tau} \rho ought to be fine for obtaining e^{- H \beta} when starting from \rho being an identity operator.
As to why you are still getting a different energy from your ED code, I don’t know the reason as it seems it is a bug you need to find, but it could be one of many different things:
- is it possible your ED code is using Pauli matrices for the spin operators whereas ITensor uses spin operators which include a factor of 1/2 ?
- perhaps you are not accounting correctly for the total amount of imaginary time evolved in one or both codes?
- I note in your code that you are printing out the energy before taking the imaginary time step. Maybe you should print it out after taking the time step?