VAMPyR—A high-level Python library for mathematical operations in a multiwavelet representation Open Access

Wavelets and multiwavelets have emerged in the past few decades as a versatile tool for computational science. Their strength derives from the combination of frequency separation with locality and a robust mathematical framework to gauge the precision of a calculation. In particular, multiwavelets have recently been employed within the field of quantum chemistry to overcome some of the known drawbacks of traditional atomic orbital (AO)-based calculations.^1–9 The code that pioneered this approach is M-A-D-N-E-S-S,¹⁰ followed by our own code, MRChem.¹¹ To date, they are the only two codes available for quantum chemistry calculations using multiwavelets.

As the development of MRChem unfolded, we found it advantageous to separate the mathematical code dealing with the multiwavelet (MW) formalism in a separate library, called MRCPP (Multiresolution Computation Program Package).¹²  MRCPP is a high-performance C++ library, which provides the tools of Multi–Resolution Analysis (MRA) with multiwavelets. It implements low-scaling algorithms for mathematical operations and convolutions and a robust mechanism for error control during numerical computations.

Once MRCPP was available, we realized that it would be useful to make its functionality easily accessible. For this purpose, we turned our attention to Python, which has become a de facto standard in the scientific community, widely taught to science students and used extensively in research. Its high-level syntax and rich ecosystem, including libraries such as NumPy,^13, SciPy,¹⁴ and Matplotlib,¹⁵ make it a powerful tool for scientific computing. Moreover, Python’s integration with Jupyter notebooks¹⁶ has facilitated interactive and reproducible research. This is in contrast to traditional quantum chemistry codes, which are implemented in lower-level languages (Fortran, C, or C++) and can pose significant barriers for non-expert programmers.

That led us to the development of VAMPyR (Very Accurate Multiresolution Python Routines),¹⁷ which builds on the MRCPP capabilities, by making them available in the high-level programming framework provided by Python, thereby simplifying the processes of code development, verification, and exploration. The purpose of VAMPyR is to allow a broader audience to make use of these tools while maintaining the original power and efficiency of the MRCPP library. This includes inheriting MRCPP’s parallelization, which employs the multi-platform shared-memory parallel programming model provided by OpenMP.^18,19

The foundations of Multiresolution Analysis (MRA) have been laid in the early 1980s, thanks to the pioneering work of Meyer, Daubechies, Strömberg, and others,^20–27 who defined for the first time the concepts of wavelet functions and nested wavelet spaces. We here provide a short survey by focusing on a particular kind of wavelets called multiwavelets.^28–32 

The matrices H⁽⁰⁾, H⁽¹⁾, G⁽⁰⁾, and G⁽¹⁾ are each of dimension (k + 1) × (k + 1). See Ref. 32 for a comprehensive derivation of these matrices. The operation illustrated in Eq. (6) is called forward wavelet transform, also known as wavelet compression. The inverse of this operation is termed backward wavelet transform or wavelet reconstruction. The operation is strictly local: the scaling function $ϕln$ and the wavelet function $ψln$ are only connected to $ϕ2ln+1$ and $ϕ2l+1n+1$⁠. This structure simplifies numerical algorithms, enabling fast implementations.

In the architecture of VAMPyR, MRCPP serves as the foundational layer, integrated into Python through Pybind11, as illustrated in Fig. 1. This setup enables seamless interoperability between Python and C++, allowing Pybind11 to automatically convert many standard C++ types to their Python equivalents and vice versa. This leads to a more Pythonic and natural interface to the MRCPP codebase.

For the development of VAMPyR, we opted for Pybind11 as our binding framework, chiefly for its robustness, maintainability, and alignment with our development objectives. Pybind11 is a lightweight, header-only library that utilizes modern C++ standards to infer type information, thus streamlining the binding process and enhancing the overall code quality.

To facilitate deployment, multiple installation options are available for VAMPyR and MRCPP. The source code for both is openly accessible and distributed under the LGPLv3 license on GitHub, specifically within the MRChemSoft organization. For users who prefer a simplified installation procedure, binary packages are also provided through the Conda package manager, which is compatible with various Linux and macOS architectures.^36,37 Each new release triggers an automated build process that uploads the binary packages to the Conda-Forge channel on anaconda.org, thereby easing the installation process and incorporating all requisite dependencies. The packages are designed to be compatible with current Python versions, ranging from 3.8 to 3.11. The installation process using Conda is straightforward:

