Automatic JWT handling in `AutoMPO` and `op` functions

ITensor Julia Questions
, ,
1

Hi all,

I’m working on a fermionic model which I create the hamiltonian using AutoMPO and siteinds("Fermion", N; conserve_qns=true). Formulation of the hamiltonian is not an issue but I want to compute some observables which includes applying an operator which has "c" or "cdag" alone. I’ve been going through the documentation and other Q&A on this topic and I’ve seen that only AutoMPO and OpSum functions support proper JWT handling in the background, however having c/cdag alone is not supported (code raises an error saying parity-odd ops are not supported). So I tried to create the operator as follows:

prod(-op(sites, "F", i) for i in 1:j-1) * op(sites, "cdag", j)

However, this makes everything very slow during computation compared to the same code I wrote for Qubit where I use MPO instead of op function e.g.

os = OpSum()
arr = vcat([("Z", qubit) for qubit in 1:j-1]...)
arr = collect(Iterators.flatten(arr))
push!(arr, "S+")
push!(arr, j)
res = add!(os, 1, arr...)
MPO(res, sites)

Note that creating JWT string with op also raises the following warning:

Contraction resulted in ITensor with 14 indices, which is greater
            than or equal to the ITensor order warning threshold 14.

I am not sure if this is related to my problem.

In this question and corresponding gist I saw that author is using op("Cdagup", sites, L) for an electron system. This makes me think that I am creating a redundant Z chain in the code. Would you mind confirming which functions take care of JWT in the background and which ones does not?

Thanks
Best regards

2

This creates a dense ITensor, with j indices (so exponentially large in j), while the MPO is more efficiently represented - which is why the code is so slow.
I would instead take the the tensors in your prod and use them in apply which will be more efficient

apply([[op(sites, "F", i) for i = 1:(j-1)]; [op(sites, "cdag", j)]], psi)

This would create c^\dagger_j|\psi\rangle, which you then could measure e.g. inner products with. You could also make an MPO of c^\dagger_j just as you did with the Z operators:

os = OpSum()
add!(os,1, ([[("F",i) for i = 1:(j-1)]; [("cdag",j)]]...)...)
Oj = MPO(os,sites)

to have a style closer to \langle \psi|O|\psi\rangle (inner(psi,Oj,psi))

1 Like
3

thanks @ryanlevy thats actually very helpful. Your first proposal

apply([[op(sites, "F", i) for i = 1:(j-1)]; [op(sites, "cdag", j)]], psi)

works but for the second one

os = OpSum()
j = 5
add!(os,1, ([[("F",i) for i = 1:(j-1)]; [("cdag",j)]]...)...)
Oj = MPO(os,sites)

I’m getting Parity-odd fermionic terms not yet supported by OpSum to MPO conversion error. Since you are able to make it work, maybe this is due to itensor version?

4

Sorry about that! You need to define your own “Cdag” that doesnt have strings

function ITensors.op!(Op::ITensor, ::OpName"Adag", ::SiteType"Fermion", s::Index)
  return Op[s' => 2, s => 1] = 1.0
end

Then you can do

os = OpSum()
add!(os,1, ([[("F",i) for i = 1:(j-1)]; [("Adag",j)]]...)...)
Oj = MPO(os,sites)
1 Like
5

Hi @ryanlevy thanks for your responses they’ve been very helpful. How would you define the “A” operator? I’m assuming it would be very similar to how you defined “Adag.” Thank you!

6

I’m just copying what is in fermion.jl
So for C → A

function ITensors.op!(Op::ITensor, ::OpName"A", ::SiteType"Fermion", s::Index)
  return Op[s' => 1, s => 2] = 1.0
end

And as a reminder this operator will not have any strings in OpSum/correlation_matrix because we have not registered that it is a fermionic operator

1 Like
7

This topic was automatically closed 10 days after the last reply. New replies are no longer allowed.