Question regarding "LoadError: setindex! not consistent with current flux"

[sujay-kazi]

I have a question regarding a particular error I am getting. Basically, I want to implement an operator that looks like Op = x(sum from j=1 to N-2)[(Sx)j (Sz)(j+1) (Sy)(j+2) - (Sy)j (Sz)(j+1) * (Sx)(j+2)] as a pseudo-momentum operator (here N and x are constants). In general, I want to able to compute <psi|Op^2|psi> for a state |psi> expressed as an MPS, which is why I am creating Op as an MPO. However, when I run the following code:

using ITensors
N = 10
x = 1
sites = siteinds("S=1/2", N, conserve_sz=true) # change to false and this code will work
ampo_Op = AutoMPO()
for j = 1:N-2
  add!(ampo_Op,x,"Sx",j,"Sz",j+1,"Sy",j+2)
  add!(ampo_Op,-x,"Sy",j,"Sz",j+1,"Sx",j+2)
end
Op = MPO(ampo_Op, sites)

I get the following error:

ERROR: setindex! not consistent with current flux
Stacktrace:
 [1] setindex!(::ITensor{2}, ::Float64, ::Int64, ::Int64) at C:\Users\user\.julia\packages\ITensors\MgdBe\src\itensor.jl:468
 [2] setindex!(::ITensor{2}, ::Float64, ::IndexVal{Index{Array{Pair{QN,Int64},1}}}, ::IndexVal{Index{Array{Pair{QN,Int64},1}}}) at C:\Users\user\.julia\packages\ITensors\MgdBe\src\itensor.jl:476
 [3] #op#265(::Base.Iterators.Pairs{Union{},Union{},Tuple{},NamedTuple{(),Tuple{}}}, ::typeof(op), ::TagType{S=1/2}, ::Index{Array{Pair{QN,Int64},1}}, ::SubString{String}) at C:\Users\user\.julia\packages\ITensors\MgdBe\src\physics\site_types\spinhalf.jl:42
 [4] op(::TagType{S=1/2}, ::Index{Array{Pair{QN,Int64},1}}, ::SubString{String}) at C:\Users\user\.julia\packages\ITensors\MgdBe\src\physics\site_types\spinhalf.jl:28
 [5] #_call_op#250(::Base.Iterators.Pairs{Union{},Union{},Tuple{},NamedTuple{(),Tuple{}}}, ::typeof(ITensors._call_op), ::Index{Array{Pair{QN,Int64},1}}, ::SubString{String}) at C:\Users\user\.julia\packages\ITensors\MgdBe\src\physics\tag_types.jl:37
 [6] _call_op at C:\Users\user\.julia\packages\ITensors\MgdBe\src\physics\tag_types.jl:22 [inlined]
 [7] #op#251(::Base.Iterators.Pairs{Union{},Union{},Tuple{},NamedTuple{(),Tuple{}}}, ::typeof(op), ::Index{Array{Pair{QN,Int64},1}}, ::String) at C:\Users\user\.julia\packages\ITensors\MgdBe\src\physics\tag_types.jl:63   
 [8] op at C:\Users\user\.julia\packages\ITensors\MgdBe\src\physics\tag_types.jl:44 [inlined]
 [9] (::ITensors.var"#calcQN#312"{Array{Index{Array{Pair{QN,Int64},1}},1}})(::Array{ITensors.SiteOp,1}) at C:\Users\user\.julia\packages\ITensors\MgdBe\src\physics\autompo.jl:393
 [10] #qn_svdMPO#300(::Base.Iterators.Pairs{Union{},Union{},Tuple{},NamedTuple{(),Tuple{}}}, ::typeof(ITensors.qn_svdMPO), ::AutoMPO, ::Array{Index{Array{Pair{QN,Int64},1}},1}) at C:\Users\user\.julia\packages\ITensors\MgdBe\src\physics\autompo.jl:399
 [11] qn_svdMPO at C:\Users\user\.julia\packages\ITensors\MgdBe\src\physics\autompo.jl:362 [inlined]
 [12] #MPO#330 at C:\Users\user\.julia\packages\ITensors\MgdBe\src\physics\autompo.jl:734 [inlined]
 [13] MPO(::AutoMPO, ::Array{Index{Array{Pair{QN,Int64},1}},1}) at C:\Users\user\.julia\packages\ITensors\MgdBe\src\physics\autompo.jl:732
 [14] top-level scope at none:0

