How to correctly do TEBD with ITensor?

ITensor Julia Questions
,
1

Hello,

I am new to Julia and I am trying learn how to use ITensor to do time evolution block decimation (TEBD) for a certain system. I already have working code written in Python that I use to see if the results of my Julia code is correct, which it is not. To do TEBD with ITensor I used this example.

To give some context, time evolution in my model is given by applying two propagators, first an even one, then an odd one

p(t+1) = U_o U_e p(t)

Each propagator is made of gates that act on three sites U_e = U_{123} U_{345} U_{567} ... U_{n-1, n} that commute (and U_o = U_{1,2} U_{234} U_{456} ... U_{n-2,n-1,n}). The boundary is special, but simply consists of adding a two site gate. These are like e^{-itH} for quantum systems, but they require no Trotterization.

So I make a random initial MPS, evolve it for a long time, and the system should decay to a steady state.

using ITensors
let
  N = 4
  cutoff = 1E-8
  tsteps = 1000

  # Maximum bond dimension
  M = 100 # No truncation
  
  # Bulk gate
  U = [1 0 0 0 0 0 0 0;
       0 0 0 1 0 0 0 0;
       0 0 1 0 0 0 0 0;
       0 1 0 0 0 0 0 0;
       0 0 0 0 0 0 1 0;
       0 0 0 0 0 0 0 1;
       0 0 0 0 1 0 0 0;
       0 0 0 0 0 1 0 0]

  # Boundary gates
  d = [1/4, 2/3, 2/5, 3/7];

  U_L = [d[1]     0    d[1]    0;
          0     d[2]    0     d[2];
         1-d[1]   0   1-d[1]   0;
          0     1-d[2]  0    1-d[2]]

  U_R = [ d[3]    d[3]    0      0;
         1-d[3]  1-d[3]   0      0;
          0       0      d[4]   d[4];
          0     0       1-d[4] 1-d[4]]

  # Make an array of 'site' indices S=1/2 
  s = siteinds("S=1/2", N)

  println("Making gates")
  gates = ITensor[]

  # First generate even, then odd gates
  for i in 1:2
    for j in i:2:(N - 2)

      s1 = s[j]
      s2 = s[j + 1]
      s3 = s[j + 2]
  
      G = op(U, s1, s2, s3)
      push!(gates, G)
      
    end

    # Boundary driving
    if i == 1
      push!(gates, op(U_R, s[N-1], s[N])) # Right boundary
      push!(gates, op(U_L, s[1], s[2]))   # Left Boundary
    end
  end

  # Initialize psi to be a random MPS
  println("Making random MPS with bond dimension " * string(M))
  p = randomMPS(s, linkdims=M);

  println("Doing time evolution.")
  # Apply gates
  for t in 0:1:tsteps
    p = apply(gates, p; cutoff)
  end
  
  println("Done.")

  return
end

After running this code I contract all my bond indices to get a vector of size 2^N and check if it’s the same vector I get from my Python code, and it isn’t. The code looks correct to me, but I’m since this is my first time using ITensor I cannot figure out what the error would be.

Thank you

2

Hi, glad you are trying out ITensor. Your code looks possibly correct, so it’s hard to say why you aren’t getting the result you expect. It could be a bug like perhaps it’s not doing the steps that you think? I would suggest altering the code to print out more information to confirm (e.g. print things like your t variable parameterizing your for loop at the end - things like that).

The other things I would suggest are to:

  1. print out some of the gates (ITensors) explicitly to see whether they are being defined properly and as you think and

  2. trying both codes for a much smaller number of layers (much smaller tsteps) and compare the results after one step, two steps, etc.

You may be able to think of some other debugging strategies too.