Partial summation of MPS coefficients

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

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])

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:

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

If you’d like to compute a sum such as:

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 your original state is an MPS, you can do that by contracting the MPS tensors, except those for sites a and b, with vectors of the form [1, 1]. Then you can contract all the ITensors together to get the resulting state which will just now have two site indices.

So you could do this with code such as:

# Say your MPS is called M
s = siteinds(M)
result = ITensor(1.)
for j=1:N
  if j == a || j == b
    result *= M[j]
    result *= (ITensor([1,1],s[j]) * M[j])
@show result
1 Like

Hi miles,

Thank you for the reply. To be clear, the M have QN conserving fermionic sites. When I tried your code, I received the following error:

Block Block(2,) has flux QN("Nf",1,-1) that is inconsistent with the desired flux QN("Nf",0,-1)

which I believed occurred during


If I do

ITensor([0, 1], s[j])

I get “QN indices must have opposite direction to contract”, and if I naively do

ITensor([0, 1], s[j])'

Then it keeps adding more indices to the resulting tensor. (I’m guessing the array passed to ITensors() here indicates the flux on the sites? I am not sure from reading the documentations)

Any ideas how to work properly with QN conserving dense blocks?

Many thanks in advance.

I’ve revisited this problem a bit,

Turns out my previous reply was not really making sense, since the site dimensions are actually the physical dimension, so in order to make an ‘identity’ I will always do for Fermions

s = siteinds(M)
ITensor([1, 1], dag(s[j])) * M[j]

Now this does not solve the incompatibility with QN conserving MPS as I hoped (I thought the dag would flip the flux on the physical indices so the contraction can proceed), resulting in error

Block Block(2,) has flux QN("Nf",-1,-1) that is inconsistent with the desired flux QN("Nf",0,-1)

There’s still something I don’t understand about the setup of QN conserving ITensor, so if you could take a look that would be great. Thanks!