Dear kwant-developers,
just a quick update. I implemented the following changes suggested by Joseph:
- instead of alternating onsite-hopping operators, allow for an arbitrary
string of operators
- instead of having the wherelist and the operator list separated, the
operators are directly initialized with their corresponding where
( e.g. sz = kwant.operator.onsite(syst, sigma_z, where=circle) )
It yields the correct results and all features needed for operators in the
kwant tutorial (e.g. bind) are working as expected.
However, it has still to be optimized in the sense that onsite operators can be
treated specially since they do not move in position space, as mentioned by
Joseph in his last mail.
(It's not yet on github. If you want to see the code, please let me know.)
As a rather short and unimportant note, I would like to mention that I was not
precise in my previous mails when talking about the hermitian conjugate terms
in the operators, e.g. the Current. Actually, we need by definition the
'reversed order term', e.g. for the Current: bra_i^\dagger H_ij ket_j -
ket_j^\dagger H_ji bra_i
In case of bra=ket and hermitian operators like the Hamiltonian, this reversed
order term is equivalent to the hermitian conjugate, but for the general case,
it makes a difference. It has been implemented so far correctly, i.e. with the
reversed term (both in kwant.operator and in my generalOperator), so it is not
of further importance, but I only wanted to clarify it here.
Coming to Joseph's question of his last mail ("'why do we need to calculate the
"paths" ahead of time?") and a previous mistake of mine: As opposed to the
sparse matrix products which are O(N), finding the right paths is of order
O(N^2), as it was expected by Joe in his first mail when he said that most of
the time is spent in finding the paths (at least for large systems).
In my reply back then, I stated that finding the paths would be O(N) but I
don't remember why I was thinking that it's linear.
Let's say we want to multiply 2 operators, op1 and op2, with different 'where',
where1 and where2. At the moment, for each hopping in where1 (O(N)), we go
through each hopping in where2 to see if they "belong to the same sites" (2nd
site of hop in where1 is the same as 1st site of hop in where2), which is in
total O(N^2). Maybe there is a way to find the paths faster by first ordering
the where-lists in a desired way and then finding the paths, similar to what
kwant.operator._bisect does for the site_ranges.
(Side note: for finding the paths, it becomes important that onsite-operators
do not couple different sites and therefore, when multiplying an
onsite-operator with a hopping-operator, we do not have to find the new paths
and the O(N^2)-step can be omitted. Again thas was already mentioned by Joseph.)
Joseph's related question in his last mail was: "Why not just evaluate the
matrix elements of the constituent operators and do a series of
sparse-matrix-dense-vector products when the wavefunction is passed in?" He
gave also a hint himself, that it might be related to the fact that the
operators that we consider cannot be expressed directly as matrices, but they
are 3- (onsite) or 4-tensors. Only when choosing a given site or hopping, we
recover N_orbs1XN_orbs2-matrices.
Although I cannot give a clear answer since I'm not sure that I completely
understand the question, I think that his hint goes in the right direction. Let
me add that for the heat current for instance, we have terms of the form:
sum_{i,j,k} Psi_k H_ki H_ij Psi_j
(in our case, i and j are sites in the lead and k is a site in the scattering
region, but that is not important for the example.)
We *cannot* just evaluate H_ki at all corresponding hoppings in where1 and H_ij
at the hoppings in where2 and then multiply them, because the site 'i' has to
be the same in the two Hamiltonians. Instead, we have to know which hoppings in
where1 have the same 'i' as the hoppings in where2. That's why we have to
compute the paths either before ( __init__) or during the __call__-methode
--- since it is O(N^2) more likely in __init__.
Again, I'm not sure if I understood the question right but hopefully my answer
helps in answering it.
Thank you for your input and helpful suggestions :-)
Best regards,
Phillipp
