How to set the orthogonality center manually?

Is there a recommended way to set the orthogonality center (OC) manually? For example, after applying a two-site gate, I will need to set the orthogonality center manually. If I follow the example code from the documentation, the resulting MPS will not have an OC.

using ITensors, ITensorMPS

psi = orthogonalize(psi, 3)

wf = (psi[3] * psi[4]) * G
wf = noprime(wf)

inds3 = uniqueinds(psi[3], psi[4])
U, S, V = svd(wf, inds3; cutoff=1E-8)
psi[3] = U
psi[4] = S * V

I found in the source code, there is a set_ortho_lims function in abstractmps.jl. However, it says the function is not exported, and thus I want to know if there is a more recommended way to do it, or should I just use set_ortho_lims.

Thank you.

Setting orthogonality limits directly can be dangerous and should be avoided when possible.
In your specific case I think

psi[3:4] = wf

should work. You can read more on it here.
In more complicated cases you may want to re-orthogonalize the state using orthogonalize[!]

One thing is you can promise that you are aren’t changing the current center with the @preserve_ortho macro

@preserve_ortho psi begin
  psi[3] = U
  psi[4] = S * V # assumes orthogonalize(psi,4) has been called
end

However if you want to apply gates you can use apply to take care of all of this for you

Thank you for the reply and the useful informaiton. I didn’t clarify my purpose clearly enough. I want to apply the gates and analyze the singular values (for example, compute the entanglement entropies) for all the bonds. Therefore, it is natural for me to assign the MPS tensors manually because I already have them anyway, but I also need to set the orthogonality center correctly. I think there are different scenarios where setting the orthogonality center could be useful. Nevertheless, I realized that I can still use the set_ortho_lims function even though it is not exported, so my problem is solved. Please let me know if there are more suggestions.