In MRCPP, C++ template classes and functions are utilized to provide abstraction over the dimensionality of the simulation box. These templates enable the generic implementation of data structures and algorithms for problems in D dimensions.³⁸ The specialization for one-, two-, and three-dimensional problems is performed at compile time, thereby eliminating any impact on runtime performance.³⁹ 

In Python, native template constructs are absent, presenting a challenge for dimension-specific specialization. In VAMPyR, we emulate MRCPP’s generic approach by implementing dimension-dependent bindings using Pybind11. Specifically, our binding code incorporates template classes and functions similar to those in MRCPP. MRCPP’s one-, two-, and three-dimensional template classes and functions are directly bound to their corresponding dimensions in VAMPyR, as demonstrated in Listing 1. Through this approach, VAMPyR keeps the flexible design of MRCPP, despite Python’s limitations in handling generic datatypes.

While the MultiResolutionAnalysis (see Sec. IV A) and FunctionTree (see Sec. IV B) classes in VAMPyR are inherited from MRCPP, there are noteworthy distinctions between vanilla C++ and the corresponding Python classes. These differences are introduced to enhance the Python version of the MRCPP classes with the so called syntactic sugar, which improves the user experience and leads to faster prototyping. Among these enhancements are the overload of various dunder (double underscore, also known as magic) methods that have been added to the classes to extend their functionalities and make them more Pythonic.

These dunder (or magic) methods enable the use of Python’s built-in operators on FunctionTree objects. For instance, the addition (__add__) and multiplication (__mul__) operators allow for the direct use of + and * operators with FunctionTree objects. This allows us to write more intuitive code, which is easier to both read and write.

As an example, consider the addition of two FunctionTree objects, which can be expressed in natural Python syntax (Listing 2), thanks to the binding code in Listing 3.

The Python code in Listing 2 can be equivalently written in the de-sugared form in Listing 4. Although Listing 4 is more involved than simply applying an arithmetic operator, it does provide more flexibility and control, which might be necessary in some specific applications.

Here, we display some of the theoretical concepts discussed in Sec. II with practical examples using VAMPyR.

The foundation of VAMPyR’s internal operations is built upon the MultiResolutionAnalysis object. This object encapsulates essential information about the physical space being modeled, along with configurations for the scaling and wavelet basis functions as defined in the theory of MRA (see Sec. II).

For users seeking a straightforward implementation, VAMPyR provides a “standard usage” option. In this approach, users are only required to specify the box size and polynomial order. The remaining parameters, such as the choice of interpolating polynomials for the scaling basis, are automatically configured by the library. Listing 5 illustrates how to create an MRA object either specifying the start and end points or the length: the former has length 2L, whereas the latter starts in 0 and ends in L.

For users requiring more control over the computational setup, VAMPyR offers an “advanced usage” approach. This method allows for the customization of various parameters, including the world box size and the specific type of scaling basis, such as Legendre polynomials. This level of customization grants the user full control over the basis configuration,

as demonstrated in Listing 6. This flexibility allows users to choose the level of interaction with the MRA object based on their specific needs and expertise.

In VAMPyR, functions are represented using a tree-based data structure: the FunctionTree object. This structure naturally arises from the MRA framework outlined in Sec. II C.

The FunctionTree object consists of interconnected nodes organized hierarchically. The root node⁴⁰ represents the entire physical domain, and each non-terminal node or branch begets 2^d child nodes while connecting to a single parent node. Terminal nodes, or leaves, possess a parent but have no children. Here, d represents the dimensionality of the system. VAMPyR currently supports d = 1, 2, 3.

Each node corresponds to a d-dimensional box within the physical domain and encapsulates information about a specific set of scaling and wavelet functions defined therein. Specifically, each node indexed by scale n and translation vector l = (l₁, l₂, …, l_d) stores (k + 1)^d scaling coefficients and (2^d − 1)(k + 1)^d wavelet coefficients. Here, the indices n and l dictate the node’s spatial size and position, respectively.

