VUMPS example, 2d heisenberg model interaction term

DMRG and Numerical Methods
1

The following is an example of creating a 2D-Heisenberg model for VUMPS calculations. However, I don’t quite understand the part where the interaction positions are specified when creating the opsum.
First, due to the use of mod, values corresponding to position 0 are included. What does the value 0 for site mean? When width=4, vertical interaction (3, 0) occurs but interaction (3,4) is absent.

Second, vertical hopping is defined only on the first line, but shouldn’t vertical hopping also be specified on the second line?

(ITensorInfiniteMPS.jl/examples/vumps/vumps_2d_heisenberg.jl at main · ITensor/ITensorInfiniteMPS.jl · GitHub)

using ITensorInfiniteMPS
using ITensors, ITensorMPS
N = width # Number of sites in the unit cell


include(
  joinpath(
    pkgdir(ITensorInfiniteMPS), "examples", "vumps", "src", "vumps_subspace_expansion.jl"
  ),
)

maxdim = 256 # Maximum bond dimension
cutoff = 1e-6 # Singular value cutoff when increasing the bond dimension
max_vumps_iters = 100 # Maximum number of iterations of the VUMPS algorithm at each bond dimension
vumps_tol = 1e-4
conserve_qns = true
solver_tol = (x -> x / 10)
outer_iters = 10 # Number of times to increase the bond dimension
width = 4

initstate(n) = isodd(n) ? "↑" : "↓"
s = infsiteinds("S=1/2", N; conserve_qns, initstate)
ψ = InfMPS(s, initstate)

function ITensorInfiniteMPS.unit_cell_terms(::Model"heisenberg2D"; width)
  opsum = OpSum()
  for i in 1:width
    # Vertical
    opsum += -0.5, "S+", i, "S-", mod(i + 1, width)
    opsum += -0.5, "S-", i, "S+", mod(i + 1, width)
    opsum += "Sz", i, "Sz", mod(i + 1, width)
    # Horizontal
    opsum += -0.5, "S+", i, "S-", i + width
    opsum += -0.5, "S-", i, "S+", i + width
    opsum += "Sz", i, "Sz", i + width
  end
  return opsum
end
2

From what I can tell, interactions with a site <= 0 are ignored, so there is a mistake in this example.
It should be

function ITensorInfiniteMPS.unit_cell_terms(::Model"heisenberg2D"; width)
  opsum = OpSum()
  for i in 1:width
    # Vertical
    opsum += -0.5, "S+", i, "S-", i + 1
    opsum += -0.5, "S-", i, "S+", i + 1
    opsum += "Sz", i, "Sz", i + 1
    # Horizontal
    opsum += -0.5, "S+", i, "S-", i + width
    opsum += -0.5, "S-", i, "S+", i + width
    opsum += "Sz", i, "Sz", i + width
  end
  return opsum
end
3

If width=4 case, why we don’t consider the vertical interactions between (5,6), (6,7), (7,8), (8,9)…

4

If you have a width=4 case and a unit cell size of N=width (4 sites), we’ll have in the unit cell from the # vertical section bonds

(1,2), (2,3), (3,4)

And from the # horizontal section

(1,5), (2,6), (3,7), (4,8)

The sites 5,6,7,8 are all in the next unit cell, so we don’t need to consider them here, only the connections within and out of our unit cell. If you want an 8+ site unit cell, then those would be created automatically from the site shifting tools within ITensorInfiniteMPS.

To see this, let’s just switch to Sz*Sz and make it the coupling larger in one direction

function ITensorInfiniteMPS.unit_cell_terms(::Model"heisenberg2D"; width)
  opsum = OpSum()
  for i in 1:width
    # Vertical
    a = 100
    opsum += a,"Sz", i, "Sz", i+1
    # Horizontal
    opsum += "Sz", i, "Sz", i + width
  end
  return opsum
end

We can see the terms of the unit cell

julia> ITensorInfiniteMPS.unit_cell_terms(model;width)
sum(
  100.0 Sz(1,) Sz(2,)
  1.0 Sz(1,) Sz(5,)
  100.0 Sz(2,) Sz(3,)
  1.0 Sz(2,) Sz(6,)
  100.0 Sz(3,) Sz(4,)
  1.0 Sz(3,) Sz(7,)
  100.0 Sz(4,) Sz(5,)
  1.0 Sz(4,) Sz(8,)
)

And if we want to make it into a finite cell of 3 unit cells (a total of 3*width=12 sites), we can see we get

julia> ITensorInfiniteMPS.opsum_finite(model,width*3; width)
sum(
  100.0 Sz(1,) Sz(2,)
  1.0 Sz(1,) Sz(5,)
  100.0 Sz(2,) Sz(3,)
  1.0 Sz(2,) Sz(6,)
  100.0 Sz(3,) Sz(4,)
  1.0 Sz(3,) Sz(7,)
  100.0 Sz(4,) Sz(5,)
  1.0 Sz(4,) Sz(8,)
  100.0 Sz(5,) Sz(6,)
  1.0 Sz(5,) Sz(9,)
  100.0 Sz(6,) Sz(7,)
  1.0 Sz(6,) Sz(10,)
  100.0 Sz(7,) Sz(8,)
  1.0 Sz(7,) Sz(11,)
  100.0 Sz(8,) Sz(9,)
  1.0 Sz(8,) Sz(12,)
  100.0 Sz(9,) Sz(10,)
  100.0 Sz(10,) Sz(11,)
  100.0 Sz(11,) Sz(12,)
)

which has the terms you’re considering

(Also note this example doesn’t have periodic boundaries for the y direction of course)

6

In the 3 unit cell output, (5,6), (6,7) and (6,8) are also included. Both these and the (9,10), (10,11), and (11,12) are generated by shifting the original unit cell over.

The important thing is that those terms on only within the unit cell, and so we dont need to include the next cell’s intra-cell terms ((5,6) etc) as they will be generated when forming the supercell (shifting the unit cell over).
In other words, we take the original unit cell

sum(
  100.0 Sz(1,) Sz(2,)
  1.0 Sz(1,) Sz(5,)
  100.0 Sz(2,) Sz(3,)
  1.0 Sz(2,) Sz(6,)
  100.0 Sz(3,) Sz(4,)
  1.0 Sz(3,) Sz(7,)
  100.0 Sz(4,) Sz(5,)
  1.0 Sz(4,) Sz(8,)
)

and add N=width here to everything to generate the next unit cell, so (1,2) → (1+N, 2+N)= (5,6), and so on. Then again for the third cell. If we included (5,6) in the original definition, it might be double counted.
If we define the unit_cell_terms properly, just by shifting the cell we should be able to generate the appropriate terms without listing them explicitly, but you can always enlarge the unit cell to list them explicitly too.

7

I understand it: In the case of a single unit cell, the interactions (1,2), (2,3), and (3,4) automatically generate the interactions (5,6), (6,7), and (7,8) when forming the next unit cell.

Thanks!

1 Like
8

Thanks again for pointing this out, the latest commit now has a fixed example that has a proper 2D model (with or without cylindrical boundaries)

1 Like
9

Hi Ryan，

Thank you for posting the example code.

In the vertical bonds, why is there a bond between width and width+1? Shoud not we limit i between 1 and width-1 for the vertical bonds?

Best,
Feng

10

Yes great point, thank you!
I think some things got lost in removing the mod, I’ll open a PR to update the code shortly

1 Like