Does ITensor support SU(3) fermion systems?

ITensor Julia Questions
1

Hi everyone,
I’m a beginner with ITensor and I’m trying to simulate an SU(3) fermionic Hubbard model. I wrote my own custom operators, but the results are clearly wrong and I don’t know where the problem is. I’ve attached my code.
Does anyone know a better way to implement SU(3) fermions in ITensor? Any advice would be really appreciated!

using ITensors, ITensorMPS

--------------------------

function ITensors.space(::SiteType"SU3F", n::Int;
conserve_qns=false,
conserve_nf=conserve_qns,
conserve_sz=conserve_qns,
qnname_nf=“Nf”,
qnname_sz=“Sz”)
if conserve_nf && conserve_sz
# (Nf, Sz)
# |empty⟩, |f1⟩, |f2⟩, |f3⟩, |f1f2⟩, |f1f3⟩, |f2f3⟩, |f1f2f3⟩
return [
QN((qnname_nf, 0, -1), (qnname_sz, 0)) => 1 # |empty⟩
QN((qnname_nf, 1, -1), (qnname_sz, +1)) => 1 #|f1⟩
QN((qnname_nf, 1, -1), (qnname_sz, 0)) => 1 # |f2⟩
QN((qnname_nf, 1, -1), (qnname_sz, -1)) => 1 # |f3⟩
QN((qnname_nf, 2, -1), (qnname_sz, +1)) => 1 # |f1f2⟩
QN((qnname_nf, 2, -1), (qnname_sz, 0)) => 1 # |f1f3⟩
QN((qnname_nf, 2, -1), (qnname_sz, -1)) => 1 # |f2f3⟩
QN((qnname_nf, 3, -1), (qnname_sz, 0)) => 1 # |f1f2f3⟩
]
elseif conserve_nf
return [
QN(qnname_nf, 0, -1) => 1
QN(qnname_nf, 1, -1) => 3
QN(qnname_nf, 2, -1) => 3
QN(qnname_nf, 3, -1) => 1
]
end
return 8
end

--------------------------

SU3

--------------------------

ITensors.state(::StateName"empty", ::SiteType"SU3F") = [1.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0]
ITensors.state(::StateName"f1", ::SiteType"SU3F") = [0.0 1.0 0.0 0.0 0.0 0.0 0.0 0.0]
ITensors.state(::StateName"f2", ::SiteType"SU3F") = [0.0 0.0 1.0 0.0 0.0 0.0 0.0 0.0]
ITensors.state(::StateName"f3", ::SiteType"SU3F") = [0.0 0.0 0.0 1.0 0.0 0.0 0.0 0.0]
ITensors.state(::StateName"f12", ::SiteType"SU3F") = [0.0 0.0 0.0 0.0 1.0 0.0 0.0 0.0]
ITensors.state(::StateName"f13", ::SiteType"SU3F") = [0.0 0.0 0.0 0.0 0.0 1.0 0.0 0.0]
ITensors.state(::StateName"f23", ::SiteType"SU3F") = [0.0 0.0 0.0 0.0 0.0 0.0 1.0 0.0]
ITensors.state(::StateName"f123", ::SiteType"SU3F") = [0.0 0.0 0.0 0.0 0.0 0.0 0.0 1.0]

--------------------------

function ITensors.op!(Op::ITensor, ::OpName"Sz", ::SiteType"SU3F", s::Index)
Op[s’=>1, s=>1] = 0.0 # empty
Op[s’=>2, s=>2] = +1.0 # f1
Op[s’=>3, s=>3] = 0.0 # f2
Op[s’=>4, s=>4] = -1.0 # f3
Op[s’=>5, s=>5] = +1.0 # f1f2
Op[s’=>6, s=>6] = 0.0 # f1f3
Op[s’=>7, s=>7] = -1.0 # f2f3
return Op[s’=>8, s=>8] = 0.0 # f1f2f3
end

N1

function ITensors.op!(Op::ITensor, ::OpName"Nf1", ::SiteType"SU3F", s::Index)
Op[s’=>2, s=>2] = 1.0
Op[s’=>5, s=>5] = 1.0
Op[s’=>6, s=>6] = 1.0
return Op[s’=>8, s=>8] = 1.0
end

function ITensors.op!(Op::ITensor, ::OpName"Nf2", ::SiteType"SU3F", s::Index)
Op[s’=>3, s=>3] = 1.0
Op[s’=>5, s=>5] = 1.0
Op[s’=>7, s=>7] = 1.0
return Op[s’=>8, s=>8] = 1.0
end

function ITensors.op!(Op::ITensor, ::OpName"Nf3", ::SiteType"SU3F", s::Index)
Op[s’=>4, s=>4] = 1.0
Op[s’=>6, s=>6] = 1.0
Op[s’=>7, s=>7] = 1.0
return Op[s’=>8, s=>8] = 1.0
end

Ntot

function ITensors.op!(Op::ITensor, ::OpName"Ntot", ::SiteType"SU3F", s::Index)
Op[s’=>2, s=>2] = 1.0
Op[s’=>3, s=>3] = 1.0
Op[s’=>4, s=>4] = 1.0
Op[s’=>5, s=>5] = 2.0
Op[s’=>6, s=>6] = 2.0
Op[s’=>7, s=>7] = 2.0
return Op[s’=>8, s=>8] = 3.0
end

f1

function ITensors.op!(Op::ITensor, ::OpName"Cf1", ::SiteType"SU3F", s::Index)
Op[s’=>1, s=>2] = 1.0
Op[s’=>3, s=>5] = 1.0
Op[s’=>4, s=>6] = 1.0
return Op[s’=>7, s=>8] = 1.0
end