To facilitate the implementation of function projections into scaling and wavelet spaces, VAMPyR provides the vp.ScalingProjector and vp.WaveletProjector classes. In Listing 7, we define a lambda function f, whose body can be any valid Python expression.⁴¹ This illustrates how to project a function using the vp.ScalingProjector and vp.WaveletProjector classes. The instances P_n and Q_n are created to perform the projection into the scaling and wavelet spaces. The resulting representations, denoted f_n and df_n, are fully described by their mathematical definitions in Eqs. (7) and (9), respectively.⁴² 

Figures 2 and 3 provide visual insights into the projection of an exponential (Slater) function, e^−a|r|, onto different scales. Each figure consists of two subfigures: the left-hand side displays the scaling function space [Fig. 2(a) and 3(a)], while the right-hand side illustrates the wavelet function space [Figs. 2(b) and 3(b)].

In Fig. 2, the projection onto the root scale is far from the target, and the large wavelet part indicates that a representation at a finer scale is mandatory. The vertical lines in both subfigures delineate the physical space spanned by k + 1 scaling and wavelet functions. In Fig. 3, the exponential (Slater) function is projected onto the fourth scale. The representation is now close to its target function: it is, however, not as sharp as the original exponential (Slater) function, and the wavelet projection indicates that in the cusp region, a finer representation would be required to attain high precision.

Building on the mathematical framework presented in Sec. II D, we cover here the practical aspects of adaptive projection within VAMPyR. Unlike the fixed-scale projection illustrated in Listing 7, adaptive projection makes use of the prec keyword to refine the function’s representation, as shown in Listing 8. The visual comparison in Fig. 4 between the adaptive and fixed-scale projections (previously shown in Figs. 2 and 3) reveals the power of adaptivity: a representation, which is at the same time more compact as well as more precise, is achieved. As Fig. 4 shows, nodes at finer and finer resolution are created only where it is necessary to represent the sharp features of the function. This strategy is also computationally more efficient as it guarantees that resources are used where and when needed, based on precision requirements.

In VAMPyR, arithmetic operations on FunctionTree objects are intuitive and flexible. Standard Python operators such as +, −, ∗, /, and ∗∗ are employed for these elementary operations. For an example, see Listing 9, where we have examples of usage for both standard arithmetic operations between FunctionTrees and between FunctionTrees and scalars. In these examples, f and g are instances of FunctionTree and a is a scalar. These operations provide an efficient and user-friendly way to manipulate FunctionTree objects in VAMPyR.

The integration between the overloaded arithmetic methods in VAMPyR and NumPy’s flexible data structures provides a range of computational advantages. Among these advantages, the support for multi-dimensional arrays, broadcasting, and linear algebra operations stands out. Although NumPy is primarily optimized for integers and floating-point numbers, its data structures can be extended to accommodate custom objects such as FunctionTrees. This capability enables the incorporation of FunctionTrees within NumPy arrays, thus allowing for the full suite of NumPy’s array-oriented computing features. For illustration, consider the example in Listing 10, where we construct a matrix of FunctionTrees within a NumPy array and execute various matrix operations.

As demonstrated in Listing 10, the integration between VAMPyR’s FunctionTree objects and NumPy arrays simplifies the implementation of complex mathematical operations. Through NumPy’s multi-dimensional arrays, it is possible to organize and manipulate arrays of function trees efficiently. Broadcasting allows for the addition of a single FunctionTree (or scalar) to an entire array of FunctionTrees. Furthermore, NumPy’s built-in linear algebra functions, such as the matrix multiplication operator @, can be used for matrix operations. This approach improves the readability of the code, as will be seen in Sec. V.

Handling derivatives with MWs is not straightforward: traditional derivative operators are not applicable. However, VAMPyR offers the Derivative(mra, type) class (see Listing 11 for a practical example), which applies the operators in the NS form as described in Subsection II E, allowing for efficient and accurate calculations. The types available are simple, center, forward, and backward, which are based on the ABGV method devised by Alpert et al. In addition, the Derivative class provides a type=b-spline option based on the method by Anderson et al.

