Hi Phillip,
Thanks again for taking the time to post to the mailing list; we really
appreciate it.
I've responded to the key parts of your proposal below. I'd also be
really keen to see what you've implemented so far; do you have a way of
sharing it (link to a code repository e.g. GitHub would be ideal)?
> So far, the kwant.operators are able to calculate:
> - Density: 'bra_i M_i ket_i'
> - Current: 'bra_i H_ij ket_j - c.c. '
> - Source: 'bra_i M_i H_ii ket_i - c.c.
> where "i" and "j" are sites, H is the Hamiltonian, M denotes an
arbitrary onsite Operator (which is a matrix with dimensions given by
the number of orbitals) and bra and ket are wave functions.
>
> For the purpose of time-dependent energy currents we need additionally:
> - offEnCur: 'bra_i H_ij H_jk ket_k - c.c.'
> - CurrHdot: 'bra_i M_i dH_ij/dt ket_j + c.c.'
> These are of a similar shape as the operators above but different
enough such that the old operators cannot be used. We already
implemented them and tested them successfully, also together with tkwant.
It's good to see some more examples of the sort of calculations that are
useful within the framework of kwant's operators.
Indeed, when we first designed the kwant operators we essentially just
wanted to provide a way for users to calculate and plot densities and
currents without having to worry about slicing out the correct parts of
the wavefunction (e.g. for spin, or for sites only in a certain region
of the system).
At the time we did think about supporting more general classes of
operators, but could not think of a convincing use-case. Given this, we
went ahead and implemented the simplest thing that did what we wanted.
-------------------------
> Since the operators above all look rather similar, we thought about
generalizing the operator module. With this general operator, one should
be able to calculate any term of the following arbitrary form:
>
> -generalOperator: 'bra_a M_a O_ab M_b O_bc M_c ....... O_yz
M_z ket_z +/-/0 c.c. '
>
> Here, 'a' to 'z' denote again sites where 'z' is used as an arbitrary
level of deepness (i.e., 'z' could be equal to 'b', if only one hopping
is desired, or 'd' if three hoppings are wished). Moreover, M_x are
still onsite Operators (matrices) and O_xy are general hopping operators
(so far, only the Hamiltonian was used) defined by the user.
> Instead of the "where" list (where to calculate e.g. the current), a
list of "where"-lists (for each hopping separately) is now needed. The
code finds all connections between the first sites "a" and the last
sites "z" given by these "where"-lists, calculates the operator along
this way and returns its value for each path.
This is the heart of the proposal. I will try to re-word what you have
written in my own terms to check that I have understood.
You essentially propose to allow the definition of arbitrary "onsite"
and "hopping" operators as functions (site → matrix) and ((site, site) →
matrix) along with a list of pairs of sites ("where") that tells us
where the operator is non-zero (we guarantee that the function will only
be evaluated with the sites in "where"). For the operator to be valid on
a particular system the sites in "where" must be a subset of the sites
in the system.
In your proposal these individual operators are not separate objects,
but rather defined implicitly by lists of "onsite" and "hopping"
functions, and the "where" specification. If we are going in the
direction of allowing arbitrary operators, then it may make sense to
expose them as concepts in their own right.
You then propose to allow the construction of "generalOperators", which
are sequences of alternating onsite/hopping operators, with the
possibility to add or subtract the hermitian conjugate of the operator
sequence.
> The code finds all connections between the first sites "a" and the
last sites "z" given by these "where"-lists, calculates the operator
along this way and returns its value for each path.
So you essentially pre-compute which blocks of the constituent operators
we will need to evaluate? If this is the case I guess the complexity is
the same as actually doing the sparse matrix-matrix products, and so
most of the time will actually be spent here. Using an appropriate
format for the constituent operators will be very important.
Also, given this, is there any reason to limit the strings of operators
to be alternating onsite/hopping? The only advantage of treating these
two on a different footing is that onsite operators cannot "move you
around" in site space, so you don't have to do the general thing when
making a product of them. If we're going to have to handle strings of
hopping operators anyway I don't think it's useful to make a distinction
anymore, and it might actually simplify the implementation.
-------------------------
> To particularly calculate a desired operator, 2 additional lists,
M_list and O_list, for the onsite and hopping operators would have to be
given.
>
> Examples:
>
> To recreate the kwant.operator.Current, we would need the lists:
M_list=[M,1] and O_list=['h'] ('h' is a shortcut for Hamiltonian)
> For kwant.operator.Density: M_list = [M], O_list = []
> For kwant.operator.Source, there are two possiblities:
> M_list = [MxH], O_list = [], where MxH is the user-defined
product of the onsite Matrix M and the Hamiltonian
> M_list = [M,'h'], O_list = ['unit'], where 'unit' is a
pre-defined function which yields the unit matrix (for both sites and
orbitals)
> for offEnCurr: M_list=[1,1,1] und O_list=['h','h']
> for CurrHDot: M_list=[M,1] und O_list=[Hdot], where Hdot would be
a user-defined function that calculates dH/dt.
>
> The working principle is the similar as for the _LocalOperators, with
the difference that we need now lists of BlockSparseMatrices, since we
need a BlockSparseMatrix for each operator.
As I said above, I think that it may be more natural to define the
constituent operators and then combine them into operator strings.
Imagine something like the following:
sz = kwant.operator.onsite(syst, sigma_z, where=circle) # onsites only
H = kwant.operator.operator(syst, syst.hamiltonian) # hoppings (and
onsites)
# define (sz @ H - H @ sz): precomputes paths
A = kwant.operator.product(sz, H, herm_conj='-')
A(psi) # compute expectation
This would be exactly equivalent to your construction of 'M_list' and
'O_list', except that the operator and the associated 'where' are kept
together, rather than separately.
-------------------------
> Since we already familiarized ourselves with the kwant.operator module
when creating the offEnCurr- and CurrHDot-operators in the "old"
fashion, we have been able to implement this general operator. It yields
exactly the same results compared to the 5 operators from above.
Moreover, it should be as flexible as the "kwant.operator" (e.g.
callables or "None" as where, binding of the operators, etc.). The only
problem that we still face at the moment is its efficiency: Due to some
python overhead, it is for large systems to 10-20 times slower than the
corresponding kwant.operator. In the next week(s) we hope to enhance the
efficiency to obtain (roughly) the same speed as the kwant.operators by
properly cythonizing our new general operator module.
I would certainly say that you should carefully profile your
implementation to see where the time is being spent. As I said at the
top, I'd be interested to see the code.
-------------------------
Regards,
Joe
