Hi Sylvain,
Thanks for your patience. I belive the iDMRG you are meaning is the code here: https://github.com/ITensor/iDMRG (just for reference of those reading).
At a high level, here is some key information about that code and what it returns. When the iDMRG finishes running, the MPS that is return describes just one unit cell, of course. A key thing to know is that the MPS tensors are in the "right orthogonal" gauge meaning each of the tensors psi(1), psi(2), ... will be right orthogonal.
There is a special tensor also returned as psi(0) which is the "center tensor" and it has no physical index, but lives on the link to the left of the first site.
So for example if you wanted to compute properties only within one unit cell of say size 2, then psi(0)*psi(1)*psi(2) would be the correct wavefunction of one unit cell. For two unit cells, conceptually psi(0)*psi(1)*psi(2)*psi(1)*psi(2) would be the wavefunction. That is, you only use the center tensor psi(0) once and keep "chaining" the other MPS tensors to make a basis for the infinite unit cells going off to the right.
From this background, you can see how the transfer matrix is schematically dag(psi(1))*psi(1)*dag(psi(2))*psi(2). I say schematically because there I am being very loose about the necessity to prime the link indices sticking out of the left and the right ends to make this contraction actually a matrix and not a scalar. Also you would never want to actually form the transfer matrix this way because it would cost the fourth power of the bond dimension of the MPS, which is a very high cost.
So finally to get the transfer matrix spectrum, which is what you want as you know for obtaining the correlation length, you can obtain it using the Arnoldi algorithm which uses iterative applications of the matrix to get only the dominant right-eigenvectors and their eigenvalues. You will need to code a custom C++ object that defines the method
void product(ITensor const& in, ITensor & out)
which, given an input ITensor "in", efficiently applies the transfer matrix to it one time and returns the resulting ITensor "out" so that it has the same collection of indices (including same prime levels) as the input ITensor. You can finally pass an instance of this custom object to the function arnoldi defined here:
https://github.com/ITensor/ITensor/blob/97f73cd7387b61420011ff36dbfd95ff87860990/itensor/iterativesolvers.h#L715
Given the right parameters, you can obtain the dominant eigenvalues of the transfer matrix efficiently that way.
Or you can use our new Julia ITensorInfiniteMPS package which includes a pre-written example code already of how to get the transfer matrix spectrum :^)
That example is located here:
https://github.com/ITensor/ITensorInfiniteMPS.jl/blob/main/examples/vumps/transfer_matrix_spectrum.jl
And we could help you get set up using that package (it has to be installed separately but only takes a second - very easy) if you are interested in that approach.
Best regards,
Miles