MultipolePreserving Quadratures for Discretization of Functions in RealSpace Electronic Structure Calculations
Abstract
Discretizing an analytic function on a uniform realspace grid is often done via a straightforward collocation method. This is ubiquitous in all areas of computational physics and quantum chemistry. An example in Density Functional Theory (DFT) is given by the external potential or the pseudopotential describing the interaction between ions and electrons. The accuracy of the collocation method used is therefore very important for the reliability of subsequent treatments like selfconsistent field solutions of the electronic structure problems. By construction, the collocation method introduces numerical artifacts typical of realspace treatments, like the socalled eggbox error, that may spoil the numerical stability of the description when the realspace grid is too coarse. As the external potential is an input of the problem, even a highly precise computational treatment cannot cope this inconvenience. We present in this paper a new quadrature scheme that is able to exactly preserve the moments of a given analytic function even for large grid spacings, while reconciling with the traditional collocation method when the grid spacing is small enough. In the context of realspace electronic structure calculations, we show that this method improves considerably the stability of the results for large grid spacings, opening the path towards reliable lowaccuracy DFT calculations with reduced number of degrees of freedom.
I Introduction
Realspace grid based techniques are very important in disciplines like Quantum Chemistry, Computational Physics and Applied Mathematics. A real space approach is mandatory in the solution of complex problem of Partial Differential Equations, as well as for the treatment of complex environments and nontrivial boundary conditions. In this framework, the collocation method is a straightforward procedure that is used to discretize a known function, to express its values in the realspace domain. This method represents the most straightforward and intuitive way to provide numerical coefficients to discretize a computational problem. In this method, a function is represented via a set of point values
(1) 
where are the sampling points of the simulation domain.
As any discretization method, function collocation introduces an error on the numerical results. This error of course decreases while increasing the number of points used for the discretization, but its convergence ratio depends of many factors. As an example, imagine that the function represents the charge density of a Poisson Equation
(2) 
discretized on the points by the collocation method. If the multipoles of the coefficients do not correspond with the ones of the original function , a numerical solution of the above equation will not produce a correct discretization of the potential , even if the adopted Poisson Solver is very accurate. This happens in a number of numerical treatment: the discrete multipoles (i.e. the momenta of the coefficients ) of the collocated functions determine the accuracy of the final results.
In this paper we present a numerical quadrature formula to obtain a set of coefficients which can be used as a generalization of the collocation method for analytic functions. This quadrature scheme is tailored to preserve the values of the discrete multipoles of the coefficients, to be fixed to the momenta of the original function, even for discretization done on grids of large spacings. Such a quadrature scheme involves the usage of the Interpolating Scaling Function (ISF) basis set.
Similar needs for quadrature formula have already been presented in Ref. MagicFilter (), in the context of the gridpoint discretization of functions expressed in the Daubechies wavelets basis, within the socalled “Magic Filter” method. Here we extend and generalize such concept to the realspace discretization of any arbitrary functions.
In the next section we will illustrate how this problem appears in a RealSpace based DFT code, and we show the importance of preservation of the monopole in this context. We then quantify the discretization errors coming from the collocation method, by explaining our need for an alternative formula, and the properties that the generalized collocation scheme should satisfy. Then we will show the improvements related on the usage of this new scheme with supporting results from the BigDFT code, showing how the behavior of the results is stabilized for large grid spacings, enabling us to perform lowaccuracy calculations with reduced number of degrees of freedom. Then we also explain how this quadrature method can be used to perform accurately and efficiently scalar product between known functions and compactly supported basis sets, showing its connections with the Magic Filter method.
Ii Collocation and Real Space Electronic Structure codes
A more elaborated example than the Poisson equation presented above is given by electronic structure calculations on realspace based codes. Density functional theory (DFT) is a widely used and recognized method to investigate material properties at the atomistic scale. In such ab initio calculations the electronic structure problem is solved by minimizing the total energy of the system which is a functional of the electronic density
(3) 
where the four terms on the right side of Eq. (3) are, respectively, the kinetic energy, the electrostatic and the exchangecorrelation energy, and the interaction energy with an external potential. Such external potential describes the interaction of the electrons with the ions, and can be considered as the input of the computational problem. In the KohnSham formalism of DFT, the potential appears in the Hamiltonian operator
(4) 
where the KohnSham potential depends selfconsistently on the charge density and it is determined during the optimization procedure. Its discretization is therefore the primary responsible of the quality of the physical and chemical information that can be extracted from the DFT calculation. As an example, let us now consider a realspace based pseudopotential code. The external potential operator usually is
(5) 
therefore separated by local and nonlocal terms. The local term, responsible for the electrostatic interaction between the ions and the valence electrons, satisfies the Poisson Equation
(6) 
where is the charge density associated to the atomic ions, screened by the core electrons. To have a reliable DFT calculation, it is therefore essential that the potential is associated to a with a correct value of the monopole, namely the sum of the ionic valence charges.
To illustrate the importance of this consideration, we have plotted in Fig. 1 the absolute DFT energy of an organic molecule, as calculated by the BigDFT code BigDFT () where is given by normconserving GTHHGH pseudopotentials HGHWilland (), possibly enhanced by a nonlinear core correction. The atomic geometry of the molecule was optimized using the LDA functional.
We show that the results become numerically unstable above the value of AU of the grid spacing. It is obvious from the figure that this effect is caused, at first approximation, by the fact that the total ionic charge is not preserved when the grid spacing increases. The discretization of based on the collocation method limits considerably the use of realspace basis sets especially for lowaccuracy or highthroughput calculations.
These results are extracted with a particular DFT code, but the message is general: each time that the collocation method is used, like for example the evaluation of the charge density on a real space mesh (that is a fundamental operation, in all DFT codes, needed to evaluate and its derivatives), attention must be paid to the accuracy of the collocation. Another example is the definition of the compensating charge in methods like the Projector Augmented Waves PAW ().
In the following we will present our solution to this crucial problem and the potential outcomes for realspace based DFT codes.
Iii The collocation method and the interpolation
For discretization on uniform grid spacings, the collocation method is welljustified when the original function can be reasonably approximated by an interpolation of its values on the grid mesh points. Let us consider a one dimensional function . Suppose we want to discretize this function on a uniform grid of spacing and coordinates . Given a family of interpolating functions , if the approximation
(7) 
is reasonably accurate, the collocation method can be applied. This fact stems from the interpolating property of the family . Indeed, an interpolating family is constituted by a set of functions , each one associated to a point of the grid, such that . From Eq. (7), and the continuous representation of may be given by .
Given the interpolating property, it is also said that an interpolating function family is dual to the Dirac deltas. In other terms, denoting the above function by the braket notations, we have
(8) 
where represents the Dirac distribution centered at point , i.e. . The above defined interpolating property implies that the duality relation holds.
iii.1 Polynomial exactness and discrete multipoles
The accuracy of the approximation (7) is of great importance for a reliable computational treatment. Clearly, such accuracy is intimately related to the family of interpolation functions chosen. The interpolating function families are normally constructed using a family of polynomial functions. An interpolating family is said to be of order if any monomial function (denoted with in the following), with is exactly expressed by the interpolating family, for all lying within a given interval . In other terms
(9) 
This is the concept of polynomial exactness. Note that the index runs over a set of grid points which might lie outside the support . We indicate by the minimum interval of grid points needed to obtain the polynomial exactness in the interval .
The collocation method is therefore meaningful for the functions for which the projector operator approaches the identity. The polynomial exactness is important in determining the accuracy of the interpolation: a smooth function can reasonably be approximated by its Taylor polynomial around a given point. The higher the order of the polynomial exactness of the functions , the better the Taylor expansion of the original function would be approximated by the function , therefore the norm of the difference will be reduced in the support by .
It is easy to understand that this condition is meaningful only when the grid spacing size is smaller than the typical oscillations of the function we want to represent. Evidence is shown in Fig. 2 for a Gaussian function of standard deviation , centered between two grid points. When , the collocation coefficients give a bad representation of the function, as their amplitudes is suppressed everywhere by decreasing . No interpolating family, even of very high order , will therefore be able to faithfully represent the original function. This situation seems unavoidable: as the expansion coefficients of the function are given in terms of the scalar products , the grid has to provide a reasonable sampling of the function such as to exhibit the convergence rate.
When the grid is small enough, one might think that also the discrete multipoles, i.e.
(10) 
follows the same convergence ratio, while approaching their exact values given in terms of . However, it is easy to see that this is not the case. The reason is that the dual function is always represented by a Dirac distribution, and it is therefore independent of the quality of the interpolating family.
As shown in Fig.2, the collocation method gives inaccurate results for the discrete multipoles when . In addition, their convergence ratio is of , as the coefficients cannot depend on the interpolating functions. At variance from the evaluation of other quantities like function derivatives, increasing the order of the interpolating function will not change the convergence ratio of the discrete multipoles of lower order.
iii.2 The need of an alternative formula for collocation
Being an analytic function, the accuracy of the multipoles is a quantitative evaluation of the accuracy of which is more severe from the one provided by the function difference, as it cannot be modified by varying the family of interpolating functions. We would like to have a different quadrature formula, such that the number of preserved discrete multipoles of the original function is of the same order of the family of interpolating functions chosen, regardless of the value of the grid spacing.
As pointed out in the previous section, we need to define an alternative set of dual functions , such that the multipoles of the original functions are preserved up to order . In other terms, we search for a family of dual functions such that
(11) 
It is easy to see that this condition can be obtained by imposing the polynomial exactness of the dual set for all lying in the support of the original function. In other terms, if the set of is such that
(12) 
then the multipole preserving property would be guaranteed.
However, we would also like to generalize the collocation method, rather than to replace it by a new quadrature formula. Firstly, we want to impose that for arbitrarily small grid spacings, . Moreover, a notable advantage of the collocation method is its closure with respect to function products. In other terms, given two collocated functions, , , with collocation coefficients , , the product function should satisfy
(13) 
This is possible only if the dual functions are such that
(14) 
for all the points where the coefficients of the function are not zero. Note that in the traditional collocation, when , Eq. (14) is always satisfied due to the interpolating property of , i.e. . That proves that the above property is a feature of the collocation method.
We therefore search for a family of biorthogonal functions , , such that all the following properties hold:
 Collocation coefficients

