[Julia] Defining custom MPOs using op vs. op!

Question

Hi,

I have defined a custom local Hilbert space given by one electron tensored with one (truncated) boson. I conserve two quantum numbers, Nf and Sz for the electronic degrees of freedom.

I was trying to use the op function to define my custom MPOs in the same manner ElectronSite MPOs are defined in ITensor:

function ITensors.op(::SiteType"MySite", s::Index, 
                 opname::AbstractString)::ITensor
    Op = ITensor(prime(s), dag(s))
    if opname == "A"
        Op[1, 2] = A12

    elseif opname == "B"
        Op[1, 3] = B12
    return Op
end

(In the above hypothetical example, A and B have different QN fluxes, which is the case for me.)
However, this threw me a QN-conservation error:

ERROR: In setindex!, the element you are trying to 
set is in a block that does not have the same flux 
as the other blocks of the ITensor. You may be 
trying to create an ITensor that does not have a 
well-defined quantum number flux.

I'm pretty sure I'm conserving QN for both operators. In particular, when I define A and B separately using op!, ITensors does not complain:

function ITensors.op!(Op::ITensor, ::OpName"A", 
                  ::SiteType"MySite", s::Index)
    Op[1, 2] = A12
end

function ITensors.op!(Op::ITensor, ::OpName"B", 
                 ::SiteType"MySite", s::Index)
    Op[1, 3] = B13
end

What might be the issue here? Why does op! work but not op? Is op not intended to be used this way for custom sites? (But then, how does this work for in-built site sets such as ElectronSite?)

Thanks!

Answer 1

Hi, so there may be a few things going on here causing the errors you are seeing.

First of all, we recommend overloading the functions called either op! or op which take an OpName"String" type argument. In your first code example, you are instead overloading a version of op which takes an AbstractString argument which is an older style we do not encourage users to adopt. So please use the newer style as you did with your later examples.

Also, I just looked at the Electron site code and I do not see it using the older-style op which takes an AbstractString anymore. You can see the latest code here:
https://github.com/ITensor/ITensors.jl/blob/master/src/physics/site_types/electron.jl

So perhaps you are looking at an older version of the Electron site code?

Finally, there is an up-to-date example of making a custom site type and custom operators with QNs here: http://itensor.org/docs.cgi?vers=julia&page=formulas/sitetype_qns

Best,
Miles

Comments

I'll add that another potential issue is that:

    Op = ITensor(prime(s), dag(s))
    Op[1, 2] = A12

is likely what is causing the flux error, since by default the constructor `ITensor(prime(s), dag(s))` makes an ITensor with zero flux (but the element you are trying to set is nonzero). Instead, you should use the constructor:

    Op = emptyITensor(prime(s), dag(s))
    Op[1, 2] = A12

which makes an empty ITensor with no specified flux, which then gets a nonzero flux when the element is set. I also agree with Miles that you should use the updated syntax.

---

Hi Miles and Matt,

Thanks for your responses, they were very helpful. Changing the constructor from ITensor to emptyITensor seemed to instantly fix the QN conservation issue - yay!  

Miles, apparently I was using an older convention, thanks for pointing that out. However, this convention happens to work very nicely for my code. This is because in my custom type, I have to define quite a few operators, and it's easier if I can overload them all in one go via conditional statements, rather than having separate op! overloads for each of them (which would involve a ton of repeat code). This was straightforward in the older convention, where I passed a string opname to the function and could then deal with it case by case.

Is there a similar way to do this using the new convention? Can one somehow overload op as op( ::OpName"opname",...) where opname is a string variable like before? I hope I am phrasing that clearly.

Thanks again.

---

Hi, so do you mean can you still call a function of the form
op(::SiteType"MySiteType", s::Index, name::AbstractString)
?

The answer is yes you can, and that older interface is still supported though we were thinking of deprecating it in the future just to simplify the code and reduce the number of options.

From our perspective, the style you like was the one that required a lot of code, because one had to make a huge if...elseif...elseif...end block with many operator definitions all cascading down the page. If you think of the function declarations in the new style as a kind of "if" statement then actually the number of lines of code is not much more in the new style, unless you count lines that are just "end" and blank lines.

But anyway it's definitely not unreasonable to like the older style. Please go ahead and use it and just be aware that (a) we might deprecate it in the future but the code will issue a warning and still work in that case and (b) the new style allows us to do some more ambitious things too like supporting multi-site operators which is another reason we prefer it.

Best,
Miles

---

