{'fraction_non_alphanumeric': 0.034544881322343825, 'fraction_numerical': 0.012610908435796965, 'mean_word_length': 4.410331384015595, 'pattern_counts': {'":': 0, '<': 0, '<?xml version=': 0, '>': 12, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 3, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Electronic structure calculations with interpolating tensor product wavelet basis

We introduce a basis set consisting of three-dimensional Deslauriers--Dubuc wavelets and numerically solve the Schr\"odinger equations of hydrogen atom, helium atom, hydrogen molecule ion, hydrogen molecule, and lithium hydride molecule with HF and DFT methods. We also compute the 2s and 2p excited states of hydrogen. The Coulomb singularity at the nucleus is handled by using a pseudopotential. Results are compared with those of CCCBDB and BigDFT. The eigenvalue problem is solved with Arnoldi and Lanczos methods, and the Poisson equation with GMRES and CGNR methods. The various matrix elements are computed using the biorthogonality relations of the interpolating wavelets.


Introduction
Standard approaches to assess properties of atoms, molecules and models of nanostructures in quantum chemistry are Hartree-Fock (HF) and Density Functional Theory (DFT). Both of these, invoke numerical solutions of the Schrödinger differential equations of the many-body system of electrons. Thus, controlled approximations are inevitable and practical numerical algorithms are necessary.
Vast majority of the algorithms are based on finding the solutions or orbitals as series expansion of basis functions, the basis set. A finite basis set leads to the Roothaan-Hall equations, a generalized matrix eigenvalue problem. Gaussian type basis functions are the most popular, due to their advantageous analytical features. Here, we introduce a different type basis set, wavelet functions, point out their advantages and drawbacks, and consider a few test cases and compare their performance with conventional approaches.
Wavelets and scaling functions are functions generated by translations and dilatations of the so called mother wavelet and mother scaling function. Interpolating wavelets use a mother scaling function satisfying the cardinal interpolating property ϕ(k) = δ k,0 where k is an integer. Orthonormal wavelets form an orthonormal basis in function space L 2 (R). Both of these wavelet types can be generalized to multivariate functions.
One-dimensional interpolating wavelets in function space C u (R) consisting of bounded and uniformly continuous functions on R are defined in Ref. [1]. One-dimensional interpolating wavelets in function space C 0 (R) consisting of functions on R vanishing at infinity are defined in Ref. [2]. Deslauriers-Dubuc wavelets have also been discussed in Refs. [3] and [4]. Compactly supported interpolating wavelets have been generalized to multiple dimensions in Refs. [5] and [6]. Fukuda, Kinoshita, and Suzuki [7] have studied unconditional convergence of wavelet expansions. They have shown that uniformly convergent wavelet expansions even for continuous functions do not always converge unconditionally in L ∞ (R). Pathak [8] has investigated translation and convolution associated with discrete wavelet transform.
Arias [9] and Engeness and Arias [10] developed formalism for electronic structure calculations with interpolating wavelets so that matrix elements of the operators are computed as usual and overlap matrices are used in the matrix form of the Schrödinger equation. Lippert et al. [11] introduced a Lagrangian based formalism for the multiresolution analysis (MRA) of electronic structure. Arias [9] uses carbon atom and N 2 molecule as examples. Engeness and Arias [10] use both calcium and aluminum atoms and molecules O 2 and H 2 O as examples.
Fischer and Defranceschi use Daubechies wavelets [12] for computation of hydrogenlike atoms [13]. They developed an iterative method based on nonstandard operator form of the Schrödinger operator. Their work shows that this method is well suited for computations of hydrogenlike atoms. Fischer and Defranceschi have presented Hartree-Fock equations in an orthonormal wavelet basis [14]. They also analyzed the Hartree-Fock method with continuous wavelet transform [15] and demonstarted it to the hydrogenic Schrödinger equation with an iterative solution scheme.
Wei and Chou [16] used orthonormal wavelets in self-consistent electronic structure calculations within the local-density approximation and demonstrated with H 2 and O 2 molecules. Tymczak and Wang [17] used orthonormal Daubechies wavelets for quantum molecular dynamic simulations and developed a wavelet selection scheme for computations. They used hydrogen atom and H 2 molecule as examples. Their method showed systematic convergence with increasing grid size. Yamaguchi and Mukoyama [18] carried out electronic structure calculations with the Hartree-Fock method and Meyer wavelets.
Our earlier contribution [19] is use of one-dimensional interpolating wavelets to solve the Hartree-Fock equations in the Central Field Approximation for orbitals of several test case atoms. We were able to derive analytic formulas for all the relevant matrix elements of Hamiltonian and Fock operator.
Iyengar and Frisch [20] studied relationship between Gaussian basis sets and wavelets. They use a time-dependent basis function set: Gaussian functions centered at nuclei of the system. When the nuclei move the basis functions move, too. Gaussian functions are an example of multiwavelets [21,20]. Comparative studies of wavelets and Gaussian functions are presented e.g. in [22] and [23].
Han, Cho, and Ihm [24] have developed an all-electron density-functional program using the Mexican hat wavelets. They analyze H 2 , CO, and H 2 O molecules and 1s core-ionized C*O and CO* molecules. Their method shows very good performance over the plane-wave based methods. Genovese et al. [25] and Mohr et al. [26] have composed a software package BigDFT that implements the DFT method for quantum physical systems using three-dimensional Daubechies wavelets as a basis function set. Here, in this study, we use the data from BigDFT as the reference for that from our interpolating tensor product Deslauriers-Dubuc wavelet basis function set.
In three dimensional space the number of basis functions may grow relatively large to give sufficient accuracy. Therefore, we do not construct the Hamiltonian matrix explicitly. Instead, we use iterative algorithms in solving the Roothaan-Hall equation, though correspondingly, the computation gets slower in self-consistent iteration of Hartree-Fock or DFT orbitals for many electron systems.
We demonstrate of the hydrogen and helium atoms, hydrogen molecule ion, hydrogen molecule, and lithium hydride molecule. Self-consistent iteration and Hartree-Fock and DFT methods are used for many electron systems.
From now on, we use atomic units throughout this article (e = m e = = 4πε 0 = 1). Thus, units for energy and length,"Hartree" and "Bohr" are used and abbreviated as Ha and B. Notations for wavelet basis functions and filters is similar to that in Ref. [5] and computation of matrix elements is similar to that in Ref. [6]. Pseudopotentials are used for two different purposes, to handle the Coulomb singularity at nuclei, and also, as frozen-core.
In this paper, we demonstrate evaluation of electronic structure with threedimensional interpolating tensor product wavelets and use of dual multiresolution analysis in computation of the matrix elements of the various operators.

General
Consider a system with n electrons and m nuclei, with atomic numbers Z i and locations R i . Within the Born-Oppenheimer approximation, dynamics of electrons and nuclei are independent, and then, the wave function separates to two factors, correspondingly. Here, we consider the electronic part, the orbitals, only, and keep the nuclear conformation {R i } fixed.
With the fixed nuclear conformation, the Coulomb potential for dynamics of electrons is and the internuclear repulsion energy For an atom we have E R = 0. Let us denote the orbital by φ, because symbols ϕ and ψ are used for scaling functions and wavelets.
If the spacing between grid points is a (usually negative) power of two we can handle this by choosing j min in equations (24) and (25) properly. Otherwise we have to make a change of variables r = ar ′ in the Schrödinger equation. Here one unit in the computation grid corresponds to a Bohrs.

Single-electron System
The Schrödinger equation of a single electron system is where E = ε and the total energy The wave equation of a single-electron system is solved by the Implicitly Restarted Arnoldi Method [27,28]. The Arnoldi method is able to find also other than the lowest eigenvalue.

Hartree-Fock Method
The HF equation for an n electron system is where the Hartree potential is given by and the charge density by The exchange potentials V i x , i = 1, . . . , n, are defined in [29, section 2.2] but it is not needed in this article. For a two electron system with both electrons in the same orbital the HF equation can be written as The Hartree potential is actually computed by solving the Poisson equation In case where there are more than one resolution level in an interpolating wavelet basis the matrix L of the Laplacian operator is not generally Hermitian so we cannot solve (9) directly with the conjugate gradient method. The nonhermiticity arises because our matrix elements are not computed as ordinary inner products between functions. In this case we use one of the following two methods: • conjugate gradient on the normal equations (CGNR): solve with the conjugate gradient method.
When the basis set consists of a single resolution level we may use ordinary conjugate gradient method to solve (9). The total energy of a two electron system is In this paper, we consider restricted Hartree-Fock approach, only.

Density Functional Theory
Suppose that we have a system with M electronic orbitals whose total wavefunction is Ψ. The Kohn-Sham equation [31,29,32] for the electronic structure is where the charge density is and the Hartree potential and V xc [ρ] is the exchange-correlation potential. We have where V x [ρ] is the exchange potential and V c [ρ] is the correlation potential. In this article we have V c [ρ] = 0. The exchange-correlation energy is defined by where E x [ρ] is the exchange energy and E c [ρ] is the correlation energy. In this article we ignore the correlation energy. Within the Local Density Approximation (LDA) we define where ε x [ρ](r) is the exchange energy per particle of a uniform electron gas at a density of ρ. It follows from the Kohn-Sham theorem [29, section 3.1] that the exchange potential is We have and The total energy of the system is The Kohn-Sham equations are solved by a similar self-consistent iteration as the HF equations.
3 Three-dimensional Wavelet Basis Set

The Basis Set
Let j min and j max be the minimum and maximum resolution levels of the point grid. Let and where j ∈ Z. Define sets Q j by The point grid G shall be some finite subset of V jmax . We define G j := G ∩ Q j for j ≥ j min . The functions ϕ j,k and ψ j,k are scaling functions and wavelets belonging to an interpolating wavelet family. Functionsφ j,k andψ j,k are dual basis functions of interpolating wavelets. Define and if j > j min and k even ψ j−1,(k−1)/2 ; if j > j min and k odd if j > j min and k eveñ ψ j−1,(k−1)/2 ; if j > j min and k odd (29) When α ∈ Q j and j ≥ j min define andζ α :=η j,k [1] ⊗η j,k [2] ⊗η j,k [3] where k = 2 j α. We also define where j ∈ Z and k ∈ Z 3 . We use notation in the definitions of basis functions sets.

Backward and Forward Wavelet Transforms
We define forward wavelet transform U and backward wavelet transform W by setting U (v) = c and W (c) = v. Mappings U and W are linear. We compute the forward wavelet transform U using an algorithm somewhat similar to [33].
where α ∈ G j and β ∈ G and matrix E (j) by where α ∈ G and β ∈ G j . Define where α ∈ G j and β ∈ G j ′ . We have For forward wavelet transform we have for j > j min we have Consequently J and we do not have to invert matrix W (j,j) . An operator representing pointwise multiplication of a function f by another in the given computation grid is where D is a diagonal matrix with values of the function f at the grid points in the diagonal. The local pseudopotentials are computed this way.

Matrix Elements of the Laplacian Operator
Laplacian operator if k x even and j > j min if k y even and j > j min t y = 0, if j = j min or j > j min and k y even 1, if j > j min and k y odd (56) if k z even and j > j min t z = 0, if j = j min or j > j min and k z even 1, if j > j min and k z odd and similar definitions for ι ′ , l ′ x , t ′ x , l ′ y , t ′ y , l ′ z , and t ′ z . Elements of matrices are computed by The filters a and s are defined by Filter a is computed with formulas a(1, 1, j, k) = a(1, 0, j + 1, 2k + 1), for j ≥ 0 (75) Filter s is computed with formulas

Constant and Interpolated Pseudopotentials
The Coulomb potential arising from a single nucleus is where Z is the charge of the nucleus. We avoid the singularity of (90) at the origin by using a pseudopotential. We define c to be the cutoff point of the pseudopotential and D to be the degree of the interpolating polynomial used in the pseudopotential. Actually we use c = 2 −jmax where j max is the highest resolution level in the wavelet basis. Parameter D has to be an odd integer and we define n = (D + 1)/2. We define and where i = 1, . . . , 2n. Let P be the interpolating polynomial of degree at most D having value t[i] at point s[i] for i = 1, . . . , 2n. Now we define the interpolated pseudopotential with The actual pseudopotential of a nucleus with charge Z and location R is then Function V interp with different values of c is plotted in figure 1. We may also use the cut pseudopotential defined by where c = 2 −jmax−1 .

HGH Pseudopotential
By using pseudopotentials the number of computed orbitals can be reduced and since the pseudo wavefunctions are smoother than all-electron wavefunctions the basis function set can also be reduced. In a pseudopotential computation only the valence electrons are actually computed and the effect of the core electrons is handled by replacing the nuclear potential with a pseudopotential. We use the Hartwigsen-Goedecker-Hutter (HGH) pseudopotential [34,25] in these computations. The HGH pseudopotential consists of a local and nonlocal part. The local part is a function of the position as the nuclear potential. The nonlocal part is a linear operator and it is not a function of position. The local pseudopotential is given by  (33).
where r is the distance from the nucleus. Note that and the local pseudopotential is defined at the origin (nucleus), too. The nonlocal pseudopotential is defined by for each nucleus. The origin of the coordinate system in (99) is located at the nucleus. The functions p l i are defined by where parameter r l is given in Bohrs. The range of values l is determined by the actual pseudopotential. The spherical harmonics Y l,m in equation (99) can be replaced by orthonormal linear combinations of Y l,m , m = −l, . . . , l. This allows us to avoid computation with complex valued functions. -0.466582 BigDFT [25,26] HGH -0.499969 exact none -0.5 Table 2: Computation results for the hydrogen atom. Quantity g is the distance between grid points in the highest resolution level. The numbers in column "basis" refer to table 1 and "TH" means this work.             iteration and Arnoldi method. When there is only one resolution level in the basis the Hamiltonian and Laplacian matrices are hermitian and the Arnoldi method reduces to a variant of the Lanczos method. When there is only one resolution level the ordinary conjugate gradient method could be used for the Poisson equation. Two methods have been used for solving the Poisson equation for the general case: CGNR and GMRES. We ran a benchmark for these methods and GMRES was 4.5 times faster when the accuracies of the results were approximately the same. All the computations use 8th order Deslauriers-Dubuc wavelets (polynomial span 7). We use interpolating polynomials of degree 7 for the pseudopotentials. In some computations the Arnoldi method did not find the desired eigenvalues. This was solved by rising the number of computed eigenvalues and the number of basis vectors in the Arnoldi method. The basis function sets (computation point grids) are presented in table 1.

Atomic and Molecular Orbitals
The data presented in the tables in this article was obtained from our own software (denoted by TH), our own computations with BigDFT [25,26], and CCCBDB [35]. BigDFT is a quantum mechanical computation package using Daubechies wavelets. CCCBDB is large database containing atomic and molecular data. The CCCBDB results in this article use the STO-3G basis set. Quantity g is the distance between grid points in the highest resolution level and quantity d is the distance between the nuclei in the result tables. E system is the total energy of the system and E binding is the binding energy. We used grid spacing 0.45 Ha (finer grid spacing 0.225 Ha) in all our BigDFT computations. The BigDFT parameters for determining the size of the basis set were crmult=10.0 and frmult=16.0. The computations were also made with values crmult=5.0 and frmult=8.0 but the results did not differ significantly.
Note that the total energy does not include the energies of the core electrons in the lithium hydride HGH computations. Computation using our own software use interpolating wavelets and BigDFT computations orthonormal Daubechies wavelets.
For the molecular computations the energy of the system as a function of the distance between the nuclei is computed in three points near the energy minimum and a second degree polynomial is fitted into these points. The distance between the nuclei is then the minimum point of the polymial and the energy of the system is computed at the minimum distance. We locate the nuclei at points (0, 0, ± d 2a ) where d is the distance between the nuclei in Bohrs. When binding energies of molecules were computed the energies of atoms were usually computed with the basis set as the molecule. When the basis was unsymmetric it was modified. For example, grid Z(4) × Z(4) × Z(10) becomes Z(4) × Z(4) × Z(4) for the atoms.
The results for hydrogen atom ground state are presented in table 2. The radially averaged ground state wavefunctions of the hydrogen atom are plotted in figures 2 and 3. A radial average of a function f : where r ∈ [0, ∞[. As the angular part of an s-type wavefunction is 1 2 √ π we estimate a radial wavefunction bȳ The computation of hydrogen excited states uses HGH pseudopotential and basis 9. Results are presented in table 3. The names of the excited states were obtained by computing inner products between the computed states and analytical states. The resulting orbitals are approximately orthonormal and the computed 2p states are approximately linear combinations of the analytical 2p states. The largest (in absolute value) inner product between different orbitals is 2p a |2p c = 5.428 × 10 −4 . The quality of the linear combinations can be measured by a quantity 1 − P f 2 2 where P is the orthogonal projection from L 2 (R 3 ) onto the space spanned by 2p x , 2p y , and 2p z . The value of this quantity is 0.1265 for all the computed orbitals 2p a , 2p b , and 2p c . Analytical expressions for hydrogenic orbitals can be found for example in Ref. [39].
Helium atom has been computed using the HGH pseudopotential. Results are presented in table 4. Grid spacing g = 0.5 B does not give sensible results with the HGH pseudopotential. We suppose that Froese Fischer's results [36] can be regarded as the Hartree-Fock limit for helium. The computed total energies of the helium atom are quite good.
We computed the hydrogen molecule with the interpolated and HGH pseudopotentials. The computation results for the hydrogen molecule are presented in table 5. The resulting dissociation curve with the HF method, HGH pseudopotential, and basis set 5 is plotted in figure 4. When the HGH pseudopotential was used the computation worked for grid spacing g = 0.25 B but did not work for spacing g = 0.5 B. The results for the hydrogen molecule ion are presented in table 6.
The computed results for lithium hydride molecule are presented in table 7. Grid spacing g = 0.5 B did not yield a physical dissociation curve for the HGH pseudopotential. Neither g = 0.5 B nor g = 0.25 B yielded a physical dissociation curve for the interpolated pseudopotential. The dissociation curve of lithium hydride computed with the Hartree-Fock method and basis set 5 is plotted in figure 5.

Conclusions
We have shown how to solve the wave equations of hydrogen and helium atoms, hydrogen molecule ion, and hydrogen and lithium hydride molecules in a threedimensional interpolating tensor product wavelet basis. As far the authors know only Arias [9] and Engeness and Arias [10] have done this before. However, they do not use the dual interpolating MRA to evaluate matrix elements. We do that and it allows us to neglect the overlap integrals of the basis functions.
It seems to require large basis sets to obtain numerically good results for quantum physical systems. Roughly, the computations require at least ten thousand basis functions. The most accurate computed bond length of the hydrogen molecule is good and the energy satisfactory. The most accurate binding energies and internuclear distances for the hydrogen molecule ion in table 6 are very accurate. HGH computations performed very well with grid spacing g = 0.25 B but not with g = 0.5 B. The same phenomenon was observed with BigDFT, too.
The BigDFT result remained approximately the same when the basis set was doubled in one dimension (i.e. the basis set was enlarged about eight times). This shows that the basis function sets in the BigDFT computations were large enough. We also found that a one level basis set can be replaced with considerably smaller two level basis set without a significant effect on the results.
It is best to use wavelet basis sets consisting of one or two resolution levels because we have then ψ s ′ ,j,k ′ , ψ s,j,k = δ s ′ ,s δ k ′ ,k for all the one-dimensional dual basis functions and basis functions, which simplifies the computation of the Laplacian operator.