For dense grid spacings, the action of should coincide with the Dirac distribution;
 Biorthogonality

Even though it is not be necessary that , for any , at least the “weak” duality relation
(15) should be verified for the points where .
 Closure wrt products

The triple product relation of Eq. (14) should hold;
 Polynomial exactness

The functions would guarantee in this way the multipolepreserving property.
In Sec. A of the Appendix, we demonstrate that all these properties are met for the quadrature formula
(16) 
where is a family of ISF of order . The coefficients defined by Eq. (16) may therefore be used at the place of the function point values . With this choice, we are guaranteed to preserve the first multipoles of during the discretization procedure. For grid spacings which are small enough, we recover the usual behaviour of the collocation method thanks to Eq. (32). Fig. 3 provides evidence of this.
When the coefficients are sensibly different from , the values of the coefficients should be rather considered as quadrature terms. However, we might interpret these coefficients as an optimal generalization of the collocation method, suitable for grid spacings that are larger than the oscillation of functions we would like to discretize. This generalization is optimal in the sense that the loosening of the accuracy in the multipoles of order , which is unavoidable for large grid spacings, still exhibits the correct convergence ratio .
It is easy to see that, by using a threedimensional separable ISF basis
(17) 
our method can be generalized straightforwardly to threedimensional grids, especially for separable functions. In the following section we will illustrate the advantage of this method for electronic structure calculations.
Iv Consequences in electronic structure calculations
We illustrate our idea using the BigDFT code BigDFT () which is based on Daubechies wavelets to express the electronic wavefunctions and on interpolating scaling functions for the electronic density and the KohnSham potential. All tests in this article will be done with the LDA functional. We have checked that these results are the same for different functionals and also for the HartreeFock approach.
The BigDFT code has an adaptive mesh with one level of refinement and the corresponding parameter hgrid specifies the grid spacing of the coarse resolution. The finer resolution which is only used near the nucleus so has a twice finer grid step by construction. As mentioned in Sec. II, normconserving GTHHGH pseudopotentials HGHWilland () are used in the BigDFT code. They are built with Gaussian functions for and . Using the collocation method, BigDFT needs a grid step of the order of the standard deviation parameter of the Gaussian function. As an example, for the case of the hydrogen atom, has a value of atomic units that obliges to use a grid spacing of the same value i.e. . In the case of BigDFT, this means that the input parameter hgrid should be of the order of AU in order to have the finer mesh of the same resolution as the Gaussian standard deviation parameter.
iv.1 Accuracy in absolute energy
In Fig. 4, we show the percentage of the difference of the total energy from the reference calculation with hgrid in function of the grid spacing. The quadrature formula (16) has been used with ISF of order . For a grid step greater than , the collocation method gives an error bigger than %, drastically increasing. On contrary, our multipole preserving method is more stable given an error of % for hgrid= with an accuracy of two orders of magnitude up to a grid step of which corresponds to five times the value of the Gaussian function.
Another intrinsic artifact of the realspace methods, especially the finite difference scheme or any methods which uses collocation technique, is the eggbox error eggbox1 (); eggbox2 (). The discretization procedure is not invariant by global translations and rigid rotations. So different relative positions of the ions with respect to the mesh of the simulation domain change the discretized values of . This eggbox effect can be considerably reduced by our method, for large grid spacings. In the figure 5, we have plotted the maximum variation of the absolute energy of a H molecule when rotating its main axis, while keeping fixed interatomic distance. As the relative positions of the centers of the pseudopotentials vary strongly with the molecule orientation, this is a good (even severe) estimator of the eggbox error. We can see that the multipole preserving method decreases considerably the eggbox error in an intermediate range between and where the eggbox energy becomes bigger than meV/atom. Not surprisingly, the traditional collocation method totally fails to capture the rotation of the H2 molecule with a very large error.
We would like to point out that the quasivariational behaviour, typical of any realspace method based on analytic potentials, is greatly improved for large grid spacings. This point is interesting because it is possible to reduce the degrees of freedom by increasing the grid spacings, therefore decreasing the accuracy, but without spoiling the physical meaning of the results. This could be an interesting alternative to tightbinding methods based on DFT which need finetuning of many parameters to represent properly a molecule with an accuracy of the order of meV per atom. In the upper x axis on the same figure 5, we indicate the number of basis set functions i.e. the number of degree of freedoms used for each grid step. These numbers vary as the inverse of the cube of the grid step and the CPU time is directly proportional to the number of the basis set functions. Thus it is very advantageous to increase the grid step parameters to save a lot of CPU time if the method is numerically stable especially when we explore the potential Energy Surface (PES) with some methods as Minima Hopping MHopping () or the Activation Relaxation Technique ART ().
iv.2 Accuracy in Geometry optimization
To illustrate this idea to have a reliable PES, we have optimized, in the figure 6, the H molecule in function of the grid spacing using the FIRE FIRE () method because this method is robust when the eggbox error becomes large. Below a value of for the grid spacing which corresponds to times the standard deviation value of the Gaussian of the pseudopotential, the collocation method still provides manageable results, giving an error less than angstroem for the bond length. For a larger value of the grid spacing, the collocation method becomes strongly unstable. On contrary, our method is stable numerically on a wider range, even though, clearly, the precision of the results decreases.
For an exploration of the PES with an accuracy of the order of angstroem a grid spacing of the order of 1 AU permits to decrease the number of basis set function by two order of magnitudes (see Fig. 5) accelerating in the same range. Such lowaccuracy setup might be sufficient for prescreening calculation in workflows like high throughput funnels. In the same spirit as high throughput screening, using the same input parameters and the same code, varying only the grid spacing i.e. the accuracy, could constitute an interesting approach to have a preliminary investigation of interesting minima and saddle points.
Finally, the multipolepreserving scheme permits to use larger grids space to calculate organic molecules containing hydrogen but also carbon, nitrogen and oxygen atoms where the coefficients of the pseudopotentials are between and . In Fig. 7 and Fig. 1, we show that the total energy of the cinchonidine molecule is numerically stable in function of the grid step.
V Collocated smoothening procedure for a compactly supported function
In this admittedly more technical section we show that our quadrature formula might be used as a reliable approach to calculate the scalar product between our original function and a compactly supported function . This constitutes an important consequence of our method when known functions have to be represented in compactlysupported basis sets.
v.1 Lagrange polynomial as the direct basis
Let us now come back to the approximated function defined in (7). On a support ranging on a interval , in Sec. A of the Appendix, we show that a reliable approximation of can be given by defined as
(18) 
where are interpolating Lagrange polynomials, defined as in Eq. (35)
These properties are interesting any time we have to perform the scalar product of the approximated function with a compactly supported function , of support in the interval . Indeed, the integral
(19) 
will be exact if is a polynomial function of order less than . This happens because also the Lagrange polynomials satisfy polynomial exactness in interval.
The above results, can be also provided in a more general sense. Given an arbitrary function of compact support, ranging from to , and its analytic moments , the coefficients
(20) 
might help us in defining a “smoothed” function
(21) 
such that
(22) 
where the coefficients are calculated by Eq. (16) and can be replaced by only when is small enough. By construction, Eq. (22) is exact for polynomial functions of order less than .
v.2 The “Magic Filter” method in BigDFT code
The above result has been already pointed out in the context of the socalled “MagicFilter” method, providing an optimal set of quadrature coefficients to convert a function expressed in the Daubechies wavelets basis into its point values on a uniform grid. In Sec. B of the Appendix we show that the usual procedure for calculating these filters leads to equation (20).
Such method is very useful in the BigDFT code to express efficiently the local potential matrix elements in Daubechies wavelet basis. In the BigDFT code, a KohnSham orbital is expressed in Daubechies wavelets
(23) 
As pointed out before, the Daubechies wavelet basis is a compact support orthogonal basis able to exactly express the polynomials up to a give order . However, the basis functions
(24) 
have the peculiar property of being less smooth than an order polynomial. Therefore, the naive collocated values would not provide an efficient discretization.
We can therefore use Eq. (21) to define an improved collocation method. We have seen that it is more accurate to express the KohnSham orbitals in real space by the following smoother function:
(25) 
whose collocated values would preserve the multipoles of , that can be derived from the original expression in terms of coefficients. This function is expressed in terms of the smoothed Daubechies scaling functions , plotted in Fig. 8. Therefore, when starting from the coefficients , the real space values of might be given by the formula
(26) 
The results of this paper might be useful to define the inverse relation. Given a set of realspace point values , these coefficients might be interpreted as generalized collocation values. With this interpretation we are able to write the piecewise polynomial expansion of the KohnSham orbital valid on a interval of size around a grid point :
(27) 
This piecewise polynomial function would have the expansion coefficients in Daubechies wavelets basis given by the equation
(28) 
This result shows that the “Magic Filter” method can be seen as the optimal passage matrix between the Daubechies wavelet basis and a realspace description in a generalized collocation scheme. As shown in MagicFilter (); BigDFT () Eqs. (26) and (28) show that this passage matrix is unitary up to . If the KohnSham orbital would be analytically known, we could apply Eq. (22), and our generalized collocation scheme would reduce to Magic Filter method. For this reason, even though it is generally applicable, the multipolepreserving quadrature presented in this paper perfectly reconciles with a Daubechies wavelets computational treatment equipped with the Magic Filter method.
Vi Conclusion
Collocation method is a universally applicable prescription for the numerical discretization of functions. However it suffers from an intrinsic limitation: highly oscillating functions cannot be well represented on a grid if the spacing is too large with respect to the typical length of the oscillations. Therefore the accuracy of the collocation is rapidly spoiled as soon as the grid spacing becomes too large. Unstable results might occur if the numerical implementation is done in such a regime. This limitation implies that there is a upper limit for the grid spacing in a real space based DFT code, and consequently a lower limit for the number of computational degrees of freedom. Results become rapidly meaningless when these limits are overcome.
With this paper, we have presented a method to generalize the collocation of arbitrary analytic function on large grid spacing, without spoiling the accuracy of the discretization. The collocation values might be replaced by the scalar products of the analytic function with the basis of Interpolating Scaling Functions. For analytic and separable functions like Gaussians, this prescription is very simple and easy to implement in three dimensions (see e.g. PSfreeBC ()), and tends to the point values when the grid is fine enough. This method has been implemented in the BigDFT code, which uses a realspace based description using Daubechies wavelets. Thanks to the inclusion of this method, the code exhibits numerical stability over a wide range of grid spacings, not accessible with traditional collocation. However, the implementation of this method is unrelated to Daubechies wavelet basis set, and can be used in other DFT codes and even in different contexts, like for example the definition of compensating charges in Fast Multipole Methods.
Having said that, we would like to point out that our method is focused in providing stabilisation of lowaccuracy results, that would be unaccessible if traditional collocation is used. If a highaccuracy calculation is needed, the grid spacing has to be adjusted in the required range. Nonetheless, as it preserves numerical stability, our multipolepreserving approach is guaranteed to be better than traditional collocation.
The outcomes of this method are very important, as it enables us to use larger grid spacings even for hard pseudopotentials and to perform coarsegrained DFT calculations. As the large grid behaviour of the code is highly stabilized with this method, the user is now able to perform lowaccuracy DFT calculations with reduced number of degrees of freedom. This is fundamental in view of rapid exploration of the energetic features of a system at DFT level, or to accelerate the convergence of iterative molecular calcualtions Tempkin (). Notable examples are the Potential Energy Surface explorations of systems at the nanoscale, as well as the recently established field of highthroughput calculations for material design.
Acknowledgements.
The authors thanks I. Duchemin, V. Perrier and S. Bertoluzza for useful discussions.Appendix A Interpolating scaling functions
Interpolating scaling functions (ISF) lazy () arise in the framework of wavelet theory daub (); ppur (). They are onedimensional functions, and their main properties are:

