Autodiff for twist?

[Confusio]

For fPEPS, I define
⊙(t1::AbstractTensorMap, t2::AbstractTensorMap)=twist(t1, filter(i -> !isdual(space(t1, i)), domainind(t1)))*t2
to take account the fermionic parity factors during the contraction. However, I get an error for the autodiff. The complete the code is following:

using MKL, TensorKit, LinearAlgebra,Zygote
⊙(t1::AbstractTensorMap, t2::AbstractTensorMap)=twist(t1, filter(i -> !isdual(space(t1, i)), domainind(t1)))*t2

V=Vect[FermionParity⊠SU2Irrep]((0, 0) => 2, (1, 1 // 2) => 1)
space_1=Vect[FermionParity⊠SU2Irrep]((0, 0) => 2, (1, 1 // 2) => 1)'
space_2=Vect[FermionParity⊠SU2Irrep]((0, 0) => 2, (1, 1 // 2) => 1)'
space_3=Vect[FermionParity⊠SU2Irrep]((0, 0) => 2, (1, 1 // 2) => 1)
space_4=Vect[FermionParity⊠SU2Irrep]((0, 0) => 2, (1, 1 // 2) => 1)
As=randn(ComplexF64,space_1⊗space_2⊗space_3⊗space_4⊗V,one(V))

f(x)=real(scalar(x'⊙x));
en, gd = withgradient(x -> f(x), As)

The output

ERROR: LoadError: Compiling Tuple{typeof(get!), TensorKit.var"#147#148"{TensorMapSpace{GradedSpace{ProductSector{Tuple{FermionParity, SU2Irrep}}, TensorKit.SortedVectorDict{ProductSector{Tuple{FermionParity, SU2Irrep}}, Int64}}, 0, 5}}, LRUCache.LRU{Any, Any}, TensorMapSpace{GradedSpace{ProductSector{Tuple{FermionParity, SU2Irrep}}, TensorKit.SortedVectorDict{ProductSector{Tuple{FermionParity, SU2Irrep}}, Int64}}, 0, 5}}: ArgumentError: array must be non-empty
Stacktrace:
  [1] macro expansion
    @ ~/.julia/packages/Zygote/3To5I/src/compiler/interface2.jl:0 [inlined]
...

Is this an error for the autodiff for twist? How to work around this in a simple way? ⊙ is heavily used in the code.

Julia and TensorKit version:
Julia Version 1.11.3 (2025-01-21)
[07d1fe3e] TensorKit v0.14.3

Many thanks.

[lkdvos]

Have you tried using @tensor to instantiate the contractions instead? This already includes these twists automatically and is compatible with AD.
Something like the following would do this:

@tensor C[a b c; d e] := A[a b c; f] * B[f; d e]

Or if you need it to be more dynamic:

using TensorOperations: tensorcontract
function mymul(A, B)
    pA = (codomainind(A), domainind(A))
    pB = (codomainind(B), domainind(B))
    pAB = (ntuple(identity, numout(A)), ntuple(i -> i + numout(A), numin(B)))
    conjA = conjB = false
    return tensorcontract(A, pA, conjA, B, pB, conjB, pAB)
end

This aside, indeed we don't currently have a rrule for twist, which might be useful so I'll look into adding it.

[Confusio]

Thanks for the dynamic solution.