Hi!

I’m writing a snippet for extracting the coefficients of an MPS with siteset Electron by following this tutorial. Here it is my attempt:

```
double get_MPS_coefficient(MPS psi,
Electron sites,
std::vector<std::string> conf)
{
// args: psi MPS representing the state psi= \sum c_gamma |gamma>
// conf is a vector of length(psi) containing the indices of the MPS |gamma> to extract: ex: configurations = [Up, Dn, Up, Dn] if length(psi) = 4
// return: coefficients c_gamma
int N = length(psi);
auto coeff = ITensor(1.);
for(auto j : range1(N))
{
// sites(j, conf[j-1]) is an IndexVal representing the index of site(j)
// corresponding to the internal state given by conf[j-1] \in {Emp, Up, Dn, UpDn}
// setElt(IndexVal) initializes a zero tensor with an element equal to 1
// at the index position given by the ival
IndexVal ival = dag(sites(j, el[j-1]));
coeff *= psi(j) * setElt(ival);
}
double val_coeff = elt(coeff);
PrintData(val_coeff);
return val_coeff;
}
```

However, I am having some troubles in understanding what is the implicit order that is assumed for the fermionic operators, because some of the coefficients that I get by an exact diagonalization have flipped signs wrt the ones I get with code above and this is not an overall minus.

My guess is that something like the following happens:

Suppose I construct a state as |\psi\rangle = c^\dag_{1\uparrow} c^\dag_{2\downarrow} |0\rangle.

If I assume that my operators are ordered with increasing indices and \uparrow < \downarrow, then the coefficient in the expansions of |\psi\rangle will be +1, but if I assume a different ordering the coefficient might be -1.

Thanks a lot for any insights/help in advance!

Best,

Davide