Matt Fishman

Very strange bug, thanks @ArtemStrashko. By the way, I would recommend updating to the latest version of ITensors (0.3.14), though I see this bug in the latest version as well...

This seems to fix the initial error: ```julia using Zygote using ITensors indices = [Index(2) for _ in 1:10] mps = randomMPS(indices; linkdims=10) state = [randomITensor(index) for index in indices]...

Actually, there's a more minimal fix. These both run without errors for me: ```julia using Zygote using ITensors indices = [Index(2) for _ in 1:10] mps = randomMPS(indices; linkdims=10) state...

So oddly enough this seems to be a Zygote bug: ```julia using FiniteDifferences using Zygote function f(x) y = [[x]', [x]] r = 0.0 o = 1.0 for n in...

This isn't particular to heterogeneous lattices, it can also happen with homogeneous lattices with certain starting states and desired bond dimensions. I think we should just make this a warning...

Interesting, looks like somehow the density matrix ends up with a nonzero flux. We should probably actually print `flux(A)` in that error message.

Look into [Broadcast.flatten](https://github.com/JuliaLang/julia/blob/v1.6.0/base/broadcast.jl#L282-L318) as a more generic way to forward the function call to NDTensors. Hopefully that wouldn't lead to type instabilities and performance issues.

Some notes on `A .+= c` and `A .= f.(B)` for QN ITensors, continued from a discussion here: https://github.com/ITensor/ITensors.jl/pull/636#discussion_r622298592. Right now there is this behavior: ```julia i = Index([QN(0) =>...

Sort of, but I think it is better to think about it in terms of the properties that we want to satisfy. I think ultimately the properties we want for...

The main corner case with defining it as `sites==siteinds(psi)` is the case when there are multiple site indices per tensor. A more general definition would be something like `all(s ->...