Add NTU tests

Author: Yue-Zhengyuan

test/bondenv/benv_ntu.jl

@@ -0,0 +1,72 @@
+using Test
+using TensorKit
+using PEPSKit
+using LinearAlgebra
+using KrylovKit
+using Random
+
+Nr, Nc = 2, 2
+Random.seed!(20)
+
+function test_ntu_env(
+        state::Union{InfinitePEPS, InfinitePEPO}, env_alg::Alg
+    ) where {Alg <: PEPSKit.NeighbourEnv}
+    @info "Testing $(typeof(env_alg))"
+    for row in 1:Nr, col in 1:Nc
+        cp1 = PEPSKit._next(col, Nc)
+        A, B = state.A[row, col], state.A[row, cp1]
+        X, a, b, Y = PEPSKit._qr_bond(A, B)
+        @tensor ab[DX DY; da db] := a[DX da D] * b[D db DY]
+        benv = PEPSKit.bondenv_ntu(row, col, X, Y, state, env_alg)
+        # this is a result of `_qr_bond`
+        @assert [isdual(space(benv, ax)) for ax in 1:numind(benv)] == [0, 0, 1, 1]
+        # NTU bond environments are exact and should be positive definite
+        @test benv' ≈ benv
+        benv = project_hermitian(benv)
+        nrm1 = PEPSKit.inner_prod(benv, ab, ab)
+        @test nrm1 ≈ real(nrm1)
+        D, U = eigh_full(benv)
+        @test minimum(D.data) >= 0
+        # gauge fixing
+        Z = PEPSKit.sdiag_pow(D, 0.5) * U'
+        @assert benv ≈ Z' * Z
+        Z2, a2, b2, (Linv, Rinv) = PEPSKit.fixgauge_benv(Z, a, b)
+        benv2 = Z2' * Z2
+        # gauge fixing should reduce condition number
+        cond = LinearAlgebra.cond(benv)
+        cond2 = LinearAlgebra.cond(benv2)
+        @test cond2 <= cond
+        @info "benv cond number: (gauge-fixed) $(cond2) ≤ $(cond) (initial)"
+        # verify gauge transformation of X, Y
+        @tensor a2b2[DX DY; da db] := a2[DX da D] * b2[D db DY]
+        nrm2 = PEPSKit.inner_prod(benv2, a2b2, a2b2)
+        X2, Y2 = PEPSKit._fixgauge_benvXY(X, Y, Linv, Rinv)
+        benv3 = PEPSKit.bondenv_ntu(row, col, X2, Y2, state, env_alg)
+        benv3 *= norm(benv2, Inf)
+        nrm3 = PEPSKit.inner_prod(benv3, a2b2, a2b2)
+        @test benv2 ≈ benv3
+        @test nrm1 ≈ nrm2 ≈ nrm3
+    end
+    return
+end
+
+@testset "NTU bond environment ($(sym))" for sym in [U1Irrep, FermionParity]
+    Pspace = Vect[sym](0 => 1, 1 => 1)
+    V2 = Vect[sym](0 => 4, 1 => 1)
+    V3 = Vect[sym](0 => 3, 1 => 2)
+    V4 = Vect[sym](0 => 2, 1 => 2)
+    V5 = Vect[sym](0 => 2, 1 => 3)
+    Pspaces = fill(Pspace, (Nr, Nc))
+    Nspaces = [V2 V2; V4 V4]
+    Espaces = [V3 V5; V5 V3]
+
+    for state in [
+            InfinitePEPS(randn, ComplexF64, Pspaces, Nspaces, Espaces),
+            InfinitePEPO(randn, ComplexF64, Pspaces, Nspaces, Espaces),
+        ]
+        normalize!.(state.A, Inf)
+        for env_alg in (NNEnv(), NNpEnv(), NNNEnv())
+            test_ntu_env(state, env_alg)
+        end
+    end
+end

test/runtests.jl

@@ -90,6 +90,9 @@ end
         @time @safetestset "Gauge fixing" begin
             include("bondenv/benv_gaugefix.jl")
         end
+        @time @safetestset "Exact NTU bond environments" begin
+            include("bondenv/benv_ntu.jl")
+        end
         @time @safetestset "Full update bond environment" begin
             include("bondenv/benv_ctm.jl")
         end
@@ -110,6 +113,9 @@ end
         @time @safetestset "J1-J2 model at finite temperature" begin
             include("timeevol/j1j2_finiteT.jl")
         end
