I’m trying to solve the following problem stemming from ED:

In exact diagonalization, if we have some basis (say of N spinless fermion) and wavefunction:

|\psi\rangle = \sum_{s_1,\cdots,s_N \in \{ 0, 1\}} c^{s_1,\cdots,s_N}|{s_1,\cdots,s_N}\rangle

I can perform the following:

Say if I’m interested in site a,b, I can sum up the exponentially many coefficients and obtain an effective state:

\phi = p_{00} |0_a0_b\rangle + p_{01} |0_a1_b\rangle + p_{10}|1_a0_b\rangle + p_{11}|1_a1_b\rangle

where

p_{00} = \sum_{s_i } c^{s_1, \cdots, s_{a-1}, 0, s_{a+1},\cdots, s_{b -1}, 0, s_{b + 1}, \cdots ,s_N } |s_1, \cdots, s_{a-1}, 0, s_{a+1},\cdots, s_{b -1}, 0, s_{b + 1}, \cdots ,s_N \rangle

and so on

Is there an equivalent way of doing this in ITensor?

I can `contract`

the MPS and obtain the coefficiently on each index combination, but that defeats the purpose of MPS

Since by definition the inner product:

\langle A|B\rangle =\sum_{s_1,\cdots, s_n} A_{s_N}\cdots A_{s_1}B_{s_1}\cdots B_{s_N}

I’m thinking to construct an all 1 MPS of bond dimension 1 (up to normalization), i.e.

```
s = siteinds("Fermion", N; conserve_qns =true)
vecs = ones( (N, 2))
psi = MPS(s, 1)
for j=1:N
psi[j] = ITensor(vecs[j,:],s[j])
end
orthogonalize!(psi,1)
normalize!(psi)
```

Then I can partially contract the two MPS (one of them being this dummy Tensor of all 1) except at sites a,b to obtain the equivalent two-site MPS.

My questions are:

- The above code does not work with QN conserving MPS and the code throws error about the block structure (flux)
- Is this the right approach to the above problem or am I oblivious to something fundamentally incorrect about this