(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.