Generalized `Operators`

[zerothi]

Describe the feature

Operators are a general way to perform operations based on
inputs.

For instance, currently we do:

es = H.eigenstate(...)
vel = es.velocity()

The velocity is actually an operator, and we might wish
to perform other operators on the eigenstate.
Would it make more sense to change operators to act on
the objects them-selves? It comes with some abstraction, but
also some more interesting automatic handling, and I think added
benefit down the road (way more complicated operators).

I think we could potentially save some things here.

For instance:

import sisl.physics.operators as si_op
op = si_op.anticommutator(si_op.SpinX, si_op.Velocity) / 2
diag_elements = op.diagonal(es[, es])
matrix_form = op.matrix(es[, es])

This would also make many of the current methods a bit more generic.
And would allow simpler usage of them.

This also brings up the question on how we should deal with bra and ket objects.

Should we do:

si_op.bra(es)
# or
es.bra
# or?

Having these one could do PDOS (unexpanded to energies) with something
like this:

op = si_op.SpinIdentity # automatic handling of spinors (1 or 2)
PDOS_charge = es.bra * op.ket(es)
PDOS_x = es.bra * si_op.SpinX @ es.ket
PDOS_y = es.bra * si_op.SpinY @ es.ket

So how do we want this to look like?

What we need to decide is how the operator should deal with these points:

1. name of diagonal method (diagonal seems obvious)
2. name of the full matrix method (diagonal + off-diagonal), matrix
3. name of the full matrix form that is similar to PDOS, i.e. state.conj() * state

How should bra, ket and matrix operations look like? It might seem weird to reuse @ and *.

Should we aim for the code to look and feel like the bra-ket notation, or what should we make it
look like?

[ialcon]

So, if I get it right, this is a way to move operators from being methods of eigenstate objects, to being actual operators that are external to eigenstates and, instead, act on them. Generally I would say this is a good idea because, as you mentioned @zerothi, it will create an organization where operators will be more general and simpler to use. Additionally, it will make the code, as you also mentioned, more in line with the fundamental way in which quantum mechanics is structured (e.g. the bra/ket notation), so it will also make more sense in the community. So, generally I would say that this approach sounds to be very reasonable.

Now, some comments/questions on your discussion. First of all, what is this? es[, es] - I'm sorry but I've never seen nor used this type of notation. Is it a way to specify the bra/ket representations of the eigenstate?

Should we do:

si_op.bra(es)
# or
es.bra
# or?

I would say, here, that es.bra is more natural, because this is just the "row" representation of the eigenstate, so I think that it makes more sense to implement it as a method of eigenstate, rather than as an external operator acting on eigenstate.

Having these one could do PDOS (unexpanded to energies) with something like this:

op = si_op.SpinIdentity # automatic handling of spinors (1 or 2)
PDOS_charge = es.bra * op.ket(es)
PDOS_x = es.bra * si_op.SpinX @ es.ket
PDOS_y = es.bra * si_op.SpinY @ es.ket

I like more doing something like si_op.SpinX(es.ket) than si_op.SpinX @ es.ket because, in your example, these are supposed to do the same thing, right? I think the first is more compact, and it is also more in line with the idea that we are giving an eigenstate to the operator that acts on it. I also think that, unless it was necessary for other things, it would be convenient to have one unique way to do all this - so, either si_op.SpinX(es.ket) or si_op.SpinX @ es.ket, but not both. Finally, couldn't one simply do:

PDOS_charge = op(es)
PDOS_x = si_op.SpinX(es)
...

as an equivalent way to do es.bra * op.ket(es), es.bra * si_op.SpinX(es.ket), ... ??

[zerothi]

Great! Thanks for your inputs!

Now, some comments/questions on your discussion. First of all, what is this? es[, es] - I'm sorry but I've never seen nor used this type of notation. Is it a way to specify the bra/ket representations of the eigenstate?

