Thank you so much!
Here are my ideas for what you're asking. Let me know if I misunderstood one of your questions.
If by convert an MPS to ITensor form, you mean obtain a single ITensor which is the entire wavefunction that an MPS represents, then I don't see any other way to do it besides multiplying all of the MPS tensors psi.A(i) together. At least there won't be any other way to do it that will be significantly faster in a generic case than just multiplying the psi.A(i) together. Of course, you might want to multiply the psi.A(i) in a different order depending on the various bond dimensions of the MPS, but in general you will end up with a single tensor with d^N entries that will be expensive to compute and will take up a lot of memory, as you know.
There could be strategies to multiply the psi.A(i) together in groups, such as multiply neighboring A tensors, then multiply the resulting N/2 tensors with their neighbors, then the remaining N/4 tensors with their neighbors etc. and this could even be done in parallel. The idea of this multiplication pattern would be to put off the most costly contractions until the very end where you finally multiply a left-half tensor with N/2 site indices times the right-half tensor with N/2 site indices.
To check time reversal symmetry, I would think you'd act on each site with the time reversal operator. Then take the overlap of the resulting MPS trpsi with the original MPS psi by doing overlap(trpsi,psi) and see if the overlap is 1 or not.
To check if your code is working correctly, construct an MPS "by hand" which you know a priori to be time reversal invariant (a chain of spin 1/2s in perfect singlet states with their neighbors) and one which is not (a product state of spin 1/2s all pointing up) and apply the above test to both.