TEBD for 2-leg Hubbard Model

[nticea]

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.

[mtfishman]

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.

[nticea]

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!