VAMPyR – A High-Level Python Library for Mathematical Operations in a Multiwavelets Representation

Wavelets and multiwavelets have lately been adopted in quantum chemistry to overcome challenges presented by the two main families of basis sets: Gaussian atomic orbitals and plane waves. In addition to their numerical advantages (high precision, locality, fast algorithms for operator application, linear scaling with respect to system size, to mention a few), they provide a framework that narrows the gap between the theoretical formalism of the fundamental equations and the practical implementation in a working code. This realization led us to the development of the Python library called VAMPyR (Very Accurate Multiresolution Python Routines). VAMPyR encodes the binding to a C++ library for multiwavelet calculations (algebra and integral and differential operator application) and exposes the required functionality to write a simple Python code to solve, among others, the Hartree-Fock equations, the generalized Poisson equation, the Dirac equation, and the time-dependent Schrödinger equation up to any predefined precision. In this study, we will outline the main features of multiresolution analysis using multiwavelets and we will describe the design of the code. A few illustrative examples will show the code capabilities and its interoperability with other software platforms.


I. INTRODUCTION
Wavelets and Multiwavelets have emerged in the last 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.2][3][4][5][6][7][8][9] The code which pioneered this approach is M-A-D-N-E-S-S, 10 followed by our own code MRChem. 11To date, they are the only two codes available for quantum chemistry calculations using multiwavelets.
As the development of MRChem unfolded, we found advantageous to separate the mathematical code dealing with the MultiWavelet (MW) formalism in a separate library, called MRCPP (Multiresolution Computation Program Package) 12 .MRCPP is a high-performance C++ library, which provides the tools of Multi-Resolution Analysis (MRA) with Multiwavelets.It implements lowscaling algorithms for mathematical operations and convolutions and a robust mechanism for error control during numerical computations.
Once MRCPP was available, we realized 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 14

and
Matplotlib 15 , make it a powerful tool for scientific computing.Moreover, Python's integration with Jupyter Notebooks 16 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 Multiwavelets Python Routines) 17 , 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,19I.MULTIWAVELETS Tha foundations of Multiresolution Analysis (MRA) have been laid in the early '80 thanks to the pioneering work of Meyer, Daubechies, Strömberg and others, [20][21][22][23][24][25][26][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][29][30][31][32] .

A. Multiresolution Analysis
The construction of a Multiresolution Analysis (MRA) starts by considering a basis of generating functions, called scaling and wavelet functions.We will here limit ourselves to the specific case of multiwavelets, referring to the literature for a wider presentation of the subject [20][21][22][23][24][25][26][27] .The most common choice is the set of k + 1 polynomials up to order k in the unit interval, such as the Legendre polynomials or interpolating Lagrange polynomials 31,32 .The original scaling function can be translated and dilated: In the equation above ϕ are the generating polynomials, j is the polynomial index ( j = 0, 1, . . .k), n is the scale index and l is the translation index.For a given scale n the functions ϕ n constitute a scaling space V n .By construction, each interval V n l in V n splits into two sub-intervals in V n+1 : V n+1 2l and V n+1 2l+1 , doubling the basis with every n along the ladder.This establishes the nested structure V n ⊂ V n+1 seen in (2).It is also possible to show that this hierarchy is dense in The wavelet space W n is defined as the orthogonal complement of V n with respect to V n+1 : The wavelet functions are also supported on the same disjoint intervals on the real line as the corresponding scaling functions.Dilation and translation relationships hold for wavelet functions as well: scaling space.This property, known as "vanishing moments", is key in achieving fast algorithms for various operations with multiwavelets.The construction of the multiwavelet functions in not unique (each W n l spans a k + 1 dimensional space and any orthonormal basis in this space is a legitimate choice).We follow the construction presented by Alpert 31 , which both maximizes the number of vanishing moments and divides the functions in symmetric and antisymmetric ones.
We finally note that Equation (3) leads to the following telescopic series:

