Valence bond states - fast implementation

Question

Hello everyone!

These days, I am writing a valence bond code for the Heisenberg model, but my problem is that it is super slow...

On paper, we write the valence bond as equation. As you can see, it is a product of dimers in different sites.

The reason my code is slow is because in order to create that state, I am literally generating all the possible slater determinants, which are 2^(N/2), where N is the number of sites.

Computing the matrix elements is fast, creating the states is slow...

Is there a way to create the valence bond MPS like we do on paper, instead of generating all the possible determinants? My concern is that each dimer corresponds to different sites, even though they have the same number of sites.

So basically what I am after is if we can code a direct product of MPS's with different sites, but same number of sites...

I would really appreciate it if someone had an idea and point me to a direction that might solve the problem.

PS. I can show my code, too, although I don't think it is helpful...

Answer 1

Hi, so if I understand correctly, you want to create a state which is like the one in the linked formula (|u1 d2> - |d1 u2>) (x) (|u3 d4> - |d3 u4>) (x) (|u5 d6> - |d5 u6>) ...

So a good way to do this is to create a two-site wavefunction (tensor with two indices) and set its elements explicitly so it is the singlet state on these two sites. It will have two non-zero elements.

Then SVD this tensor (really a matrix) to obtain a two-site MPS.

Next, create a longer MPS on N sites. You can then copy this two-site MPS you made above into the MPS over and over again in a loop, copying the tensors into the sites spanning each odd-numbered bond (sites 1 & 2, 3 & 4, 5 & 6, etc.).

Finally, update all of the physical indices of the MPS so that they are distinct, rather than repeating with a 2-site periodicity. You can do this by using a delta tensor:

psi.ref(n) *= delta(s,sites(n));

where s is one of the indices of your original 2-site tensor and sites(n) is the n'th index obtained from a site set on N sites. (To see how to make a site set see some of the example codes included with ITensor.)

Alternatively you can do the two-site-dimer and SVD procedure from scratch each time on each bond, remaking it again in a loop.

For an actual example code that does this, you can see the following code:
https://github.com/emstoudenmire/finiteTMPS/blob/master/mpo_ancilla.cc
around line 96.
Please note this code is using ITensor version 2 notation, and we would recommend using the ITensor version 3 interface. Mainly this means you change "Aref" to just "ref".

Best regards,
Miles

Comments

Thank you very much for your reply.

I made it work perfectly, but I also would like to try dimers that are not next-nearest neighbors. What I tried is the following:

    int N = 8;
    vector <int> dimer;
    dimer.push_back(1);
    dimer.push_back(3);
    dimer.push_back(5);
    dimer.push_back(7);
    dimer.push_back(2);
    dimer.push_back(4);
    dimer.push_back(6);
    dimer.push_back(8);
    dimer.push_back(9);
    dimer.push_back(11);
    dimer.push_back(13);
    dimer.push_back(15);
    dimer.push_back(10);
    dimer.push_back(12);
    dimer.push_back(14);
    dimer.push_back(16);

    for(int n = 1; n <= 2*N; n += 2)
        {

        auto s1 = sites(dimer[n-1]);
        auto s2 = sites(dimer[n+1-1]);
        auto wf = TensorT(s1,s2);
        wf.set(s1(1),s2(2), ISqrt2);
        wf.set(s1(2),s2(1), -ISqrt2);
        TensorT D;
        psi.ref(n) = TensorT(s1);
        psi.ref(n+1) = TensorT(s2);
        svd(wf,psi.ref(n),D,psi.ref(n+1));
        psi.ref(n) *= D;
      }

But when I calculate the expectation value of the Hamiltonian I get the error:
In findIndex: more than one Index found

What I might be doing wrong?

---

Hi, so it's tricky to define an MPS such that the dimers are able to be non-nearest-neighbors.

It's not as simple as just changing the sites in the code; one has to preserve the MPS structure, which is that a tensor on site j has two indices: one connecting to the tensor on site j-1 and one connecting to site j+1 (except at the edges). Without this structure, we can't make any guarantees about how our algorithms will behave and you are likely to encounter errors.

There is a way to make an MPS that does connect further away sites (as long as each site is only connected to at most one other site). But the construction is more technical and also might lead to very high costs if the pattern is chosen to be arbitrary.

One route you may want to pursue instead of constructing the state directly would be to define operators which construct dimer states and act these on an MPS, starting with some trivial "reference" MPS such as all spins up.

