ctmrg for fermions

src/algorithms/ctmrg.jl

@@ -120,17 +120,19 @@ function left_move(
             end
             #@ignore_derivatives @show norm(Q1*Q2)
 
-            (U, S, V) = tsvd(Q1 * Q2; trunc=trscheme, alg=SVD())
+            #(U, S, V) = tsvd(Q1 * Q2; trunc=trscheme, alg=SVD())
+            @tensor QQ[-1 -2 -3;-4 -5 -6] := Q1[-1 -2 -3;1 2 3] * Q2[1 2 3;-4 -5 -6]
+            (U, S, V) = tsvd(QQ; trunc=trscheme, alg=SVD())
 
-            @ignore_derivatives n0 = norm(Q1 * Q2)^2
+            @ignore_derivatives n0 = norm(QQ)^2
             @ignore_derivatives n1 = norm(U * S * V)^2
             @ignore_derivatives ϵ = max(ϵ, (n0 - n1) / n0)
 
             isqS = sdiag_inv_sqrt(S)
-            #Q = isqS*U'*Q1;
-            #P = Q2*V'*isqS;
-            @tensor Q[-1; -2 -3 -4] := isqS[-1; 1] * conj(U[2 3 4; 1]) * Q1[2 3 4; -2 -3 -4]
-            @tensor P[-1 -2 -3; -4] := Q2[-1 -2 -3; 1 2 3] * conj(V[4; 1 2 3]) * isqS[4; -4]
+            Q = isqS*U'*Q1;
+            P = Q2*V'*isqS;
+            #@tensor Q[-1; -2 -3 -4] := isqS[-1; 1] * conj(U[2 3 4; 1]) * Q1[2 3 4; -2 -3 -4]
+            #@tensor P[-1 -2 -3; -4] := Q2[-1 -2 -3; 1 2 3] * conj(V[4; 1 2 3]) * isqS[4; -4]
 
             @diffset above_projs[row] = Q
             @diffset below_projs[row] = P
@@ -296,12 +298,11 @@ function contract_ctrmg(
             envs.corners[SOUTHWEST, r, c][16; 1] *
             peps_above[r, c][17; 6 10 14 2] *
             conj(peps_below[r, c][17; 7 11 15 3])
-        total *= tr(
-            envs.corners[NORTHWEST, r, c] *
-            envs.corners[NORTHEAST, r, mod1(c - 1, end)] *
-            envs.corners[SOUTHEAST, mod1(r - 1, end), mod1(c - 1, end)] *
-            envs.corners[SOUTHWEST, mod1(r - 1, end), c],
-        )
+
+        total *=  @tensor envs.corners[NORTHWEST, r, c][1; 2] *
+            envs.corners[NORTHEAST, r, mod1(c - 1, end)][2; 3] *
+            envs.corners[SOUTHEAST, mod1(r - 1, end), mod1(c - 1, end)][3; 4] *
+            envs.corners[SOUTHWEST, mod1(r - 1, end), c][4; 1]
 
         total /= @tensor envs.edges[WEST, r, c][1 10 11; 4] *
             envs.corners[NORTHWEST, r, c][4; 5] *

src/states/abstractpeps.jl

@@ -62,4 +62,5 @@ Abstract supertype for a 2D projected entangled pairs operator.
 abstract type AbstractPEPO end
 
 Base.rotl90(t::PEPSTensor) = permute(t, ((1,), (3, 4, 5, 2)))
+Base.rotr90(t::PEPSTensor) = permute(t, ((1,), (5, 2, 3, 4)))
 Base.rotl90(t::PEPOTensor) = permute(t, ((1, 2), (4, 5, 6, 3)));

src/states/infinitepeps.jl

@@ -96,3 +96,4 @@ Base.axes(T::InfinitePEPS, args...) = axes(T.A, args...)
 TensorKit.space(t::InfinitePEPS, i, j) = space(t[i, j], 1)
 
 Base.rotl90(t::InfinitePEPS) = InfinitePEPS(rotl90(rotl90.(t.A)));
+Base.rotr90(t::InfinitePEPS) = InfinitePEPS(rotr90(rotr90.(t.A)));

src/utility/util.jl

@@ -44,6 +44,21 @@ function ChainRulesCore.rrule(::typeof(rotl90), a::AbstractMatrix)
     end
     return rotl90(a), pb
 end
+function ChainRulesCore.rrule(::typeof(rotr90), a::AbstractMatrix)
+    function pb(x)
+        if !iszero(x)
+            x = if x isa Tangent
+                ChainRulesCore.construct(typeof(a), ChainRulesCore.backing(x))
+            else
+                x
+            end
+            x = rotl90(x)
+        end
+
+        return (ZeroTangent(), x)
+    end
+    return rotr90(a), pb
+end
 
 structure(t) = codomain(t) ← domain(t);