function ITensors.op!(Op::ITensor, ::OpName"Cdagf1", ::SiteType"SU3F", s::Index)
Op[s’=>2, s=>1] = 1.0
Op[s’=>5, s=>3] = 1.0
Op[s’=>6, s=>4] = 1.0
return Op[s’=>8, s=>7] = 1.0
end

f2

-------------------

function ITensors.op!(Op::ITensor, ::OpName"Cf2", ::SiteType"SU3F", s::Index)
Op[s’=>1, s=>3] = 1.0
Op[s’=>2, s=>5] = -1.0
Op[s’=>4, s=>7] = 1.0
return Op[s’=>6, s=>8] = -1.0
end

function ITensors.op!(Op::ITensor, ::OpName"Cdagf2", ::SiteType"SU3F", s::Index)
Op[s’=>3, s=>1] = 1.0
Op[s’=>5, s=>2] = -1.0
Op[s’=>7, s=>4] = 1.0
return Op[s’=>8, s=>6] = -1.0
end

-------------------

#f3

-------------------

function ITensors.op!(Op::ITensor, ::OpName"Cf3", ::SiteType"SU3F", s::Index)
Op[s’=>1, s=>4] = 1.0
Op[s’=>2, s=>6] = -1.0
Op[s’=>3, s=>7] = -1.0
return Op[s’=>5, s=>8] = 1.0
end

function ITensors.op!(Op::ITensor, ::OpName"Cdagf3", ::SiteType"SU3F", s::Index)
Op[s’=>4, s=>1] = 1.0
Op[s’=>6, s=>2] = -1.0
Op[s’=>7, s=>3] = -1.0
return Op[s’=>8, s=>5] = 1.0
end

--------------------------

--------------------------

ITensors.has_fermion_string(::OpName"Cf1", ::SiteType"SU3F") = true
ITensors.has_fermion_string(::OpName"Cdagf1", ::SiteType"SU3F") = true
ITensors.has_fermion_string(::OpName"Cf2", ::SiteType"SU3F") = true
ITensors.has_fermion_string(::OpName"Cdagf2", ::SiteType"SU3F") = true
ITensors.has_fermion_string(::OpName"Cf3", ::SiteType"SU3F") = true
ITensors.has_fermion_string(::OpName"Cdagf3", ::SiteType"SU3F") = true

N = 6
sites = siteinds(“SU3F”, N; conserve_nf=true,conserve_sz=false)

--------------------------

t=1.0 # hopping
U=4.0 # on-site interaction

ampo = AutoMPO()
for j in 1:N-1
# hopping: f1
add!(ampo, -t, “Cdagf1”, j, “Cf1”, j+1)
add!(ampo, -t, “Cdagf1”, j+1, “Cf1”, j)
# hopping: f2
add!(ampo, -t, “Cdagf2”, j, “Cf2”, j+1)
add!(ampo, -t, “Cdagf2”, j+1, “Cf2”, j)
# hopping: f3
add!(ampo, -t, “Cdagf3”, j, “Cf3”, j+1)
add!(ampo, -t, “Cdagf3”, j+1, “Cf3”, j)
end

on-site interaction:

for j in 1:N
add!(ampo, U, “Nf1”, j, “Nf2”, j)
add!(ampo, U, “Nf2”, j, “Nf3”, j)
add!(ampo, U, “Nf1”, j, “Nf3”, j)
end

H = MPO(ampo, sites)

--------------------------

--------------------------

state = [“f123”, “f13”, “empty”, “empty”, “empty”, “f123”]
psi0 = productMPS(sites, state)

--------------------------

4. DMRG

--------------------------

sweeps = Sweeps(30)
maxdim!(sweeps, 50,100,200,400,400)
cutoff!(sweeps, 1E-6)

energy, psi = dmrg(H, psi0, sweeps;)

println(“Ground state energy = $energy”)

N1_vals = expect(psi, “Nf1”)
println(N1_vals)
println(sum(N1_vals))

N2_vals = expect(psi, “Nf2”)
println(N2_vals)
println(sum(N2_vals))

N3_vals = expect(psi, “Nf3”)
println(N3_vals)
println(sum(N3_vals))

println(sum(N1_vals)+sum(N2_vals)+sum(N3_vals))

2

This all looks reasonable at a glance. Sorry to hear it’s not working. I was specifically looking to see if you overloaded the has_fermion_string method for your “C” and “Cdag” operators and you did. That is a crucial step.

Some questions back to you:

  1. When you say SU(3), do you mean just an Abelian subgroup of SU(3)? I.e. are you needing to conserve the entire non-Abelian SU(3) symmetry to get the kind of results you’re expecting or is an Abelian subgroup enough? The reason is that ITensors.jl and ITensorMPS.jl just supports Abelian symmetries at this time.
  2. Your system may be challenging to converge. Have you tried including the noise parameter in your call to DMRG?

As a minor thing, I’d also recommend the newer syntax for DMRG, where you pass nsweeps, maxdim, and cutoff (and noise) as keyword arguments:
https://docs.itensor.org/ITensorMPS/stable/tutorials/DMRG.html

3

First, thank you very much for your reply. I believe this issue indeed stems from ITensor’s inability to handle non-Abelian symmetries.

  1. I need to conserve the full non-Abelian SU(3) symmetry because I require each flavor to be completely equivalent and the operators “C” and “Cdag” to anticommute pairwise. While I manually inserted the Jordan–Wigner signs, I still could not achieve the expected behavior. Currently, ITensor only enforces conservation of the total particle number and the particle number for each flavor, but the real-space charge density distributions of the three flavors are still different.
  2. I also tried including noise in the DMRG procedure. Even for very small system sizes with many sweeps, the expected results were not obtained, indicating that this is not a limitation of the DMRG algorithm itself.