B. Multiwavelet transform
The ladder of spaces structure which is summarized in Eq. ( 3) and Eq. ( 5) implies that it is possible to span V n+1 either through the scaling functions ϕ n+1 or through the scaling functions ϕ n and the wavelet functions ψ n .The basis set transformation between these two bases is known as the multiwavelet transform or two-scale filter relationship.The unitary transformation matrix is assembled by making use of the four wavelet filters G (0) , G (1) , H (0) and H (1)   The matrices H (0) , H (1) , G (0) , and G (1) 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 φ n l and wavelet function ψ n l are only connected to φ n+1 2l and φ n+1 2l+1 .This structure simplifies numerical algorithms, enabling fast implementations.

C. Function projection onto the multiwavelet space
The most basic operation in MRA is the projection of a function.It generates the representation of an arbitrary function in a MW basis.Let us consider the scaling space V n , and the associated projector P n .The MW projection of a given function f (x) can be obtained as: Here, s n, f j,l are the scaling coefficients, given by: Similarly, a wavelet projector Q n is associated with the wavelet space W n : The wavelet coefficients w n, f j,l can also be obtained as:

D. Adaptive Projection
Keeping in mind the projections defined in the previous section, the complete representation of a function in a given MRA can be written as follows: This equation gains special significance when considering that the wavelet coefficients, w n l , approach zero for smooth functions.The precision of the representation can thus be assessed by inspecting the wavelet coefficients at a specific scale n and translation l: This consideration provides the foundation for an adaptive projection strategy: rather than fixing the projection to a predefined scale N, one can incrementally increase the scale from 0 only in those intervals where the wavelet norm is larger than the requested precision.Utilizing the twoscale filter relations, Eq. ( 6), we can compute the wavelet coefficients and assess their magnitude.
If these coefficients meet predefined precision criteria at a specific translation l, that branch of the function representation can be truncated.This truncation effectively focuses computational and data resources only where high precision is necessary.Should the desired precision level not be reached, the refinement process continues in a recursive manner.This recursion can in principle be carried out indefinitely until precision requirements are met.In practice, a maximum allowed scale is set.

E. Operator Application in Multiwavelet Framework -Non-Standard Form
The most convenient strategy to apply an operator using MWs is to construct its non-standard (NS) form.Similarly to functions, one can construct the operator projection T N = P N T P N .By recalling that for every scale n, P n+1 = P n + Q n one obtains: where for each n the following definitions are used: This structure is the operator analogue of the two-scale relationship for functions and it can be extended telescopically to obtain: A n , B n , and C n components exhibit sparsity properties, similar to what is observed for the wavelet coefficients of function representations, and are therefore leading to fast, and often linearly scaling algorithms, for any arbitrary predefined precision.The purely scaling part T 0 of the operator is only required at the coarsest scale, where only a handful of grid points are present.Another important feature of the NS form is the absence of coupling between different scales, which allows to preserve the adaptive precision of the representation on the one hand, and independent or asynchronous operator application across different scales, which boosts computational efficiency, on the other hand.For additional details about the application of operators the non-standard form 33 , the reader is referred to the literature on the subject 32,34 .

III. VAMPYR
In the architecture of VAMPyR, MRCPP serves as the foundational layer, integrated into Python through Pybind11, as illustrated in Figure 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 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 35 .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: conda install vampyr -c conda-forge

A. Implementation
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. 38The specialization for 1-, 2-, and 3dimensional problems is performed at compile-time, thereby eliminating any impact on runtime performance. 39 Python, native template constructs are absent, presenting a challenge for dimension-specific specialization.In VAMPyR, we emulate MRCPP's generic approach by implementing dimensiondependent bindings using Pybind11.Specifically, our binding code incorporates template classes and functions similar to those in MRCPP.MRCPP's 1-, 2-, and 3-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.
from vampyr import vampyr1d from vampyr import vampyr2d from vampyr import vampyr3d While the MultiResolutionAnalysis (see section IV A) and FunctionTree (see section IV B) classes in VAMPyR are inherited from MRCPP, there are noteworthy distinctions between the vanilla C++ and the corresponding Python classes.These differences are introduced to enhance the Python version of the MRCPP classes with so call 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 to write more intuitive code, which is easier to both read and write.
As an example, consider the addition of two FunctionTree objects, which be expressed in natural Python syntax (Listing 2), thanks to the binding code in Listing 3.
Listing 2. Syntax sugar for FunctionTree addition in action.The Python code in Listing 2 can be equivalently written in the de-sugared form in Listing 4.Although the 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.

