Dear Thomas,
to answer your question :
"Maybe I missed the point, but do you have a special (fast) implementation of
the energy current and the time-dependent hopping operator that uses only
the 'old' kwant.operator._LocalOperator?"
Yes, I do have the specialized versions. This is what I was talking about in part I. of my rather long mail. For example, you could compare in the file 'operatorsHeatCurrent.pyx' on my github repository ( https://github.com/PhiReck/heatCurrent ) the class 'CurrentWithArbitHop', which inherits from the kwant.operator._LocalOperator, to the kwant.operator.Current class. They are the same except that an arbitrary hopping function can be used instead of the Hamiltonian. For the energy current, there is the offEnergyCurrent class in the same file with two Hamiltonians instead of one.
Concerning the question of whether they should be integrated to kwant or to tkwant: I guess I agree that a good place for the specialized currents is tkwant whereas the generalOperator would fit better into kwant. The only problem and/or inconsistency, that I see, is the fact that there would be then two different parent classes for operators in kwant, the _localOperator and the generalOperator. The former is faster since specialized, the latter simplifies the generation of new operators since there is no need for an _operate method in each new operator, as opposed to the _localOperator.
So I don't know what the best place for the generalOperator would be.
Best regards,
Phillipp
________________________________
Von: Thomas Kloss [kloss@itp.uni-frankfurt.de]
Gesendet: Freitag, 22. Februar 2019 17:04
An: RECK Phillipp
Cc: kwant-discuss@kwant-project.org
Betreff: Re: Heat and energy currents using tkwant
Dear Phillip,
thank you very much for your contribution.
I agree with Christoph that the two operators for the energy current and the time-dependent hopping elements
are specific to time-dependent problems and should be integrated into tkwant.
I also like the idea of the general operator module to create more complex
objects. However, since this is not limited to time-dependent systems only and could be of greater interest,
I think that kwant would be a better place for this module.
This would not exclude the general operator module from being used for time-dependent problems,
as one can always cast the time argument in the calling signature of the operator
in a similar way as we do it right now with the three standard kwant operators.
Maybe I missed the point, but do you have a special (fast) implementation of
the energy current and the time-dependent hopping operator that uses only
the 'old' kwant.operator._LocalOperator?
Best,
Thomas
On Wednesday, February 20, 2019 12:41 CET, RECK Phillipp wrote:
Dear Thomas, dear kwant community,
the purpose of this mail is to discuss the technical details to calculate time-dependent heat and energies currents using tkwant (extension of kwant which allows to create a time-dependent Hamiltonian and compute the time-dependent observables like densities and currents) and if/how our additional modules to do so can be integrated into tkwant. As far as I know, tkwant is not yet officially released but publicly available.
Upon request of Christoph Groth, we will share our discussion with the kwant community and invite anyone to participate. As far as I know, the kwant mailing list was not (only) meant to be a Q&A platform for the community to ask questions to the developers but also for the whole community to discuss topics. Moreover, with the recent retiring of the developer mailing list, it was recommended to use the general mailing list also for development discussions. Still, this mail is mainly addressed to Thomas, which is why I apologize in advance for code details which are not explained and thus not understandable for everyone.
That being said, let's come to the physics.
I. Heat and energy currents in time-dependent systems
( https://github.com/PhiReck/heatCurrent )
We want to calculate time-dependent heat currents into a lead. The heat current is defined as I_h = I_E - μ I_N, with I_E being the energy current, μ is the chemical potential in the given lead and I_N is the particle current. Since μ is assumed to be known and I_N(t) can be calculated directly by (t)kwant with the help of the kwant.operator.Current class, the only missing part is the energy current.
In a static setup, the energy current is obtained easily by an additional factor of the energy in the integrand of the Landauer-Büttiker formula (e.g. for two leads: I_E =e/h \int dE E T(E) [f_L(E)-f_R(E)] ). However for time-dependent systems, the energy is not conserved anymore which makes things more complicated.
I do not want to explain how tkwant works here, since this would take too long. The basic idea is that the lesser Green's function can be expressed in terms of the time-dependent scattering states. The integration over the energy in a Landauer-Büttiker-type formula is done by an instance of the kwant.manybody.Solve class. The provided operator, e.g. kwant.operator.Current or kwant.operator.Density, tells this solver which observable is to be calculated.
So far, the kwant.operators available are:
- 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 instance, for spin densities and spin currents, M would be the appropriate Pauli matrix (at least if spin is the only additional degree of freedom).
For the purpose of time-dependent heat and energy currents we need additionally:
- Energy Current: 'bra_i H_ij H_jk ket_k - c.c.'
- Explicit time dependence of hopping: '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 ( https://github.com/PhiReck/heatCurrent ) and tested them successfully in comparison to the analytical Keldysh-Green's function approach in the resonant level model. As far as we see, the calculation of the time-dependent wave function takes distinctly more time than the computation of the energy current (factor of >10) such that the calculation time is not really affected by the energy current.
(Note that we did not only implement the operators above, but also utility functions that makes it as easy for the user to calculate particle currents in tkwant as the heat current into/from a lead.)
II. A general operator module
( https://github.com/PhiReck/generalOperator )
Moreover, since the operators above all look rather similar, we generalized the operator module and discussed it already on the development mailing list. With this general operator, one is able to calculate any term of the following arbitrary form:
- generalOperator: 'bra_a O1_ab O2_bc ... O25_yz 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, On_xy are user-defined operators (if x==y, it is also an onsite operator). '+/-/0' is meant to imply that the complex conjugate can either be added subtracted or not taken into account at all, depending on what the user wants to do.
In practice, we would first initialize each operator as
O1 = generalOperator.Operator(syst, O1-FUNCTION, where=where_O1)
O2 = generalOperator.Operator(syst, O2-FUNCTION, where=where_O2)
...
with On-FUNCTION being the function that defines the operator, e.g. syst.hamiltonian or any other function that takes 1-2 sites and returns a matrix with dimensions given by the additional degrees of freedom (spin, orbitals, ...) at the sites at hand. Moreover, each operator has to know where it is to be considered/computed, specified normally as a list of hoppings (where_On).
The product is achieved by
prodOp = generalOperator.Op_Product(O1, O2, ..., complex_conj, const_fac)
'complex+conj' takes care of the +/-/0 option above and 'const_fac' is a constant factor which is multiplied to the result (in case of the current, it would be the imaginary unit 1j).
To tell tkwant to calculate the matrix product with the wave functions as shown above, it is used the same way as for the kwant.operators:
result = prodOp(bra,ket)
which returns a list of results for all possible site connections.
Speaking of the site connections (or as I call it: the path finding), note that now, the different where_On have to be matched: since we want to have something like "O1_ab O2_bc", we have to make sure that only hoppings in where_O1 are connected to hoppings in where_O2 which have the same site 'b' at the according position. This path finding is done automatically in the initialization of the class generalOperator.Op_Product.
The path finding is the only new conceptual thing, but the user does not have to worry about it which is why I don't want to go into details here now. (There is a documentation in the git repository if someone is interested.)
How to create a new operator from the generalOperator at the example of the Energy Current :
This is a copy of the source code, where I omitted the check_hermiticity option and the sum option for better readability.
class offEnergyCurrent:
def __init__(self, syst, where1, where2):
self.hamil1 = Operator(syst, syst.hamiltonian, where1)
self.hamil2 = Operator(syst, syst.hamiltonian, where2)
self.offEn = Op_Product(self.hamil1, self.hamil2, complex_conj=-1, const_fac=1j)
def __call__(self, bra, ket=None, args=(), *, params=None):
return self.offEn(bra, ket, args, params=params)
Long story short: With the help of the generalOperator, it is very easy to create new operators (just define the individual operators and take the product via generalOperator.Op_Product). Before, for each operator inheriting from the kwant.operator._Localoperator, an individual "_operate" method was needed, although they are very similar in all 5 cases from above. At the moment, I am not sure if anyone ever needs to create new operators. So far, for kwant, it seems that the 3 kwant.operators have been enough (and would have been enough for our purpose of thermoelectric properties). With t-kwant however, I have the feeling that certain operators will be more complicated as compared to their time-independent counterparts as we have seen in the case of the energy current. Therefore, this general operator might be of use in the future.
However, the larger universality of the general Operator comes at a price: Calculations are slower, which is mostly due to the fact that recursive functions are called for each single site combination. For the three examples of the kwant.operators, we see that the __call__ of the same operators initialized with the general operator takes roughly 3-5 times longer. This means that the general operator is not (yet?) suited to replace the 'specialized' ones, but at least it will help potential users to easily define and use new operators: they won't have to write and test an operate method in the 'old' kwant.operator._LocalOperator way.
So this is the current state of progress. We hope to merge both the 'specialized' heat/energy current as well as the general operator with tkwant soon.
@ Thomas: Do you see any problems or do you have any concerns about merging it with tkwant? Maybe I should have added that all the features of the kwant operators should be able to use, e.g. bind or functions for where. For instance, all the kwant examples from the tutorial also work fine with the general operator.
Best regards,
Phillipp Reck