Refactor operator types

docs/src/man/models.md

@@ -28,7 +28,7 @@ function MPSKitModels.transverse_field_ising(
     return LocalOperator(
         spaces,
         (neighbor => ZZ for neighbor in nearest_neighbours(lattice))...,
-        ((idx,) => X for idx in vertices(lattice))...,
+        ([idx,] => X for idx in vertices(lattice))...,
     )
 end
 ```

src/algorithms/contractions/bp_contractions.jl

@@ -60,13 +60,51 @@ function MPSKit.expectation_value(peps::InfinitePEPS, O::LocalOperator, env::BPE
 end
 
 function contract_local_operator(
-        inds::NTuple{1, CartesianIndex{2}},
+        inds::Vector{CartesianIndex{2}},
+        O::AbstractTensorMap,
+        ket::InfinitePEPS,
+        bra::InfinitePEPS,
+        env::BPEnv,
+    )
+    length(inds) == 1 && return contract_local_operator1x1(only(inds), O, ket, bra, env)
+
+    if length(inds) == 2
+        ind_relative = inds[2] - inds[1]
+        if ind_relative == CartesianIndex(1, 0)
+            return contract_local_operator2x1(inds[1], O, ket, bra, env)
+        elseif ind_relative == CartesianIndex(0, 1)
+            return contract_local_operator1x2(inds[1], O, ket, bra, env)
+        end
+    end
+    error("No implementation for contractions for BP environments with $inds")
+end
+function contract_local_norm(
+        inds::Vector{CartesianIndex{2}},
+        ket::InfinitePEPS,
+        bra::InfinitePEPS,
+        env::BPEnv,
+    )
+    length(inds) == 1 && return contract_local_norm1x1(only(inds), ket, bra, env)
+
+    if length(inds) == 2
+        ind_relative = inds[2] - inds[1]
+        if ind_relative == CartesianIndex(1, 0)
+            return contract_local_norm2x1(inds[1], ket, bra, env)
+        elseif ind_relative == CartesianIndex(0, 1)
+            return contract_local_norm1x2(inds[1], ket, bra, env)
+        end
+    end
+    error("No implementation for contractions for BP environments with $inds")
+end
+
+function contract_local_operator1x1(
+        ind::CartesianIndex{2},
         O::AbstractTensorMap{<:Any, <:Any, 1, 1},
         ket::InfinitePEPS,
         bra::InfinitePEPS,
         env::BPEnv,
     )
-    row, col = Tuple(only(inds))
+    row, col = Tuple(ind)
     M_north = env.messages[NORTH, _prev(row, end), mod1(col, end)]
     M_east = env.messages[EAST, mod1(row, end), _next(col, end)]
     M_south = env.messages[SOUTH, _next(row, end), mod1(col, end)]
@@ -83,10 +121,10 @@ function contract_local_operator(
     end
 end
 
-function contract_local_norm(
-        inds::NTuple{1, CartesianIndex{2}}, ket::InfinitePEPS, bra::InfinitePEPS, env::BPEnv
+function contract_local_norm1x1(
+        ind::CartesianIndex{2}, ket::InfinitePEPS, bra::InfinitePEPS, env::BPEnv
     )
-    row, col = Tuple(only(inds))
+    row, col = Tuple(ind)
     M_north = env.messages[NORTH, _prev(row, end), mod1(col, end)]
     M_east = env.messages[EAST, mod1(row, end), _next(col, end)]
     M_south = env.messages[SOUTH, _next(row, end), mod1(col, end)]
@@ -102,24 +140,7 @@ function contract_local_norm(
     end
 end
 
-function contract_local_operator(
-        inds::NTuple{2, CartesianIndex{2}},
-        O::AbstractTensorMap{<:Any, <:Any, 2, 2},
-        ket::InfinitePEPS,
-        bra::InfinitePEPS,
-        env::BPEnv,
-    )
-    ind_relative = inds[2] - inds[1]
-    return if ind_relative == CartesianIndex(1, 0)
-        contract_vertical_operator(inds[1], O, ket, bra, env)
-    elseif ind_relative == CartesianIndex(0, 1)
-        contract_horizontal_operator(inds[1], O, ket, bra, env)
-    else
-        error("Only contractions for nearest neighbor bonds are implemented.")
-    end
-end
-
-function contract_vertical_operator(
+function contract_local_operator2x1(
         coord::CartesianIndex{2},
         O::AbstractTensorMap{<:Any, <:Any, 2, 2},
         ket::InfinitePEPS,
@@ -147,7 +168,7 @@ function contract_vertical_operator(
         O[dNb dSb; dNt dSt]
 end
 
-function contract_horizontal_operator(
+function contract_local_operator1x2(
         coord::CartesianIndex{2},
         O::AbstractTensorMap{<:Any, <:Any, 2, 2},
         ket::InfinitePEPS,
@@ -181,20 +202,7 @@ function contract_horizontal_operator(
     end
 end
 
-function contract_local_norm(
-        inds::NTuple{2, CartesianIndex{2}}, ket::InfinitePEPS, bra::InfinitePEPS, env::BPEnv
-    )
-    ind_relative = inds[2] - inds[1]
-    return if ind_relative == CartesianIndex(1, 0)
-        contract_vertical_norm(inds[1], ket, bra, env)
-    elseif ind_relative == CartesianIndex(0, 1)
-        contract_horizontal_norm(inds[1], ket, bra, env)
-    else
-        error("Only contractions for nearest neighbor bonds are implemented.")
-    end
-end
-
-function contract_vertical_norm(
+function contract_local_norm2x1(
         coord::CartesianIndex{2}, ket::InfinitePEPS, bra::InfinitePEPS, env::BPEnv
     )
     row, col = Tuple(coord)
@@ -217,7 +225,7 @@ function contract_vertical_norm(
         M_northwest[DNWb; DNWt]
 end
 
-function contract_horizontal_norm(
+function contract_local_norm1x2(
         coord::CartesianIndex{2}, ket::InfinitePEPS, bra::InfinitePEPS, env::BPEnv
     )
     row, col = Tuple(coord)

src/algorithms/contractions/localoperator.jl

@@ -17,23 +17,24 @@ specified by `inds`. The `peps` is contracted with `O` from above and below, and
 sandwich is surrounded with the appropriate environment tensors.
 """
 function contract_local_operator(
-        inds::NTuple{N, CartesianIndex{2}},
-        O::AbstractTensorMap{T, S, N, N},
-        ket::InfinitePEPS, bra::InfinitePEPS,
-        env::CTMRGEnv,
-    ) where {T, S, N}
-    static_inds = Val.(inds)
+        inds::Vector{CartesianIndex{2}}, O::AbstractTensorMap,
+        ket::InfinitePEPS, bra::InfinitePEPS, env::CTMRGEnv,
+    )
+    static_inds = Tuple(Val.(inds))
     return _contract_local_operator(static_inds, O, (ket, bra), env)
 end
 function contract_local_operator(
-        inds::NTuple{N, Tuple{Int, Int}},
-        O::AbstractTensorMap{T, S, N, N},
-        ket::InfinitePEPS, bra::InfinitePEPS,
-        env::CTMRGEnv,
-    ) where {T, S, N}
+        inds::Vector{Tuple{Int, Int}}, O::AbstractTensorMap,
+        ket::InfinitePEPS, bra::InfinitePEPS, env::CTMRGEnv,
+    )
     return contract_local_operator(CartesianIndex.(inds), O, ket, bra, env)
 end
 
+Base.@deprecate(
+    contract_local_operator(inds::NTuple, ket::InfinitePEPS, bra::InfinitePEPS, env::CTMRGEnv),
+    contract_local_operator(collect(inds), ket, bra, env)
+)
+
 # This implements the contraction of an operator acting on sites `inds`.
 # The generated function ensures that we can use @tensor to write dynamic contractions (and maximize performance).
 
@@ -264,14 +265,14 @@ on a rectangular patch based on `inds` but replacing the operator with an identi
 that the PEPS norm is computed. (Note that this is not the physical norm of the state.)
 """
 function contract_local_norm(
-        inds::NTuple{N, CartesianIndex{2}}, ket::InfinitePEPS, bra::InfinitePEPS, env::CTMRGEnv
-    ) where {N}
-    static_inds = Val.(inds)
+        inds::Vector{CartesianIndex{2}}, ket::InfinitePEPS, bra::InfinitePEPS, env::CTMRGEnv
+    )
+    static_inds = Tuple(Val.(inds))
     return _contract_local_norm(static_inds, (ket, bra), env)
 end
 function contract_local_norm(
-        inds::NTuple{N, Tuple{Int, Int}}, ket::InfinitePEPS, bra::InfinitePEPS, env::CTMRGEnv
-    ) where {N}
+        inds::Vector{Tuple{Int, Int}}, ket::InfinitePEPS, bra::InfinitePEPS, env::CTMRGEnv
+    )
     return contract_local_norm(CartesianIndex.(inds), ket, bra, env)
 end
 @generated function _contract_local_norm(
@@ -299,6 +300,11 @@ end
     return macroexpand(@__MODULE__, returnex)
 end
 
+Base.@deprecate(
+    contract_local_norm(inds::NTuple, ket::InfinitePEPS, bra::InfinitePEPS, env::CTMRGEnv),
+    contract_local_norm(collect(inds), ket, bra, env)
+)
+
 @doc """
 $(SIGNATURES)
 
@@ -310,52 +316,84 @@ specified by `inds`. The result is normalized such that `tr(ρ) = 1`.
 """ reduced_densitymatrix
 
 function reduced_densitymatrix(
-        inds::NTuple{N, CartesianIndex{2}}, ket::InfinitePEPS, bra::InfinitePEPS, env::CTMRGEnv
-    ) where {N}
-    static_inds = Val.(inds)
+        inds::Vector{CartesianIndex{2}}, ket::InfinitePEPS, bra::InfinitePEPS, env::CTMRGEnv
+    )
+    length(inds) == 1 && return reduced_densitymatrix1x1(only(inds), ket, bra, env)
+
+    if length(inds) == 2
+        if inds[2] - inds[1] == CartesianIndex(1, 0)
+            return reduced_densitymatrix2x1(inds[1], ket, bra, env)
+        elseif inds[2] - inds[1] == CartesianIndex(0, 1)
+            return reduced_densitymatrix1x2(inds[1], ket, bra, env)
+        end
+    end
+
+    static_inds = Tuple(Val.(inds))
     return _contract_densitymatrix(static_inds, (ket, bra), env)
 end
 function reduced_densitymatrix(
-        inds::NTuple{N, CartesianIndex{2}}, state::InfinitePEPO, env::CTMRGEnv
-    ) where {N}
+        inds::Vector{CartesianIndex{2}}, state::InfinitePEPO, env::CTMRGEnv
+    )
     size(state, 3) == 1 || throw(DimensionMismatch("only single-layer densitymatrices are supported"))
-    static_inds = Val.(inds)
+    static_inds = Tuple(Val.(inds))
     return _contract_densitymatrix(static_inds, (state,), env)
 end
 function reduced_densitymatrix(
-        inds::NTuple{N, CartesianIndex{2}}, ket::InfinitePEPO, bra::InfinitePEPO, env::CTMRGEnv
-    ) where {N}
+        inds::Vector{CartesianIndex{2}}, ket::InfinitePEPO, bra::InfinitePEPO, env::CTMRGEnv
+    )
     size(ket) == size(bra) || throw(DimensionMismatch("incompatible bra and ket dimensions"))
     size(ket, 3) == 1 || throw(DimensionMismatch("only single-layer densitymatrices are supported"))
-    static_inds = Val.(inds)
+    static_inds = Tuple(Val.(inds))
     return _contract_densitymatrix(static_inds, (ket, bra), env)
 end
-function reduced_densitymatrix(
+reduced_densitymatrix(inds, ket::InfinitePEPS, env::CTMRGEnv) =
+    reduced_densitymatrix(inds, ket, ket, env)
+
+Base.@deprecate(
+    reduced_densitymatrix(
+        inds::NTuple{N, CartesianIndex{2}}, ket::InfinitePEPS, bra::InfinitePEPS, env::CTMRGEnv
+    ) where {N},
+    reduced_densitymatrix(collect(inds), ket, bra, env)
+)
+Base.@deprecate(
+    reduced_densitymatrix(
+        inds::NTuple{N, CartesianIndex{2}}, state::InfinitePEPO, env::CTMRGEnv
+    ) where {N},
+    reduced_densitymatrix(collect(inds), state, env)
+)
+Base.@deprecate(
+    reduced_densitymatrix(
+        inds::NTuple{N, CartesianIndex{2}}, ket::InfinitePEPO, bra::InfinitePEPO, env::CTMRGEnv
+    ) where {N},
+    reduced_densitymatrix(collect(inds), ket, bra, env)
+)
+
+Base.@deprecate(
+    reduced_densitymatrix(
         inds::NTuple{N, Tuple{Int, Int}}, ket::InfinitePEPS, bra::InfinitePEPS, env::CTMRGEnv
-    ) where {N}
-    return reduced_densitymatrix(CartesianIndex.(inds), ket, bra, env)
-end
-function reduced_densitymatrix(
+    ) where {N},
+    reduced_densitymatrix(collect(CartesianIndex.(inds)), ket, bra, env)
+)
+Base.@deprecate(
+    reduced_densitymatrix(
         inds::NTuple{N, Tuple{Int, Int}}, ket::InfinitePEPO, bra::InfinitePEPO, env::CTMRGEnv
-    ) where {N}
-    return reduced_densitymatrix(CartesianIndex.(inds), ket, bra, env)
-end
-function reduced_densitymatrix(inds, ket::InfinitePEPS, env::CTMRGEnv)
-    return reduced_densitymatrix(inds, ket, ket, env)
-end
-function reduced_densitymatrix(
+    ) where {N},
+    reduced_densitymatrix(collect(CartesianIndex.(inds)), ket, bra, env)
+)
+Base.@deprecate(
+    reduced_densitymatrix(
         inds::NTuple{N, Tuple{Int, Int}}, state::InfinitePEPO, env::CTMRGEnv
-    ) where {N}
-    return reduced_densitymatrix(CartesianIndex.(inds), state, env)
-end
+    ) where {N},
+    reduced_densitymatrix(collect(CartesianIndex.(inds)), state, env)
+)
 
 # Special case 1x1 density matrix:
 # Keep contraction order but try to optimize intermediate permutations:
 # EE_SWA is largest object so keep largest legs to the front there
-function reduced_densitymatrix(
-        inds::Tuple{CartesianIndex{2}}, ket::InfinitePEPS, bra::InfinitePEPS, env::CTMRGEnv
+function reduced_densitymatrix1x1(
+        inds::CartesianIndex{2}, ket::InfinitePEPS, bra::InfinitePEPS, env::CTMRGEnv
     )
-    row, col = Tuple(inds[1])
+    row, col = Tuple(inds)
 
     # Unpack variables and absorb corners
     A = ket[mod1(row, end), mod1(col, end)]

src/algorithms/time_evolution/gaugefix_su.jl

@@ -28,7 +28,7 @@ function _trivial_gates(elt::Type{<:Number}, lattice::Matrix{S}) where {S <: Ele
         h = TensorKit.id(elt, V1 ⊗ V2)
         return [site1, site2] => h
     end
-    return TrotterGates(lattice, vec(gates))
+    return LocalCircuit(lattice, vec(gates))
 end
 
 """

src/algorithms/time_evolution/simpleupdate.jl

@@ -108,6 +108,7 @@ function _su_iter!(
     # rotate back
     A = _bond_rotation(A, bond[1], rev; inv = true)
     B = _bond_rotation(B, bond[1], rev; inv = true)
+    rev && (s = transpose(s))
     # remove environment weights
     siteA, siteB = map(sites) do site
         return CartesianIndex(mod1(site[1], Nr), mod1(site[2], Nc))
@@ -127,12 +128,12 @@ end
 One iteration of simple update
 """
 function su_iter(
-        state::InfiniteState, gates::TrotterGates,
+        state::InfiniteState, circuit::LocalCircuit,
         alg::SimpleUpdate, env::SUWeight
     )
     Nr, Nc, = size(state)
     state2, env2, ϵ = deepcopy(state), deepcopy(env), 0.0
-    for (sites, gate) in gates.terms
+    for (sites, gate) in circuit.gates
         if length(sites) == 1
             # 1-site gate
             # TODO: special treatment for bipartite state