PEPSKit.jl/test/timeevol/tf_ising_finiteT.jl at 885faaf16a9b2938df7438b7d345832f8fe3c534 · QuantumKitHub/PEPSKit.jl · GitHub

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
using Test
using LinearAlgebra
using TensorKit
import MPSKitModels: σˣ, σᶻ
using PEPSKit

# Benchmark data of [σx, σz] from HOTRG
# Physical Review B 86, 045139 (2012) Fig. 15-16
bm_β = [0.5632, 0.0]
bm_2β = [0.5297, 0.8265]

function converge_env(state, χ::Int)
    env0 = initialize_ctmrg_environment(state, ProductStateInitialization())
    trunc = truncrank(χ) & truncerror(; atol = 1.0e-12)
    env, = leading_boundary(env0, state; alg = :sequential, trunc, tol = 1.0e-10)
    return env
end

function measure_mag(pepo::InfinitePEPO, env::CTMRGEnv; purified::Bool = false)
    r, c = 1, 1
    lattice = physicalspace(pepo)
    Mx = LocalOperator(lattice, ((r, c),) => σˣ(Float64, Trivial))
    Mz = LocalOperator(lattice, ((r, c),) => σᶻ(Float64, Trivial))
    if purified
        magx = expectation_value(pepo, Mx, pepo, env)
        magz = expectation_value(pepo, Mz, pepo, env)
    else
        magx = expectation_value(pepo, Mx, env)
        magz = expectation_value(pepo, Mz, env)
    end
    return [magx, magz]
end

Nr, Nc = 2, 2
ham = transverse_field_ising(Float64, Trivial, InfiniteSquare(Nr, Nc); J = 1.0, g = 2.0)
pepo0 = PEPSKit.infinite_temperature_density_matrix(ham)
wts0 = SUWeight(pepo0)

trunc_pepo = truncrank(8) & truncerror(; atol = 1.0e-12)

dt, nstep = 1.0e-3, 400
β = dt * nstep

# when g = 2, β = 0.4 and 2β = 0.8 belong to two phases (without and with nonzero σᶻ)
@testset "Finite-T SU (force_mpo = $(force_mpo))" for force_mpo in (false, true)
    # use second order Trotter decomposition
    symmetrize_gates = true
    bipartite = true

    # PEPO approach: results at β, or T = 2.5
    alg = SimpleUpdate(; trunc = trunc_pepo, purified = false, bipartite, force_mpo)
    pepo, wts, info = time_evolve(pepo0, ham, dt, nstep, alg, wts0; symmetrize_gates)

    ## BP gauge fixing
    bp_alg = BeliefPropagation(; maxiter = 100, tol = 1.0e-9, bipartite)
    bp_env, = leading_boundary(BPEnv(ones, Float64, pepo), pepo, bp_alg)
    pepo, = gauge_fix(pepo, BPGauge(), bp_env)

    env = converge_env(InfinitePartitionFunction(pepo), 16)
    result_β = measure_mag(pepo, env)
    @info "tr(σ(x,z)ρ) at T = $(1 / β): $(result_β)."
    @test β ≈ info.t
    @test isapprox(abs.(result_β), bm_β, rtol = 1.0e-2)

    # continue to get results at 2β, or T = 1.25
    pepo, wts, info = time_evolve(pepo, ham, dt, nstep, alg, wts; t0 = β, symmetrize_gates)
    env = converge_env(InfinitePartitionFunction(pepo), 16)
    result_2β = measure_mag(pepo, env)
    @info "tr(σ(x,z)ρ) at T = $(1 / (2β)): $(result_2β)."
    @test 2 * β ≈ info.t
    @test isapprox(abs.(result_2β), bm_2β, rtol = 1.0e-4)

    # Purification approach: results at 2β, or T = 1.25
    alg = SimpleUpdate(; trunc = trunc_pepo, purified = true, bipartite, force_mpo)
    pepo, wts, info = time_evolve(pepo0, ham, dt, 2 * nstep, alg, wts0; symmetrize_gates)
    env = converge_env(InfinitePEPS(pepo), 8)
    result_2β′ = measure_mag(pepo, env; purified = true)
    @info "⟨ρ|σ(x,z)|ρ⟩ at T = $(1 / (2β)): $(result_2β′)."
    @test 2 * β ≈ info.t
    @test isapprox(abs.(result_2β′), bm_2β, rtol = 1.0e-2)
end