Ising-Heisenberg spin models

DMRG and Numerical Methods
1

Hi miles,

I an interested to simulate on-site magnetization of a typical hybrid Ising-Heisenberg spin 1/2 model on a 2D frustrated lattice.
I reproduced the corresponding Hamiltonian and applied the dmrg to do so.
Finally, I get unconventional results that seem not to be consistence with the exact analytical method.
Can you please let me know where is the problem.
Please find in below the Hamiltonian and the dmrg code.

using ITensors
using DelimitedFiles

# Hamiltonian:
function H_IH_martini(Nx,Ny,JH,JI,Delta,sites)
  ampo = AutoMPO()
  N = Nx*Ny
  for n in 1:N
    x = div(n - 1, Ny) + 1
    y = mod(n - 1, Ny) + 1
    a = 0.5
    if x < Nx && y < Ny
       if mod(y, 2) == 1 && mod(x, 2) == 1  # Red balls
          ampo += a*JH,"S+",n,"S-",n+1
          ampo += a*JH,"S-",n,"S+",n+1
          ampo += Delta,"Sz",n,"Sz",n+1
          # latt[bH += 1] = LatticeBond(n, n + 1, x, y, x, y + 1)  
          #
          ampo += a*JH,"S+",n,"S-",n+Ny+1
          ampo += a*JH,"S-",n,"S+",n+Ny+1
          ampo += Delta,"Sz",n,"Sz",n+Ny+1

          # latt[bH += 1] = LatticeBond(n, n + Ny + 1, x, y, x + 1, y + 1)
       end
    end
    if x < Nx && y < Ny
       if mod(y, 2) == 1 && mod(x, 2) == 1  # Red balls
          ampo += JI,"Sz",n,"Sz",n+Ny
          # latt[bI += 1] = LatticeBond(n, n + Ny, x, y, x + 1, y)  
       end
    end
    if x < Nx && 1 < y <= Ny
       if mod(y, 2) == 0 && mod(x, 2) == 1 # Blue balls
          ampo += a*JH,"S+",n,"S-",n+Ny
          ampo += a*JH,"S-",n,"S+",n+Ny
          ampo += Delta,"Sz",n,"Sz",n+Ny
          # latt[bH += 1] = LatticeBond(n, n + Ny, x, y, x + 1, y)
       end
    end
    if 1 < x <= Nx && 1 < y < Ny  
       if mod(y, 2) == 0 && mod(x, 2) == 0 # Yellow balls
          ampo += JI,"Sz",n,"Sz",n+1
          # latt[bI += 1] = LatticeBond(n, n + 1, x, y, x, y + 1)
       end
    end
    if 1 < x < Nx && 1 < y < Ny 
       if mod(y, 2) == 1 && mod(x, 2) == 0  # White balls
          ampo += JI,"Sz",n,"Sz",n+Ny-1
           # latt[bI += 1] = LatticeBond(n, n + Ny - 1, x, y, x + 1, y - 1)
       end
    end
    
     # periodic bonds
    if 1 < x <= Nx && y == 1
      if mod(x, 2) == 0 # White-Yellow balls
         ampo += JI,"Sz",n,"Sz",n+Ny-1
         # latt[bI += 1] = LatticeBond(n, n + Ny - 1, x, y, x, y + Ny)
      end
    end
    if 1 < x < Nx && y == 1 # White-Blue balls
      if mod(x, 2) == 0
         ampo += JI,"Sz",n,"Sz",n+2*Ny-1
         # latt[bH += 1] = LatticeBond(n, n + 2*Ny - 1, x, y, x + 1, y + Ny)
      end
    end
    
    if x == 1 # White-Blue balls in x dir
      if mod(y, 2) == 0
         ampo += JI,"Sz",n,"Sz",Nx*Ny-n+1
         # latt[bH += 1] = LatticeBond(n, n + Nx*Ny -n+1, x, y, x + 1, y + Ny)
      end
    end
    
  end
    
    H = MPO(ampo,sites)

    return(H)
end

and the dmrg part of the code is as following:

let
@time begin

# physical parameters
Nx = 10
Ny = 8
N = Nx*Ny
JH = 0.5
JI = 1
Delta = JH*1

linkdim = 2000
sweeps = Sweeps(15)
maxdim!(sweeps, 500, 800, 1000, 1000, 1200)
cutoff!(sweeps,1E-10,1E-12,1E-16,1E-18,1E-20,1E-24,1E-26)
noise!(sweeps,1E-10,1E-12,1E-16,1E-18,1E-20,1E-24,1E-26,0)

# system
sites = siteinds("S=1/2",N; conserve_sz=true)

sz_min = 0
sz_max = 40
sz_max_psiAF = Int(N/2)
S2 = S2_op(N,sites)

Szin = [0. for m in 1:N]
SzT = [0. for m in 1:sz_max-sz_min+1]
E_GS = [0. for m in 1:sz_max-sz_min+1]
h = [0. for m in 1:sz_max-sz_min+1]
Szih = [0 for m in 1:sz_max-sz_min+1]

statei = ["Dn" for n in 1:N]
for j in 1:Int(N/2)
     statei[j] = "Up"
end
@show statei
mz = 0
for j in sz_min:sz_max
     mz += 1
     Szih[mz] = j
end
@show Szih
m_up = 0
for i in sz_max_psiAF+sz_min:sz_max_psiAF + sz_max  
    # if i <= sz_max_psiAF + sz_max 
    @show m_up + sz_min
    m_up += 1
    # end
    for j in 1:i
        statei[j] = "Up"
    end
    @show statei
    ψi = randomMPS(sites,statei,linkdim)
    Szi = [Szi_op(szi,sites) for szi in 1:N]
    H = H_IH_martini(Nx,Ny,JH,JI,Delta,sites)
    E0,ψ0 = dmrg(H,ψi,sweeps)
    println()
    # E_GS[n]
    E_GS[m_up] = E0
    # Szi(i)
    for k in 1:N
        Szin[m_up] = inner(ψ0,Szi[k],ψ0)
    end
    # SzT(i)
    for n in 1:N
        SzT[m_up] += Szin[n]
    end
    @show SzT[m_up]]
end

# Local field h(i)
h[1] = 0
for i in 2:length(E_GS)
    h[i] = E_GS[i] - E_GS[i-1]
end

As you can see, I increased linkdim upto 2000 and decreased the cutoff down to 1E-26. But still the results of the magnetization is not correct.
The running time is too long, while the model is an Ising-Heisenberg model where the run time should be much more shorter than full Heisenberg model. I used other packages such as TeNeS and also dmrgpy both of them get reasonable results. Why I am interested in using the ITensor, because it allows to conserve the quantum number to fix the magnetization magnitude.
Please know that the code works well for the full Heisenberg model.

Thank you for your support.

2

If your question is about how to get good results with DMRG, I’d refer you to the following FAQ section of our website:
https://itensor.github.io/ITensors.jl/dev/faq/DMRG.html

It’s interesting that you got good results with other packages. There’s a certain parameter called eigsolve_krylovdim that is set by default to a higher value in other packages, such as TeNPy, but a lower default value in ITensor. So one thing you could try first is to set eigsolve_krylovdim to a larger value, such as 10 to see if that helps. For example,

dmrg(H,psi; nsweeps, maxdim, cutoff, eigsolve_krylovdim=10)

Lastly, please do not post long codes and request that we read through them and find bugs for you. While it’s helpful to post minimal example codes to illustrate what you’re doing, please don’t post codes longer than a few lines. Thank you.