The full basis set can be obtained from all the translations by a certain grid spacing of the mother function centered at the origin.We indicate the basis set with
(29) 
The mother function is symmetric, with compact support from to . It satisfies the interpolating property .

They satisfy the refinement relation
(30) where the ’s are the elements of a filter that characterizes the wavelet family, equal to in the case of ISF, and is the order of the scaling function. Eq. (30) establishes a relation between the scaling functions on a grid with grid spacing and the ones defined higher resolution level with spacing .

The filters in Eq. (30) are defined such that (as proven in Ref. beyl ()) that the lowest moments of the scaling function are all vanishing but the first, i.e.:
(31) This enables us to show that (see e.g. PSfreeBC ())
The ISF families exhibit polynomial exactness: indeed, the socalled lifting procedure allows to define a set of functions – the lifted wavelets – that are both orthogonal to and have (at most) vanishing moments. As the multiresolution basis formed by the at lowest resolution and the lifted wavelets at all the resolution levels forms a complete set, this proves that the basis of the exhibits polynomial exactness up to order . Fig. 9 shows an interpolating scaling function , together with the corresponding lower resolution lifted wavelets.
The polynomial exactness of the basis, together with their compact support, allows us to demonstrate the collocation property: we can demonstrate that
(32) 
where we have expressed in terms of the dimensionless unit . To show this, let us suppose that the function can be well approximated by its Taylor polynomial close to the point :
(33) 
which would lead, thanks to Eq. (31), to
(34) 
which proves Eq. (32).
We have seen that for a grid spacing which is small enough, the coefficients approach the collocation coefficients . Therefore in this case any set of interpolating functions with polynomial exactness of order would be a good choice for , as the corresponding would be similar up to . As a matter of fact, ISF basis can be treated as an orthogonal basis, even this is not true: nonorthogonality effects are visible only up to .
However, even in this case, it might be useful to identify the direct basis which better generalizes the collocation approach. The Lagrange polynomials
(35) 
might constitute a basis of interest. The matrix is the inverse of the Vandermonde matrix, i.e.
(36) 
In the following we show by using this basis set as the direct interpolating basis, we recover the biorthogonality and the triple product property (Eq. (14)).
Let us consider a set of discretization points ranging from to . The properties of ISF imply
(37)  
(38) 
Therefore, for all the lying in the interval , the basis set constitutes a biorthogonal basis generalizing the collocation method to ^{a}^{a}aThe triple product formula of Eq. (38) would require a ISF family of order 2n to be exact, however the accuracy is still limited by the given by the Lagrange polynomial..
Appendix B Magic Filter, the original idea
Given a set of family Daubechies scaling functions centered on a uniform mesh of spacing , the expansion coefficients of a given function in this set are defined as
(39) 
This discretization of the function in this basis set has an algebraic convergence rate. In other terms
(40) 
A wavelet quadrature in this context is based on the idea of approximating the expression above by a collocation formula. In other words, we should define some coefficients , where such that
(41) 
This can be possible only if the above formula would give the exact result when , . By comparing (39) and (41) in the case of polynomials we thus find the equation defining the Magic Filter:
(42) 
which is solved if and only if
(43) 
the Magic Filters are given by Eq.(10) of Neelov and Goedecker paper MagicFilter ():
(44) 
This equation shows that the magic filters can be viewed as the expansion coefficients of the Lagrange polynomials in the Daubechies scaling functions basis, and is identical to Eq. (20). Also the paper from Johnson JohnsonMF () can be used as a reference in this regard.
References
 (1) A. Neelov, S. Goedecker, J. Comput. Phys. 217, 312 (2006)
 (2) L. Genovese, A. Neelov, S. Goedecker, T. Deutsch, S. A. Ghasemi, A. Willand, D. Caliste, O. Zilberberg, M. Rayson, A. Bergman, R. Schneider, J. Chem. Phys. 129, 014109 (2008)
 (3) A. Willand, Y. O. Kvashnin, Á. VázquezMayagoitia, A. K. Deb A. Sadeghi, T; Deutsch, S. Goedecker, J. Chem. Phys. 138, 104109 (2013)
 (4) G. Deslauriers and S. Dubuc, Constr. Approx. 5, 49 (1989).
 (5) I. Daubechies, “Ten Lectures on Wavelets”, SIAM, Philadelphia (1992)
 (6) ) S. Goedecker, “Wavelets and their application for the solution of partial differential equations”, Presses Polytechniques Universitaires Romandes, Lausanne, Switzerland 1998, (ISBN 2880743982)
 (7) N. Saito, G. Beylkin, G., Signal Processing, IEEE Transactions on [see also Acoustics, Speech, and Signal Processing, IEEE Transactions on] Vol. 41 (12), (1993) 3584–3590
 (8) L. Genovese, T. Deutsch, A. Neelov, S. Goedecker, G. Beylkin, J. Chem. Phys. 125, 074105 (2006)
 (9) B. R. Johnson, J. P. Modisette, P; J. Nordlander, J. L. Kinsey, J. Chem. Phys., 110, 8309 (1999)
 (10) F. Jollet, M. Torrent, N. Holzwarth, Comp. Phys. Commun., 185, 4, (2014), 12461254
 (11) E. Artacho, A. Anglada, O. Diéguez, J.D. Gale, A. Garcia, J. Junquera, R. M. Martin, P. Ordejón, J. M. Pruneda, D. SánchezPortal, J. M. Soler J. Phys.: Condens. Matter 20, 064208 (2008)
 (12) T. Ono, M. Heide, N. Atodiresei, P. Baumeister, S. Tsukamoto, S. BlÃ¼gel Phys. Rev. B 82, 205115 (2010)
 (13) J. O. B. Tempkin, B. Qi, M. G. Saunders, B. Roux, A. R. Dinner, J, Weare, J. Chem. Phys. 140, 2184114 (2014)
 (14) B. Schaefer, S. Mohr, M. Amsler, S. Goedecker, J. Chem. Phys. 140, 2014102 (2014)
 (15) E. MachadoCharry, L. K. Béland, D. Caliste, L. Genovese, T. Deutsch, N. Mousseau, P. Pochet J. Chem. Phys. 135, 034102 (2011)
 (16) E. Bitzek, P. koskinen, F. Gähler, M. Moseler, P. Gumbsch, Phys. Rev. Lett. 97, 170201 (2006)
 (17) S. Kotochigova, Z. H. Levine, E. L. Shirley, M. D. Stiles, C.. W. Clark Atomic reference data for electronic structure calculations, NIST, http://www.nist.gov/pml/data/dftdata