Add support for fermionic systems

Project.toml

@@ -36,6 +36,8 @@ julia = "1.6"
 [extras]
 SafeTestsets = "1bc83da4-3b8d-516f-aca4-4fe02f6d838f"
 Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
+ChainRulesTestUtils = "cdddcdb0-9152-4a09-a978-84456f9df70a"
+FiniteDifferences = "26cc04aa-876d-5657-8c51-4c34ba976000"
 
 [targets]
-test = ["Test", "SafeTestsets"]
+test = ["Test", "SafeTestsets", "ChainRulesTestUtils", "FiniteDifferences"]

src/PEPSKit.jl

@@ -58,7 +58,7 @@ module Defaults
 end
 
 export CTMRG, CTMRGEnv
-export NLocalOperator, OnSite, NearestNeighbor
+export NLocalOperator, AnisotropicNNOperator, OnSite, NearestNeighbor
 export expectation_value, costfun
 export leading_boundary
 export PEPSOptimize, GeomSum, ManualIter, LinSolve

src/algorithms/ctmrg.jl

@@ -39,30 +39,30 @@ function MPSKit.leading_boundary(envinit, state, alg::CTMRG)
 
     for i in 1:(alg.maxiter)
         env, ϵ = ctmrg_iter(state, env, alg)  # Grow and renormalize in all 4 directions
+        conv_condition, normold, CSold, TSold, ϵ = ignore_derivatives() do
+            # Compute convergence criteria and take max (TODO: How should we handle logging all of this?)
+            Δϵ = abs((ϵold - ϵ) / ϵold)
+            normnew = norm(state, env)
+            Δnorm = abs(normold - normnew) / abs(normold)
+            CSnew = map(c -> tsvd(c; alg=TensorKit.SVD())[2], env.corners)
+            ΔCS = maximum(zip(CSold, CSnew)) do (c_old, c_new)
+                # only compute the difference on the smallest part of the spaces
+                smallest = infimum(MPSKit._firstspace(c_old), MPSKit._firstspace(c_new))
+                e_old = isometry(MPSKit._firstspace(c_old), smallest)
+                e_new = isometry(MPSKit._firstspace(c_new), smallest)
+                return norm(e_new' * c_new * e_new - e_old' * c_old * e_old)
+            end
+            TSnew = map(t -> tsvd(t; alg=TensorKit.SVD())[2], env.edges)
 
-        # Compute convergence criteria and take max (TODO: How should we handle logging all of this?)
-        Δϵ = abs((ϵold - ϵ) / ϵold)
-        normnew = norm(state, env)
-        Δnorm = abs(normold - normnew) / abs(normold)
-        CSnew = map(c -> tsvd(c; alg=TensorKit.SVD())[2], env.corners)
-        ΔCS = maximum(zip(CSold, CSnew)) do (c_old, c_new)
-            # only compute the difference on the smallest part of the spaces
-            smallest = infimum(MPSKit._firstspace(c_old), MPSKit._firstspace(c_new))
-            e_old = isometry(MPSKit._firstspace(c_old), smallest)
-            e_new = isometry(MPSKit._firstspace(c_new), smallest)
-            return norm(e_new' * c_new * e_new - e_old' * c_old * e_old)
-        end
-        TSnew = map(t -> tsvd(t; alg=TensorKit.SVD())[2], env.edges)
+            ΔTS = maximum(zip(TSold, TSnew)) do (t_old, t_new)
+                MPSKit._firstspace(t_old) == MPSKit._firstspace(t_new) ||
+                    return scalartype(t_old)(Inf)
+                # TODO: implement when spaces aren't the same
+                return norm(t_new - t_old)
+            end
 
-        ΔTS = maximum(zip(TSold, TSnew)) do (t_old, t_new)
-            MPSKit._firstspace(t_old) == MPSKit._firstspace(t_new) ||
-                return scalartype(t_old)(Inf)
-            # TODO: implement when spaces aren't the same
-            return norm(t_new - t_old)
-        end
+            conv_condition = max(Δnorm, ΔCS, ΔTS) < alg.tol && i > alg.miniter
 
-        conv_condition = max(Δnorm, ΔCS, ΔTS) < alg.tol && i > alg.miniter
-        ignore_derivatives() do # Print verbose info
             if alg.verbosity > 1 || (alg.verbosity == 1 && (i == 1 || conv_condition))
                 @printf(
                     "CTMRG iter: %3d   norm: %.2e   Δnorm: %.2e   ΔCS: %.2e   ΔTS: %.2e   ϵ: %.2e   Δϵ: %.2e\n",
@@ -80,14 +80,9 @@ function MPSKit.leading_boundary(envinit, state, alg::CTMRG)
                 @warn(
                     "CTMRG reached maximal number of iterations at (Δnorm=$Δnorm, ΔCS=$ΔCS, ΔTS=$ΔTS)"
                 )
+            return conv_condition, normnew, CSnew, TSnew, ϵ
         end
         conv_condition && break  # Converge if maximal Δ falls below tolerance
-
-        # Update convergence criteria
-        normold = normnew
-        CSold = CSnew
-        TSold = TSnew
-        ϵold = ϵ
     end
 
     # Do one final iteration that does not change the spaces
@@ -96,11 +91,22 @@ function MPSKit.leading_boundary(envinit, state, alg::CTMRG)
     )
     env′, = ctmrg_iter(state, env, alg_fixed)
     envfix = gauge_fix(env, env′)
