Initialising an MPS from a reduced density matrix input

I’m trying to implement the following. Given a tensor with indices on sites q (vector of integers), I want to initialise it as part of a larger MPS where the sites to the left and right of sites q are alternating between empty and filled (no entanglement between the three sections). I’m using the code below, which works without quantum numbers but doesn’t with them. The reason is because the fluxes of the link dimensions are incorrect. The fluxes of the linkinds of sys_MPS only consider the QNs of
of sys_MPS, so don’t include the flux from the left or the right (obviously only one is needed depending on the direction of orthogonalization). Is there a way to manually override the fluxes, so if I put the orthogonalization centre at the first site and there are N particles to the left of the sites q, I can add N to the flux of the linkinds of sites q?
The way sys_Tensor is constructed it will have the same fluxes as alternating sites so the fluxes of the sites to the right of sites q should be correct.

I also realise my implementation may not be the right way to do this, I would welcome any suggestions. Also let me know if anything is unclear.

"""

This code initialises an MPS with 3 regions, the first 1:Nb sites

are alternating empty and filled (for spinless and spinful sites).

The sites q=Nb+1:Nb+Ns are initialised from a tensor over those sites.

The remaining sites are again alternating empty and filled.

Inputs:

-sys_Tensor

-s, the siteinds.

-q, the sites which are filled from the MPS created from sys_Tensor.

"""

##To initialise the tensor as a subset of an MPS, buffer sites

##are attached so they have link indices at the boundaries.

buff_q = (q[1]-1):(q[end]+1)

sys_ITensor = diagITensor(1,s[buff_q[1]])*sys_ITensor

sys_ITensor = sys_ITensor*diagITensor(1,s[buff_q[end]])

sys_MPS = MPS(sys_ITensor,s[buff_q])

##Initialise the MPS as alternating empty and filled sites

Empty = "Emp"

Full = "UpDn"

states = [isodd(n) ? Full : Empty for n in 1:N]

therm = MPS(ComplexF64, s, states)

sys_lk = linkinds(sys_MPS)

therm_lk = linkinds(therm)

##Change the sites of the sys_MPS so they have the correct link indices.

delta_L = delta(sys_lk[1],dag(therm_lk[2*N_L]))

sys_MPS[2] = sys_MPS[2]*delta_L

delta_R = delta(dag(sys_lk[end]),therm_lk[2*N_L+2*Ns])

sys_MPS[end-1] = sys_MPS[end-1]*delta_R

##Filling in therm with sys_MPS, excluding the buffer sites.

b = 1

for i in q

b += 1

therm[i] = sys_MPS[b]

end