+        @time @safetestset "Transverse Field Ising model real-time evolution" begin
+            include("timeevol/tf_ising_realtime.jl")
+        end
     end
     if GROUP == "ALL" || GROUP == "TOOLBOX"
         @time @safetestset "Density matrix from double-layer PEPO" begin

test/timeevol/j1j2_finiteT.jl

@@ -5,9 +5,8 @@ import MPSKitModels: σˣ, σᶻ
 using PEPSKit
 
 # Benchmark energy from high-temperature expansion
-# at β = 0.3, 0.6
-# Physical Review B 86, 045139 (2012) Fig. 15-16
-bm = [-0.1235, -0.213]
+const βs = [0.2, 0.4, 0.6]
+const bm = [-0.08624893, -0.15688984, -0.21300888]
 
 function converge_env(state, χ::Int)
     trunc1 = truncrank(χ) & truncerror(; atol = 1.0e-12)
@@ -23,37 +22,36 @@ ham = j1_j2_model(
 )
 pepo0 = PEPSKit.infinite_temperature_density_matrix(ham)
 wts0 = SUWeight(pepo0)
-# 7 = 1 (spin-0) + 2 x 3 (spin-1)
-trunc_pepo = truncrank(7) & truncerror(; atol = 1.0e-12)
-check_interval = 100
-dt, nstep = 1.0e-3, 600
+dt, nstep, check_interval = 5.0e-3, 40, 40
 
-# PEPO approach
-alg = SimpleUpdate(; trunc = trunc_pepo, purified = false)
-evolver = TimeEvolver(pepo0, ham, dt, nstep, alg, wts0)
-pepo, wts, info = time_evolve(evolver; check_interval)
-env = converge_env(InfinitePartitionFunction(pepo), 16)
-energy = expectation_value(pepo, ham, env) / (Nr * Nc)
-@info "β = $(dt * nstep): tr(ρH) = $(energy)"
-@test dt * nstep ≈ info.t
-@test energy ≈ bm[2] atol = 5.0e-3
-
-# PEPS (purified PEPO) approach
-alg = SimpleUpdate(; trunc = trunc_pepo, purified = true)
-evolver = TimeEvolver(pepo0, ham, dt, nstep, alg, wts0)
-pepo, wts, info = time_evolve(evolver; check_interval)
-env = converge_env(InfinitePartitionFunction(pepo), 16)
-energy = expectation_value(pepo, ham, env) / (Nr * Nc)
-@info "β = $(dt * nstep) / 2: tr(ρH) = $(energy)"
-@test energy ≈ bm[1] atol = 5.0e-3
-
-# test BP gauge fixing for purified iPEPO
-bp_alg = BeliefPropagation(; maxiter = 100, tol = 1.0e-9)
-bp_env, = leading_boundary(BPEnv(ones, Float64, pepo), pepo, bp_alg)
-pepo, = gauge_fix(pepo, BPGauge(), bp_env)
+@testset "Simple update" begin
+    # 7 = 1 (spin-0) + 2 x 3 (spin-1)
+    trunc_pepo = truncrank(7) & truncerror(; atol = 1.0e-12)
+    alg = SimpleUpdate(; trunc = trunc_pepo, purified = true)
+    pepo, wts = deepcopy(pepo0), deepcopy(wts0)
+    for (β, bme) in zip(βs, bm)
+        t0 = β - βs[1]
+        pepo, wts, info = time_evolve(pepo, ham, dt, nstep, alg, wts; t0, check_interval)
+        # measure energy
+        env = converge_env(InfinitePEPS(pepo), 16)
+        energy = expectation_value(pepo, ham, pepo, env) / (Nr * Nc)
+        @info "β = $(info.t): ⟨ρ|H|ρ⟩ = $(energy)"
+        @test energy ≈ bme atol = 5.0e-3
+    end
+end
 
-env = converge_env(InfinitePEPS(pepo), 16)
-energy = expectation_value(pepo, ham, pepo, env) / (Nr * Nc)
-@info "β = $(dt * nstep): ⟨ρ|H|ρ⟩ = $(energy)"
-@test dt * nstep ≈ info.t
-@test energy ≈ bm[2] atol = 5.0e-3
+@testset "Neighbourhood tensor update" begin
+    trunc_pepo = truncrank(4) & truncerror(; atol = 1.0e-12)
+    opt_alg = ALSTruncation(; trunc = trunc_pepo, tol = 1.0e-10)
+    alg = NeighbourUpdate(; opt_alg, bondenv_alg = NNEnv())
+    pepo = deepcopy(pepo0)
+    for (β, bme) in zip(βs, bm)
+        t0 = β - βs[1]
+        pepo, info = time_evolve(pepo, ham, dt, nstep, alg; t0, check_interval)
+        # measure energy
+        env = converge_env(InfinitePEPS(pepo), 16)
+        energy = expectation_value(pepo, ham, pepo, env) / (Nr * Nc)
+        @info "β = $(info.t): ⟨ρ|H|ρ⟩ = $(energy)"
+        @test energy ≈ bme atol = 2.0e-2
+    end
+end