-    check_elementwise_convergence(env, envfix; atol=alg.tol^(3 / 4)) ||
+    check_elementwise_convergence(env, envfix; atol=alg.tol^(1 / 2)) ||
         @warn "CTMRG did not converge elementwise."
     return envfix
 end
 
+#very ugly fix to avoid Zygoting @plansor during fermionic gauge fix
+@generated function MPSKit.transfer_right(
+    v::AbstractTensorMap{S,1,N₁}, A::GenericMPSTensor{S,N₂}, Ā::GenericMPSTensor{S,N₂}
+) where {S,N₁,N₂}
+    t_out = MPSKit.tensorexpr(:v, -1, -(2:(N₁ + 1)))
+    t_top = MPSKit.tensorexpr(:A, (-1, reverse(3:(N₂ + 1))...), 1)
+    t_bot = MPSKit.tensorexpr(:Ā, (-(N₁ + 1), reverse(3:(N₂ + 1))...), 2)
+    t_in = MPSKit.tensorexpr(:v, 1, (-(2:N₁)..., 2))
+    return :(return @tensor $t_out := $t_top * conj($t_bot) * $t_in)
+end
+
 """
     gauge_fix(envprev::CTMRGEnv{C,T}, envfinal::CTMRGEnv{C,T}) where {C,T}
 
@@ -145,6 +151,7 @@ function gauge_fix(envprev::CTMRGEnv{C,T}, envfinal::CTMRGEnv{C,T}) where {C,T}
             scalartype(T),
             MPSKit._lastspace(Tsfinal[end])' ← MPSKit._lastspace(M[end])',
         )
+
         ρprev = eigsolve(TransferMatrix(Tsprev, M), ρinit, 1, :LM)[2][1]
         ρfinal = eigsolve(TransferMatrix(Tsfinal, M), ρinit, 1, :LM)[2][1]
 
@@ -329,7 +336,9 @@ function left_move(state, env::CTMRGEnv{C,T}, alg::CTMRG) where {C,T}
             else
                 alg.trscheme
             end
-            U, S, V, ϵ_local = tsvd!(Q_sw * Q_nw; trunc=trscheme, alg=TensorKit.SVD())  # TODO: Add field in CTMRG to choose SVD function
+            @tensor QQ[-1 -2 -3; -4 -5 -6] := Q_sw[-1 -2 -3; 1 2 3] * Q_nw[1 2 3; -4 -5 -6]
+            U, S, V, ϵ_local = tsvd!(QQ; trunc=trscheme, alg=TensorKit.SVD())
+            #U, S, V, ϵ_local = tsvd!(Q_sw * Q_nw; trunc=trscheme, alg=TensorKit.SVD())  # TODO: Add field in CTMRG to choose SVD function
             ϵ = max(ϵ, ϵ_local / norm(S))
             # TODO: check if we can just normalize enlarged corners s.t. trunc behaves a bit better
 
@@ -409,8 +418,10 @@ function build_projectors(
     U::AbstractTensorMap{E,3,1}, S, V::AbstractTensorMap{E,1,3}, Q_sw, Q_nw
 ) where {E<:ElementarySpace}
     isqS = sdiag_inv_sqrt(S)
-    @tensor Pl[-1 -2 -3; -4] := Q_nw[-1 -2 -3; 1 2 3] * conj(V[4; 1 2 3]) * isqS[4; -4]
-    @tensor Pr[-1; -2 -3 -4] := isqS[-1; 1] * conj(U[2 3 4; 1]) * Q_sw[2 3 4; -2 -3 -4]
+    Pl = Q_nw * V' * isqS
+    Pr = isqS * U' * Q_sw
+    #@tensor Pl[-1 -2 -3; -4] := Q_nw[-1 -2 -3; 1 2 3] * conj(V[4; 1 2 3]) * isqS[4; -4]
+    #@tensor Pr[-1; -2 -3 -4] := isqS[-1; 1] * conj(U[2 3 4; 1]) * Q_sw[2 3 4; -2 -3 -4]
     return Pl, Pr
 end
 
@@ -448,12 +459,10 @@ function LinearAlgebra.norm(peps::InfinitePEPS, env::CTMRGEnv)
             peps[r, c][17; 6 10 14 2] *
             conj(peps[r, c][17; 7 11 15 3])
 
-        total *= tr(
-            env.corners[NORTHWEST, r, c] *
-            env.corners[NORTHEAST, r, mod1(c - 1, end)] *
-            env.corners[SOUTHEAST, mod1(r - 1, end), mod1(c - 1, end)] *
-            env.corners[SOUTHWEST, mod1(r - 1, end), c],
-        )
+        total *= @tensor env.corners[NORTHWEST, r, c][1; 2] *
+            env.corners[NORTHEAST, r, mod1(c - 1, end)][2; 3] *
+            env.corners[SOUTHEAST, mod1(r - 1, end), mod1(c - 1, end)][3; 4] *
+            env.corners[SOUTHWEST, mod1(r - 1, end), c][4; 1]
 
         total /= @tensor env.edges[WEST, r, c][1 10 11; 4] *
             env.corners[NORTHWEST, r, c][4; 5] *

src/algorithms/peps_opt.jl

@@ -112,6 +112,7 @@ function _rrule(
     alg::CTMRG,
 )
     envs = leading_boundary(envinit, state, alg)
+    #TODO: fixed space for unit cells
 
     function leading_boundary_pullback(Δenvs′)
         Δenvs = unthunk(Δenvs′)

src/operators/localoperator.jl

@@ -24,6 +24,18 @@ struct NLocalOperator{I<:AbstractInteraction}
     op::AbstractTensorMap
 end
 
+"""
+    struct NLocalOperator{I<:AbstractInteraction}
+
+Operator in form of a `AbstractTensorMap` which is parametrized by an interaction type.
+Mostly, this is used to define Hamiltonian terms and observables.
+"""
+struct AnisotropicNNOperator
+    h0::NLocalOperator{OnSite}
+    hx::NLocalOperator{NearestNeighbor}
+    hy::NLocalOperator{NearestNeighbor}
+end
+
 @doc """
     operator_env(peps::InfinitePEPS, env::CTMRGEnv, ::AbstractInteraction)
 
