METTS seems not to work

DMRG and Numerical Methods
1

Hello everybody,

I tried to code up METTS to estimate thermal expectation values for the Ising mode, by modifying the example code given in the ITensorMPS repository. But I couldn’t get my code to agree with exact diagonalization results. I attach my code here. Could someone take a look and let me know where I screw up?

Thanks!
Xiaojun

using ITensors, ITensorMPS
using Printf

N = 8
Nmid = ceil(Int, N/2)
Nm1 = N-1
J = -1
hx = -1
hz = -1
dtau = 0.01
βmax = 0.1
βmax_half = βmax/2
bond = 150
cutoff_set = 1e-12

function avg_err(v::Vector)
    N = length(v)
    avg = v[1] / N
    avg2 = v[1]^2 / N
    for j in 2:N
        avg += v[j] / N
        avg2 += v[j]^2 / N
    end
    return avg, √((avg2 - avg^2) / N)
end

function ITensors.op(::OpName"exptauZZ", ::SiteType"Qubit", s1::Index, s2::Index; tau)
    Hzz = J * op("Z", s1) * op("Z", s2)
    return exp(tau * Hzz)
end

function ITensors.op(::OpName"exptauZ", ::SiteType"Qubit", s::Index; tau)
    Hz = hz * op("Z", s)
    return exp(tau * Hz)
end

function ITensors.op(::OpName"exptauX", ::SiteType"Qubit", s::Index; tau)
    Hx = hx * op("X", s)
    return exp(tau * Hx)
end

s = siteinds("Qubit", N; conserve_qns = false)
### Make gates for thermal state cooling
gates_expX = [("exptauX", (n,), (tau = -dtau/2,)) for n in 1:N]
gates_expZ = [("exptauZ", (n,), (tau = -dtau/2,)) for n in 1:N]
gates_expZZ = [("exptauZZ", (2n, 2n + 1), (tau = -dtau/2,)) for n in 1:(div(N,2))-1]
append!(gates_expZZ, [("exptauZZ", (N, 1), (tau = -dtau/2,))])
append!(gates_expZZ, [("exptauZZ", (2n-1, 2n), (tau = -dtau/2,)) for n in 1:(div(N,2))])
gates_T = ops(gates_expX, s)
append!(gates_T, ops(gates_expZ, s))
append!(gates_T, ops(gates_expZZ, s))
# Include gates in reverse order to complete Trotter formula
append!(gates_T, reverse(gates_T))

### Hamiltonian
function Hamil(N_)
    os = OpSum()
    for n in 1:N_
        if n == N_
            os += J, "Z", N_, "Z", 1
            os += hx, "X", N_
            os += hz, "Z", N_
        else
            os += J, "Z", n, "Z", n+1
            os += hx, "X", n
            os += hz, "Z", n
        end
    end
    return os
end
H = MPO(Hamil(N), s)

### Make y-rotation gates to use in METTS collapses
Ry_gates = ops([("Ry", n, (θ = π / 2,)) for n in 1:N], s)

Nwarm = 50
Nmetts = 1000
Ntot = Nwarm+Nmetts

function run_metts(gates_T, Ry_gates, s, cutoff_set, bond)
    energies = Float64[]
    init_state = rand(["Z+", "Z-"], N)
    psi = MPS(s, init_state)
    for step in 1:Ntot
        for β in dtau:dtau:βmax_half
            psi = apply(gates_T, psi; cutoff=cutoff_set, maxdim = bond)
            normalize!(psi)
        end
        if step > Nwarm
            energy = inner(psi', H, psi)
            push!(energies, energy)
            println(avg_err(energies))
        end
        # Measure in X or Z basis on alternating steps
        if step % 2 == 1
            psi = apply(Ry_gates, psi)
            samp = sample!(psi)
            new_state = [samp[j] == 1 ? "X+" : "X-" for j in 1:N]
        else
            samp = sample!(psi)
            new_state = [samp[j] == 1 ? "Z+" : "Z-" for j in 1:N]
        end
        psi = MPS(s, new_state)
    end
end