Is there a OpProd function like OpSum to have products of operators?
I need to create an operator in a bosonic space that is like:
(a + a\dagger)^N + a\dagger a
I tried with: (freq is just a constant and m is the site number)
ϕ = OpSum()
ϕ += freq, “a”,m
ϕ += freq, “a†”,m
cosine = OpSum()
cosine += c_p, ϕ, [1:N_sites]
But it is not working
BoundsError: attempt to access 0-element Vector{Int64} at index [1]
No hints?
Could you be more specific? It is easier to help if you give more information. I don’t really understand the equation or the code you wrote. Please write a self contained equation defining all of the variables, and a full code that you hope would work with all of the variables defined (and the code formatted, see Welcome to ITensor Discourse).
I am sorry.
I want to write a bosonic hamiltonian that is:
I was wondering if there is a built in method to build a product (or in this case the power) of an operator.
To recap, I cannot build (a_s + a_s^\dagger)^n and sum it to other operator.
using ITensors
N_sites = 2
N = 5
sites = siteinds("Boson",N_sites,dim=10)
op = OpSum()
for m in 1:N_sites
op += 1, "a",m
op += 1, "a†",m
end
pow_of_op = OpSum()
for n in 1:N
pow_of_op += 1, op^n, sites
end
end
I fully understand my first post was a mess, if you prefer I can create a new topic
One way to repeat an operator multiple times on a site would be:
julia> os = OpSum() + ("Id", 1)
sum(
1.0 Id(1,)
)
julia> for j in 1:5
os *= ("X", 1)
end
julia> os
sum(
1.0 Id(1,) X(1,) X(1,) X(1,) X(1,) X(1,)
)
Alternatively, you could define your own operator type (see Physics (SiteType) System Examples · ITensors.jl) for "op^n", though it may be subtle to do that for any n.
It does not work with my problem though.
I would want to define an operator as an OpSum()
my_op = OpSum()
for m in 1:N_sites
op += 1, "a",m
op += 1, "a†",m
end
and then apply this self defined operator multiple times
pow_of_op = OpSum()
for n in 1:N
pow_of_op *= my_op
end
In such a case pow_of_op is no more an OpSum() but its type is prod(), and I cannot use it to build an MPO
You could try:
julia> using ITensors
julia> os = OpSum() + ("X", 1) + ("X", 2)
sum(
1.0 X(1,)
1.0 X(2,)
)
julia> os2 = os * os
prod(
sum(
1.0 X(1,)
1.0 X(2,)
)
sum(
1.0 X(1,)
1.0 X(2,)
)
)
julia> Ops.expand(os2)
sum(
1.0 X(1,) X(1,)
1.0 X(2,) X(1,)
1.0 X(1,) X(2,)
1.0 X(2,) X(2,)
)
It works!
Thank you so much!
Is it in the documentation? Otherwise it seems to be something possibly useful to general users
Glad that worked!
It’s a bit of a “secret” feature. It came about because I rewrote the design of the OpSum object in terms of a general lazy operator algebra system (actually more generally a full lazy algebra system, similar to the design of Julia packages like GitHub - JuliaArrays/LazyArrays.jl: Lazy arrays and linear algebra in Julia), with the motivation of performing more general operations on OpSum, such as implementing a system for automatically creating the Trotterization of a Hamiltonian. I wasn’t quite ready to advertise that system since I wanted to see it tested out “in the wild” and be free to change some of the interface if needed. I think I would like to see it tested out more before fully documenting and advertizing it, but specific functions like Ops.expand are pretty self-explanatory and the functionality is quite clear and unlikely to change, so I would be open to documenting that function.
It is working great and we are using it a lot!
Thank you so much!
@EmilioRui - Thanks for asking this question back in March. I encountered a similar situation with a product of sum of operators (specifically a square of an OpSum() defined at each site), which led me to this thread. Thanks a lot to Matt for posting its resolution.
However, I am facing a new problem with this - I am unable to create an MPO representation of the Hamiltonian with this. Specifically, from Matt’s example above, if I define at first a HAM = OpSum(), and then starting a loop over sites ‘j’, where at each ‘j’ I have something like the aforementioned os2, then how are you doing HAM += os2 in a manner such that H=MPO(HAM,sites) works (where ‘sites’ is whatever your siteinds are) ? I am being shown either a method matching error for MPO or “Overload of “op” or “op!” functions not found for operator name os2…” depending on a few things that I tried around.
@mtfishman - any suggestion ? I can post a code snippet if needed in case what I said above is not clear.
I apologise if the resolution is something embarrassingly simple and straightforward - I am new to ITensors (and even Julia) and this is just the beginning of my first MPS project 
Hi adityab999, Does this code do something similar to what you need?
using ITensors
let
N=5
os2=OpSum()
for i in 1:N-1
osi=OpSum() + ("X", i) + ("X", i+1)
os2+=Ops.expand(osi*osi)
end
@show os2
s=siteinds("S=1/2",N;conserve_qns=false)
HAM=MPO(os2,s)
end
Kind Regards
Jan
Hi adityab999,
I was having similar issues and Jan solution worked for me.
As a note you can also create the square an operator by applying it twice.
using ITensors
let
N = 5
s = siteinds("S=1/2", N; conserve_qns=false)
H_vec = []
for i in 1:N-1
osi=OpSum() + ("X", i) + ("X", i+1)
H_i = MPO(osi, s)
H_i_square = apply(H_i, H_i)
push!(H_vec, H_i_square)
end
HAM = sum(H_vec)
end
@mtfishman Is there any reason to prefer one over the other?
@JanReimers - Hi Jan, yes that’s very much what my code is essentially trying to do. Thanks a lot really, your suggestion works without any issue.
@EmilioRui - Hi Emilio, thanks very much, your suggestion is also working equally well as Jan’s suggestion.
On small systems both are showing similar running times after starting the DMRG engine with the same parameters for both cases, perhaps if there’s any performance difference between the two it shows up for large enough systems (to check later!)
ITensor Website
·
Twitter @ITensorLib
© 2022 ITensor Collaboration