@@ -178,3 +190,17 @@ function costfun(peps::InfinitePEPS, env, op::NLocalOperator{NearestNeighbor})
     ov = sum(expectation_value(rotl90(peps), rotl90(env), op))
     return real(oh + ov)
 end
+
+"""
+    costfun(peps::InfinitePEPS, env, op::AnisotropicNNOperator)
+
+Compute the expectation value of an on-site and an anisotropic nearest-neighbor operator.
+This is used to evaluate and differentiate the energy in ground-state PEPS optimizations.
+"""
+function costfun(peps::InfinitePEPS, env, op::AnisotropicNNOperator)
+    oos = expectation_value(peps, env, op.h0)[1]
+    oh = sum(expectation_value(peps, env, op.hx))
+    ov = sum(expectation_value(rotr90(peps), rotate_north(env, WEST), op.hy))
+    #ov = sum(expectation_value(rotl90(peps), rotl90(env), op.hy))
+    return real(oos + oh + ov)
+end

src/states/abstractpeps.jl

@@ -62,4 +62,5 @@ Abstract supertype for a 2D projected entangled-pair 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

@@ -5,7 +5,7 @@ Represents an infinite projected entangled-pair state on a 2D square lattice.
 """
 struct InfinitePEPS{T<:PEPSTensor} <: AbstractPEPS
     A::Matrix{T}
-
+    InfinitePEPS{T}(A::Matrix{T}) where {T<:PEPSTensor} = new{T}(A)
     function InfinitePEPS(A::Array{T,2}) where {T<:PEPSTensor}
         for (d, w) in Tuple.(CartesianIndices(A))
             space(A[d, w], 2) == space(A[_prev(d, end), w], 4)' || throw(
@@ -160,3 +160,11 @@ function ChainRulesCore.rrule(::typeof(rotl90), peps::InfinitePEPS)
     end
     return peps′, rotl90_pullback
 end
+
+function ChainRulesCore.rrule(::typeof(rotr90), peps::InfinitePEPS)
+    peps′ = rotr90(peps)
+    function rotr90_pullback(Δpeps)
+        return NoTangent(), rotl90(Δpeps)
+    end
+    return peps′, rotr90_pullback
+end

src/utility/symmetrization.jl

@@ -61,7 +61,7 @@ end
 
 # rotation invariance
 
-Base.rotr90(t::PEPSTensor) = permute(t, ((1,), (5, 2, 3, 4)))
+#Base.rotr90(t::PEPSTensor) = permute(t, ((1,), (5, 2, 3, 4)))
 Base.rot180(t::PEPSTensor) = permute(t, ((1,), (4, 5, 2, 3)))
 
 function rot_inv(x)