Should QDense constructor error if zero blocks?

[emstoudenmire]

The QDense storage constructor taking an IndexSet and a QN sometimes results in a storage with zero blocks, in the case that the requested QN is not compatible with any settings of the indices. This is ok, but the current behavior is that the constructor just makes no blocks (empty list of blocks) and goes on. As a result, calling flux on the resulting ITensor gives QN() which is interpreted as a zero QN. Then calling something like T.generate inside of a function like randomITensor now will generate blocks corresponding to a zero-flux ITensor. So the final QN ends up different from the one requested.

Where should we throw an error to prevent this? Is it reasonable for the QDense constructor to error? It's probably an easy change, but we just need to be sure of the right design.

Example to reproduce:

auto sites = SpinHalf(4,{"ConserveQNs=",true});
auto T = randomITensor(QN({"Sz",1}),sites(1),sites(2),sites(3),sites(4));
Print(flux(T) != QN({"Sz",1}));
Print(flux(T) == QN({"Sz",0}));

[mtfishman]

Good question. The way this is handled in the Julia version right now is:

using ITensors
s = siteinds("S=1/2", 4)
T = randomITensor(QN("Sz", 1), s...)
@assert store(T) isa NDTensors.BlockSparse
@assert iszero(nnzblocks(T))
@assert isnothing(flux(T))

so it handles it in a similar way (makes an ITensor with no blocks), however for that case flux returns nothing.

Perhaps trying to make that ITensor should error (or at least give a warning), since likely it is a mistake to try to make it in the first place. At least in Julia, the "correct" way to make an ITensor with no blocks is:

T = emptyITensor(s...)
@assert store(T) isa NDTensors.Empty

but I forget if there is an equivalent thing in the C++ version.

On a related note, to explicitly make an ITensor with zero blocks but with an actual BlockSparse storage (instead of Empty storage), I was considering a notation like:

T = ITensor(Block[], s...)

which would generalize to being able to specify a list of flux-consistent nonzero blocks, i.e.:

T = ITensor([Block(2, 2, 1, 1), Block(1, 1, 2, 2)], s...)

[mtfishman]

So I would lean towards it being an error to try to make the tensor randomITensor(QN("Sz", 1), s...) in the first place, since probably it is a mistake to do that and if you want an ITensor with empty blocks, you should use a different constructor.

[emstoudenmire]

I agree. There is the small issue that randomITensor works by calling a specific, fairly general-purpose QDense constructor so it could interfere with other functions that call that constructor and rely on it to just silently create a QDense with no blocks. But I think we'll just have to run the tests we have and if those pass then assume it's not breaking any other code and/or wait for a user to raise an issue if it errors in theirs.

[mtfishman]

Yeah. I think in terms of it being a breaking change, I would say that right now it is ill-defined behavior so it is ok to break it. We can also have the error message point users to the "correct" way to make a QN ITensor with no blocks (so in Julia, emptyITensor, and in C++, I guess the plain ITensor constructor with no QNs specified?).

[mtfishman]

Oh I see what you are saying. I guess it depends on where we implement it. We could just implement this at the level of randomITensor and ITensor (where we just run the QDense constructor, and if it ends up with no blocks then we throw an error), as opposed to actually in the QDense constructor itself.

[emstoudenmire]

Ok that's the right solution of course!