TEBD for 2-leg Hubbard Model

Hi there! I am trying to implement time evolving block decimation for the Hubbard model. I believe my DMRG works just fine (the state convergences, my equilibrium correlations look good, and the ground state energy matches up with benchmark values), but there is a bug in how I am constructing the Trotter gates for my TEBD evolution.

Specifically: when I calculate the inner product of the time-evolved ground state with its original self, the overlap decreases with time. I encountered this previously in the case of the 1D TEBD, and ultimately realized that it was due to a mistake in how I was constructing my Trotter gates. The only difference between the 1D case and the 2-leg ladder, as far as I am aware, is the insertion of fSWAP gates before and after the application of a two-site operator between non-adjacent sites. Here is my code for constructing the Trotter gates:

  gates = ITensor[]
  for b in lattice
      s1_idx = min(b.s1, b.s2)
      s2_idx = max(b.s1, b.s2)

      # Check if we need to apply any swaps (sites are not adjacent) 
      if abs(s2_idx-s1_idx) > 1
          # Apply SWAP to get to new configuration
          SWAP_forward = []
          s_last = s1_idx
          for s_next in s1_idx+1:(s2_idx-1)
              push!(SWAP_forward, op("fSWAP", sites[s_last], sites[s_next]))
              s_last = s_next
          end
          append!(gates, SWAP_forward)

          # Apply gate
          hj_twosite = make_twosite(sites[s_last], sites[s2_idx])
          push!(gates,hj_twosite)

          # Apply SWAP to get back 
          SWAP_backward = []
          for s_next in s_last-1:s1_idx
              push!(SWAP_backward, op("fSWAP", sites[s_last], sites[s_next]))
              s_last = s_next
          end
          append!(gates, SWAP_backward)
          #append!(gates, reverse(SWAP_forward))

      else # no need to apply swaps 
          hj_twosite = make_twosite(sites[s1_idx], sites[s2_idx])
          push!(gates,hj_twosite)
      end

  end

  # One-site terms (no swap needed)
  for j=1:N-1
      hj_onesite = make_onesite(sites[j])     
      push!(gates,hj_onesite)
  end

  # End site (no swap needed)
  hn = make_endsite(sites[N])
  push!(gates,hn)

  # Add all the gates in reverse 
  append!(gates,reverse(gates))

I am fairly certain that my make_twosite(...), make_onesite(...), and make_endsite(...) gates are correct, because when I run this code with only one leg in the ladder, I am getting good overlap between the initial ground state and the time-evolved ground state. I’ve also attached a screenshot of my fSWAP operator in (what I believe to be) the correct basis.

Thank you so much for all the continued support! Hope you have a great week.

You should not need to implement your own (fermionic) swap gates, the apply function inserts them automatically. You can just input nonlocal gates into apply and it should just work.

1 Like

Hi Matt,

Thanks so much for the help! Your solution ended up fixing my problem. When I tried adding phonons (I have my own Hubbard Holstein site type) I got an out-of-bounds indexing error that I traced back to the _fermionic_swap(...) function. I think this is due to the fact that there is a bug in that function at the part where we fill in the diagonals of the blocks of the fSWAP operator. Here is what I did to fix the problem:

function _fermionic_swap(s1::Index, s2::Index)
  T = ITensor(QN(), s1', s2', dag(s1), dag(s2))
  for b in nzblocks(T) 
    dval = 1.0 
    # Must be a diagonal block
    ((b[1] ≠ b[3]) || (b[2] ≠ b[4])) && continue
    n1, n2 = b[1], b[2]
    if _isodd_fermionic_parity(s1, n1) && _isodd_fermionic_parity(s2, n2)
      dval = -1.0 
    end
    Tb = ITensors.blockview(ITensors.tensor(T), b) 
    mat_dim = prod(dims(Tb)[1:2])  
    Tbr = reshape(Tb, mat_dim, mat_dim)
    for i in diagind(Tbr)
      # NDTensors.setdiagindex!(Tbr, dval, i) <-- OLD CODE
      setindex!(Tbr, dval, i) # <-- MY MODIFICATION  
    end
  end
  return T
end 

Not sure if this is actually a bug, but changing that line avoids the out-of-bounds error and (seemingly) results in the correct time evolution, as measured by the fact that the overlap between |\psi_{GS}\rangle and |\psi_{GS}(t)\rangle remains close to 1.

Thanks again for all the help!