# METTS for Majorana fermion

**URL:** <https://itensor.discourse.group/t/metts-for-majorana-fermion/2426>\
**Category:** ITensor Julia Questions\
**Tags:** julia, mps, mpo, fermions, tdvp\
**Created:** [June 21, 2025, 11:07am UTC](https://itensor.discourse.group/t/metts-for-majorana-fermion/2426 "2025-06-21T11:07:54Z")\
**Posts on this page:** 3\
**Page:** 1

<div class="post-metadata">

**Author:** ![prange](https://avatars.discourse-cdn.com/v4/letter/p/c4cdca/32.png) [@prange](https://itensor.discourse.group/u/prange)\
**Post date:** [June 21, 2025, 11:07am UTC](https://itensor.discourse.group/t/metts-for-majorana-fermion/2426/1 "2025-06-21T11:07:54Z")

</div>

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

---

<div class="post-metadata">

**Author:** ![miles](https://yyz2.discourse-cdn.com/free1/user_avatar/itensor.discourse.group/miles/32/6_2.png) [@miles](https://itensor.discourse.group/u/miles)\
**Post date:** [June 21, 2025, 10:41pm UTC](https://itensor.discourse.group/t/metts-for-majorana-fermion/2426/2 "2025-06-21T22:41:10Z")

</div>

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.

---

<div class="post-metadata">

**Author:** ![prange](https://avatars.discourse-cdn.com/v4/letter/p/c4cdca/32.png) [@prange](https://itensor.discourse.group/u/prange)\
**Post date:** [June 22, 2025, 11:13am UTC](https://itensor.discourse.group/t/metts-for-majorana-fermion/2426/3 "2025-06-22T11:13:53Z")

</div>

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.