Meaning of Nup parameter for spin_basis_1d initializer unclear when S >= 1

[tomohiro-soejima]

Hello,

I've been playing around with S = 1 AKLT model with your package. I tried ED for all symmetry sectors worked great, and I was able to find the four degenerate ground states at once. However, once I tried using symmetries, I encountered some puzzling situations I could not resolve.

1. Nup is ambiguous According to the documentation, it looked like the way to specify a symmetry sector based on Sz magnetization is something like

   L = 12
   basis = spin_basis_1d(L, pauli = False, S = "1", Nup = 6)

   However, when S = 1 or larger, specifying the number of up spins does not specify magnetization of the system, and I could not figure out what the convention is from the documentation.

2. Meaning of Nup seems to be internally inconsistent. Next, I tried using pzblock = 1 to speed up my calculations as follows:

   L = 12
   basis = spin_basis_1d(L, pauli = False, S = "1", Nup = 6, pzblock = 1)

   However, this resulted in the following error

   
   ---------------------------------------------------------------------------
   ValueError                                Traceback (most recent call last)
   <ipython-input-210-98143b7a24a1> in <module>
   ----> 1 basis = spin_basis_1d(L, pauli = False, S = "1", Nup = 6, pzblock=1)

~/miniconda3/envs/py3/lib/python3.7/site-packages/quspin/basis/basis_1d/spin.py in init(self, L, Nup, m, S, pauli, *blocks) 193 raise ValueError("spin inversion/particle-hole symmetry with particle/magnetization conservation must be used with chains with 0 magnetization sector or at half filling") 194 if Nup != L(self.sps-1)//2: --> 195 raise ValueError("spin inversion/particle-hole symmetry only reduces the 0 magnetization or half filled particle sector") 196 197 if (Nup_list is not None) and ((type(zAblock) is int) or (type(zBblock) is int)):

ValueError: spin inversion/particle-hole symmetry only reduces the 0 magnetization or half filled particle sector

After playing around with it, I realized `basis.sps == 3` for ` S = "1"`, so I need `Nup == L` (!!) for this system. Accordingly, I tried the following:
```python
L = 12
basis = spin_basis_1d(L, pauli = False, S = "1", Nup = 12, pzblock = 1)

This indeed ran without any error. However, this does not make physical sense to me. If I have `Nup == L``, my magnetization is clearly not 0, and spin inversion takes me out of this magnetization sector.

I would love to be able to use both Sz conservation and spin inversion symmetry. Do you guys mind looking into it?

[mgbukov]

Indeed, using Nup for higher spin is not optimal. There is however, the basis optional argument m which sets the magnetization density, e.g.

spin_basis_1d(L,S="1",m=0,kblock=0,pblock=1,zblock=1) # 

Does this help?

[tomohiro-soejima]

okay, that line worked. Thanks! I will report back if I can observe expected physics.

If I could offer a humble, novice user's opinion, I would deprecate or introduce a warning for Nup for large spins. The fact Nup can be larger than the system size (e.g. basis = spin_basis_1d(12, pauli = False, S = "1", Nup = 14) works ) is a pretty confusing behavior.

[tomohiro-soejima]

Using the argument m worked! I was able to look at energy vs. entanglement entropy collapse in AKLT. Here is a plot for your entertainment. I'm very glad it works beautifully for spin-1 systems.

Thanks for your help!

[mgbukov]

re-opening this, so we don't forget to consider deprecating N_up for higher spins before te next release.

[weinbe58]

We could also add in a total Sz projection argument instead of Nup. I think a useful feature of this argument is that it can take both float and string arguments to avoid rounding errors.

[mgbukov]

consider adding input S_tot and Sz

[HamidArianZad]

How it is possible to use conserved S_tot and Sz for a S=1 chain? Indeed I want to diagonalize the Hamiltonian of a spin-1 chain, when Sz total is conserved. I tried to use the coservation of Nup and Ndn with following code:

sz_min = -L
sz_max = L
for k in range(sz_min, sz_max + 1):
       basis_up = spin_basis_1d(L, pauli=False, S="1", Nup=sz_max + k)
       basis_dn = spin_basis_1d(L, pauli=False, S="1", Ndn=sz_min + k)

unfortunately, I got an error that says ValueError: unexpected optional argument(s): Ndn . Hence, I could not get any result and the code failed.