test/timeevol/tf_ising_realtime.jl

@@ -0,0 +1,90 @@
+using Test
+using TensorKit
+import MPSKitModels: S_zz, σˣ
+using PEPSKit
+using Printf
+using Accessors: @set
+
+const hc = 3.044382
+const formatter = Printf.Format("t = %.2f, ⟨σˣ⟩ = %.7e + %.7e im.")
+# real time evolution of ⟨σx⟩
+# benchmark data from Physical Review B 104, 094411 (2021)
+# Figure 6(a) calculated with D = 8 and χ = 32
+const data = [
+    # 0.01 9.9920027e-1
+    0.06 9.7274912e-1
+    0.11 9.1973182e-1
+    0.16 8.6230618e-1
+    0.21 8.1894325e-1
+    0.26 8.0003708e-1
+    0.31 8.0081082e-1
+    0.36 8.0979257e-1
+    # 0.41 8.1559623e-1
+    # 0.46 8.1541661e-1
+    # 0.51 8.1274128e-1
+]
+
+# the fully polarized state
+peps0 = InfinitePEPS(zeros, ComplexF64, ℂ^2, ℂ^1; unitcell = (2, 2))
+for t in peps0.A
+    t[1, 1, 1, 1, 1] = 1.0
+    t[2, 1, 1, 1, 1] = 1.0
+end
+lattice = collect(space(t, 1) for t in peps0.A)
+
+# Hamiltonian
+op = LocalOperator(lattice, ((1, 1),) => σˣ())
+ham = transverse_field_ising(ComplexF64, Trivial, InfiniteSquare(2, 2); J = 1.0, g = hc)
+
+# truncation strategy
+Dcut, chi = 4, 16
+trunc_peps = truncerror(; atol = 1.0e-10) & truncrank(Dcut)
+trunc_env = truncerror(; atol = 1.0e-10) & truncrank(chi)
+
+ctm_alg = SequentialCTMRG(;
+    tol = 1.0e-8, maxiter = 50, verbosity = 2,
+    trunc = trunc_env, projector_alg = :fullinfinite
+)
+
+interval = 5
+ntu_alg = NeighbourUpdate(;
+    opt_alg = FullEnvTruncation(; trunc = trunc_peps, tol = 1.0e-10),
+    bondenv_alg = NNEnv(), imaginary_time = false
+)
+
+# do one step of NTU to match benchmark data
+peps0, = time_evolve(peps0, ham, 0.01, 6, ntu_alg)
+@info "Space of `peps0[1, 1]` = $(space(peps0[1, 1]))."
+env0 = CTMRGEnv(ones, ComplexF64, peps0, ℂ^1)
+env0, = leading_boundary(env0, peps0, ctm_alg)
+# measure magnetization
+magx = expectation_value(peps0, op, env0)
+@info Printf.format(formatter, 0.06, real(magx), imag(magx))
+@test isapprox(magx, data[1, 2]; atol = 0.005)
+
+@testset "Neigborhood tensor update" begin
+    peps, env = deepcopy(peps0), deepcopy(env0)
+    count = 2
+    evolver = TimeEvolver(peps, ham, 0.01, 30, ntu_alg; t0 = 0.06)
+    spaces0 = collect(space(t) for t in peps.A)
+    for (peps, info) in evolver
+        !(evolver.state.iter % interval == 0) && continue
+        spaces = collect(space(t) for t in peps.A)
+        if spaces0 == spaces
+            env, = leading_boundary(env, peps, ctm_alg)
+        else
+            env = complex(CTMRGEnv(info.wts))
+            env, = leading_boundary(env, peps, ctm_alg)
+        end
+        # monitor the growth of env dimension
+        corner = env.corners[1, 1, 1]
+        corner_dim = dim.(space(corner, ax) for ax in 1:numind(corner))
+        @info "Dimension of env.corner[1, 1, 1] = $(corner_dim)."
+        # measure magnetization
+        magx = expectation_value(peps, op, env)
+        @info Printf.format(formatter, info.t, real(magx), imag(magx))
+        @test isapprox(magx, data[count, 2]; atol = 0.005)
+        count += 1
+        spaces0 = spaces
+    end
+end