Dear there,
I’ve a quick question. Just wondering is there a way or function in julia version Itensor to do exact diagonalization? I mean after I have a Hamiltonian mpo, for a small system, and just want to know all energy eigenvalues of it.
Thank you so much!
Hi Kristin, The short answer is yes. Here is a code sample:
using ITensors
let
N=6
s=siteinds("S=1/2",N;conserve_qns=true)
os = OpSum()
for j in 1:(N - 1)
os += "Sz", j, "Sz", j + 1
os += 0.5, "S+", j, "S-", j + 1
os += 0.5, "S-", j, "S+", j + 1
end
cb=combiner(dag(s)...)
cbp=combiner(s'...)
H = prod(MPO(os, s))*cb*cbp
E,V=eigen(H,combinedind(cb),combinedind(cbp))
@show sort(array(real(diag(E))))
end
Efficiency can be improved with ishermitian=true in the eigen call.
I hope this helps.
JR
Also there is an exact diagonalization example code one can find in our examples folder here:
Thanks for your code sample! I just tried but got a message saying:
Contraction resulted in ITensor with 15 indices, which is greater
than or equal to the ITensor order warning threshold 14.
You can modify the threshold with macros like @set_warn_order N,
@reset_warn_order, and @disable_warn_order or functions like
ITensors.set_warn_order(N::Int), ITensors.reset_warn_order(), and
ITensors.disable_warn_order().
Do you happen to know what does this mean…? I’m not very sure if it’s becasue my H is a little unusual, my site type is composite, I have fermion type on odd sites (1,3,5,…) and an effective spin1/2 type on even sites (2,4,6…). Here is what I just tried:
H = MPO(ampo, sites)
cb=combiner(dag(sites)...)
cbp=combiner(sites'...)
H = prod(MPO(ampo, sites))*cb*cbp
E,V=eigen(H,combinedind(cb),combinedind(cbp))
@show sort(array(real(diag(E))))
Thank you!
Thank you! I just tried but I can’t Pkg add MKL in the cluster I’m using for some reason. I need contact cluster people tomorrow before I can really try this…
But may I ask what do “fuse” and “tree” do here? I’m not very familiar with these, maybe a naive question…
if fuse
if binary
println("Fuse the indices using a binary tree")
T = fusion_tree_binary(s)
H_full = @time fuse_inds_binary(H, T)
ψ0_full = @time fuse_inds_binary(ψ0, T)
else
println("Fuse the indices using an unbalances tree")
T = fusion_tree(s)
H_full = @time fuse_inds(H, T)
ψ0_full = @time fuse_inds(ψ0, T)
end
else
println("Don't fuse the indices")
@disable_warn_order begin
H_full = @time contract(H)
ψ0_full = @time contract(ψ0)
end
end
end
Working with tensors that have a very large number of indices can use a lot of memory, and also a common type of bug in ITensor codes can be mistakenly not contracting an Index, so we put that warning by default to catch such issues. The meaning of the warning is exactly what the text says, which is that somewhere in the code an ITensor was made with 15 or more indices, and you can follow one of the steps listed in the warning to disable the warning.
I would recommend adding this line:
ITensors.disable_warn_order()
to your code.
The part about “fuse” are basically different optimizations of the code, trying to perform the tensor operations in different ways by merging indices together or not. You can safely ignore most of that code if you are just trying to do a quick exact diagonalization (ED) calculation. On the other hand, if you are trying to push to a very large calculation, you may want to experiment with enabling different branches of that if statement to see which setting (such as trying fuse=true or fuse=false) gives you better performance.
I see! I just tried with a small testing case, a tight-binding model for L=6 sites. But the result looks a little strange… I got lots of eigenvalues while I only expect 6… Is there anything I did wrong in the code?
using ITensors
let
N=6
s=siteinds("Electron",N;conserve_qns=true)
os = OpSum()
for j in 1:(N - 1)
os += -1, "Cdagup", j, "Cup", j + 1
os += -1, "Cdagup", j + 1, "Cup", j
end
cb=combiner(dag(s)...)
cbp=combiner(s'...)
H = prod(MPO(os, s))*cb*cbp
E,V=eigen(H,combinedind(cb),combinedind(cbp))
@show sort(array(real(diag(E))))
end
Sorry, but we cannot take requests from users to debug their codes for them. Please try running your code on different examples, printing out various quantities, etc to check what the code is doing and what might be the issue.
For a system of size L=6, you would expect to get 2^6 eigenvalues/eigenstates.
For Electrons it should be 4^6
Right, thanks for the correction!
@Kristin perhaps your confusion comes from the fact that there are O(N) (i.e. 6 for each spin, in your case) single particle energy levels. ITensor works with the full many-body Hamiltonian, and the Hilbert space of N electron (i.e. spinful fermion) sites is 4^N. You would have to map your Hamiltonian to a O(N) \times O(N) hopping Hamiltonian to get the single particle energies, perhaps https://github.com/ITensor/ITensors.jl/tree/main/ITensorGaussianMPS might help (see for example https://github.com/ITensor/ITensors.jl/blob/main/ITensorGaussianMPS/examples/hubbard_1d_no_spin_conservation.jl, where you can see how to convert an OpSum of quadratic fermionic terms into a hopping Hamiltonian).
Yes, sorry I was thinking of a single particle spectrum, then fill in noninteracting electrons to the energy level etc., thanks for remindering me and the link!
Thanks for reminding me this! I’ve one more quick question, just wondering what’s the purpose of cb and cbp? Because in my previous calculation, I usually just define sites, then set H=MPO(ampo, sites), then doing dmrg, I haven’t used combiner of dag(sites) or site’ (also not sure what’s difference between dag(s) and s’…), and the prod of H etc, I’m curious if these are needed for doing ED or I’ve missed some regular setting up steps in dmrg also? Thank you!
cb=combiner(dag(s)...)
cbp=combiner(s'...)
H = prod(MPO(os, s))*cb*cbp
Jan can give his reasons for the cb and cbp combiners in his code, which is really nice and simple, but there could be multiple reasons for introducing these in the context of an ED calculation.
The simplest one would be to map the Hamiltonian to an actual matrix (2 index tensor) so that matrix operations like eigen are easier to do and to reason about, even though ITensor itself lets you call eigen on whole tensors.
The deeper reason is that it could make the code perform better (run faster) because Hamiltonians typically are very sparse and by combining the indices, the various small, sparse blocks can be combined together into larger dense blocks, making certain low-level matrix operations perform better. More specialized ED codes (ITensor itself is not so specialized for ED) would directly take advantage of the sparse structure of the Hamiltonian without needing to combine blocks togther. I think in the end one would have to try both approaches (sparse oriented approach versus blocking approach) to see which one performs the best in which regimes (i.e. for large and small problems with different amounts of symmetry etc.).
Yes the combiners are not essential. Under the hood H is stored in symmetry adapted blocks. For the 4^6 Hilbert space in the example, with combiners you get 35 non-zero blocks, and the eigen routine takes 0.25 sec on my machine. Without the combiners you get 10240 non-zero blocks and the eigen routine takes 0.81 sec.
It is interesting to see that the combiners do help! (Even if both versions are relatively fast.)
I see! Thank you!
Amazing Thank you!
ITensor Website
·
Twitter @ITensorLib
© 2022 ITensor Collaboration