@mingruyang, these new changes to the density matrix algorithm make the code look more complicated than before. What problems were you seeing before that this fixes? (I mean besides the write to disk, which seems like a good thing to add.)

The problem before this fix is that it uses huge memory when the MPO's bond dimension is large even if the bond dimension of the MPS is small, because there is a vector of N edge tensor which has 2 bond indices of the MPS and 2 bond indices of the MPO, which make it much bigger than in the DMRG case. When I have a MPS of bond dimension 100 and MPO of bond dimension 100, the density-matrix-apply function uses about 50GB. In addition to write-to-disk, I also change the code to only make a temporary conjugate of each site tensor of the MPS and MPO when constructing the edge tensor at each site, instead of constructing a conjugate of the whole MPS and MPO at a time which is only for one-off use to construct the edge tensor K but exist in the memory all the time until the function returns.