I’ve seen a few posts alluding to time evolution of a Hubbard Hamiltonian (or more generally a Hamiltonian with more than one term), but I’m not sure how it would work in practice. I’m getting an error about how the network can’t be contracted (Indices to not have opposite directions), but I feel like my approach might be incorrect in general. Here’s my code:
Sorry if this question is trivial, I’m definitely still getting used to ITensor and navigating this discourse group. Thanks for all of your help, and for the amazing package!
Could you share a full runnable code? See Please Read: Make It Easier to Help You where it says to post a minimal working example of your code (i.e. something we can run as-is, including the definition of all of the objects/variables needed and the packages you loaded). It should be something we could copy and paste and be able to run without modification.
Of course, sorry for the oversight. The idea of this code is to compute the time evolution of the operator given in equation 3 of this paper: https://arxiv.org/pdf/1503.02010
The paper uses iTebd and therefore does not define the length of the lattice, but uses L = 10 when comparing to a solution through exact diagonalization, so thats what I used in my code.
using ITensors, ITensorMPS, Plots
plt = plot()
let
N = 40
cutoff = 1E-10
maxdim = 1024
tau = 0.1
ttotal = 5
s = siteinds("Electron", N; conserve_qns=true)
gates = ITensor[]
U = 8
t = 1
for j in 1:(N-1)
s1 = s[j]
s2 = s[j + 1]
H_t = -t * (op("Cdagup", s2) * op("Cup", s1) +
op("Cdagup", s1) * op("Cup", s2) +
op("Cdagdn", s2) * op("Cdn", s1) +
op("Cdagdn", s1) * op("Cdn", s2))
H_U = U * op("Nupdn", s1)
Gj = exp(-im * tau / 2 * H_t) * exp(-im * tau / 2 * H_U)
push!(gates, Gj)
end
# Include gates in reverse order
append!(gates, reverse(gates))
psi = MPS(s, n -> isodd(n) ? "Up" : "Dn")
os = OpSum()
for i in 1:(10)
os += "Nup", i, "Ndn", i
end
dt = MPO((1/10)*(os), s)
Counting = []
for t in 0.0:tau:ttotal
D = inner(psi', dt, psi)
push!(Counting, real.(D))
println(t)
t ≈ ttotal && break
psi = apply(gates, psi; cutoff, maxdim)
normalize!(psi)
end
plot!((t*(0.0:tau:ttotal)), Counting)
end
display(plt)
As for my version and package info: Julia 1.10, ITensors v0.5.2. Thanks so much!
I haven’t tested it but I think this line is the one giving you trouble:
You might want to change that to:
Gj = product(exp(-im * tau / 2 * H_t), exp(-im * tau / 2 * H_U))
The reason is that those two are ITensors with some matching primed and unprimed indices, so just using * contracts them over the corresponding primed and unprimed indices, which is more like an inner product, but they have matching arrows so the arrow directions are not correct in that contraction. product acts more like a matrix multiplication of the two gates by manipulating the prime levels correctly.
Glad to hear it! It’s one of the trickier parts of how ITensor works, * has a particular definition in ITensor as tensor contraction over shared indices but is often used for matrix multiplication in other contexts. I sometimes wish we could use * for matrix multiplication-like contractions like you were expecting in your original code but that would be too confusing and magical to have it have both meanings. We plan to define “ITensor maps” which would be ITensors that are maps between two sets of indices, and those could act more like matrices between those two sets of indices with a matrix algebra interface, but functions like product work well for the time being.
That makes a lot of sense, and ITensor is already like magic as is. Looking forward to ITensor Maps, sounds like they might compute matrix operations more quickly than ‘’’ product ‘’'. Thanks so much for the help and thorough explanation!