Add gauge transformation for iPEPS and its CTMRGEnv

Author: Yue-Zhengyuan

src/PEPSKit.jl

@@ -71,6 +71,7 @@ include("algorithms/ctmrg/projectors.jl")
 include("algorithms/ctmrg/simultaneous.jl")
 include("algorithms/ctmrg/sequential.jl")
 include("algorithms/ctmrg/gaugefix.jl")
+include("algorithms/ctmrg/gaugetrans.jl")
 
 include("algorithms/truncation/truncationschemes.jl")
 include("algorithms/truncation/fullenv_truncation.jl")
@@ -123,7 +124,7 @@ export InfinitePEPO, InfiniteTransferPEPO
 
 export BPEnv, BeliefPropagation
 export BPGauge, SUGauge
-export gauge_fix
+export gauge_fix, gauge_transform
 
 export initialize_mps, initializePEPS
 export ReflectDepth, ReflectWidth, Rotate, RotateReflect

src/algorithms/ctmrg/gaugetrans.jl

@@ -0,0 +1,41 @@
+"""
+    gauge_transform(env::CTMRGEnv, XXinv)
+
+Transform the CTMRG environment `env` of an InfinitePEPS after it is
+transformed by gauge transformations `XXinv` on its virtual bonds.
+
+`XXinv` consists of `(X, X⁻¹)` pairs on each virtual bond of the InfinitePEPS.
+```
+    T[r-1,c]
+    |
+    X⁻¹
+    |   [NORTH,r,c]
+    X
+    |
+    T[r,c]----X---X⁻¹----T[r,c+1]
+            [EAST,r,c]
+```
+"""
+function gauge_transform(env::CTMRGEnv, XXinv)
+    edges = map(eachcoordinate(env, 1:4)) do (d, r, c)
+        if d == NORTH
+            X⁻¹ = XXinv[NORTH, _next(r, end), c][2]
+            return @tensor edge[χ1 d0 d1; χ2] :=
+                env.edges[d, r, c][χ1 d0′ d1′; χ2] * X⁻¹[d0; d0′] * conj(X⁻¹[d1; d1′])
+        elseif d == EAST
+            X⁻¹ = XXinv[EAST, r, _prev(c, end)][2]
+            return @tensor edge[χ1 d0 d1; χ2] :=
+                env.edges[d, r, c][χ1 d0′ d1′; χ2] * X⁻¹[d0; d0′] * conj(X⁻¹[d1; d1′])
+        elseif d == SOUTH
+            X = XXinv[NORTH, r, c][1]
+            return @tensor edge[χ1 d0 d1; χ2] :=
+                env.edges[d, r, c][χ1 d0′ d1′; χ2] * X[d0′; d0] * conj(X[d1′; d1])
+        else # d == WEST
+            X = XXinv[EAST, r, c][1]
+            return @tensor edge[χ1 d0 d1; χ2] :=
+                env.edges[d, r, c][χ1 d0′ d1′; χ2] * X[d0′; d0] * conj(X[d1′; d1])
+        end
+    end
+    # corners are unaffected
+    return CTMRGEnv(env.corners, edges)
+end

src/states/infinitepeps.jl

@@ -270,3 +270,34 @@ function ChainRulesCore.rrule(
     end
     return network, InfiniteSquareNetwork_pullback
 end
+
+## Gauge transform
+
+"""
+    gauge_transform(ψ::InfinitePEPS, XXinv)
+
+Transform the InfinitePEPS `ψ` with gauge transformations `XXinv` on its virtual bonds.
+
+`XXinv` consists of `(X, X⁻¹)` pairs on each virtual bond of `ψ`.
+```
+    T[r-1,c]
+    |
+    X⁻¹
+    |   [NORTH,r,c]
+    X
+    |
+    T[r,c]----X---X⁻¹----T[r,c+1]
+            [EAST,r,c]
+```
+"""
+function gauge_transform(ψ::InfinitePEPS, XXinv)
+    A2 = map(eachcoordinate(ψ)) do (r, c)
+        Xn = XXinv[NORTH, r, c][1]
+        Xe = XXinv[EAST, r, c][1]
+        Xs = XXinv[NORTH, _next(r, end), c][2]
+        Xw = XXinv[EAST, r, _prev(c, end)][2]
+        return @tensor t[p; n e s w] := ψ[r, c][p; n′ e′ s′ w′] *
+            Xn[n′; n] * Xe[e′; e] * Xs[s; s′] * Xw[w; w′]
+    end
+    return InfinitePEPS(A2)
+end

test/ctmrg/gaugetrans.jl

@@ -0,0 +1,42 @@
+using Test
+using Random
+using TensorKit
+using PEPSKit
+using PEPSKit: random_dual!, eachcoordinate, NORTH
+
+ds = Dict(
+    U1Irrep => U1Space(i => d for (i, d) in zip(-1:1, (1, 1, 2))),
+    FermionParity => Vect[FermionParity](0 => 2, 1 => 1)
+)
+Ds = Dict(
+    U1Irrep => U1Space(i => D for (i, D) in zip(-1:1, (1, 3, 2))),
+    FermionParity => Vect[FermionParity](0 => 3, 1 => 2)
+)
+χs = Dict(
+    U1Irrep => U1Space(i => D for (i, D) in zip(-1:1, (2, 4, 2))),
+    FermionParity => Vect[FermionParity](0 => 4, 1 => 4)
+)
+
+@testset "CTMRGEnv of InfinitePEPS ($S)" for S in keys(ds)
+    d, D, χ, uc = ds[S], Ds[S], χs[S], (2, 3)
+    ψds = fill(d, uc)
+    ψDNs = random_dual!(fill(D, uc))
+    ψDEs = random_dual!(fill(D, uc))
+    ψ0 = InfinitePEPS(ψds, ψDNs, ψDEs)
+    env0 = CTMRGEnv(randn, ComplexF64, ψ0, χ)
+    # create a random set of gauge transformation on each bond
+    XXinv = map(eachcoordinate(ψ0, 1:2)) do (d, r, c)
+        V = (d == NORTH) ? ψDNs[r, c] : ψDEs[r, c]
+        X = randn(ComplexF64, V → V)
+        return (X, inv(X))
+    end
+    # gauge transform ψ and env
+    ψ1 = gauge_transform(ψ0, XXinv)
+    env1 = gauge_transform(env0, XXinv)
+    # reduced density matrices should remain unchanged
+    for (r, c) in eachcoordinate(ψ0)
+        ρ0 = reduced_densitymatrix(((r, c),), ψ0, env0)
+        ρ1 = reduced_densitymatrix(((r, c),), ψ1, env1)
+        @test ρ0 ≈ ρ1
+    end
+end

test/runtests.jl

@@ -28,6 +28,9 @@ end
         @time @safetestset "Gauge Fixing" begin
             include("ctmrg/gaugefix.jl")
         end
+        @time @safetestset "Gauge Transform" begin
+            include("ctmrg/gaugetrans.jl")
+        end
         @time @safetestset "Unit cell" begin
             include("ctmrg/unitcell.jl")
         end