# MPO for Anyon Hubbard model

**URL:** <https://itensor.discourse.group/t/mpo-for-anyon-hubbard-model/2275>\
**Category:** DMRG and Numerical Methods\
**Tags:** julia, dmrg, mpo\
**Created:** [March 4, 2025, 6:11am UTC](https://itensor.discourse.group/t/mpo-for-anyon-hubbard-model/2275 "2025-03-04T06:11:51Z")\
**Posts on this page:** 13\
**Page:** 1

<div class="post-metadata">

**Author:** ![Ashinxinl](https://yyz2.discourse-cdn.com/free1/user_avatar/itensor.discourse.group/ashinxinl/32/265_2.png) [@Ashinxinl](https://itensor.discourse.group/u/Ashinxinl)\
**Post date:** [March 4, 2025, 6:11am UTC](https://itensor.discourse.group/t/mpo-for-anyon-hubbard-model/2275/1 "2025-03-04T06:11:51Z")

</div>

Hello everyone,  
I am trying to reproduce the results from a paper on the Anyon Hubbard model.

 ![image](https://global.discourse-cdn.com/free1/uploads/itensor/original/1X/e632599575f4c7a34e401d185fc87df40cead039.png)  
(ImageSource:DOI: 10.1103/PhysRevB.97.115126)

After defining the operators, I encountered an issue: when summing terms in the Hamiltonian using the following expressions:  
`os += -J,"adag",j,"a",j+1,"ExpIθn",j`  
`os += -J,"adag",j+1,"a",j,"ExpIθndag",j`  
my program identifies the Hamiltonian as a non-Hermitian Hamiltonian. However, in theory, the Hamiltonian should be Hermitian.

 ![image](https://global.discourse-cdn.com/free1/uploads/itensor/original/1X/bc2574b554889595da2549dbee21df0bd8c7b22a.png)

Due to my limited programming experience, I couldn’t find a solution on the forum. Additionally, I am not sure if my `ExpIθn` function is implemented correctly.  
Any suggestions or guidance would be greatly appreciated! Thank you all for your help.  
this is a code demo:

```julia
using ITensors, ITensorMPS
let
function ITensors.op(::OpName"ExpIθn", ::SiteType"Qudit", d::Int)
    mat = zeros(ComplexF64, d, d)
    for k in 1:d
        mat[k, k] = exp(im *1*pi * (k - 1))
    end
    return mat
end
function ITensors.op(::OpName"ExpIθndag", ::SiteType"Qudit", d::Int)
  mat = zeros(ComplexF64, d, d)
  for k in 1:d
      mat[k, k] = exp(im *1*pi * (k - 1))
  end
  return mat'
end

#----------------------------parameter---------------------
N=60
L_list = collect(Int64, 1:1:N)
U=10
J=1.0
V=10
alpha=1/3
delta=2*pi/3
#----------------------------parameter dmrg----------------
    maxdim = [10,10,20,20,30,30,40,40,80,80,200,200,400,400,600,600]
    nsweeps=40
    cutoff = [1E-9]
    noise=[1E-6,1E-6,1E-6,1E-7,1E-7,1E-7,1E-7,1E-7,1E-7,1E-8,1E-8,1E-8,1E-9,1E-9,1E-9,1E-9,1E-9,1E-9,1E-10,1E-10,1E-10,1E-11,1E-11,1E-11,1E-11,1E-11,1E-11,1E-11,1E-11,1E-11,1E-11,1E-11,1E-11,1E-11,1E-11,1E-11,1E-11,0.0]
    etol=1E-7
#---------------------------- dmrg----------------
    state = [n % 3 == 0 ? "1" : "0" for n=1:N]
    sites = siteinds("Boson",N,conserve_qns=true,dim=5)
    psi0 = random_mps(sites,state)
    @show flux(psi0)
    os = OpSum()
    for j=1:N-1
        os += -J,"adag",j,"a",j+1,"ExpIθn",j
        os += -J,"adag",j+1,"a",j,"ExpIθndag",j
    end
    for j=1:N
        os += U/2,"n",j, "n",j
        os -= U/2,"n",j  
    end
    for j=1:N
        os += V*cos(2*pi*alpha*j+delta),"n",j
    end
#pbc
    os += -J,"adag",N,"a",1,"ExpIθn",N
    os += -J,"ExpIθndag",N,"adag",1,"a",N
    H = MPO(os,sites)
#----------------------------strat calculate---------------
    energy,psi = dmrg(H,psi0;nsweeps,noise,maxdim,cutoff)
end

```

---

<div class="post-metadata">

**Author:** ![VinceNeede](https://yyz2.discourse-cdn.com/free1/user_avatar/itensor.discourse.group/vinceneede/32/731_2.png) [@VinceNeede](https://itensor.discourse.group/u/VinceNeede)\
**Post date:** [March 4, 2025, 10:40am UTC](https://itensor.discourse.group/t/mpo-for-anyon-hubbard-model/2275/2 "2025-03-04T10:40:50Z")

</div>

Hi @Ashinxinl  
On the top of my head, it looks like there is an error on how you are writing the hermitian conjugate part, in particular that should be

\text{e}^{-i\theta n\_j} b\_{j+1}^\dagger b\_j

The term \text{e}^{-i\theta n\_j} commutes with b\_{j+1}^\dagger since they act on different sites, but it does not commute with b\_j. So try to correct the line

```auto
        os += -J,"adag",j+1,"a",j,"ExpIθndag",j

```

to

```auto
        os += -J,"ExpIθndag",j,"adag",j+1,"a",j

```

and let us know if the problem persists.

EDIT: I found some time to test it, strangely enough I didn’t get any warning, but by explicitly testing if the MPO was hermitian, it was not. By changing the line as said previously, the mpo is now hermitian so your code should now work

---

<div class="post-metadata">

**Author:** ![VinceNeede](https://yyz2.discourse-cdn.com/free1/user_avatar/itensor.discourse.group/vinceneede/32/731_2.png) [@VinceNeede](https://itensor.discourse.group/u/VinceNeede)\
**Post date:** [March 4, 2025, 12:05pm UTC](https://itensor.discourse.group/t/mpo-for-anyon-hubbard-model/2275/3 "2025-03-04T12:05:03Z")

</div>

I’m a bit concerned about the absence of the warning in my case. Can you please check what versions do you have of `ITensorMPS` and of `KrylovKit`?

---

<div class="post-metadata">

**Author:** ![Ashinxinl](https://yyz2.discourse-cdn.com/free1/user_avatar/itensor.discourse.group/ashinxinl/32/265_2.png) [@Ashinxinl](https://itensor.discourse.group/u/Ashinxinl)\
**Post date:** [March 4, 2025, 12:43pm UTC](https://itensor.discourse.group/t/mpo-for-anyon-hubbard-model/2275/4 "2025-03-04T12:43:00Z")

</div>

Hi VinceNeede  
Apologies for the delayed response. Your suggestion was excellent, and my code is now working correctly. However, I don’t quite understand why the order of operators in the `opsum()` function matters. As far as I know, such an issue usually only arises in the case of fermions. Could this be due to the **Jordan-Wigner transformation**?

For reference, my package versions are:

- **ITensorMPS** : v0.2.6
- **KrylovKit** : v0.8.3

In my original code, warnings started appearing **after 16 sweeps**.

 ![image](https://global.discourse-cdn.com/free1/uploads/itensor/original/1X/dcc2c152b4eb69fa6f0d9b580d3b2acd02fe5a9b.png)

---

<div class="post-metadata">

**Author:** ![VinceNeede](https://yyz2.discourse-cdn.com/free1/user_avatar/itensor.discourse.group/vinceneede/32/731_2.png) [@VinceNeede](https://itensor.discourse.group/u/VinceNeede)\
**Post date:** [March 4, 2025, 1:16pm UTC](https://itensor.discourse.group/t/mpo-for-anyon-hubbard-model/2275/5 "2025-03-04T13:16:59Z")

</div>

@Ashinxinl It isn’t related to Jorgan Wigner transformations at all, it’s just how creation and annihilation operators are defined for bosons, i.e. such that

[a,a^\dagger]=1 

So the order matters, cause a^\dagger a \neq a a^\dagger but a a^\dagger=1+a^\dagger a.  
The difference in the ordering that you are mentioning is related to operators that acts on different sites, that is because the commutation rule for bosons on a lattice is

[a\_i,a^\dagger\_j]=\delta\_{i,\,j} 

i.e. it is one if they act on the same site, otherwise it is zero, that means you can write either a\_j a\_i^\dagger or a\_i^\dagger a\_j as long as i\neq j. For fermions it is different cause on different sites they anticommute, therefore you would have f\_i^\dagger f\_j=-f\_j f\_i^\dagger for i\neq j, but if they are acting on the same site, order matters both for fermions and for bosons cause the commutator is not zero.

---

<div class="post-metadata">

**Author:** ![VinceNeede](https://yyz2.discourse-cdn.com/free1/user_avatar/itensor.discourse.group/vinceneede/32/731_2.png) [@VinceNeede](https://itensor.discourse.group/u/VinceNeede)\
**Post date:** [March 4, 2025, 4:06pm UTC](https://itensor.discourse.group/t/mpo-for-anyon-hubbard-model/2275/6 "2025-03-04T16:06:37Z")

</div>

Anyway, I would recommend you to upgrade `julia` and the other packages to the latest versions.  
You can download the last julia release using `juliaup` if you already have it installed, otherwise you can download it using the instructions on the official Julia page [Download Julia](https://julialang.org/downloads/) .  
For the packages it should be easier, once upgraded julia, from the REPL just digit the following commands:

```julia
julia>]
pkg> up

```

---

<div class="post-metadata">

**Author:** ![Ashinxinl](https://yyz2.discourse-cdn.com/free1/user_avatar/itensor.discourse.group/ashinxinl/32/265_2.png) [@Ashinxinl](https://itensor.discourse.group/u/Ashinxinl)\
**Post date:** [March 5, 2025, 8:26am UTC](https://itensor.discourse.group/t/mpo-for-anyon-hubbard-model/2275/7 "2025-03-05T08:26:51Z")

</div>

@ [VinceNeede](https://itensor.discourse.group/u/VinceNeede)Thank you for your patience and detailed response. I have already updated my version.  
But When attempting to reproduce the image from Figure 2 in the paper,

 ![image](https://global.discourse-cdn.com/free1/uploads/itensor/original/1X/af60fd05570945c87e24fc4abcf8e9389630c58e.png)  
(ImageSource:DOI: 10.1103/PhysRevB.97.115126)  
I found that the results I obtained for ⟨ψₙ₊₁|nj|ψₙ₊₁⟩ - ⟨ψₙ|nj|ψₙ⟩ did not match those in the figure when setting θ to π or any other angle.  
⟨ψₙ₊₁|nj|ψₙ₊₁⟩  
 ![image](https://global.discourse-cdn.com/free1/uploads/itensor/original/1X/c441d32a4a1bd4af43864d4aa5c2459f2029e97b.png)  
⟨ψₙ|nj|ψₙ⟩  
 ![image](https://global.discourse-cdn.com/free1/uploads/itensor/original/1X/34737b143346312dc512a80166c9f45f652aabc3.png)  
Δnj  
 ![image](https://global.discourse-cdn.com/free1/uploads/itensor/original/1X/f8811047ee132409a7f661a77d071703a1bcfa7b.png)  
My conditions are consistent with those in the paper, so I suspect that the issue might still lie in the construction of the MPO.

---

<div class="post-metadata">

**Author:** ![VinceNeede](https://yyz2.discourse-cdn.com/free1/user_avatar/itensor.discourse.group/vinceneede/32/731_2.png) [@VinceNeede](https://itensor.discourse.group/u/VinceNeede)\
**Post date:** [March 5, 2025, 8:37am UTC](https://itensor.discourse.group/t/mpo-for-anyon-hubbard-model/2275/8 "2025-03-05T08:37:58Z")

</div>

The graph you’re showing is obtained using Open Boundary Condition (OBC), while you are using PBC, could that be the problem? Otherwise, can you show how you are taking the expectation values?

---

<div class="post-metadata">

**Author:** ![Ashinxinl](https://yyz2.discourse-cdn.com/free1/user_avatar/itensor.discourse.group/ashinxinl/32/265_2.png) [@Ashinxinl](https://itensor.discourse.group/u/Ashinxinl)\
**Post date:** [March 5, 2025, 8:52am UTC](https://itensor.discourse.group/t/mpo-for-anyon-hubbard-model/2275/9 "2025-03-05T08:52:39Z")

</div>

Yes, I am using OBC conditions, and I obtained the expectation values using the following code.  
 ![image](https://global.discourse-cdn.com/free1/uploads/itensor/original/1X/ad67735ff4c25317bed11e9f7cddf3856b06bbd5.png)

---

<div class="post-metadata">

**Author:** ![Ashinxinl](https://yyz2.discourse-cdn.com/free1/user_avatar/itensor.discourse.group/ashinxinl/32/265_2.png) [@Ashinxinl](https://itensor.discourse.group/u/Ashinxinl)\
**Post date:** [March 5, 2025, 8:55am UTC](https://itensor.discourse.group/t/mpo-for-anyon-hubbard-model/2275/10 "2025-03-05T08:55:39Z")

</div>

@ [VinceNeede](https://itensor.discourse.group/t/mpo-for-anyon-hubbard-model/2275/8)This is my whole code:

```auto
using ITensors, ITensorMPS
using CairoMakie

let

function ITensors.op(::OpName"ExpIθn", ::SiteType"Qudit", d::Int)
    mat = zeros(ComplexF64, d, d)
    for k in 1:d
        mat[k, k] = exp(im *pi * (k - 1))
    end
    return mat
end
function ITensors.op(::OpName"ExpIθndag", ::SiteType"Qudit", d::Int)
  mat = zeros(ComplexF64, d, d)
  for k in 1:d
      mat[k, k] = exp(im *pi * (k - 1))
  end
  return mat'
end

#----------------------------parameter---------------------
N=90
L_list = collect(Int64, 1:1:N)
f=30
U=10
J=1
V=10
alpha=1/3
delta=2*pi/3
#----------------------------parameter dmrg----------------
    maxdim = [10,10,20,20,30,30,40,40,80,80,200,200,400,400,600,600]
    nsweeps=60
    cutoff = [1E-10]
    noise=[1E-6,1E-6,1E-6,1E-7,1E-7,1E-7,1E-7,1E-7,1E-7,1E-8,1E-8,1E-8,1E-9,1E-9,1E-9,1E-9,1E-9,1E-9,1E-10,1E-10,1E-10,1E-11,1E-11,1E-11,1E-11,1E-11,1E-11,1E-11,1E-11,1E-11,1E-11,1E-11,1E-11,1E-11,1E-11,1E-11,1E-11,0.0]
    etol=1E-7
#----------------------------parameter dmrg----------------
    state = [n <= f ? "1" : "0" for n=1:N]
   # state[Int64(N/3+1)] = "1" 
    sites = siteinds("Boson",N,conserve_qns=true,dim=5)
    psi0 = random_mps(sites,state)
    @show flux(psi0)
    os = OpSum()
    for j=1:N-1
        os += -J,"adag",j,"a",j+1,"ExpIθn",j
        os += -J,"ExpIθndag",j,"adag",j+1,"a",j
    end
    for j=1:N
        os += U/2,"n",j, "n",j
        os += -U/2,"n",j  
    end

    for j=1:N
        os += V*cos(2*pi*alpha*j+delta),"n",j
    end
    #pbc
   # os += -J,"adag",N,"a",1,"ExpIθn",N
   # os += -J,"ExpIθndag",N,"adag",1,"a",N
    H = MPO(os,sites)

    H = MPO(os,sites)
#----------------------------strat calculate---------------
    energy,psi = dmrg(H,psi0;nsweeps,eigsolve_krylovdim=6,noise,maxdim,cutoff)

density= expect(psi, "n")
print(density)

f1 = Figure(backgroundcolor = :white) 
ax = Axis(f1[1,1],#xticks = [1,L/2,L],#yticks = [-0.25,0,0.25],
#title = "L=($L)N=($(L/2))h=($h)beta=($beta)V=($V)U=($U)t=($t)",
title = "alpha=(1/3)delta=(2*pi/3)L=($N)N=($(f))V=($V)U=($U)_θ=pi",
xlabel = "L",
)
#ylims!(ax,0.0, 0.7)
resize_to_layout!(f1)
scatterlines!(L_list,density,label="ni")
axislegend(position = :rb)
display(f1)

  end

```

---

<div class="post-metadata">

**Author:** ![VinceNeede](https://yyz2.discourse-cdn.com/free1/user_avatar/itensor.discourse.group/vinceneede/32/731_2.png) [@VinceNeede](https://itensor.discourse.group/u/VinceNeede)\
**Post date:** [March 5, 2025, 3:47pm UTC](https://itensor.discourse.group/t/mpo-for-anyon-hubbard-model/2275/11 "2025-03-05T15:47:47Z")

</div>

@Ashinxinl have you actually checked whether the dmrg converged? i.e. incrementing the number of sweeps or the maxdim does not change the result? Cause this is the result I’ve obtained:

 ![display](https://global.discourse-cdn.com/free1/uploads/itensor/original/1X/f4e1fe44bb4dce57238f69e94f4089bd0569d3da.png)

The mode is localized on the last site instead of the first one but it shouldn’t really matter, what is important is that it is localized on the extrema.

Anyway, I took the freedom to rewrite your code in a more readable way. Please take a look at ([Style Guide · The Julia Language](https://docs.julialang.org/en/v1/manual/style-guide/)), in particular use functions when possible, add logging messages, and use variable names that actually represent what the variable means (i.e. try to avoid naming variables that are used for the whole program with a single letter, like `N`, `L` or `f`)

> **Rewriting of the code**
>
> ```auto
> using ITensors, ITensorMPS
> using LinearAlgebra
> using CairoMakie
> 
> let
> function ITensors.op(::OpName"ExpIθn", ::SiteType"Boson", dim::Int64; kwargs...)
> return Diagonal([exp(im * pi * (k - 1)) for k in 1:dim])
> end
> function ITensors.op(::OpName"ExpIθndag", ::SiteType"Boson", dim::Int64; kwargs...)
> return Diagonal([exp(-im * pi * (k - 1)) for k in 1:dim])
> end
> #----------------------------parameter---------------------
> chain_length = 90
> Nbosons = 30
> U = 10.
> J = 1.
> V = 10.
> alpha = 1 // 3
> delta = 2 * pi / 3
> #----------------------------parameter dmrg----------------
> maxdim = [10, 10, 20, 20, 30, 30, 40, 40, 80, 80, 200, 200, 400, 400, 600, 600]
> nsweeps = 120
> cutoff = 1.e-12
> #----------------------------parameter dmrg----------------
> """
> function constructBoseHubbardHamiltonian(J::T, U::T, V::T, alpha::T, delta::T, chain_length::Int, pbc::Bool)::OpSum where T<:Float64
> 
> Construct the Bose-Hubbard Hamiltonian from the parameters, returns it as a sum of operators.
> """
> function constructBoseHubbardHamiltonian(J::T, U::T, V::T, alpha::T, delta::T, chain_length::Int, pbc::Bool)::OpSum where T<:Float64
> os = OpSum()
> for j = 1:chain_length-1
> os += -J, "adag", j, "a", j + 1, "ExpIθn", j
> os += -J, "ExpIθndag", j, "adag", j + 1, "a", j
> end
> for j = 1:chain_length
> os += U / 2, "n", j, "n", j
> os += -U / 2, "n", j
> end
> 
> for j = 1:chain_length
> os += V * cos(2 * pi * alpha * j + delta), "n", j
> end
> if pbc
> os += -J, "adag", chain_length, "a", 1, "ExpIθn", chain_length
> os += -J, "ExpIθndag", chain_length, "adag", 1, "a", chain_length
> end
> return os
> end
> 
> """
> function constructBoseHubbardHamiltonian(J::T, U::T, V::T, alpha::T, delta::T, chain_length::Int, pbc::Bool, sites::Vector{<:Index})::MPO where T<:Float64
> 
> Construct the Bose-Hubbard Hamiltonian from the parameters, returns it as an MPO acting on the sites.
> """
> function constructBoseHubbardHamiltonian(J::T, U::T, V::T, alpha::T, delta::T, chain_length::Int, pbc::Bool, sites::Vector{<:Index})::MPO where T<:Float64
> os = constructBoseHubbardHamiltonian(J, U, V, alpha, delta, chain_length, pbc)
> return MPO(os, sites)
> end
> 
> """
> function constructInitialState(chain_length::Int, Nbosons::Int, sites::Vector{<:Index})::Vector{String}
> 
> Construct an initial random state with fixed number of bosons.
> """
> function constructInitialState(chain_length::Int, Nbosons::Int, sites::Vector{<:Index})
> state = [n <= Nbosons ? "1" : "0" for n = 1:chain_length]
> psi0 = random_mps(sites, state)
> return psi0
> end
> 
> #----------------------------strat calculate---------------
> sites = siteinds("Boson", chain_length, conserve_qns=true, dim=5)
> 
> global densities = Dict{Int,Vector{Float64}}()
> global energies = Dict{Int,Float64}()
> if !(@isdefined ψs)
> global ψs = Dict{Int,MPS}()
> end
> for N in Nbosons-1:Nbosons+1
> tmp = get(ψs, N, nothing)
> if !isnothing(tmp)
> @info "Taking previous result"
> psi0 = tmp
> maxdim = 2 * maxlinkdim(psi0)
> sites = siteinds(psi0)
> else
> psi0 = constructInitialState(chain_length, N, sites)
> end
> @info "Constructing Hamiltonian for $N bosons"
> H = constructBoseHubbardHamiltonian(J, U, V, float(alpha), delta, chain_length, false, sites)
> @info "Starting DMRG for $N bosons"
> energy, psi = dmrg(H, psi0; nsweeps=nsweeps, maxdim=maxdim, cutoff=cutoff)
> ψs[N] = psi
> energies[N] = energy
> densities[N] = expect(psi, "n")
> end
> 
> quasiparticle_density = densities[Nbosons+1] - densities[Nbosons]
> 
> f1 = Figure(backgroundcolor=:white)
> ax = Axis(f1[1, 1],
> title = "α = $alpha δ = $(rationalize(delta/π))π L = $chain_length N = $(Nbosons) V = $V U = $U ",
> xlabel = "Lattice Site i",
> ylabel = "Quasiparticle Density"
> )
> resize_to_layout!(f1)
> scatterlines!(collect(1:chain_length), quasiparticle_density, label="ni")
> axislegend(position=:rb)
> display(f1)
> 
> end
> 
> ```
> 
> Something more complicated that I’ve added is that the eigenvector obtained is saved in a global variable, and if this variable is already defined it doesn’t get redefined. In this way when working from the REPL and calling `include("file_name.jl")` more than one time, the previous ground state is used as starting state, this allows to simply increment the number of sweeps and check for the convergence in a more easy way. If you have more questions on the code feel free to ask 😄.

---

<div class="post-metadata">

**Author:** ![Ashinxinl](https://yyz2.discourse-cdn.com/free1/user_avatar/itensor.discourse.group/ashinxinl/32/265_2.png) [@Ashinxinl](https://itensor.discourse.group/u/Ashinxinl)\
**Post date:** [March 6, 2025, 6:12am UTC](https://itensor.discourse.group/t/mpo-for-anyon-hubbard-model/2275/12 "2025-03-06T06:12:06Z")

</div>

@VinceNeede  
Thank you for your patient and detailed response. I believe I previously encountered the issue of DMRG getting stuck in a local minimum during my calculations. After following your advice and increasing the precision, I was able to achieve the results shown in the figure. Wishing you a happy life and success in your work!

---

<div class="post-metadata">

**Author:** ![system](https://global.discourse-cdn.com/free1/uploads/itensor/original/1X/d3072b13e047cd06df7f3981547c0917d940dfa4.png) [@system](https://itensor.discourse.group/u/system)\
**Post date:** [March 16, 2025, 6:12am UTC](https://itensor.discourse.group/t/mpo-for-anyon-hubbard-model/2275/13 "2025-03-16T06:12:59Z")

</div>

This topic was automatically closed 10 days after the last reply. New replies are no longer allowed.