• ABGV Types: Designed for piecewise continuous functions, these operators provide a weak formulation for first-order derivatives.

• B-Spline Type: This type is more versatile but works best for smooth continuous functions. The implementation covers up to second-order derivatives by transforming the function onto a B-spline basis before differentiation.

A more general method for constructing convolution operators based on Gaussian expansions is also available in VAMPyR. This approach is documented with an example in Listing 13, and an unusual application is illustrated in Fig. 5 by blurring an image with a Gaussian convolution kernel, although VAMPyR is not generally optimized for such tasks.

One of the main technical challenges of quantum chemistry is the large gap between the equations describing the physical nature of the system and their practical realization in a working code. This gap arises because the equations in their original form are, in general, too complicated to be solved. To achieve equations that have a manageable computational cost, it has so far been necessary to make use of representations, which on the one hand lead to practical numerical implementation, but on the other hand obfuscate the original physical equations.

We will here present four little pieces of quantum chemistry where we show how the MW representation of functions and operators in MRCPP and the simple Python interface provided by VAMPyR can close this semantic gap.

For each of the four examples presented, we will highlight how VAMPyR has been employed to obtain working implementations that closely resembles the theoretical framework. The underlying equations will here be presented briefly to highlight their correspondence with the code. We refer to separate Appendixes, or previously published work, for a more in-depth exposition of the theory.

All examples are available as Jupyter notebooks in a separate GitHub repository,⁴⁴ openly accessible under the CC-BY-4.0 license. We have also put in Zenodo: VAMPyR: https://doi.org/10.5281/zenodo.10290360 VAMPyR-coven: https://doi.org/10.5281/zenodo.10887800.

To achieve an even more compact structure, the set of orbitals can be collected in a NumPy vector. This requires that the implementation of operators in the code accepts a vector of orbitals as an input. The operators V_nuc, J_n, and K_n [defined in Eqs. (A2)–(A4), respectively] are implemented in a vectorized manner using NumPy. For example, the CoulombOperator class is presented in Listing 15.

These equations are easily represented in a MW framework, and each of them can be written as a single line of code with VAMPyR. This is shown in Listing 16. The model is extensively discussed in a previous paper⁴⁷ where a thorough study of theory and implementation was presented. The example in Listing 16 shows how continuum solvation has been implemented using VAMPyR. Here, we have assumed that we already have a projected permittivity function and the total solute density function together with a suitable mra and a function that computes γ called computeGamma.

The above integral equation can be solved iteratively as in the non-relativistic case, and the algorithm is thus very similar to the non-relativistic case. The main difference is that the infrastructure required to deal with four-component spinors needs to be implemented. Our design is built essentially on two Python classes: one to deal with complex functions (a single component of a spinor is a complex function) and another one to deal with the four components. Thanks to VAMPyR, these classes are easy to implement and use, for the most part overloading dunder operators.

Listing 17 shows a portion of the class implementation for the four-component Spinor. The objects are made callable and work as NumPy arrays. Each component is itself a complex function object (defined in a separate class). Some dunder methods are shown, and the methods to perform the operation c(α ⋅p)ψ are also reported.

First, the energy is computed in order to obtain the parameter μ, then comes the convolution of Vψ with the Helmholtz kernel $Ĝμ$⁠, and finally the Dirac Hamiltonian plus the energy h_D + ɛ is applied to the convolution. The iteration is repeated until the norm of the spinor update is below the requested threshold (Listing 18). At the end of the iteration, the energy must be computed once more (not reported in the listing).

In Listing 21, we define all the necessary operators and initialize Scheme (33). The numerical integration of the time evolution is performed in Listing 22. Figure 6 shows the result of the simulation: the oscillation movement of the density $Ψ(t)2$ and the difference between the numerical and exact solution at the time moment t = τ.

One of the main benefits of introducing a Python interface is the seamless interoperability it offers with a number of other packages. This not only simplifies the usage of the package but also enhances its capabilities by allowing the integration of features provided by other packages. In this section, we will explore some potential applications of this interoperability, showcasing how VAMPyR can interact with other packages to perform tasks that may otherwise be computationally demanding or require a more complex implementation.