It is a programming way of saying arg1[, arg2] that arg1 is mandatory (call like func(arg1)), and arg2is optional (call likefunc(arg1, arg2)`).

Should we do:

si_op.bra(es)
# or
es.bra
# or?

I would say, here, that es.bra is more natural, because this is just the "row" representation of the eigenstate, so I think that it makes more sense to implement it as a method of eigenstate, rather than as an external operator acting on eigenstate.

Ok, I that makes sense. :)

Having these one could do PDOS (unexpanded to energies) with something like this:

op = si_op.SpinIdentity # automatic handling of spinors (1 or 2)
PDOS_charge = es.bra * op.ket(es)
PDOS_x = es.bra * si_op.SpinX @ es.ket
PDOS_y = es.bra * si_op.SpinY @ es.ket

I like more doing something like si_op.SpinX(es.ket) than si_op.SpinX @ es.ket because, in your example, these are supposed to do the same thing, right? I think the first is more compact, and it is also more in line with the idea that we are giving an eigenstate to the operator that acts on it. I also think that, unless it was necessary for other things, it would be convenient to have one unique way to do all this - so, either si_op.SpinX(es.ket) or si_op.SpinX @ es.ket, but not both.

I agree, I like the first one, better, having thought some more, however, see further down! ;)

Finally, couldn't one simply do:

PDOS_charge = op(es)
PDOS_x = si_op.SpinX(es)
...

as an equivalent way to do es.bra * op.ket(es), es.bra * si_op.SpinX(es.ket), ... ??

No, because you might want to do something different with si_op.SpinX(es). And here, perhaps users can be confused about si_op.SpinX(es.ket) being different than si_op.Spin(es)? I think we should probably leave the interface open for what happens when passing a state object, since you don't know if a ket or a bra is asked for.

I have also looked more closely into what sympy is doing:
And there they have this:

1. es.bra * es.ket == inner-product
2. es.ket * es.bra == outer-product
3. es.bra * si_op.SpinX * es.ket is the way to perform operators. It is a bit weird (going againts the way numpy works) i.e. * for scalar multiplication and @ for matrix-multiplication.

[zerothi]

Otherwise we should have specific methods to call the necessary methods, in this way arguments can be added as needed. Using * and @ will not allow one to pass arguments.

So then it would be something like:

import sisl.operator as si_op

op = si_op.SpinX()
es = H.eigenstate()
value = si_op.innerprod(es.bra, op, es.ket)
value = si_op.outerprod(es.ket, op, es.bra)
# since the method is generic, and has order specific bra-ket notation, one should also be able to do
value = si_op.innerprod(es op, es)
value = si_op.outerprod(es, op, es)

This is a bit opposite of what sympy does, they actually implement OuterProduct as an operator. However, they only allow 2 arguments.

I just don't know how we should call the es.bra.T * op @ es.ket which in the current codes we name as project=True.

There are a couple of things that is necessary:

1. should we have a specific innerprod operator (we might have it anyways),
2. or should we have methods on operators,
   • inner
   • outer
   • inner_matrix, or inner(..., matrix=True)
   • inner_?, or inner(..., ?=True)

[ialcon]

Personally, I think there shouldn't be an innerprod operator... So, to get inner/outer products of operators one could simply do:

import sisl.operator as si_op

op = si_op.SpinX()
es = H.eigenstate()
value_IN = op.inner(es)
value_OUT = op.outer(es)

The use of iner/outer methods already specifies that we are internally doing es.bra * op(es.ket), so there is no need to specify anything else than es in the .inner/.outer calls. I think in this way it is quite unambiguous what one is doing, and pretty simple as well. Finally, to get the inner product between two vectors we could simply do:

density = es.bra * es.ket
overlap = es_1.bra * es_2.ket

I think this is also simple enough to, instead, require an independent inner/outer operators... Unless you would like to be able to provide more arguments to the bra-ket multiplication, which then may require an actual operator. However, as of now, I do not see which additional arguments may be necessary for such type of bra-ket inner/outer multiplications beyond eigenstate-properties (e.g. a specific spin-channel) which may already be specified within each eigenstate object (e.g. es(spin=1).ket).

[ialcon]

ah, and I guess that in the above proposition we could also have:

value_IN = op.inner(es_1, es_2)  # --> es_1.bra * op(es_2.ket)
value_OUT = op.outer(es_1, es_2) # --> es_1.ket * op(es_2.bra)

[zerothi]

Personally, I think there shouldn't be an innerprod operator... So, to get inner/outer products of operators one could simply do:

import sisl.operator as si_op

op = si_op.SpinX()
es = H.eigenstate()
value_IN = op.inner(es)
value_OUT = op.outer(es)

The use of iner/outer methods already specifies that we are internally doing es.bra * op(es.ket), so there is no need to specify anything else than es in the .inner/.outer calls. I think in this way it is quite unambiguous what one is doing, and pretty simple as well. Finally, to get the inner product between two vectors we could simply do:

density = es.bra * es.ket
overlap = es_1.bra * es_2.ket

I think this is also simple enough to, instead, require an independent inner/outer operators... Unless you would like to be able to provide more arguments to the bra-ket multiplication, which then may require an actual operator. However, as of now, I do not see which additional arguments may be necessary for such type of bra-ket inner/outer multiplications beyond eigenstate-properties (e.g. a specific spin-channel) which may already be specified within each eigenstate object (e.g. es(spin=1).ket).

One thing is that we will necessarily have the si_op.inner operator because in some cases you really want additional arguments. As I mentioned above:

   * `inner`
   * `outer`
   * `inner_matrix`, or `inner(..., matrix=True)`
   * `inner_?`, or `inner(..., ?=True)`

here the matrix forms returns orbital-resolved inner products (in numpy terms this would be: bra * ket.T). We will quite often need this. Actually, come to think of it, I think this way of using the operators may be confusing because $\langle\psi|\psi\rangle=s$ where $s$ is a scalar.
In terms of what happens in bra-ket notation then: $\langle \psi | A$, $A|\psi\rangle$ and $|\psi\rangle\langle\psi|$ are matrix multiplications.
While the inner product is not a matrix multiplication, but rather the trace of the matrix-multiplication.

Perhaps it would be smartest to start with only operators and methods (no * or @ overloading) until we find an intuitive way of dealing with these things.
My main concern is that users will be confused as to when we use * and when we use @ because es.bra * es.ket is really a trace of a matrix multiplication.

ah, and I guess that in the above proposition we could also have:

value_IN = op.inner(es_1, es_2)  # --> es_1.bra * op(es_2.ket)
value_OUT = op.outer(es_1, es_2) # --> es_1.ket * op(es_2.bra)

Yes, but I think it rather should be: es_1.ket.outer(es_2.bra) for clarity.
However, I do believe that reading si_op.outer(es_1, es_2) is more inline with how you read bra-ket notation. But I think it would be ok to support both, even though you'll then have two ways of doing things. ;)

I think I have a pretty clear idea of how this should look, thanks for the inputs! (feel free to add more if you want to!)