To answer your question I thought it would be easier to show you some diagrams, which I uploaded to this link:
Can you please let me know if those notes answer your question and give a clear enough idea of what the code there is doing?
The purpose of this ProjMPS overall is to make a Hamiltonian term that is equal to |M><M| for some state M (represented as an MPS) then include this Hamiltonian term in an efficient way within a DMRG calculation that's optimizing an MPS psi.
I think the code there already partially works for n-site DMRG, because you can see that makeL! and makeR! (which are called by position!) query a parameter P.nsite which right now is always set to 2 but could be adjusted. Then the code in
product would need to be extended to multiply on the center "n" tensors of M in the lower part of the diagram in the notes, rather than just the two center tensors right now. These two tensors right now are getting included by the lines:
Lpm = dag(prime(P.M[P.lpos+1],"Link"))
Rpm = dag(prime(P.M[P.rpos-1],"Link"))