Multiple exchange interactions in 2D lattices

ITensor Julia Questions
,
1

Dear miles,

I would like to define two different latticebonds in a square lattice.
I made below changes in the lattices.jl file:

function square_lattice(Nx::Int, Ny::Int; kwargs...)::Lattice
  yperiodic = get(kwargs, :yperiodic, false)
  yperiodic = yperiodic && (Ny > 2)
  N = Nx * Ny
  Nbond = 2N - Ny + (yperiodic ? 0 : -Nx)
  
  latt = Lattice(undef, Nbond)
  b1 = 0
  b2 = 0
  for n in 1:N
    x = div(n - 1, Ny) + 1
    y = mod(n - 1, Ny) + 1
    if x < Nx
      latt[b1 += 1] = LatticeBond(n, n + Ny, x, y, x + 1, y)
    end
    if Ny > 1
      if y < Ny
        latt[b2 += 1] = LatticeBond(n, n + 1, x, y, x, y + 1)
      end
      if yperiodic && y == 1
        latt[b1 += 1] = LatticeBond(n, n + Ny - 1, x, y, x, y + Ny)
      end
    end
  end
  return latt
end

and also I make the Hamiltonian function as:

function H_fst_2D_square(Nx,Ny,JH,JI,Delta,sites)
    
    lattice = square_lattice(Nx, Ny; yperiodic = false)

    # Define the Heisenberg spin Hamiltonian on this lattice
      ampo = OpSum()
      for b1 in lattice
          ampo .+= 0.5*J1, "S+", b1.s1, "S-", b1.s2
          ampo .+= 0.5*J1, "S-", b1.s1, "S+", b1.s2
          ampo .+= J1, "Sz", b1.s1, "Sz", b1.s2
      end
      for b2 in lattice
          ampo .+= 0.5*J2, "S+", b2.s1, "S-", b2.s2
          ampo .+= 0.5*J2, "S-", b2.s1, "S+", b2.s2
          ampo .+= J2, "Sz", b2.s1, "Sz", b2.s2
      end
      H = MPO(ampo,sites)

    return(H)
end

after running the code I get below error:

ERROR: LoadError: UndefRefError: access to undefined reference
Stacktrace:
 [1] getindex
   @ ./array.jl:924 [inlined]
 [2] iterate
   @ ./array.jl:898 [inlined]
 [3] H_fst_2D_square(Nx::Int64, Ny::Int64, JH::Int64, JI::Float64, Delta::Int64, sites::Vector{Index{Vector{Pair{QN, Int64}}}})

Should I do something else?

For the uniform exchange coupling I can run the code without problem.

QUESTION:

  1. I think your pre-defined periodic boundary conditions along the y-direction is not true (what I see in lattices.jl file after installing ITensor).
    In fact, in line:
if yperiodic && y == 1
        latt[b += 1] = LatticeBond(n, n + Ny - 1, x, y, x, y + Ny)
end

it should be like:

if yperiodic && y == 1
        latt[b += 1] = LatticeBond(n, n + Ny - 1, x, y, x, y + Ny - 1)
end

am I right?

  1. I also obtained the magnetization per site for a square lattice with uniform exchange coupling by using ED method under full periodic boundary conditions. In ITensor I considered periodic boundary conditions along the y-direction. For 24 number of spins I got large discrepancy in magnetic field position of the true magnetization plateaus and jumps with the ED results.
    I increased the number of sweeps up to 20 and linkdim=200 with cutoff!(sweeps, 1E-5,1E-5,1E-5,1E-8,1E-8).
    Do I encounter this discrepancy (magnetic field position of the magnetization plateaus) for large-size lattice, let say for N=144?

  2. When I compose my own lattice function with an arbitrary name, i.e., function triangle_hexagon_lattice(Nx::Int, Ny::Int; kwargs...)::Lattice, I get the error: ERROR: LoadError: UndefVarError: triangle_hexagon_lattice not defined.
    Should I define the lattice name triangle_hexagon_lattice in another file?

Thanks!

2

Hi Hamid,
Glad you’re trying out a 2D calculation with ITensor.

Regarding your first question, it looks like there is a bug in your code where an array is accessed out of bounds (beyond the size of that array). What is the line of your H_fst_2D_square function on which the error occurs? Basically it looks like a bug in your own code and I encourage you to look through the stack trace to understand the source of the error.

