Implement evolve of an MPS with an MPO

src/Ansatz/Chain.jl

@@ -652,6 +652,63 @@ function unpack_2sitewf!(ψ::Chain, bond, left_inds, right_inds, virtualind)
     return ψ
 end
 
+"""
+    evolve!(ψ::Chain, mpo::Chain)
+
+Applies a Matrix Product Operator (MPO) `mpo` to the [`Chain`](@ref).
+"""
+function evolve(ψ::Chain, mpo::Chain)
+    updated_tensors = Tensor[]
+    Λ = Tensor[]
+    L = nsites(ψ)
+
+    for i in 1:L
+        t = contract(select(ψ, :tensor, Site(i)), select(mpo, :tensor, Site(i)); dims = (select(ψ, :index, Site(i)),))
+        physicalind = select(mpo, :index, Site(i))
+
+        # Fuse the two right legs of t into one
+        if i == 1
+            wanted_inds = (physicalind, rightindex(ψ, Site(i)), rightindex(mpo, Site(i)))
+            new_inds = (physicalind, rightindex(ψ, Site(i)))
+        elseif i < L
+            wanted_inds = (
+                physicalind,
+                leftindex(ψ, Site(i)),
+                leftindex(mpo, Site(i)),
+                rightindex(ψ, Site(i)),
+                rightindex(mpo, Site(i)),
+            )
+            new_inds = (physicalind, leftindex(ψ, Site(i)), rightindex(ψ, Site(i)))
+        else
+            wanted_inds = (physicalind, leftindex(ψ, Site(i)), leftindex(mpo, Site(i)))
+            new_inds = (physicalind, leftindex(ψ, Site(i)))
+        end
+
+        perm = Tenet.__find_index_permutation(wanted_inds, inds(t))
+        t = PermutedDimsArray(parent(t), perm)
+
+        t = Tensor(
+            reshape(t, tuple(size(t, 1), [size(t, k) * size(t, k + 1) for k in 2:2:length(wanted_inds)]...)),
+            new_inds,
+        )
+        push!(updated_tensors, t)
+
+        if i < L
+            d = size(TensorNetwork(mpo), rightindex(mpo, Site(i)))
+            Λᵢ = select(ψ, :between, Site(i), Site(i + 1))
+            Λᵢ = Tensor(diag(kron(Matrix(LinearAlgebra.I, d, d), diagm(parent(Λᵢ)))), inds(Λᵢ))
+            push!(Λ, Λᵢ)
+        end
+    end
+
+    ψ_ev = MPS(updated_tensors)
+    for i in 1:L-1
+        push!(TensorNetwork(ψ_ev), Tensor(parent(Λ[i]), (rightindex(ψ_ev, Site(i)),)))
+    end
+
+    return ψ_ev
+end
+
 function expect(ψ::Chain, observables)
     # contract observable with TN
     ϕ = copy(ψ)