Follow up on this question: https://itensor.org/support/2072/products-of-2-site-and-4-site-wavefunctions

The tetramers that I used to have are:

Now my purpose is instead of having only 6 states for each tetramer (the Sz=0 part) is having all 2^4 = 16 states (adding the 3 ups-1 down, 1 up-3 downs, 4 ups, 4 downs cases).

I tried:

```
for(int n = 1; n <= N; n += 4)
{
auto s1 = sites(n);
auto s2 = sites(n+1);
auto s3 = sites(n+2);
auto s4 = sites(n+3);
normalization = 0;
int pos = n/4*ms;
for (int i=pos; i<=pos+ms-1; i++)
normalization += pow(gsl_vector_get(cicoeff,i),2);
auto wf = ITensor(s1,s2,s3,s4);
wf.set(s1(1),s2(1),s3(2),s4(2), gsl_vector_get(cicoeff, pos)/sqrt(normalization));
wf.set(s1(1),s2(2),s3(1),s4(2), gsl_vector_get(cicoeff, pos+1)/sqrt(normalization));
wf.set(s1(1),s2(2),s3(2),s4(1), gsl_vector_get(cicoeff, pos+2)/sqrt(normalization));
wf.set(s1(2),s2(1),s3(1),s4(2), gsl_vector_get(cicoeff, pos+3)/sqrt(normalization));
wf.set(s1(2),s2(1),s3(2),s4(1), gsl_vector_get(cicoeff, pos+4)/sqrt(normalization));
wf.set(s1(2),s2(2),s3(1),s4(1), gsl_vector_get(cicoeff, pos+5)/sqrt(normalization));
wf.set(s1(2),s2(1),s3(1),s4(1), gsl_vector_get(cicoeff, pos+6)/sqrt(normalization));
wf.set(s1(1),s2(2),s3(1),s4(1), gsl_vector_get(cicoeff, pos+7)/sqrt(normalization));
wf.set(s1(1),s2(1),s3(2),s4(1), gsl_vector_get(cicoeff, pos+8)/sqrt(normalization));
wf.set(s1(1),s2(1),s3(1),s4(2), gsl_vector_get(cicoeff, pos+9)/sqrt(normalization));
wf.set(s1(2),s2(2),s3(2),s4(1), gsl_vector_get(cicoeff, pos+10)/sqrt(normalization));
wf.set(s1(2),s2(2),s3(1),s4(2), gsl_vector_get(cicoeff, pos+11)/sqrt(normalization));
wf.set(s1(2),s2(1),s3(2),s4(2), gsl_vector_get(cicoeff, pos+12)/sqrt(normalization));
wf.set(s1(1),s2(2),s3(2),s4(2), gsl_vector_get(cicoeff, pos+13)/sqrt(normalization));
wf.set(s1(2),s2(2),s3(2),s4(2), gsl_vector_get(cicoeff, pos+14)/sqrt(normalization));
wf.set(s1(1),s2(1),s3(1),s4(1), gsl_vector_get(cicoeff, pos+15)/sqrt(normalization));
psi.ref(n) = ITensor(s1);
psi.ref(n+1) = ITensor(s2);
auto [U,S,V] = svd(wf,{inds(psi.ref(n))},{"Cutoff=",1E-8});
psi.ref(n) = U;
psi.ref(n+1) = S*V;
psi.ref(n+2) = ITensor(s3);
ITensor U2(commonIndex(psi.ref(n), psi.ref(n+1)),inds(ITensor(s2))),S2,V2;
svd(S*V,U2,S2,V2);
psi.ref(n+1) = U2;
psi.ref(n+2) = S2*V2;
psi.ref(n+3) = ITensor(s4);
ITensor U3(commonIndex(psi.ref(n+1), psi.ref(n+2)),inds(ITensor(s3))),S3,V3;
svd(S2*V2,U3,S3,V3);
psi.ref(n+2) = U3;
psi.ref(n+3) = S3*V3;
}
```

but I am getting a "Setting Tensor element non-zero would violate its symmetry." error. The code above works perfectly if I consider only the Sz=0 part, so I am pretty sure that the problem is that I consider more Sz sectors now.

Is there a way to code that? Reading the Itensor's documentation, I found the mixedIQTensor() function, but I didn't quite understand if it does exactly what I want to do.

Any thoughts on that?

Thank you very much,

Thanos