As an example, the SCF solver implemented in VAMPyR, though fully general, encounters practical limitations with more complex molecular systems. One primary challenge is the selection of an appropriate initial guess for the orbitals; a poor choice can severely impact the convergence rate and even result in a failure to converge. The first step beyond the HF approximation might be to incorporate exchange and correlation density functionals, which, while beneficial, can be tedious to implement.

Here is where interoperability steps in. We can leverage the capabilities of different Python packages to address these issues. We will walk through examples demonstrating how VAMPyR can:

• Use VeloxChem⁵³ to generate an initial guess for a SCF, speeding up the process and increasing the chance of successful convergence.

• Use VAMPyR and VeloxChem to compute a potential energy surface grid, effectively utilizing computational resources.

• Perform stability analysis on Density Functional Theory (DFT) functionals from Libxc.^54,55

Here, we illustrate how to generate an initial guess using VeloxChem that can be imported into VAMPyR. Reconstructing the Molecular Orbitals (MOs) (or electron density) from the AO basis set using the MO (or density) matrix is a rather complex task, especially if the basis contains high angular momentum functions. We thus want to make use of VeloxChem own internal evaluator for these objects and wrap a simple function around it which can be projected onto the MW basis using a ScalingProjector. In principle, any $R3→R$ function can be projected in this way, so we just have to define a function that takes as argument a point in real space, runs it through VeloxChem internal AO evaluator code, and returns the function value of a given MO at that point.

First, we set up and run the SCF calculation with VeloxChem:

We can define a function that returns the value of the MO at a given point in space. This function will be used for projecting onto an MRA in VAMPyR:

Finally, we use a projector from VAMPyR to project the MOs from VeloxChem onto an MRA, generating a list of FunctionTrees for the initial guess:

In the last step, we loop over the desired number of orbitals, project each one onto the MRA, and store the result in Phi_0. This list of FunctionTrees represents an approximation to the molecular orbitals of the system and serves as an excellent initial guess for the SCF procedure in VAMPyR.

A Potential Energy Surface (PES) can be computed as a series of single point energy calculations while varying one (or more)⁵⁵ bond distance(s) in a molecule. We will here use VAMPyR’s adaptive function projector as a driver for computing the PES of the hydrogen molecule using VeloxChem’s single-point Hartree–Fock evaluator. We define a Python function, pes(r), that computes the energy for two hydrogen atoms given the bond distance as an input:

The function above computes the energy of a hydrogen molecule as a function of the bond distance (the 0.2 offset is to avoid the singularity at r = 0), by performing an SCF calculation at the given geometry. Being a $R→R$ mapping, it can be projected on a 1D MRA using VAMPyR.

The adaptive projector will automatically sample the pes(r) function in appropriate points in order to produce a smooth surface. The potential energy surface obtained in this manner can be visualized, as shown in Fig. 7.

The vertical lines in Fig. 7 represent the adaptive MRA grid, and the potential energy surface depicted is interpolated from this grid rather than computed directly at each point. The adaptive algorithm focuses computational resources where needed to achieve the desired precision. Further efficiency gains could be achieved realizing that high precision is only required for the region of low potential energy. This has not been exploited in the above example.

Inspired by the paper by Lehtola and Marques,⁵⁷ we investigate the numerical stability of density functional approximations (DFAs) lda-c-vwn and lda-c-gk72, see Fig. 8, with VAMPyR and Pylibxc.⁵⁴ The choice of these DFAs is due to their contrasting numerical stability, making them interesting subjects for our analysis. The lda-c-gk72, developed by Gordon and Kim in 1972, was known for its problematic convergence, while the lda-c-vwn (Vosko, Wilk, and Nusair functional) was noted for its numerical stability. Our objective is to verify these claims, illustrating the potential of VAMPyR as a numerical analysis tool.

First, we initialize a neon atom and calculate the self-consistent field (SCF) using VeloxChem:

Next, we calculate the density and project it onto a Multiresolution Analysis (MRA) grid:

Subsequently, we use Pylibxc to define the chosen DFA functionals:

Following this, an exchange potential is generated. We iteratively refine the grid based on precision, one scale at a time, until no new nodes are created:

