Hello!

I'd like to ask a question about constructing the MPO operator being the product of local operators. In my problem I consider two types of bosons "a" and "b" with total number @@N = N_a + N_b@@ conserved, but @@N_a - N_b@@ is not conserved. Total number of sites is @@M = N@@.

I model my system by creating an artificial lattice with @@M = 2N@@ sites, and each two lattice sites form a unit cell where the first site corresponds to boson @@a@@ and second to boson @@b@@.

On such a lattice I can define spin operators @@S^x_i, S^y_i, S^z_i@@ as follows:

```
auto tmp_Sx = AutoMPO(sites);
auto tmp_Sy = AutoMPO(sites);
auto tmp_Sz = AutoMPO(sites);
for (int j = 1; j < M_sites; j+=2){
tmp_Sx += 0.5, "Adag",j,"A", j+1;
tmp_Sx += 0.5, "Adag",j+1,"A", j;
tmp_Sy += 0.5/Cplx_i, "Adag",j,"A", j+1;
tmp_Sy += -0.5/Cplx_i, "Adag",j+1,"A", j;
tmp_Sz += 0.5, "N", j;
tmp_Sz += -0.5, "N", j+1;
}
auto Sx = toMPO(tmp_Sx);
auto Sy = toMPO(tmp_Sy);
auto Sz = toMPO(tmp_Sz);
```

I'm interested in calculating expectation value of the @@D = \prod_j d_j@@

operator, where @@d_j = S^y_{j} + i S^z_{j}@@. I've tried the following construction (let's assume N = 2 particles):

```
vector<MPO> d_vec(N);
int j = 0;
for(int i = 1; i<M; i+=2){
auto tmp_d_i = AutoMPO(sites);
tmp_d_i += 0.5/Cplx_i, "Adag",i, "A", i+1;
tmp_d_i += -0.5/Cplx_i, "Adag",i+1,"A", i;
tmp_d_i += 0.5*Cplx_i, "N", i;
tmp_d_i += -0.5*Cplx_i, "N", i+1;
d_vec[j] = toMPO(tmp_d_i);
j = j+1;
}
auto D = nmultMPO(d_vec[1],prime(d_vec[0]), "Method=",nmultMPO_method,"MaxDim",nmultMPO_MaxDim,"Cutoff",nmultMPO_cutoff});
```

However, it doesn't work (I have exact-diagonalization results to compare).

Can someone point me to what I'm missing? Thanks in advance!

Best,

Marcin