Best,
Miles

---

Hi,

I see...

 Can you please maybe elaborate a bit more on your last point? Doesn't this approach still require an exponential number of MPOs or have I misunderstood something?

Thank you.

PS. The swapIndices() still causes a problem I guess, right?

---

Hi, so to elaborate, I wasn't necessarily suggesting MPOs. The operators could also just be local operators acting on two neighboring sites (4 index tensors). You could make an operator which maps an initial state such as two spins up to a dimer state.

Then to create a dimer which is not nearest-neighbor, you could use swap gates (crossed identity operators) to move one of the ends of the dimer further and further away from its original position.

In fact, a good strategy would be to just make a product state of dimers, then act with a pattern of swap gates which moves the ends of the dimers around until they are at the locations you want.

Miles

---

The T.swapInds(is1,is2) function should work ideally; if it doesn't please let us know here or file a bug report.

However, usually we suggest people swap indices using delta tensors. These are like "copy" tensors with all diagonal elements equal to 1 and otherwise zero. You can use one to replace an index like this:

T *= delta(i,j);

which would replace Index "i" with Index "j".

---

I tried to do something like that, but no luck yet... I am getting an error about mismatched IQIndex arrows.


My code:

for(int n = 1; n <= N; n += 2)
        {

        auto s1 = sites(n);
        auto s2 = sites(n+1);
        auto wf = TensorT(s1,s2);
        wf.set(s1(1),s2(2), ISqrt2);
        wf.set(s1(2),s2(1), -ISqrt2);
        TensorT D;
        psi.ref(n) = TensorT(s1);
        psi.ref(n+1) = TensorT(s2);
        svd(wf,psi.ref(n),D,psi.ref(n+1));
        psi.ref(n) *= D;
      }
      // Change index 1 with index 3
      TensorT D = delta(psi.ref(1).index(1),prime(psi.ref(3).index(1))); // What I thought you mean
      // The code below is from http://www.itensor.org/docs.cgi?vers=cppv2&page=formulas/gate
      psi.position(1);
      auto wf = psi.A(1)*psi.A(3);
      wf *= D;
      wf.noPrime();
      ITensor S,V;
      ITensor U = psi.A(1);
      svd(wf,U,S,V,{"Cutoff=",1E-8});
      psi.setA(1,U);
      psi.setA(1,S*V);

---

Hi, so you are on the right track. But one important comment is that, very generally, you should try to avoid referring to ITensor indices positionally, that is, by what order they appear on the tensor. In fact, the major design goal of ITensor is to permit writing code that does not depend on index ordering. We still offer that interface just for very special cases, and when writing developer-level, advanced code.

So for the specific case of getting the 'site' index of an MPS tensor (as you are doing, I believe, when making the delta tensor in your code), please try using the "siteIndex" function which is described on this page:
http://itensor.org/docs.cgi?vers=cppv3&page=classes/mps

What that function does in more detail is that it analyzes the connectivity of the MPS, and sees which index of a given MPS tensor is *not* shared with its neighboring tensors, and returns that one.

Best,
Miles

---

Hi again,

thank you so much, I finally managed to do it and it works like a charm.

---

Great - thanks for updating us

---

Hi Miles,

I am not very clear about the final point in your answer, which is :
"
Finally, update all of the physical indices of the MPS so that they are distinct, rather than repeating with a 2-site periodicity. You can do this by using a delta tensor:

psi.ref(n) *= delta(s,sites(n));
"
Neither, I don't see this line of code or something similar in https://github.com/emstoudenmire/finiteTMPS/blob/master/mpo_ancilla.cc .

I think the corresponding part in mpo_ancilla.cc should be:
auto psi = MPST(sites);
for(int n = 1; n <= 2*N; n += 2)
        {
        auto s1 = sites(n);
        auto s2 = sites(n+1);
        auto wf = TensorT(s1,s2);
        wf.set(s1(1),s2(2), ISqrt2);
        wf.set(s1(2),s2(1), -ISqrt2);
        TensorT D;
        psi.Aref(n) = TensorT(s1);
        psi.Aref(n+1) = TensorT(s2);
        svd(wf,psi.Aref(n),D,psi.Aref(n+1));
        psi.Aref(n) *= D;
        }

Can you please explain the reason why I should do this? I also get some confusions of the comments below the answers.

Besides, I  am also concerned about the normalization. After doing the first two steps, the ITensors are not normalized to 1.

Thanks!

Eric