(a) (b)
Figure 2. Schematics describing the algorithm to construct the
hierarchy of ADMs and compute their mappings. In this example,\(d=3\). a) Each ADM of a given tier is coupled to d ADMs of the
following tier. Indices in red are duplicates of the already added ADM
indices and are discarded. b) List of enumerated d-dimensional index
vector for tiers 0, 1, and 2. Note that vectors of a given tier follow
before the vectors of the following tier and there are not duplicate
vectors. The matrices \(n^{+}\) and \(n^{-}\) contain the indices of the
ADMs that are coupled to a given ADM via raising or lowering one of its
components by one.
The tree traversal returns an enumerated list of d-dimensional vectors.
In computing time-derivatives of the ADMs, Eqs. 5 and 13, one also needs
to quickly access the \(\rho_{\mathbf{n}^{(+)}}\) and\(\rho_{\mathbf{n}^{(-)}}\) ADMs knowing just the index of the\(d\)-dimensional vector \(\mathbf{n}\). For quick access, the
“gen_hierarchy” function also computes the mapping tables for such
quick access (Figure 2b). The mapping \(\mathbf{n}^{\mathbf{+}}\), for
instance, can be thought as a \(N_{\text{adm}}\times d\) matrix, such
that its element \((i,\ j)\) contains a sequential index (scalar
integer) of the \(\rho_{\mathbf{n}_{j}^{+}}\) , if the analogous
sequential (scalar) index of the vector \(\mathbf{n}\) is \(i\).
Likewise, the element \((i,\ j)\) of the matrix\(\mathbf{n}^{\mathbf{-}}\) contains the sequential index of the\(\rho_{\mathbf{n}_{j}^{-}}\) such that the index of the vector\(\mathbf{n}\) is \(i\). In other words, if we start with the i-th
vector, \(\mathbf{n}\), and increment its j-th element by one, the
element \(n^{+}\left[i\right][j]\) will
return the index of such a multi-component vector in the overall set of
vectors. If there is no such vector (e.g. due to the truncation of the
hierarchy), the matrix elements are set to -1. This concept is better
understood from a careful examination of Figure 2b.
The second group of C++/Python functions includes: “compute_deriv_n”
which computes the time-derivative of a given ADM, n, according to Eqs.
5 or 13, “compute_heom_derivatives” that performs the computation of
time-derivatives of all ADMs (this is where the openMP parallelization
is currently involved), and “filter” function that updates the flags
which determine whether each of the ADMs should be used in computing the
right-hand side of the ADM time-derivatives in Eqs. 5 and 13. In
addition, the Python-exposed function “compute_heom_derivatives” has
a particular importance in constructing the solver of the HEOM.
The integration of the HEOM is done via a general-purpose 4-th order
Runge-Kutta integrator (RK4) method, implemented in Libra in the
“libintegrators” sub-library. It solves a set of differential
equations that can be represented in a matrix format as:
\(\dot{X}=f(X;params)\). (16)
One of the parameters the RK4 function accepts is a function itself, the
one that computes the time-derivatives of the ADMs. Passing functions as
arguments is one of the convenient features of Python, since the
functions are treated as objects and are similar to other types of
arguments in this regard. Unlike C++, where one needs to specify the
signature of the function that is being passed, Python doesn’t require
such extra care. Of course, the RK4 function that calls the argument
function would need to know the signature of the latter as well as how
to unpack the returned results, so the developer needs to ensure the
consistency of the expected function call invocation with the function’s
signature and expected output. However, these details can be hidden in
the function being called, whereas the interface of the RK4 function
won’t change.
Since every function passed to the RK4 general-purpose solver may have a
varied number of parameters, the RK4 function also expects a dictionary
of parameters. These parameters would then be passed to the RK4
argument-function. In the present context, the
“compute_heom_derivatives” is the function that is being passed to
the “RK4” solver. The parameters of “compute_heom_derivatives” are
therefore passed to the “RK4” function a dictionary of arguments. This
dictionary contains parameters such as Matsubara frequencies, Matsubara
coefficients, electron-phonon coupling matrices, and so on. The
parameters dictionary can have extra key-value pairs, not needed
directly by the “compute_heom_derivatives”, but used to control the
execution of the RK4 algorithm. This way, one may define a common
dictionary of all parameters and pass it to many functions with the idea
that only the parameters that are relevant to the calculations will be
utilized.
The feature of passing functions as arguments of other functions is also
utilized in other Libra modules, such as in the trajectory surface
hopping approaches. There, one can define a Python function that would
compute some essential properties of the system (Hamiltonian,
derivatives, etc.) for a given input of coordinates. Such a Python
function can be passed to the dynamical algorithms (e.g. propagation)
using a common interface. The user is then free to redefine the
arguments functions to suit their needs – to either define a model
Hamiltonian or to setup a call of external electronic structure packages
and extracting the needed information from their output. In this mode,
the user doesn’t need to know the details of the implementation of the
dynamical procedures and only needs to care about the input/output
signatures and the internal computations of the function passed as an
argument.
Finally, the third group of C++/Python functions comprises the auxiliary
functions “scale_rho” and “inv_scale_rho” to convert between
pristine and scaled ADMs, Eqs. 12, and the functions “pack_mtx” and
“unpack_mtx” to transform between two representations of the ADMs.
Namely, the set of all ADMs used in HEOM calculations,\(\left\{\rho_{n},\ n=0,\ldots,N_{\text{adm}}-1\right\}\), can be
represented as a list of \(M\times M\) matrices. However, the RK4
function implemented in the “libintegrators” library expects a single
matrix of the function arguments. Because of this expected format, the
list of the \(M\times\ M\) ADMs is packed into a single\(N_{\text{adm}}M\times M\) matrix \(\rho\), which can then be passed
into the RK4 function. The function \(f\left(\ldots\right)\) appearing
in Eq. 16 is the “compute_heom_derivatives”. Internally, the input
ADMs matrix is unpacked inside of this function, and the derivatives are
computed for each ADM independently. In fact, this is where we have
added the OpenMP parallelization. The computed time-derivatives,\(\left\{\frac{d}{\text{dt}}{\tilde{\rho}}_{\mathbf{n}}\right\}\), are
then packed into a single matrix format to be used within the
general-purpose RK4 integrator.
The “scaled” ADMs and their time-derivatives are considered
fundamental in our implementation: the integration and the computation
of the time-derivatives, as well as filtering and truncation, are done
within the framework of the “scaled” variable, with the “unscaled”
ones being a special case. The unscaled variables are only updated
before storing them to the results (memory or on disk), but do not
directly participate in the dynamics.