As far as I can tell, this has to do with the QN conservation. I get this error if I have either conserve_qns or conserve_sz true, whereas if they are both false (as they are by default), then I don't get any error.

What exactly is the motivation for this restriction? It seems to be that there would be many instances (such as mine) where one would want to apply operators involving Sx and Sy even if Sz (summed over all sites) is the conserved quantity, and it seems that the code just doesn't let you do this.

[emstoudenmire]

Hi Sujay, so this error is most likely due to the use of the “Sx” and “Sy” operators. For total-Sz-conserving indices, please use the “S+” and “S-“ operators instead. The reason is that we have a strict convention in ITensor that when conserving quantum numbers, all tensors must have a well-defined QN “flux”, that is they must change the QNs which are being conserved by a definite amount. The Sx and Sy operators don’t have this property (they sometimes increase Sz, or sometimes decrease Sz, or a superposition of both).

Please let me know if you rewrite your terms using S+ and S- and it still doesn’t work.

Best, Miles

[mtfishman]

We could expand on that error message (add something like "You may be trying to create an ITensor that does not have a well defined quantum number flux.")

[emstoudenmire]

@mtfishman that’s a good idea to have an improved error message

@sujay-kazi regarding why we have the rule I mentioned above about all QN conserving ITensors needing to have a well-defined flux, thus ruling out the use of Sx and Sy with QNs, the reason right now is that although we could allow Sx and Sy etc., in many cases it could quickly lead to nominally block-sparse ITensors where all of the blocks are actually present and non-zero. In such cases, calculations would generally be slower than if block-sparsity wasn’t used at all. But as far as I can see, these cases would be very hard to detect on a tensor-by-tensor level without doing some analysis of like the typical fraction of non-zero blocks according to some heuristic. So we’d end up in a situation where users could easily get slower code, or incorrect code (changing of QNs) without any warning or reliable way to detect that it’s happening. Whereas in contrast, if one just uses S+ and S-, which is not too hard to do, then we can always be 100% sure that codes are as fast as possible and can’t superpose different QN sectors.

[sujay-kazi]

Thank you for the clarification! It works now with the updated version of the code.

If you don't mind, I have another question about computing <psi|Op^2|psi>. What's the logical way of doing this? I couldn't seem to find a built-in function that does this, so I was tempted to try to compute |psi_2> = Op|psi> as an ITensor and then contract |psi_2> with itself (this is equivalent because Op is Hermitian). I wanted to make |psi_2> be an MPS, the same way |psi> is, but I didn't see how to do this, so I was forced to multiply the tensors together to make one big ITensor.

There are two problems I have. First, when I try to do the contraction as psi_2 * psi_2, I end up getting an error like "ERROR: LoadError: Indices must have opposite directions to contract." Second, for any decently sized chain, this would result in warnings about the number of the indices in the ITensor (for instance, if N = 100, then |psi_2> will have 100 indices, and the program starts giving warnings if any ITensor has at least 15 indices).

In general, I assume there is a better way of handling this, but I have been pretty confused about this for a while now. In general, I want to be able to apply an MPO to an MPS to get a new MPS (that way I can take an inner product to find an expectation value). Any help on this would be appreciated.

[emstoudenmire]

Hi Sujay, I see you are asking two different questions, kind of:

1. how to multiply an MPS times an MPO to get a new MPS
2. how to compute an expectation value of the form <psi|Op^2|psi> where |psi> is an MPS and Op an MPO

The answer to #2 is that you can use the function inner(B::MPO, y::MPS, A::MPO, x::MPS) as inner(Op,psi,Op,psi) to compute <psi|Op*Op|psi>.

The answer to #1 is that you can use the function contract(A::MPO,psi::MPS) to multiply an MPO times an MPS using various approximate algorithms. However, there are some issues we need to work out with the contract method for MPO*MPS currently, such as the default algorithm having a bug for QN-conserving tensors (issue #320).

Finally, please consider posting support questions that would be of general interest, and which are not bug reports, on the ITensor Support Forum here: http://itensor.org/support/ If you need an account, please just email support@itensor.org with the username you'd like - thanks!

And thanks again for all the useful bug reports so far; it's really helping us to catch a lot of issues before we register the package soon.

Best, Miles

[mtfishman]

I've improved the error message. I'll close this since I don't think there is an issue to fix, but feel free to continue the discussion on the support forum.