Could you give an example of what you are doing where the older style of overloading `op`/`op!` allows more code reuse? Perhaps we can give you some pointers for how to use it in a better way. For example, the initial example you showed could be rewritten more simply by overloading `op!` with:

    # This makes it so you can overload without writing `ITensors.op!`
    import ITensors: op!

    function op!(Op::ITensor, ::OpName"A", ::SiteType"MySite", s::Index)
        Op[s' => 1, s => 2] = A12
    end

    function op!(Op::ITensor, ::OpName"B", ::SiteType"MySite", s::Index)
        Op[s' => 1, s => 3] = B12
    end

which I think is pretty minimal code (and like Miles said, in general we have found this style to be a lot cleaner and easier to maintain compared to the older style where you need a long list of `if ... elseif ...` statements). If you have a common operation you are performing across multiple operators, you could include a helper function:

    function mysite_op_helper!(Op::ITensor, s::Index)
        # Perform some common operation that works for multiple op definitions
        Op[s' => 1, s => 1] = sqrt(2)
    end

    function op!(Op::ITensor, ::OpName"A", ::SiteType"MySite", s::Index)
        mysite_op_helper!(Op, s)
        Op[s' => 1, s => 2] = A12
    end

    function op!(Op::ITensor, ::OpName"B", ::SiteType"MySite", s::Index)
        mysite_op_helper!(Op, s)
        Op[s' => 1, s => 3] = B12
    end

This example is a bit contrived, but I think it helps illustrate the added flexibility of working with function overloads.

---

A less contrived example would be if one `op!` definition could be derived from another one, for example you could define a `sigma_z` operator from an identity operator by flipping the sign of the last entry:

    function op!(Op::ITensor, ::OpName"Id", ::SiteType"MySite", s::Index)
        for n in 1:dim(s)
          Op[n, n] = 1.0
        end
    end

    function op!(Op::ITensor, on::OpName"sigma_z", st::SiteType"MySite", s::Index)
        op!(Op, OpName("Id"), st, s)
        d = dim(s)
        Op[d, d] = -Op[d, d]
    end

This is also slightly contrived, but in practice this is useful for more complicated situations such as defining gates with parameters, where many gates can be defined in terms of setting specific parameters of a more general gate (such as the phase shift gate: https://en.wikipedia.org/wiki/Quantum_logic_gate#Phase_shift_(%7F'%22%60UNIQ--postMath-00000037-QINU%60%22'%7F)_gates).

Anyway, these situations may not be applicable for your use case, but I'm just bringing them up for general reference and to point out that the newer function overload style has many potential advantages over the older style.

---

Hi Miles and Matt,

Sorry for the slow response and thanks for your detailed answers. I guess it made sense for me to use op because I had several operators whose entries were best filled by a triple-for loop. This is because my custom site consists of two fermions and a boson, and the matrix entries depend on the occupation numbers of all three.
So, it seemed better to have a single triple-for loop and put the if-else conditionals within that. However, I might just switch to the op! convention as you recommend.

Actually, a different issue has come up with my code and I was wondering if you could help with that. As I mentioned, I have two spinless fermions and one boson at each site. I am having trouble implementing anti-commuting operators on this set.

I followed the format given in, e.g., https://github.com/ITensor/ITensors.jl/blob/master/src/physics/site_types/electron.jl.

I have defined a parity string operator and labeled it "F". This is nothing but the product (1-2*n1)(1-2*n2) where n1 and n2 are the number operators for the two fermions at a site.
I have also specified which operators are fermionic using the has_fermion_string function. However, when I try making an MPO using autoMPO, I'm pretty sure it isn't correctly making these operators fermionic. Is there anything I am missing? I can provide details if that would be helpful. Thanks.

---

Hi, I'd have to think about it a bit more since this seems like a fairly complicated setup, but is there a reason you put the fermions and bosons together on a single site, or can you just have a repeated 3 site pattern like:

-- Fermion1 -- Fermion2 -- Boson -- Fermion1 -- Fermion2 -- Boson --

and use the built-in `"Fermion"` type, and then implement a `"Boson"` type? That seems like it would be a lot simpler, and it would handle the fermionic statistics automatically (if you use the built-in `"Fermion"` operator definitions and AutoMPO()).

---

Thanks, Matt. Some of my primitive operators would be hard to break down into single-site operators if I implemented the sites as f1 - f2  - b1 - ...

As an example, the operators might involve matrix elements such as sqrt(nf+nb) where nf and nb are the fermionic and bosonic numbers at some site. I'm don't think I can write this as a product of single-site observables on the fermions and bosons.

The only way I can see is to use projectors onto nf=0 and nf=1 and define each observable as a sum over these projectors, namely, sqrt(nb)*Proj(nf=0) + sqrt(nb+1)*Proj(nf=1). In theory this is possible but since I have many such operators, it would probably increase the size of my code considerably.

For a custom set, is it not sufficient to define the parity string "F" and specify which operators are fermionic? Or is there an additional hidden assumption?

---

I see. I guess I don't particularly see why it would be much more complicated to rewrite it in terms of factorized operators with projectors. It seems like at least for defining the operators/site types themselves it should be much simpler (just using standard "Fermion" and "Boson" site types with pretty simple operator definitions), but I guess you are saying that the MPOs themselves become considerably more complicated? Also keep in mind the DMRG may benefit from keeping separate sites since you may be able to do more truncation. With 2 fermions and 1 boson per site, the local physical Hilbert space could easily get quite large, like with 4 bosons per site you already have a local physical dimension of 16.

It's my understanding that you should only need to define the "F" operator and specify which operators are fermionic with `has_fermion_string`, though perhaps Miles can confirm that.

Maybe you could try debugging by implementing certain simpler operators in two ways, one with composite sites and one with separate sites `-f1-f2-b-`. Then you could compare generating simple fermionic operators like `Cdag_i C_j` for both cases and see which operators may be failing.

---

Yes, the MPOs would become more complicated, although your second point about more efficient code with smaller sites makes me reconsider my choice, so I might try separate sites to see if it works out better. Thank you!

(As for the issue of fermionic anti-commutation, I have for the time being explicitly included Jordan-Wigner strings which seems to be working fine.)