METTS for Majorana fermion

ITensor Julia Questions
, , , ,
1

Hi,

I am trying to implement the METTS ( Phys. Rev. Lett. 102, 190601 (2009); New J. Phys. 12 055026 (2010)) for interacting Kitaev chain (Grand canonical ensemble). But with a problem on the choice/implementation of changing basis to reduce autocorrelation time, e.g., result does not match the known result for the free kitaev chain. – In particular, my concern is about the problem of fermonic string, so I have replaced the TEBD method with TDVP, while the rest part remains unclear to me.

It would be greatly appreciate if you can provide help, e.g., exemplary code, or suggestion.

Below, I attached my code,

using ITensors, ITensorMPS
using Printf
using LinearAlgebra, Random
using Random, RandomMatrices
using Plots

ITensors.op(::OpName"Ue",::SiteType"Fermion") =
[1/√2 1/√2;
1/√2 -1/√2]

ITensors.state(::StateName"f+“, ::SiteType"Fermion”) = [1/sqrt(2), 1/sqrt(2)]
ITensors.state(::StateName"f-“, ::SiteType"Fermion”) = [1/sqrt(2), -1/sqrt(2)]

function H_evo(sites, N, t_hop, Δ_pair, bf, bc)
##t_hop, Δ_pair for parameter in the kitaev chain
##bf = 0, 1 for Z2 flux
##bc = 0, 1 obc pbc,
os = OpSum()

for j in 1:N
next_site_j = mod1(j + 1, N)
δ_bd = Int(j == N)
bf_pos = (-1)^(bf * δ_bd) * (bc * δ_bd + (1 - δ_bd))

os += - t_hop * bf_pos , "Cdag", j, "C", next_site_j 
os +=  t_hop * bf_pos, "C", j, "Cdag", next_site_j
    
os +=  Δ_pair * bf_pos, "C", j, "C", next_site_j 
os +=  -Δ_pair * bf_pos, "Cdag", j, "Cdag", next_site_j

end
H = MPO(os, sites)
return H

end

“”"
Given a Vector of numbers, returns
the average and the standard error
(= the width of distribution of the numbers)
“”"
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, √abs(((avg2 - avg^2)) / N)
end

function main(; N=10, cutoff=1E-10, δτ=0.01, β=1.0, NMETTS=3000, Nwarm=10, t_hop, Δ_pair, bf, bc)

Make an array of ‘site’ indices

s = siteinds("Fermion", N; conserve_nfparity=false)

Make gates (1,2),(2,3),(3,4),…

next_site(i) = mod1(i, N)


# Arbitrary initial state
ψ = random_mps(ComplexF64, s)
Ry_gates = ops([("Ue", n) for n in 1:N], s)

P = MPO(op.("F", s))


fps = ComplexF64[]

if norm(Int(β/δτ) * δτ - β) > 1E-10
error("Time step δτ=$δτ not commensurate with β=$(β)")
end

for step in 1:(Nwarm + NMETTS)
    if step <= Nwarm
      println("Making warmup METTS number $step")
    else
      println("Making METTS number $(step-Nwarm)")
    end

    # Do the time evolution via TDVP
    H = H_evo(s, N, t_hop, Δ_pair, bf, bc)
    ψ =  tdvp(H,
                -β,
                ψ;
                time_step= -δτ,
                maxdim=100,
                cutoff=cutoff,
                normalize=true,
                reverse_step=false,
                outputlevel=0,)
    normalize!(ψ)


# Measure properties after >= Nwarm
# METTS have been made
 if step > Nwarm
   fp = inner(ψ', P, ψ)
   push!(fps, fp)
   a_fp = mean(fps)
   println("current FP = $fp, average FP = $a_fp")
 end

# Measure in X or Z basis on alternating steps
if isodd(step)
    
  ψ = apply(Ry_gates, ψ)
  samp = sample!(ψ)
  new_state = [samp[j] == 1 ? "f+" : "f-" for j in 1:N]
else

        
  samp = sample!(ψ)
  new_state = [samp[j] == 1 ? "Emp" : "Occ" for j in 1:N]

  
  
end
ψ = MPS(ComplexF64, s, new_state) #MPS(s, new_state)
normalize!(ψ)

end

return nothing
end

2

Do you have a specific question you could ask about? It’s a little hard to know how to help with the information provided.

Also do you have a colleague or supervisor who you could discuss this with? Devising a strategy for METTS of a new system is mostly about (1) time evolving systems in a correct way and (2) finding good bases for collapsing the states which ideally generate samples which change a lot from the previous ones (i.e. are ergodic updates). Most of these issues are not specific to ITensor and are about testing your code on examples or working out examples on small systems mathematically.

3

Hi, Miles.

Thanks for your reply!

I have attached my code above.

ps: I’m monitoring the expectation value of the total fermion parity as a check on the symmetry sector. In the collapsed states I see a pattern (that does not align with the free Kitaev chain): Odd collapses, \langle \text{Parity}\rangle \sim \pm 1 (looks fine); Even collapses: \langle \text{Parity}\rangle \sim 0. ==> |\langle \text{Parity}\rangle| \leq 1/2.