Regarding your other questions:

  1. yes I think you’re correct that the formula should be y+Ny-1. I’ll file an issue or a fix for that right now, thanks.

  2. the DMRG calculation you are doing is using open boundary conditions in the x direction, while your ED calculation is using fully periodic boundary conditions. So in general one would just expect the results to be different. One can do a careful finite-size analysis and extrapolation to estimate bulk properties from both calculations and that comparison might match if both calculations have small finite-size effects. But in general there’s no reason to expect that the results should match. Most of all, the fully periodic calculation might not break the spin symmetry whereas the DMRG one definitely could.

  3. an error of this type is less about whether you should or should not put something into a file, and more about whether the file you put it in is included into the code you are running (either by doing using ModuleName if you put the code into a module named ModuleName or by doing include("filename") if you put the code into a file named "filename").

3

Dear miles,

Thank you very much for your reply.
I tried to introduce a simple square lattice with two different latticebonds (b1 and b2) or in other words two exchange couplings J1 and J2. But I get the below error.
Can you please specify that how can I introduce different latticebonds for 2D lattices. I changed the function square_lattice in lattices.jl file as:

function square_lattice(Nx::Int, Ny::Int; kwargs...)::Lattice
  yperiodic = get(kwargs, :yperiodic, false)
  yperiodic = yperiodic && (Ny > 2)
  N = Nx * Ny
  Nbond = 2N - Ny + (yperiodic ? 0 : -Nx)
  
  latt = Lattice(undef, Nbond)
  b1 = 0
  b2 = 0
  for n in 1:N
    x = div(n - 1, Ny) + 1
    y = mod(n - 1, Ny) + 1
    if x < Nx
      latt[b1 += 1] = LatticeBond(n, n + Ny, x, y, x + 1, y)
    end
    if Ny > 1
      if y < Ny
        latt[b2 += 1] = LatticeBond(n, n + 1, x, y, x, y + 1)
      end
      if yperiodic && y == 1
        latt[b1 += 1] = LatticeBond(n, n + Ny - 1, x, y, x, y + Ny - 1)
      end
    end
  end
  return latt
end

and the Hamiltonian is as following form:

using ITensors

let
  Ny = 6
  Nx = 4

  N = Nx*Ny

  sites = siteinds("S=1/2", N; conserve_qns = true)

  # Obtain an array of LatticeBond structs
  # which define nearest-neighbor site pairs
  # on the 2D square lattice (wrapped on a cylinder)
  lattice = square_lattice(Nx, Ny; yperiodic = false)

  # Define the Heisenberg spin Hamiltonian on this lattice
  os = OpSum()
  for b1 in lattice
    os .+= 0.5, "S+", b1.s1, "S-", b1.s2
    os .+= 0.5, "S-", b1.s1, "S+", b1.s2
    os .+=      "Sz", b1.s1, "Sz", b1.s2
  end
  
  for b2 in lattice
    os .+= 0.5, "S+", b2.s1, "S-", b2.s2
    os .+= 0.5, "S-", b2.s1, "S+", b2.s2
    os .+=      "Sz", b2.s1, "Sz", b2.s2
  end
  
  H = MPO(os,sites)

  state = [isodd(n) ? "Up" : "Dn" for n=1:N]
  # Initialize wavefunction to a random MPS
  # of bond-dimension 10 with same quantum
  # numbers as `state`
  psi0 = randomMPS(sites,state,10)

  nsweeps = 5
  maxdim = [20,60,100]
  cutoff = [1E-4]

  energy,psi = dmrg(H,psi0; nsweeps, maxdim, cutoff)

  return
end

The error is:

ERROR: LoadError: UndefRefError: access to undefined reference
Stacktrace:
 [1] getindex
   @ ./array.jl:924 [inlined]
 [2] iterate(A::Vector{LatticeBond}, i::Int64)
   @ Base ./array.jl:898
 [3] top-level scope
4

Dear miles,

Could you please address this problem, too?
Thanks!

5

Hi Hamid,
It seems like you encountered a bug in your code. Judging from the error message, a part of the code tried to access the vector of LatticeBonds at an index that was outside of its length (so the index was either < 1 or > length of that vector).

If you look at the rest of the error message stack trace, it ought to tell you exactly which function or line of code led to this error, then I would encourage you to print out or otherwise look at the indices going into the vector to make sure they are within range.