Dear all,

First, the background is that I am currently trying performing some DMRG calculation for the 2-Dimensional Spin-1 BBH model "on a cylinder" (periodic in y direction) with only nearest neighbor interaction, which goes as follows.

$$

\hat{H} = \sum_{i, j} \cos\theta \vec{S}_i\cdot \vec{S}_j + \sin\theta (\vec{S}_i\cdot \vec{S}_j)^2

$$

Now, I want to start with different types of initial states, such as Neel-like state, dimer(spin-singlet) state, or other initial states of interests. There are already a number of great discussions about how to set a product initial state in 1D, as shown in the following links,

http://itensor.org/support/1357/setting-random-initial-product-state

http://itensor.org/support/1359/how-to-create-a-particular-initial-state-as-a-mps

Also, as previously pointed by Miles, there is a very good example in the "finiteT/ancilla.cc". https://github.com/ITensor/ITensor/blob/v3/tutorial/finiteT/ancilla.cc

But after reading these answers, now I want to construct a dimer state as follow.

$$

\left| Initial\right>=\otimes_{1\leq n\leq N_x, 1\leq m\leq N_y} (2n-1, 2n; m)

$$

here my notation is that I use @@(2n-1, 2n; m)@@ to indicate a spin-singlet (dimer) state between site (x=2n-1, y=m) and site (x=2n, y=m). A spin-singlet state is like @@1/\sqrt{3} (\left|+1\right>_i\otimes\left|-1\right>_j+\left|-1\right>_i\otimes\left|+1\right>_j-\left|0\right>_i\otimes\left|0\right>_j)@@, if you wish.

And a naive thought of mine is just to do something similar in 1D case that I product all the pairs of site and make a product state. Here're my codes in c++ v3.

```
auto psi0 = MPS(sites);
for (int n = 1; n < N; n += 2 * Ny)
{
for (int m = 0; m < Ny; m++)
{
auto s1 = sites(n + m);
auto s2 = sites(n + m + Ny);
auto wf = ITensor(s1, s2);
wf.set(s1(1), s2(3), ISqrt3); // 1- Up; 2- Z0; 3 - Dn; ISqrt=1/\sqrt{3}
wf.set(s1(3), s2(1), ISqrt3);
wf.set(s1(2), s2(2), -ISqrt3);
ITensor D;
psi0.ref(n + m) = ITensor(s1);
psi0.ref(n + m + Ny) = ITensor(s2);
svd(wf, psi0.ref(n + m), D, psi0.ref(n + m + Ny));
psi0.ref(n + m) *= D;
}
}
```

But it returns some error as follow.

```
davidson: size of initial vector should match linear matrix size
```

I guess that since in ITensor, a 2D model is simply mapped into a snake-like 1D chain, I should have not made such naive generalization from 1D to 2D. So, I hope someone can help me figure out why this does not work from the perspective of building a product state in ITensor, as well as how to realize the setup of a product state in 2D. Btw, I think this question can also be generalized to how to set an general initial product state in the case of 2D models.

Thank you very much!

Best,

Yi