Observing the output for each DFA, we note that the number of nodes converges to a finite value for the chosen precision in case of the lda_c_vwn functional. For the lda_c_gk72 functional, the number of nodes grows exponentially before the requested precision is reached. Looking at the wavelet norm, we notice that we are still far from a converged result. In this example, we employed a rather high precision (ɛ = 10⁻⁸).

Multiresolution analysis is a relatively new tool in quantum chemistry, compared to traditional methods based on atomic orbitals. Its main advantages are the possibility to reach any predefined arbitrary precision (the only limits are dictated by the machine precision and the available memory), and a numerical implementation that stays close to the theoretical formalism. VAMPyR has been designed to capitalize on these two aspects and, at the same time, enable fast code development for prototype applications.

We have employed VAMPyR to test new ideas and methods and to create interfaces with existing codes. In the first category, we have shown here a simple SCF optimization, a continuum solvation model, the solution of the Dirac equation for one electron, and the application of the time evolution operator. In the second category, we illustrated how to import a starting guess for the orbitals, the plotting of a PES, a numerical analysis of two density functionals. For all the above examples, we have provided extracts of the code, which show how VAMPyR is used. The complete set of working examples illustrated in Sec. V has been collected in an openly accessible GitHub repository,⁴⁴ under the CC-BY-4.0 license.

This work was supported by the Research Council of Norway throughits Centres of Excellence scheme (262695), through the FRIPRO grant ReMRChem (324590), and from NOTUR – The Norwegian Metacenter for Computational Sciencethrough grant of computer time (nn14654k).

The authors have no conflicts to disclose.

Magnar Bjørgve: Conceptualization (lead); Methodology (lead); Project administration (lead); Software (lead); Writing – original draft (equal). Christian Tantardini: Conceptualization (equal); Supervision (equal); Writing – original draft (equal); Writing – review & editing (equal). Stig Rune Jensen: Conceptualization (equal); Software (equal). Gabriel A. Gerez S.: Conceptualization (equal); Software (equal); Writing – original draft (equal). Peter Wind: Conceptualization (equal); Software (equal). Roberto Di Remigio Eikås: Conceptualization (equal); Software (equal); Supervision (equal); Writing – review & editing (equal). Evgueni Dinvay: Conceptualization (equal); Software (equal); Writing – original draft (equal). Luca Frediani: Conceptualization (equal); Project administration (equal); Supervision (lead); Writing – review & editing (lead).

The complete set of working examples illustrated in Sec. V has been collected in an openly accessible GitHub repository https://github.com/MRChemSoft/vampyr-coven, under the CC-BY-4.0 license.

The problem is then cast in matrix form, and the eigenvalues (energies) and eigenvectors (orbital expansion coefficients) are obtained by standard linear algebra techniques. The immediate advantage is a representation closely related to the physics of the system (atomic orbitals), but this comes at a price: the implementation deals with the representation of such functions, their integrals, their overlaps, and seemingly simple operations, such as applying an operator or multiplying two functions, become technically complicated obfuscating the physical significance. Moreover, the initial scaling of the problem becomes n⁴, where n is the number of basis functions, due to the two-electron integrals, and a lot of effort is necessary to recover the original n² scaling, which is, instead, straightforward for a MW implementation.

In these equations, we have first split the one-electron operator $ĥ$ in the kinetic energy $T̂$ and the potential energy $V̂N$⁠, then rearranged the terms, including the diagonal element F_ii of the Fock matrix, and finally applied the BSH Green’s function to obtain the final expression. Here, there is no need to transform the problem to a different basis. The last term, which represents the Lagrangian multipliers, can be omitted if the canonical HF equations are used. The final equation can, moreover, be interpreted as a preconditioned steepest descent, where the gradient is the expression in the square bracket (it can be obtained by taking the differential of the expectation value of the energy of a Slater determinant with respect to a generic orbital variation) and the preconditioner is the BSH kernel.

The VAMPyR implementation of the time evolution operator $expit∂x2$ is under optimization currently, although it is already available in VAMPyR (see Listing 23).

Finally, the tree structure of the potential semigroup $exp−itV$ is introduced in Listing 24.

Together with a splitting scheme, for example (33), this completes the algorithm description for the time evolution simulations.