IV. MATHEMATICS WITH VAMPYR
In this section we display some of the theoretical concepts discussed in Section II with practical examples using VAMPyR.

A. The MultiResolutionAnalysis Object in 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 Section II).
a. Standard Usage 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 illustrate 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.
b. Advanced Usage 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, the MRA object based on their specific needs and expertise.

B. Function Projectors
In VAMPyR, functions are represented using a tree-based data structure: the FunctionTree

C. Scaling and Wavelet projectors in VAMPyR
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. 41This illustrates how to project a function using the vp.ScalingProjector and vp.WaveletProjector classes.
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 Equations ( 7) and ( 9), respectively.In Figure 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 Figure 3 4 between the adaptive and fixed-scale projections (previously shown in Figures 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 Figure 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.

E. Arithmetic Operations in FunctionTrees
In VAMPyR, arithmetic operations on FunctionTree objects are intuitive and flexible.Standard Python operators like +, −, * , /, 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.

F. Vectorizing FunctionTrees with NumPy
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 stand out.Although NumPy is primarily optimized for integers and floating-point numbers, its data structures 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 Section V.

G. Derivative Operators in VAMPyR
Handling derivatives with MWs is not straightforward: traditional derivative operators are not applicable.However, VAMPyR offers two specialized operators to address this issue, each tailored for specific requirements concerning the continuity of the function and the desired order of the derivative.These methods are based on the works by Alpert et al., giving the ABGVOperator, 32 and Anderson et al., leading to the BSOperator. 42th the ABGVOperator and the BSOperator apply the operators following the NS form as described in Subsection II E, allowing for efficient and accurate calculations: ABGVOperator: Designed for piece-wise continuous functions, this operator provides a weak formulation for first-order derivatives.It is primarily suitable for handling first-order derivatives and offers a balance between accuracy and computational efficiency.Usage is demon- The two kernels are numerically implemented in the MRCPP library and available in VAMPyR.They are constructed as a sum of Gaussians which represents the numerical quadrature of the integral of a superexponentially decaying function 43 which depends on the distance d = |x − y|.
In the expression above, µ = 0 yields the Poisson kernel, and µ > 0 corresponds to the bound-state Helmholtz kernel.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 14 and an unusual application is illustrated in Figure 5 by blurring an image with a Gaussian convolution kernel, although VAMPyR is not generally optimized for such tasks.

V. FOUR LITTLE PIECES OF QUANTUM CHEMISTRY
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 which 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 appendices, 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 44 , openly accessible under the CC-BY-4.0 license.

A. The Hartree-Fock Equations
The Hartree-Fock equations 45 can be concisely written as: where F is the Fock operator, ϕ i is an occupied orbital and F i j are the Fock matrix elements in the basis of the occupied orbitals.
Traditionally, a basis-set representation is introduced and most of the successive discussion regards methods to achieve a good performance and a sufficient precision.A MW representation avoids this complication and solves the equations directly.As shown in Appendix A, this is achieved by recasting the problem as an integral equation: In the above equation, Ĝµ n i is the bound-state Helmholtz operator, V N , J and K are the nuclear, Coulomb and exchange operator, respectively and we have restricted ourselves to closed-shell systems.i is the orbital index and n represents the iteration index and is displayed for all quantities that are iteration-dependent.The resulting orbitals φi are not normalized (here indicated with a ˜superscript) and an orthonormalization step must then be performed.The complete algorithm can be found in the accompanying Python notebook 44 and includes a Löwdin orthogonalization of the orbitals, ensuring the desired orthogonality condition.We base our implementation on these equations.The Self-Consistent Field (SCF) equations ( 22 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 input.The operators V_nuc, J_n, and K_n (defined in equations (A2), (A3), and (A4), respectively) are implemented in a vectorized manner using NumPy.For example, the CouloumbOperator class is presented in Listing 16.

B. Continuum solvation
A very common approach to including solvent effects in a quantum chemistry calculation is to consider the molecule immersed in a dielectric continuum.The governing equation is the generalized Poisson equation (GPE): where V is the electrostatic potential, ρ is the charge distribution of the molecule and ε is the (position-dependent) permittivity of the dielectric.The practical solution of the equation depends on how ε is parametrized.Traditionally, a cavity boundary is defined and the permittivity is a step function of the cavity: ε i = 1 inside of it and ε o = ε s outside, where ε s is the bulk permittivity of the solvent.This parametrization leads to an equivalent boundary integral equation, solved with a suitable boundary element method. 46The main advantage of the approach is transferring the problem from the (infinite) 3-dimensional space to the (finite) cavity surface.The main disadvantage is a less physical parametrization of the model, due to the presence of the sharp boundary.
Using MWs, Equation 23can instead be solved directly, without assumptions on the functional form of ε.The main working equation is the application of the Poisson operator P to the effective charge density ρ eff (x) and the surface polarization γ(x): where the effective density ρ eff and the surface polarization γ(x) are defined respectively as: The surface polarization depends on the potential, which implies that the equation must be solved iteratively.
The equations above 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 17.The model is extensively discussed in a previous paper 47 where a thorough study of theory and implementation was presented.The example in Listing 17 shows how continuum solvation has been implemented using VAMPyR.Here we have assumed 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 solvation energy E R is finally computed from the reaction potential V R = V − dr ′ ρ(r ′ ) |r−r ′ | , and the solute charge density ρ: Dirac equation. 48,49It is well known that the level of complexity required to deal with the Dirac equation increases significantly: the scalar non-relativistic eigenvalue problem, dealing generally with real-valued functions, becomes a 4-component problem involving complex functions: where V is the external potential, ε is the energy eigenvalue, c is the speed of light, p is the momentum operator and σ is the vector collecting the three Pauli matrices.As shown by Blackledge and Babajanov 50 and later exploited by Anderson et al. 51 , the Dirac equation can also be reformulated as an integral equation as: where Ĝµ is the BSH kernel as in the non relativistic case (see Section IV H), and µ = m 2 c 4 −ε mc 2 .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 4-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.With the Spinor class implementation at hand, the iterative scheme to converge the Hydrogen atom is easily implemented: 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.At the end of the iteration the energy must be computed once more (not reported in the listing).

D. Time-dependent Schrödinger equation
The MWs framework can successfully be employed for real-time simulations of a wavepacket by directly integrating the time-dependent Schrödinger equation: where the potential V may be time-dependent.For simplicity we consider here a time-independent potential and a one-dimensional problem: V = V (x) and Ψ = Ψ(x,t) with x,t ∈ R, and we assume that the initial conditions are given as Ψ(x, 0) = Ψ 0 (x).The standard numerical treatment of (30)   is to choose a small time step t > 0 and construct the time evolution operator on the interval [0,t].
By applying it iteratively, the solution at any finite time can in principle be obtained.The wave propagation can be expressed as where we used the fact that the potential is time-independent.We proceed by splitting the propagator in the kinetic exp −it T and potential exp −it V parts.The former is a multiplicative operator in momentum space, whereas the latter is a multiplicative operator in real space.This separation is not exact, because T and V do not commute.The resulting operator has an error of This simple splitting is too rough for practical applications, therefore we make use of the following fourth order scheme 52 where With A = −i T and B = −i V , B becomes a multiplicative operator containing the potential gradient Remarkably, this high order scheme requires only two applications of the free-particle semigroup operator exp −it T /2 per time step.For an example of the implementation see Listing where self-adjoint H stands for either the kinetic energy −∂ 2 x , or the potential V (x) or the full Hamiltonian −∂ 2 x +V (x).The implementation is shown in Listing 21.As a worked-out example,  we show a simulation of the time evolution of a Gaussian wave packet Ψ(x,t) in harmonic potential V (x): It is well known that the density |Ψ(t)| 2 oscillates in the harmonic potential with the period τ = π/ √ V 0 .More precisely, Ψ (τ) = −Ψ 0 .This can immediately be seen taking into account that the eigenvalues for the Hamiltonian are √ V 0 (2n + 1).We take x 0 = 0.3, σ = 0.04 and x 1 = 0.5, V 0 = 25000.The parameters are chosen in such a way that the solution stays localized mainly on the [0, 1] interval.Note that the B operator simplifies to .
In Listing 22 we define all the necessary operators and initialize Scheme (33).The numerical integration of the time evolution is performed in Listing 23. 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 = τ.

VI. INTEROPERABILITY WITH OTHER PACKAGES
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 53 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.

A. Generating an Initial Guess with VeloxChem
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 R 3 → 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.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: 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.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 region of low potential energy.This has not been exploited in the above example.

C. Numerical Stability Analysis using VAMPyR and Pylibxc
Inspired by the paper by Lehtola and Marques, 56 we investigate the numerical stability of density functional approximations (DFAs) lda-c-vwn and lda-c-gk72, 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 numerical analysis tool.
Firstly, we initialize a Neon atom and calculate the self-consistent field (SCF) using VeloxChem: 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 requated precision is reached.Looking at the wavelet norm, we notice that we are i still fa r from a converged result.In this example we employed a rather high precision (ε = 10 −8 ). and B j m obey the same recurrence for j ⩾ 1.
The VAMPyR implementation of the time evolution operator exp it∂ 2 x is under optimization currently, though it is already available in VAMPyR (see Listing 24).
Finally, the tree structure of the potential semigroup exp (−itV ) is introduced Listing 25.Together with a splitting scheme, for example (33), this completes the algorithm description for the time evolution simulations.

Figure 1 .
Figure 1.Connections between MRChem, MRCPP and VAMPyR.MRCPP implements a high-performance MRA framework with MWs.It implements the representation of functions and operators as well as fast algorithms for operator application.MRChem implements electronic structure methods and uses MRCPP for the underlying mathematical operations.VAMPyR imports features from MRCPP using Pybind11 and brings intuitive design and easy prototyping through Python.It is used as a Python library can therefore be interfaced other quantum chemistry software.VAMPyR doesn't include functionality from MRChem, at the current stage.All software packages are available on GitHub, under the MRChemSoft organization 35 .

Listing 5 .
Standard usage of the MRA object.# For a unit cell centered at the origin MRA = vp.MultiResolutionAnalysis(box=[-L,L], order=k) # For a unit cell with left corner at the origin MRA = vp.MultiResolutionAnalysis(box=[L], order=k) as demonstrated in Listing 6.This flexibility allows users to choose the level of interaction with Listing 6. Advanced usage of the MRA object.# Construct scaling basis scaling_basis = InterpolatingBasis(order) # or LegendreBasis(order) # Define world box world_box = BoundingBox(scale=0, corner=[x0, y0, z0], nboxes=[Nx, Ny, Nz], scaling=[a, b, c]) # in 3D, for 1D and 2D replace list of size 3 with list of size 1 or 2 # Construct MRA MRA = vp.MultiResolutionAnalysis(box=world_box, basis=scaling_basis) object.This structure naturally arises from the MRA framework outlined in Section II C.a.The Tree The FunctionTree object consists of interconnected nodes organized hierarchically.The root node40 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. b.The Node 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 1 , l 2 , . . ., 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.

Listing 7 .
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 (Figures 2(a) and 3(a)), while the right-hand side illustrates the wavelet function space (Figures 2(b) and 3(b)).

Figure 2 .
Figure 2. Projection of a Slater function at root scale.Left: Projection onto scaling function space V 0 4 with vertical lines marking the physical space spanned by k + 1 scaling functions.Right: Projection onto wavelet function space W 0 4 , where vertical lines indicate the physical space spanned by k + 1 wavelet functions.

Figure 3 .
Figure 3. Projection of a Slater function at the 4th scale.Left: Projection onto scaling function space V 4 , with vertical lines delineating the physical space spanned by k + 1 scaling functions.Right: Projection onto wavelet function space W 4 , where vertical lines mark the physical space spanned by k + 1 wavelet functions.

Figure 4 .
Figure 4. Adaptive projection of the Slater function with prec = 1.0 × 10 −1 .The non-uniform vertical lines indicate the physical space spanned by the adaptive basis functions, underscoring the method's efficiency and precision.

Listing 9 .
Basic Arithmetic Operations in VAMPyR.# Arithmetics with scalars a_mult_f = a * f # Scalar on the left f_mult_a = f * a # Scalar on the right f_div_a = f / a # Division by scalar f_pow_a = f ** a # Power operator # Arithmetics between FunctionTrees f_mult_g = f * g # Multiplication operator f_add_g = f + g # Addition operator f_sub_g = f -g # Subtraction operator # In-place arithmetics f *= a # Multiplication with scalar f *= g # Multiplication with function f += g # Addition with function f -= g # Subtraction with function

Listing 10 .# 2 )
Performing matrix operations using function trees in a NumPy array.# Creation of a matrix containing function trees matrix_a = np.array([[a11_tree,a12_tree], [a21_tree, a22_tree]]) matrix_b = np.array([[b11_tree,b12_tree], [b21_tree, b22_tree]]) # Performing matrix multiplication matrix_mul = matrix_a @ matrix_b # Performing addition and subtraction matrix_sum = matrix_a + matrix_b matrix_neg = matrix_a -matrix_b # Multiplication by a scalar matrix_mul_scalar = 2 * matrix_a * matrix_b # Demonstrating broadcasting by adding a single function tree to all elements matrix_a_bcast_add = matrix_a + single_tree strated in Listing 11.BSOperator: This operator is more versatile but it works best for smooth continuous functions.It accommodates higher-order derivatives by transforming the function onto a B-spline basis before differentiation.See Listing 12 for a practical example.H. Convolution Operators in VAMPyR Convolution operators are essential tools for solving integral forms of key equations in Quantum Chemistry, such as the integral formulation of the Hartree-Fock (HF) and Kohn-Sham (KS) equations or the Poisson equation used in self-consistent reaction field (SCRF) models.Both equations can be recast from the more familiar differential form, to an equivalent integral form Listing 11.First-order derivatives using ABGVOperator.D = vp.ABGVOperator(MRA) # x, y, z derivatives using ABGVOperator f_x = D(f, 0) f_y = D(f, 1) f_z = D(f, 2) Listing 12. First-and second-order derivatives using BSOperator.D1 = vp.BSOperator(MRA, 1) # First-order derivative operator D2 = vp.BSOperator(MRA, 2) # Second-order derivative operator by making use of the Bound-State Helmholtz (BSH) kernel (see Appendix A) and the Poisson kernel.The key step in both problems is the convolution of the corresponding Green's function kernel with a test function:

Listing 14 .
Construction of a Custom Convolution Operator in VAMPyR.# Create a 1D GaussExp object, serving as a container for Gaussians kernel = vp1.GaussExp() # Append a Gaussian function to the kernel kernel.append(vp1.GaussFunc(beta, alpha)) # Construct the convolution operator T = vp.ConvolutionOperator(mra, kernel, prec) # Apply the operator to a function tree g_tree = T(f_tree) (a) Original image projected onto the Multi-Resolution Analysis (MRA).(b) Resulting image after performing convolution with a Gaussian Kernel, leading to smoothing or blurring of the original image.

Figure 5 .
Figure 5. Illustration of convolution process in VAMPyR.

The
Listing 18 shows a portion of the class implementation for the 4-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.

20 .cos
Complex algebra is not supported natively in VAMPyR, therefore we represent complex exponential operators as follows exp −i Ht =   Ht sin Ht − sin Ht cos Ht   , operating on vector-functions (a) Density evolution in the harmonic potential.(b) The difference between the numerical and exact solution.

#
Define a function to get the MOs, this can be projected onto an MRA by vampyr def mo_i(mol_orbs, r): R = np.array([r])return vis_drv.get_mo(R,molecule, basis, mol_orbs, i, "alpha")[0] 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: # This can now be imported into vampyr as: P_eps = vp.ScalingProjector(mra, prec) Phi_0 = [P_eps(mo_i) for i in range(nr_orbs)]

Figure 7 .
Figure 7. Potential energy surface calculated using VeloxChem and VAMPyR.The grid lines represent the adaptive MRA grid, from which the PES is interpolated.

Figure 8 .
Figure 8.Comparison of node growth between the lda-c-vwn and lda-c-gk72 DFAs.The plot illustrates the sharp contrast in the convergence behavior, with lda-c-vwn showing a stable, linear growth, and lda-c-gk72 displaying an exponential growth.