Help: Custom Majorana Operators

I attempted to compute the ground-state energy of a critical Kitaev chain using DMRG with the “Fermion” already implemented in ITensor. In the process, I custom-defined the Majorana operators β and γ. The code is as follows:

L = 20  
sites = siteinds("Fermion", L; conserve_qns=false)

function ITensors.op(::OpName"Beta", ::SiteType"Fermion", s::Index)
  c = op("C", s)
  cdag = op("Cdag", s)
  return c + cdag

function ITensors.op(::OpName"Gamma", ::SiteType"Fermion", s::Index)
  c = op("C", s)
  cdag = op("Cdag", s)
  return -im * (c - cdag)

os = OpSum()
for j=1:L
  os += -1* im, "Beta", j, "Gamma", j
for j=1:L-1
  os += 2*im, "Beta", j+1, "Gamma", j

H = MPO(os, sites)
psi0 = randomMPS(sites)

# Numerical accuracy parameters
nsweeps = 50
maxdim = [10, 20, 100, 100, 200, 300,1000]
cutoff = 1E-10

# Find the ground state
energy, psi = dmrg(H, psi0; nsweeps, maxdim, cutoff,outputlevel=1)

println("ground-state energy: ", energy)

However, when I wrote the Hamiltonian using the original “C” and “Cdag”, the efficiency of the DMRG was much better, and it required significantly fewer bond dimensions. Additionally, the energy appears to be incorrect. Why might this be?

It’s an interesting question, but the short answer is “I don’t know” until you explore a bit more what is happening. Some possibilities:

(1) is your Hamiltonian Hermitian? At a glance, since Beta and Gamma look Hermitian I would say overall your Hamiltonian looks to be anti-Hermitian. Sometimes when people accidentally put non-Hermitian Hamiltonians in our DMRG code the result is very unstable behavior accompanied by large bond dimensions. (We should probably add an automated check for that.)

(2) assuming your Hamiltonian is Hermitian, then another possibility is that in this Majorana basis it is finding a “Schrodinger cat” state which is a superposition of the ground states one would find in the regular electron basis. But please start by checking #1 above.

Also, when you say the energy appears to be incorrect, you mentioned two different calculations. In which one did you find the energy to be incorrect?

Dear miles-san,

Thank you for your response. Theoretically, since β and γ are anticommutative, if implemented correctly, the Hamiltonian should be Hermitian. Indeed, the energy output by DMRG is real. However, I would like to verify its Hermitian nature just to be sure. Would this be a valid way to check?

function is_hermitian(M::MPO, tol=1e-10)
    Mdagger = adjoint(M)

    if M==Mdagger
        return true
        return false

the issue with incorrect energy occurs with the implementation using Majorana operators. For reference, here is an implementation example using “C” and “Cdag” (the correct one):

sites = siteinds("Fermion", L; conserve_qns=false)
os = OpSum()
for j = 1:L-1
    os += 1, "C", j+1, "C", j
    os += -1, "Cdag", j+1, "Cdag", j
    os += 1, "Cdag", j+1, "C", j
    os += -1, "C", j+1, "Cdag", j

H = MPO(os, sites)
psi0 = randomMPS(sites)

# Parameters that determine the numerical accuracy
nsweeps = 50
maxdim = [10, 20, 100, 100, 200, 300]
cutoff = [1E-10]

# Setting outputlevel to 0 to suppress detailed output
energy, psi = dmrg(H, psi0; nsweeps, maxdim, cutoff,outputlevel=1)
println("ground state energy: ", energy)

Notably, when the ground state is a product state and the energy of the ground state becomes -L , then for
H=\sum_j i \beta_j \gamma_j
or H=\sum_j i \beta_j+1 \gamma_j
the former proceeds correctly with DMRG, but the latter was problematic. Therefore, there might be an issue with implementing the anticommutativity between different sites.

As a note about formatting your code in the forum, please use three “back ticks” which are the tick marks that lean to the left. I have edited your post above to use these so it is formatted correctly.

I’ll think about your issue and get back to you soon.

I just remembered something pretty fundamental about this system that you might have missed. For fermionic operators to be handled correctly by the OpSum system, you have to overload the function has_fermion_string. For example:

function ITensors.has_fermion_string(on::OpName"cdag", st::SiteType"Fermion")
  return has_fermion_string(alias(on), st)

is the overload for the “Fermion” “cdag” operator.
Here is a link to some code inside the library where more of these are defined:

Please try defining these for your custom operators if they are anticommuting & let me know if it works after that or not.

Thank you for the guidance on implementing the function has_fermion_string for the fermionic operators. Following your instructions, I successfully integrated the Majorana operators, and they are now functioning correctly in the OpSum system. I really appreciate your help with this!

Glad it is working now!

1 Like