# On basis set optimisation in quantum chemistry

Eric Cancès, Geneviève Dusson, Gaspard Kemlin, Laurent Vidal

► To cite this version:

Eric Cancès, Geneviève Dusson, Gaspard Kemlin, Laurent Vidal. On basis set optimisation in quantum chemistry. ESAIM: Proceedings and Surveys, 2023, 73, pp.107-129. 10.1051/proc/202373107. hal-03736412v2

HAL Id: hal-03736412

https://inria.hal.science/hal-03736412v2

Submitted on 1 Dec 2022

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons CC BY 4.0 - Attribution - International License
ON BASIS SET OPTIMISATION IN QUANTUM CHEMISTRY

ERIC CANCÈS<fnref n="1" />, GENEVIÈVE DUSSON<fnref n="2" />, GASPARD KEMLIN<fnref n="3" /> AND LAURENT VIDAL<fnref n="4" />

**Abstract.** In this article, we propose general criteria to construct optimal atomic centered basis sets in quantum chemistry. We focus in particular on two criteria, one based on the ground-state one-body density matrix of the system and the other based on the ground-state energy. The performance of these two criteria are then numerically tested and compared on a parametrized eigenvalue problem, which corresponds to a one-dimensional toy version of the ground-state dissociation of a diatomic molecule.

**Résumé.** Nous proposons dans cet article des critères généraux pour construire des bases atomiques localisées optimales en chimie quantique. Nous nous intéressons en particulier à deux critères: l'un basé sur la matrice densité à un corps du système, l'autre basé sur l'énergie de l'état fondamental du système. Les performances de ces deux critères sont évaluées et comparées sur un problème aux valeurs propres paramétré, correspondant à une version simplifiée du problème de dissociation de l'état fondamental d'une molécule diatomique.

# 1. INTRODUCTION

In quantum chemistry, a central problem is the computation of the electronic ground-state (GS) of a given molecular system. For many-electron systems, it is not possible to solve the N-body Schrödinger equations and most calculations are thus based on variational (e.g. Hartree-Fock) or non-variational (e.g. coupled cluster) approximations of the latter, or on Kohn-Sham density functional theory (DFT). For all these models, the continuous equations (e.g. a nonlinear elliptic eigenvalue problem in the Hartree-Fock or Kohn-Sham settings) are discretized into a finite-dimensional approximation space. Approximation spaces constructed from atomic orbitals (AO) basis sets [15,23] have many advantages and are therefore the most common choice in the quantum chemistry community. An AO basis set consists of a collection of functions $\chi = (\chi_\mu^z)_{z \in \text{CE}, 1 \le \mu \le n_z}$ where CE is a set of atomic numbers (e.g. CE = {1, ..., 92} for the natural chemical elements of the periodic table), $n_z$ a positive integer depending on the electronic shell-structure of the chemical element with atomic number z, and $\chi_\mu^z \in H^1(\mathbb{R}^3)$ a fast decaying function centered at the origin called an atomic orbital. Consider an atomic configuration $\omega$ consisting of M nuclei with nuclear charges $z_1, ..., z_M$ (in atomic units) and positions $\mathbf{R}_1, ..., \mathbf{R}_M$ in the three dimensional physical space. If the AO basis set $\chi$ is chosen by the user, the (spatial component of the) one-electron finite-dimensional space in which the chosen electronic structure model of a molecular system with atomic configuration $\omega$ is discretized is

$$ \mathcal{X}_{\omega} := \operatorname{span}(\chi_{1}^{z_{1}}(\cdot - \mathbf{R}_{1}), \dots, \chi_{n_{z_{1}}}^{z_{1}}(\cdot - \mathbf{R}_{1}), \dots, \chi_{1}^{z_{M}}(\cdot - \mathbf{R}_{M}), \dots, \chi_{n_{z_{M}}}^{z_{M}}(\cdot - \mathbf{R}_{M})). $$

<footnote n="1">¹ CERMICS, Ecole des Ponts and Inria Paris, 6 & 8 avenue Blaise Pascal, 77455 Marne-la-Vallée, France, e-mail: eric.cances@enpc.fr</footnote>

<footnote n="2">² Laboratoire de Mathématiques de Besançon, UMR CNRS 6623, Université Bourgogne Franche-Comté, 16 route de Gray, 25030 Besançon, France, e-mail: genevieve.dusson@math.cnrs.fr</footnote>

<footnote n="3">³ CERMICS, Ecole des Ponts and Inria Paris, 6 & 8 avenue Blaise Pascal, 77455 Marne-la-Vallée, France, e-mail: gaspard.kemlin@enpc.fr</footnote>

<footnote n="4">⁴ CERMICS, Ecole des Ponts and Inria Paris, 6 & 8 avenue Blaise Pascal, 77455 Marne-la-Vallée, France, e-mail: laurent.vidal@enpc.fr</footnote>
The accuracy of the approximation therefore crucially depends on the quality of the AO basis set. The main advantage of AO basis sets is that only a small number of AO per atoms (typically a dozen) are necessary to obtain a relatively accurate result on most quantities of interest. This is in sharp contrast with standard discretization methods used in the simulation of partial differential equations such as finite-element methods. To make connection with discretization methods used in mechanical and electrical engineering, AO basis set discretization methods can be considered as spectral methods [9], and share common features with the modal synthesis method [11, Chapter 7], [10]. A drawback of AO basis sets is that conditioning quickly blows up when increasing the size of the basis by including polarization and diffuse basis functions, a problem known as overcompleteness [18]. The numerical errors due to this large condition number can deteriorate the accuracy of the computed solutions and/or significantly increase computational times. AO basis sets can therefore not be systematically improved in a straightforward way.

In the early days, AOs were Slater functions [29], with exponential decay and a cusp at the origin. It was then realized by Boys [7] in 1950 that it was much more efficient from a computational viewpoint to use Gaussian-type orbitals (GTO), that are linear combinations of polynomials times Gaussian functions. Indeed the multi-center overlap, kinetic and Coulomb integrals

$$ \int_{\mathbb{R}^3} \chi_i^{z_a}(\mathbf{r} - \mathbf{R}_a) \chi_j^{z_b}(\mathbf{r} - \mathbf{R}_b) \mathrm{d}\mathbf{r}, $$

$$ \int_{\mathbb{R}^3} \frac{\chi_i^{z_a}(\mathbf{r} - \mathbf{R}_a) \chi_j^{z_b}(\mathbf{r} - \mathbf{R}_b)}{|\mathbf{r} - \mathbf{R}_k|} \mathrm{d}\mathbf{r}, $$

$$ \int_{\mathbb{R}^3} \nabla \chi_i^{z_a}(\mathbf{r} - \mathbf{R}_a) \cdot \nabla \chi_j^{z_b}(\mathbf{r} - \mathbf{R}_b) \mathrm{d}\mathbf{r}, $$

$$ \int_{\mathbb{R}^3 \times \mathbb{R}^3} \frac{\chi_i^{z_a}(\mathbf{r} - \mathbf{R}_a) \chi_j^{z_b}(\mathbf{r} - \mathbf{R}_b) \chi_k^{z_c}(\mathbf{r}' - \mathbf{R}_c) \chi_l^{z_d}(\mathbf{r}' - \mathbf{R}_d)}{|\mathbf{r} - \mathbf{r}'|} \mathrm{d}\mathbf{r} \mathrm{d}\mathbf{r}', $$

arising in discretized electronic structure models can then be computed analytically by means of explicit calculations and recursion formulas.

However, individual Gaussian function poorly describes the cusps of the bound states electronic wavefunctions at nuclear positions. Contracted Gaussians [20], that are linear combinations of Gaussians with different variances, were quickly introduced as they overcome this deficiency. Several classes of GTO basis sets have been proposed since the 50's: STO-ng basis sets [14] were built as the contraction of n Gaussians that fit Clementi STO SCF AOs in an $L^2$ least-squares sense [30]. It was quickly realized that better GTO basis sets could be obtained by minimizing atomic Hartree-Fock ground state energy. This approach led to the split-valence basis sets (e.g. 6-31G) with core and valence orbitals being approximated differently, developed by Pople et al. [5]. Basis sets better suited for correlated electronic structure methods were then introduced, notably Atomic Natural Orbitals (ANO) [2] and Dunning basis sets [13]. ANO basis sets are built by selecting occupied and virtual orbitals from Hartree-Fock and natural orbitals from correlated computations of atomic systems. Dunning bases provide a (finite) hierarchy of bases obtained by consistently increasing the number of basis functions corresponding to different angular momenta. This optimization strategy yields the so-called correlation consistent cc-pVXZ basis sets, which are, with their augmented version, still commonly used nowadays.

Mathematical studies proving convergence rates or proposing systematic enrichment of GTO basis sets are so far quite limited. The approximability of solutions to electronic structure problems by Gaussian functions was studied in [16], and later on in [27, 28]. An a priori error estimate on the approximation of Slater-type functions by Hermite and even-tempered Gaussian functions was derived in [3]. A construction of Gaussian bases combined with wavelets was proposed on a one-dimensional toy model in [25].

Commonly used Pople and Dunning GTO basis sets were optimized from atomic configuration energies and Hartree-Fock (and/or natural) atomic orbitals. Let us also mention [12, 26] where system specific optimization of AO bases has been investigated, however focusing on specific models (e.g. one electron periodic Hamiltonian) or optimization criteria. In this article, we propose a different approach, which is adaptable to any criterion one might be interested in, and involves molecular configurations. In Section 2, we introduce an abstract mathematical framework for the construction of optimal AO basis sets, based on the choices of

(1) a set of admissible atomic configurations $\Omega$;
(2) a probability measure $\mathbb{P}$ on $\Omega$;

(3) a set of admissible AO basis sets $\mathcal{B}$;

(4) a criterion $j(\chi, \omega)$ quantifying the error between the exact values of the quantities of interest when the system has atomic configuration $\omega \in \Omega$ – for the continuous model under consideration – and the ones obtained after discretization in the basis set $\chi \in \mathcal{B}$.

We also provide examples of possible choices of $\Omega$, $\mathbb{P}$, $\mathcal{B}$, and $j$. As a proof of concept (Section 3), we apply this strategy to a simple toy model of a 1D homonuclear diatomic “molecule” with two 1D non-interacting spinless “electrons”, which allows for extremely accurate reference calculations. Finally, we present numerical results in Section 4, where we compare the efficiency of the so-optimized AO bases compared to AO basis constructed from the occupied and unoccupied orbitals of the isolated “atom”.

## 2. OPTIMIZATION CRITERIA

### 2.1. Abstract framework

We start by formulating the problem of basis set optimization in an abstract setting. The procedure can be divided into four steps.

First, we select the set $\Omega$ of all possible atomic configurations we are a priori interested in. For instance, depending on the foreseen applications, one can consider the set of all possible finite atomic configurations containing only hydrogen, nitrogen, carbon, and oxygen atoms, or the set of all possible periodic arrangements of chemical elements with less than 20 atoms per unit cell.

Second, we equip $\Omega$ with a probability measure $\mathbb{P}$ in order to allow for different configurations to have different weights in the optimization procedure. We will see later that our method requires the computation of very accurate reference solutions for all $\omega$'s in the support of $\mathbb{P}$. For practical reasons we therefore need to choose $\mathbb{P}$ of the form

$$ \mathbb{P} = \sum_{n=1}^{N_c} \beta_n \delta_{\omega_n}, \qquad (1) $$

where $\{\omega_1, \dots, \omega_{N_c}\}$ is a finite (not too large) subset of $\Omega$, $\delta_{\omega_n}$ the Dirac mass at $\omega_n$, and $\{\beta_1, \dots, \beta_{N_c}\}$ are positive weights such that $\sum_{n=1}^{N_c} \beta_n = 1$. Assume that we are solely interested in reproducing accurately the dissociation curve of the HF (Hydrogen Fluoride) diatomic molecule. Then the set $\Omega$ should be identified with the interval $(0, +\infty)$, and a configuration $\omega \in \Omega$ with the H–F interatomic distance $R \in (0, +\infty)$, and $\mathbb{P}$ should be a probability measure on the interval $(0, +\infty)$. The selection of the $\omega_n$'s and $\beta_n$'s can be done in various ways. An option is to

i) choose a continuous probability distribution $\mathbb{P}$ on $(0, +\infty)$ on the basis of chemical arguments, putting little weight on usually unimportant very small interatomic distances, more weight on interatomic distances close to the equilibrium distance ($d \approx 0.92$ Å), sufficient cumulated weight on very large interatomic distance to correctly reproduce the dissociation energy, and more or less weight on intermediate interatomic distances in the range 2–8 Å, depending on its importance for the targeted application;

ii) fix the number $N_c$ according to the available computational means;

iii) compute the $\omega_n$'s and $\beta_n$'s using e.g. quantization algorithms [21] possibly based on optimal transport or clustering algorithms [24].

Third, we select the set $\mathcal{B}$ of admissible AO basis sets. Restricting ourselves to the framework of GTOs, this can be done by choosing, for each chemical element arising in $\Omega$, the number, symmetries, and contraction patterns of the Gaussian polynomials of the AO associated with this particular element. In this case, $\mathcal{B}$ has the geometry of a convex polyedron of $\mathbb{R}^d$.

Given an atomic configuration $\omega \in \Omega$ and an AO basis set $\chi \in \mathcal{B}$, we denote by $\chi_{\omega}$ the one-electron finite-dimensional space obtained by using the AO basis set $\chi$ to describe the electronic structure of a molecular system with atomic configuration $\omega$ and an arbitrary number $N$ of electrons.
The fourth and final step consists in choosing a criterion $j(\chi, \omega)$ quantifying the quality of the results obtained when using the basis set $\chi \in \mathcal{B}$ to compute the electronic structure of a molecular system with atomic configuration $\omega$. The choice of the function

$$j : \mathcal{B} \times \Omega \to \mathbb{R}_+$$

depends on the quantity of interest (QoI) to the user, and on the respective weights of these quantities in the case of multicriteria analyses. For instance, if one focuses on the ground-state energy of the electrically neutral molecular system, a natural criterion is

$$j_E(\chi, \omega) := |E_\omega - E_\omega^\chi|^2, \qquad (2)$$

where $E_\omega$ is the exact ground-state energy of the neutral system with atomic configuration $\omega$ for the chosen continuous model (e.g. Hartree-Fock, MCSCF, Kohn-Sham B3LYP...) and $E_\omega^\chi$ the ground-state energy obtained with the model discretized in the AO basis set $\chi$. Note that the absolute value of the difference is squared to make $j_E$ differentiable. Another possible choice is to use a criterion based on the one-body reduced density matrices (1-RDM), for instance

$$j_A(\chi, \omega) := -\mathrm{Tr}(\Pi_{\chi\omega}^A \gamma_\omega \Pi_{\chi\omega}^A A), \qquad (3)$$

where $A$ is a given self-adjoint, positive, definite operator on the one-particle state space $\mathcal{H}$ with form domain $Q(A)$, $\gamma_\omega$ the exact ground-state 1-RDM of the neutral system with atomic configuration $\omega$ for the chosen continuous model, and $\Pi_{\chi\omega}^A : Q(A) \to \chi_\omega \subset \mathcal{H}$ the orthogonal projector on $\chi_\omega$ for the inner product $A$ on $Q(A)$. If $A = I_\mathcal{H}$, then the $Q(A) = \mathcal{H}$ and $\Pi_{\chi\omega}^A$ is the orthogonal projector on $\chi_\omega$ for the inner product of $\mathcal{H}$. If $A = (1 - \Delta)$, then $Q(A)$ is the Sobolev space $H^1(\mathbb{R}^3)$, and $\Pi_{\chi\omega}^A$ is the orthogonal projector on $\chi_\omega$ for the $H^1$-inner product. Diagonalizing $\gamma_\omega$ as

$$\gamma_\omega = \sum_j n_{\omega,j} |\psi_{\omega,j}\rangle\langle\psi_{\omega,j}|, \quad 0 \le n_{\omega,j} \le 1, \quad \langle\psi_{\omega,j}|\psi_{\omega,j'}\rangle = \delta_{jj'}$$

where the $n_{\omega,j}$'s are the natural occupation numbers (NON) and $\psi_{\omega,j}$'s the natural orbitals (NO) for the chosen continuous model of the neutral system with atomic configuration $\omega$, it holds

$$j_A(\chi, \omega) = -\sum_j n_{\omega,j} \|\Pi_{\chi\omega}^A \psi_{\omega,j}\|_{Q(A)}^2.$$

Minimizing $j_A(\chi, \omega)$ thus amounts to maximizing the NON-weighted sum of the $Q(A)$-norms of $Q(A)$-orthogonal projections of the NON on the discretization space $\chi_\omega$. Other criteria may include errors on molecular properties, or a weighted sum of several elementary criteria, each of them targeting a specific property. The criterion should be chosen according to the intended application.

The aggregated criterion to be optimized then reads as an integral over the configuration space $\Omega$ with respect to the probability measure $\mathbb{P}$

$$J(\chi) := \int_{\Omega} j(\chi, \omega) d\mathbb{P}(\omega), \qquad (4)$$

and the problem of basis set optimization can be formulated as

$$\text{find } \chi_0 \in \underset{\chi \in \mathcal{B}}{\operatorname{argmin}} J(\chi)$$

In the following, $J_E$ and $J_A$ denote the evaluation of the criterion (4) with $j = j_E$ and $j = j_A$ respectively.
**Remark 1** (Reference solutions). The evaluation of criteria $J_E$ and $J_A$ hinges on the knowledge of exact GS energy $E_\omega$ or 1-RDM $\gamma_\omega$ for $\omega$ in the support of $\mathbb{P}$. In practice, these data can be approximated by very accurate off-line reference computations for a small, wisely chosen, sample of configurations $\omega$. This is the reason why the probability measure $\mathbb{P}$ can only be a finite weighted sum of Dirac distributions, as defined in (1).

### 3. APPLICATION TO 1D TOY MODEL

In this section, we focus on a linear one-dimensional toy model, mimicking a homonuclear diatomic molecule.

#### 3.1. Description of the model

Let us consider a system of two 1D point-like “nuclei” and two 1D spinless non-interacting quantum “electrons”. The one-particle state space is then $\mathcal{H} = L^2(\mathbb{R})$ and the configuration space $\Omega = (0, +\infty)$. In this section, a configuration of $\Omega$ will be labelled by the positive real number $a > 0$ such that the nuclei are located at $-a$ and $a$. The one-particle Hamiltonian at configuration $a$ then is

$$H_a = -\frac{1}{2} \frac{d^2}{dx^2} + V_a, \qquad (5)$$

where $V_a$ models the nuclei-electron interaction. We choose $V_a$ to be a double-well potential with minima at $-a$ and $+a$, defined by

$$\forall x \in \mathbb{R}, \quad V_a(x) = \frac{1}{8a^2 + 4}(x - a)^2(x + a)^2. \qquad (6)$$

Several considerations led us to define the potential as such. First, $V_a$ is designed so that i) each $H_a$ admits a non-degenerate ground-state of energy $E_a$, and ii) the function $a \mapsto E_a$ has the shape of the ground-state dissociation curve of a homonuclear diatomic molecule with atoms at $-a$ and $+a$. Since the two “electrons” do not interact, the ground-state energy $E_a$ and density matrices $\gamma_a \in \mathcal{G}_2$ are given by

$$E_a = \mathrm{Tr}(H_a \gamma_a) = \min_{\gamma \in \mathcal{G}_2} \mathrm{Tr}(H_a \gamma), \qquad (7)$$

where

$$\mathcal{G}_2 := \{ \gamma \in \mathcal{L}(L^2(\mathbb{R})), \gamma^2 = \gamma = \gamma^*, \mathrm{Tr}(\gamma) = 2 \},$$

$\mathcal{L}(L^2(\mathbb{R}))$ denoting the space of bounded linear operators on $L^2(\mathbb{R})$. The existence and uniqueness of the solution to problem (7) can be shown by elementary arguments of functional analysis and spectral theory that we do not detail here. Second, $V_0(x) = \frac{1}{4}x^4$ so that (5) corresponds to the quartic oscillator, for which we have reference numerical solutions (e.g. [6]). Third, $V_a$ behaves like $x^2/2$ around $\pm a$ for large values of $a$ and $V_a(0) \sim a^4/8 \to +\infty$ when $a \to +\infty$. Therefore, in the limit $a \to +\infty$, problem (7) is similar to two decoupled quantum harmonic oscillators centered in $-a$ and $+a$ whose bound states are all explicitly known. For the sake of illustration, we display in Figure 1 the potential $V_a$ for two different values of $a$.

In practice, it is convenient to compute $\gamma_a$ and $E_a$ from the lowest two eigenvalues $\lambda_{a,1} < \lambda_{a,2}$ of $H_a$ and an associated pair $(\varphi_{a,1}, \varphi_{a,2})$ of orthonormal eigenvectors

$$\left\{ \begin{array}{ll} H_a \varphi_{a,i} = \lambda_{a,i} \varphi_{a,i}, & i = 1, 2 \\ \langle \varphi_{a,i} | \varphi_{a,j} \rangle = \delta_{ij}, & i, j = 1, 2, \end{array} \right. \qquad (8)$$

$\langle \cdot | \cdot \rangle$ denoting the $L^2$ inner product. We indeed have

$$E_a = \lambda_{a,1} + \lambda_{a,2} \quad \text{and} \quad \gamma_a = |\varphi_{a,1}\rangle\langle\varphi_{a,1}| + |\varphi_{a,2}\rangle\langle\varphi_{a,2}|. \qquad (9)$$
<figure><img src="image_2.png" /><figcaption>FIGURE 1. $x \mapsto V_a(x)$ for $a = 1$ and $a = 5$.</figcaption></figure>

The evaluation of criteria $J_A$ and $J_E$ requires the computation of reference ground-state density matrices or energies, which amounts to find very accurate solutions of (8) for the configurations $a_k$ in the support of the chosen atomic probability measure

$$ \mathbb{P} = \sum_{n=1}^{N_c} \beta_n \delta_{a_n}, \quad 0 < a_1 < a_2 < \dots < a_{N_c}, \quad \beta_n > 0, \quad \sum_{n=1}^{N_c} \beta_n = 1. \qquad (10) $$

We chose to compute these reference data using a 3-point finite-difference (FD) scheme on a large enough interval $[-x_{\max}, x_{\max}]$ discretized into a uniform grid with $N_g$ grid points:

$$ x_j = -x_{\max} + j\delta x, \quad 1 \le j \le N_g, \quad \delta x = \frac{2x_{\max}}{N_g + 1}. $$

We then impose homogeneous Dirichlet boundary conditions at $-x_{\max}$ and $x_{\max}$. The parameter $x_{\max}$ is chosen such that $x_{\max} = a_{\max} + r_{\max}$, where $a_{\max} = \max(\text{supp}(\mathbb{P}))$ and $r_{\max} > 0$ is the radius beyond which atomic densities are zero at machine (double) precision. Note that this numerical scheme is independent of the configuration $a$. The FD discretization of problem (13) gives rise to the eigenvalue problem

$$ \left\{ \begin{array}{l} H_a^{\text{FD}} \varphi_{a,i}^{\text{FD}} = \lambda_{a,i}^{\text{FD}} \varphi_{a,i}^{\text{FD}} \quad i = 1, 2 \\ \delta x (\varphi_{a,i}^{\text{FD}})^T \varphi_{a,j}^{\text{FD}} = \delta_{ij}, \end{array} \right. \qquad (11) $$

where $H_a^{\text{FD}} \in \mathbb{R}_{\text{sym}}^{N_g \times N_g}$ is a real symmetric matrix of size $N_g \times N_g$, and the reference data are obtained as

$$ E_a^{\text{FD}} = \lambda_{a,1}^{\text{FD}} + \lambda_{a,2}^{\text{FD}} \quad \text{and} \quad P_a^{\text{FD}} = \varphi_{a,1}^{\text{FD}} (\varphi_{a,1}^{\text{FD}})^T + \varphi_{a,2}^{\text{FD}} (\varphi_{a,2}^{\text{FD}})^T \in \mathbb{R}_{\text{sym}}^{N_g \times N_g}, \qquad (12) $$

where $P_a^{\text{FD}}$ can be interpreted as an approximation of the matrix $\gamma_a(x_j, x_{j'})$ containing the values of the (integral kernel of the) density matrix $\gamma_a$ at the grid points.

## 3.2. Variational approximation in AO basis sets

For any given configuration $a \in \mathbb{R}_+$ and basis $\chi = \{\chi_\mu\}_{1 \le \mu \le N_b} \in \mathcal{B}$, problem (8) is solved using a Galerkin method with the basis $\chi_a = \{\chi_{a,\mu}\}_{1 \le \mu \le 2N_b}$ composed of two copies of the basis $\chi$, the first one translated to $a$,
and the second one to $-a$:

$$ \chi_{a,1} = \chi_1(\cdot - a), \dots, \chi_{a,N_b} = \chi_{N_b}(\cdot - a), \chi_{a,N_b+1} = \chi_1(\cdot + a), \dots, \chi_{a,2N_b} = \chi_{N_b}(\cdot + a). $$

Defining the Hamiltonian matrix

$$ H_a^\chi = \left( \left\langle \chi_{a,\mu} \left| \left( -\frac{1}{2} \frac{d^2}{dx^2} + V_a \right) \right| \chi_{a,\nu} \right\rangle \right)_{1 \le \mu, \nu \le 2N_b} $$

and the overlap matrix

$$ S_a^\chi = (\langle \chi_{a,\mu} | \chi_{a,\nu} \rangle)_{1 \le \mu, \nu \le 2N_b}, $$

the discretization of problem (8) in the AO basis set $\chi$ then reads as the generalized eigenvalue problem: find $(C_{a,i}^\chi, \lambda_{a,i}^\chi) \in \mathbb{R}^{2N_b} \times \mathbb{R}, i = 1, 2$ such that

$$ \left\{ \begin{array}{ll} H_a^\chi C_{a,i}^\chi = \lambda_{a,i}^\chi S_a^\chi C_{a,i}^\chi & i = 1, 2 \\ (C_{a,i}^\chi)^T S_a^\chi C_{a,j}^\chi = \delta_{ij}. & \end{array} \right. \qquad (13) $$

The approximation $\varphi_{a,i}^\chi$ of $\varphi_{a,i}$ in the AO basis set $\chi$ can then be recovered as the linear combination of atomic orbitals (LCAO)

$$ \forall x \in \mathbb{R}, \quad \varphi_{a,i}^\chi(x) = \sum_{\mu=1}^{2N_b} [C_{a,i}^\chi]_\mu \chi_{a,\mu}(x). \qquad (14) $$

One way to compare the LCAO ground-state 1-RDM to the reference FD solution $P_a^{\text{FD}}$ is to simply evaluate the former at the grid points $x_j$, which gives rise to the matrix $P_a^\chi \in \mathbb{R}_{\text{sym}}^{N_g \times N_g}$ with entries

$$ [P_a^\chi]_{jj'} = \sum_{i=1}^{2} \varphi_{a,i}^\chi(x_j) \varphi_{a,i}^\chi(x_{j'}). $$

Due to numerical errors, the matrix $P_a^\chi$ is however not a rank-2 orthogonal projector. We therefore chose to follow a slightly different route (leading to very similar results). The finite difference grid gives a reference discrete setting in which any quantity of interest for any configuration and AO basis set can be expressed. For all $a$'s, the basis $\chi_a$ is represented by a matrix $X_a \in \mathbb{R}^{N_g \times 2N_b}$. For any vectors $Y_1, Y_2 \in \mathbb{R}^{N_g}$, the discrete $A$ inner product simply reads $\delta x Y_1^T A Y_2$ where the notation $A$ stands for both the continuous inner product and its finite-difference discretization matrix. We denote by $\| \cdot \|_A$ the associated norm on $\mathbb{R}^{N_g}$. Solutions to (13) are then obtained by approximating respectively the Hamiltonian and overlap matrix by

$$ H_a^\chi \approx H_a^X := (\delta x X_{a,\mu}^T H_a^{\text{FD}} X_{a,\nu})_{1 \le \mu, \nu \le 2N_b}, \quad S_a^\chi \approx S_a^X := (\delta x X_{a,\mu}^T X_{a,\nu})_{1 \le \mu, \nu \le 2N_b}, $$

and finding $(C_{a,i}^X, \lambda_{a,i}^X) \in \mathbb{R}^{2N_b} \times \mathbb{R}, i = 1, 2$, such that

$$ \left\{ \begin{array}{ll} H_a^X C_{a,i}^X = \lambda_{a,i}^X S_a^X C_{a,i}^X, & i = 1, 2 \\ (C_{a,i}^X)^T S_a^X C_{a,j}^X = \delta_{ij}, & i, j = 1, 2, \end{array} \right. \qquad (15) $$

from which we get the discrete approximations

$$ \varphi_{a,i}^X = X_a C_{a,i}^X, \quad i = 1, 2. \qquad (16) $$
Let us gather the coefficients $C_a^X$ into the $2N_b \times 2$ matrix $C_a^X = (C_{a,1}^X | C_{a,2}^X)$. The ground-state density matrix in the basis $\chi_a$ is approximated by

$$P_a^X = \varphi_{a,1}^X(\varphi_{a,1}^X)^T + \varphi_{a,2}^X(\varphi_{a,2}^X)^T = (X_a C_a^X) (X_a C_a^X)^T \in \mathbb{R}_{\text{sym}}^{N_g \times N_g}. \quad (17)$$

### 3.3. Overcompleteness of Hermite Basis Sets

Before getting into basis set optimization, we introduce the following standard Hermite Basis Set (HBS), constructed from eigenfunctions of the quantum harmonic oscillator. Those functions are solutions to the eigenvalue problem $(-\frac{1}{2}\frac{d^2}{dx^2} + \frac{1}{2}x^2)h_n = \epsilon_n h_n$ and are explicitly given by

$$h_n(x) = c_n p_n(x) \exp\left(-\frac{x^2}{2}\right), \quad \epsilon_n = n + \frac{1}{2}, \quad n \in \mathbb{N}, \quad (18)$$

where $p_n$ is the Hermite polynomial of degree $n$ (with the same parity as $n$) and $c_n$ a normalization constant such that $(h_n)_{n \in \mathbb{N}}$ forms an orthonormal basis of $L^2(\mathbb{R})$. The $h_n$'s are the analogues of the standard atomic orbitals obtained by solving atomic electronic structure problems. Let us first consider the AO basis set made of the first $N_b$ Hermite functions

$$\chi_{\mu}^{\text{HBS}} = \{\chi_{\mu}^{\text{HBS}}\}_{1 \le \mu \le N_b} = \{h_n\}_{0 \le n \le N_b - 1}.$$

The overlap matrix for the configuration $a$ then is of the form

$$S_a^{\chi_{\mu}^{\text{HBS}}} = \begin{pmatrix} I_{N_b} & \Sigma_a \\ \Sigma_a^T & I_{N_b} \end{pmatrix} \quad \text{where} \quad \Sigma_a := (\langle h_n(\cdot - a) | h_m(\cdot + a) \rangle)_{0 \le n, m \le N_b - 1}.$$

The matrix $\Sigma_a$ corresponds to the overlap of functions that are localized at different atomic positions. It satisfies $\Sigma_a \simeq 0$ when $a$ is large and $\Sigma_a \simeq I_{N_b}$ when $a$ is close to 0, therefore causing conditioning issues on the overlap matrix $S_a^{\chi_{\mu}^{\text{HBS}}}$, a phenomenon known as overcompleteness: when $a$ is too small, the basis functions centered at $\pm a$ are almost equal, hence almost linearly dependent in the basis set. We illustrate this problem by plotting the condition number of the overlap matrix $S_a^{\chi_{\mu}^{\text{HBS}}}$ for different values of $a$ in Figure 2, which indeed blows up for small values of $a$. This is a well-known issue, and several methods have been proposed in the literature to cure this phenomenon, such as the standard canonical orthonormalization procedure [19] or more recent work based on a Cholesky decomposition of the overlap matrix [17]. Such methods are however not directly related to the optimization procedure presented in this paper.

### 3.4. Practical computation of the criterion $J_A$ and $J_E$

The rest of this section is dedicated to the rewriting and the computation of criteria $J_A$ and $J_E$ for our 1D model in the discrete setting.

#### 3.4.1. Reference orthonormal basis

In order to avoid potential numerical stability issues, each of the $N_b$ atomic orbital $\chi_\mu$ is decomposed on a given truncated orthonormal basis of $L^2(\mathbb{R})$ of size $N$ such that $N_b \ll N \ll N_g$. We choose here the orthonormal basis introduced in (18). Hence, the matrix $X_a \in \mathbb{R}^{N_g \times 2N_b}$ is written as

$$X_a = B_a I_R, \quad (19)$$

with

$$B_a = (h_0(x. - a) | \cdots | h_{N-1}(x. - a) | h_0(x. + a) | \cdots | h_{N-1}(x. + a)) \in \mathbb{R}^{N_g \times 2N},$$
<figure><img src="image_3.png" /><figcaption>FIGURE 2. Condition number of the HBS overlap matrix $S_a^{\text{HBS}}$ for different values of a in log-log scale. The larger the basis set, the faster the condition number blows up for small values of a.</figcaption></figure>

and

$$I_R = \begin{pmatrix} R & 0 \\ 0 & R \end{pmatrix} \in \mathbb{R}^{(2N) \times (2N_b)}, \qquad (20)$$

where $R \in \mathbb{R}^{N \times N_b}$ gathers the coefficients of the atomic orbitals $\chi_\mu$ in the truncated HBS orthonormal basis. Note that we have duplicated $R$ in $I_R$ as we consider the same basis at each position $\pm a$, but everything that follows can be easily adapted to the case where we would like to optimize the bases at each position separately (to deal with heteronuclear molecular systems for instance). We moreover impose that $R^T R = I_{N_b}$, so that the overlap matrix of $X_a$, denoted by $S(X_a)$, has the same form as in Section 3.3, that is

$$S(X_a) := \delta x X_a^T X_a = \begin{pmatrix} I_{N_b} & \Sigma_a \\ \Sigma_a^T & I_{N_b} \end{pmatrix}, \qquad (21)$$

where $\Sigma_a$ is the overlap between functions localized at $+a$ and functions localized at $-a$. To avoid any issues arising from the conditioning of $S(X_a)$, the minimal sampled distance $a_{\min}$ should not be taken too small.

In the following, we detail the computation of each of the two criteria using the matrix $R$ as the main variable. We will subsequently optimize the criteria $J_A$ and $J_E$ with respect to $R$ to obtain optimal AO basis sets. In order to ease the reading of the following computations, every vector of $\mathbb{R}^{N_g}$ is rescaled by a factor $\sqrt{\delta x}$ so that for any given $Y_1, Y_2 \in \mathbb{R}^{N_g}$ the discrete A inner product simply reads $Y_1^T A Y_2$. The same holds for overlap matrices: with this convention, $S(X_a) = X_a^T X_a$. The output of the optimization is then scaled back to its former state by a factor $1/\sqrt{\delta x}$ to recover the original normalization.

### 3.4.2. Criterion $J_A$

Let $a \in \mathbb{R}_+$ be fixed and let $S^A(Y) = Y^T A Y$ denote the overlap matrix for the A-inner product of any rectangular matrix $Y \in \mathbb{R}^{N_g \times d}$. Since the columns of $X_a[S^A(X_a)]^{-\frac{1}{2}}$ are orthonormal for the A inner product, that is

$$(X_a[S^A(X_a)]^{-\frac{1}{2}})^T A (X_a[S^A(X_a)]^{-\frac{1}{2}}) = I,$$
the projection $\Pi_{X_a}^A$ takes the simple form

$$ \Pi_{X_a}^A = (X_a[S^A(X_a)]^{-\frac{1}{2}})(X_a[S^A(X_a)]^{-\frac{1}{2}})^T A = X_a[S^A(X_a)]^{-1}X_a^T A. \quad (22) $$

Hence, using the cyclicity of the trace and definitions (3), (12) and (22), one has

$$ \begin{aligned} j_A(\chi, a) &\approx -\mathrm{Tr}(P_a^{\mathrm{FD}} \Pi_{X_a}^A A \Pi_{X_a}^A) \\ &= -\mathrm{Tr}(P_a^{\mathrm{FD}} \times (AB_a I_R)[S^A(B_a I_R)]^{-1}(AB_a I_R)^T) \\ &= -\mathrm{Tr}(M_A^{\mathrm{offline}}(a) I_R [S^A(B_a I_R)]^{-1} I_R^T), \end{aligned} $$

where we have collected in the last expression all matrices independent of $R$ into the matrix

$$ M_A^{\mathrm{offline}}(a) = (AB_a)^T P_a^{\mathrm{FD}} AB_a \in \mathbb{R}^{2N \times 2N}. \quad (23) $$

Then, using the probability measure $\mathbb{P}$ in (10), we get

$$ J_A(R) = - \int_{\Omega} \mathrm{Tr} (M_A^{\mathrm{offline}}(a) I_R [S^A(B_a I_R)]^{-1} I_R^T) \, d\mathbb{P}(a) = - \sum_{n=1}^{N_c} \beta_n \mathrm{Tr} (M_A^{\mathrm{offline}}(a_n) I_R [S^A(B_{a_n} I_R)]^{-1} I_R^T) $$

and the optimization problem finally writes, with unknown $R \in \mathbb{R}^{N \times N_b}$ and for a given inner product $A$

$$ \text{Find } R_{\text{opt}} \in \underset{R \in \mathbb{R}^{N \times N_b}, R^T R = I_{N_b}}{\operatorname{argmin}} J_A(R) \quad (24) $$

### 3.4.3. Criterion $J_E$

Let again $a \in \mathbb{R}_+$ be fixed. We denote by

$$ G(N_g, 2) := \{ P \in \mathbb{R}^{N_g \times N_g} | P^2 = P = P^T, \mathrm{Tr}(P) = 2 \} $$

the discrete counterpart of the Grassmann manifold $\mathcal{G}_2$, and write $E_a^R$ (resp. $H_a^R$) instead of $E_a^\chi$ (resp. $H_a^\chi$), so that the dependence in the matrix $R$ appears explicitly. Equation (7) reads in the discrete setting

$$ \begin{aligned} E_a^R &= \min_{P \in G(N_g, 2)} \mathrm{Tr}(P H_a^R) = \min_{\substack{C \in \mathbb{R}^{2N_b \times 2} \\ (C)^T S(B_a I_R) C = I_2}} \mathrm{Tr}(C C^T \times (B_a I_R)^T H_a^{\mathrm{FD}}(B_a I_R)) \\ &= \mathrm{Tr}(C_a^R (C_a^R)^T \times I_R^T M_E^{\mathrm{offline}}(a) I_R) \end{aligned} \quad (25) $$

where, as for the previous case, all matrices independent of $R$ have been gathered in the matrix

$$ M_E^{\mathrm{offline}}(a) = B_a^T H_a^{\mathrm{FD}} B_a, \quad (26) $$

and the matrix $C_a^R$ is solution to the minimization problem

$$ \min_{\substack{C^R \in \mathbb{R}^{2N_b \times 2} \\ (C^R)^T S(B_a I_R) C^R = I_2}} \mathrm{Tr} (C^R (C^R)^T \times I_R^T M_E^{\mathrm{offline}}(a) I_R) \quad (27) $$

and is given in practice by $C_a^R = [S(B_a I_R)]^{-\frac{1}{2}} (u_{a,1} | u_{a,2})$ where $u_{a,1}$ and $u_{a,2}$ are orthonormal eigenvectors associated to the lowest two eigenvalues of

$$ [S(B_a I_R)]^{-\frac{1}{2}} I_R M_E^{\mathrm{offline}}(a) I_R^T [S(B_a I_R)]^{-\frac{1}{2}}. $$
From (10) and (25), one can compute

$$J_E(R) = \int_{\Omega} |E_a^{\text{FD}} - E_a^R|^2 d\mathbb{P}(a) = \sum_{n=1}^{N_c} \beta_n |E_{a_n}^{\text{FD}} - E_{a_n}^R|^2$$

and the optimization problem reads

$$\text{Find } R_{\text{opt}} \in \underset{R \in \mathbb{R}^{N \times N_b}, R^T R = I_{N_b}}{\operatorname{argmin}} J_E(R) \qquad (28)$$

# 4. NUMERICAL RESULTS

## 4.1. Numerical setting and first results

Problems (24) and (28) are solved by direct minimization algorithms over the Stiefel manifold [1]

$$\text{St}(\mathcal{N}, N_b) = \{R \in \mathbb{R}^{\mathcal{N} \times N_b} | R^T R = I_{N_b}\}.$$

The explicit computation of the gradients of $J_A$ and $J_E$ with respect to $R$ is detailed in the Appendix. We used a L-BFGS algorithm (with tolerance $10^{-7}$ on the norm of the projected gradient), as implemented in the Optim.jl package [22] in the Julia language [4]. As initial guess, we picked the first $N_b$ Hermite functions introduced in Section 3.3.

In this subsection, we choose a probability distribution $P$ supported in the interval $\mathcal{I} = [1.5, 5]$ so as to retain the physics of interest that takes place around the equilibrium configuration $a_0 \simeq 1.925$ and all the way to dissociation. In particular $a_{\min} = 1.5$ is taken sufficiently large to avoid the conditioning issues on the overlap matrices described in Section 3.3. More precisely, all the results in this subsection are obtained with the probability

$$P = \frac{1}{10} \sum_{n=1}^{10} \delta_{a_n} \quad \text{with } a_n = 1.5 + (n-1) \frac{3.5}{9}. \qquad (29)$$

The quantities $M_A^{\text{offline}}(a_n)$ and $M_E^{\text{offline}}(a_n)$ are computed offline beforehand. We will discuss this choice and consider other probability measures $P$ in Sections 4.2.2 and 4.2.3.

The finite-difference grid is a uniform grid on the interval $[-20, 20]$ discretized into $N_g = 1999$ points ($\delta x = 0.02$). Finally, we decompose the basis functions to be optimized in the HBS $\{h_n\}_{0 \le n \le N-1}$ of $L^2(\mathbb{R})$ of size $\mathcal{N} = 10$. Regarding the choice of the inner product for the first criterion $J_A$, we used the standard $L^2(\mathbb{R})$ and the $H^1(\mathbb{R})$ inner products, and denoted by $J_{L^2}$ and $J_{H^1}$ the corresponding. This translates at the discrete level by choosing $A = I_{N_g}$ for $J_{L^2}$ and $A = I_{N_g} - \Delta$ for $J_{H^1}$ where $\Delta$ is the 3-point finite-difference discretization matrix of the 1D Laplace operator. Once obtained, the optimal bases are used to solve the variational problem (15) on a much finer sampling of $\mathcal{I}$ and their accuracy is compared to the HBS. The code performing the simulations and plotting the results is available online<fnref n="1" />. Also, for the sake of clarity in the plots, $\tilde{E}_a$ (resp. $\tilde{\rho}_a$) denotes the GS energy (resp. the density) in the configuration $a$ with a given basis (specified by the context) and $E_a$ (resp. $\rho_a$) stands for the reference energy (resp. density) on the finite difference grid. Note that we write HBS for the (nonoptimized) Hermite basis set, and $L^2$-OBS, $H^1$-OBS or $E$-OBS for optimized basis sets with respect to the criterion $J_{L^2}$, $J_{H^1}$, or $J_E$.

Figure 3 displays the dissociation curve and the energy difference on the interval $\mathcal{I}$ for different values of $N_b$, the size of the AO basis set. For $N_b = 1$, i.e. only one basis function at $\pm a$, criterion $J_E$ shows better performance than the criterion $J_A$, regardless of the choice of norm to perform the projections. It also very closely matches the accuracy of the standard HBS. When $N_b$ becomes larger however, the different criteria

<footnote n="1">¹https://github.com/gkemlin/1D_basis_optimization</footnote>
behave in a similar fashion and we observe that they approach the dissociation curve better than the Hermite basis. Comparing the values of criterion $J_E$ for all bases, which directly measures the distance to the dissociation curve, we see in Table 1 that all optimized bases give an increased accuracy of roughly four orders of magnitude over the interval $\mathcal{I}$ for $N_b = 4$.

In Figure 4, we plot the density for a given value of $a$ and the error on the density for different norms, with varying values of $N_b$. The error is plotted with respect to three different distances: the $L^1$-norm, which corresponds to the $L^2$-norm on eigenvectors, the $H^1$-norm of the error on the density, as it is common to compute the forces $\int_{\mathbb{R}} \rho \nabla_a V_a$ with good estimates on the $H^{-1}$-norm of $\nabla_a V_a$ (see e.g. [8]), and the distance

$$ \|\nabla\sqrt{\rho_1} - \nabla\sqrt{\rho_2}\|_{L^2} $$

(recall that the von-Weizsäcker kinetic energy reads $\frac{1}{2}\int_{\mathbb{R}} |\nabla\sqrt{\rho}|^2$). We observe similar behaviors between these different distances. For $N_b = 1$, both bases obtained with the first criterion behave slightly better than the standard Hermite basis and the basis computed with the second criterion. For $N_b = 3$, we observe again that all optimal bases yield better accuracy than the Hermite basis. Table 1 gives the confirmation that each basis for a given criterion indeed performs better than the other bases for that particular criterion. As for dissociation curves, we read from the values of $J_{L^2}$ and $J_{H^1}$ that the optimized bases yield similar results for large $N_b$, all of them giving lower values than the HBS. Note that the optimal bases for criterion $J_{L^2}$ and $J_{H^1}$ give similar results for any number of basis functions $N_b$, so that the $L^2$ and $H^1$ norm optimizations seem equivalent.

In terms of computational time, first note that criterion $J_{H^1}$ is always more expensive to compute than $J_{L^2}$ as it requires additional matrix-vector products with the matrix $A$, this having noticeable impact on the computational time. Second, criterion $J_E$ requires less off-line data as it only needs to be given the reference eigenvalues while criterion $J_A$ requires the reference GS eigenvectors (or density matrices). In addition, the use of orthonormality constraints as detailed in appendix allows one to compute the gradient of $J_E$ at very low cost. In turn, criterion $J_E$ is more than twice faster to minimize than criterion $J_{L^2}$ in our implementation.

Finally, for the sake of completeness, we plot in Figure 5 the different basis functions built with each criterion for different values of $N_b$, confirming again the previous observations that the optimal basis functions are quite close to the standard Hermite basis functions.

The main conclusion of these observations is that, for $N_b$ large enough, there is no real difference between the proposed criteria. Still, if the bases we built do not seem to be very different from the standard Hermite basis (Figure 5), building optimal bases allows to increase accuracy on the quantities of interest we focused on by one order of magnitude in average.
<table><caption>Value of $J_{L^2}$ for the different basis sets</caption><thead><tr><th>Basis</th><th>N<sub>b</sub> = 1</th><th>N<sub>b</sub> = 2</th><th>N<sub>b</sub> = 3</th><th>N<sub>b</sub> = 4</th></tr></thead><tbody><tr><td>HBS</td><td>-7.40829</td><td>-7.70051</td><td>-7.74312</td><td>-7.77138</td></tr><tr><td>L<sup>2</sup>-OBS</td><td>-7.43954</td><td>-7.76479</td><td>-7.77725</td><td>-7.77773</td></tr><tr><td>H<sup>1</sup>-OBS</td><td>-7.43928</td><td>-7.76466</td><td>-7.77724</td><td>-7.77772</td></tr><tr><td>E-OBS</td><td>-7.39410</td><td>-7.76425</td><td>-7.77720</td><td>-7.77772</td></tr></tbody></table>

<table><caption>Value of $J_{H^1}$ for the different basis sets</caption><thead><tr><th>Basis</th><th>N<sub>b</sub> = 1</th><th>N<sub>b</sub> = 2</th><th>N<sub>b</sub> = 3</th><th>N<sub>b</sub> = 4</th></tr></thead><tbody><tr><td>HBS</td><td>-10.5613</td><td>-11.0566</td><td>-11.1451</td><td>-11.2402</td></tr><tr><td>L<sup>2</sup>-OBS</td><td>-10.6256</td><td>-11.2338</td><td>-11.2630</td><td>-11.2650</td></tr><tr><td>H<sup>1</sup>-OBS</td><td>-10.6265</td><td>-11.2342</td><td>-11.2630</td><td>-11.2651</td></tr><tr><td>E-OBS</td><td>-10.5334</td><td>-11.2313</td><td>-11.2626</td><td>-11.2650</td></tr></tbody></table>

<table><caption>Value of $J_E$ for the different basis sets</caption><thead><tr><th>Basis</th><th>N<sub>b</sub> = 1</th><th>N<sub>b</sub> = 2</th><th>N<sub>b</sub> = 3</th><th>N<sub>b</sub> = 4</th></tr></thead><tbody><tr><td>HBS</td><td>3.77956×10<sup>-2</sup></td><td>3.98301×10<sup>-3</sup></td><td>1.86537×10<sup>-3</sup></td><td>1.35309×10<sup>-4</sup></td></tr><tr><td>L<sup>2</sup>-OBS</td><td>6.52016×10<sup>-2</sup></td><td>2.18282×10<sup>-4</sup></td><td>1.01365×10<sup>-6</sup></td><td>3.22260×10<sup>-8</sup></td></tr><tr><td>H<sup>1</sup>-OBS</td><td>6.83537×10<sup>-2</sup></td><td>2.40548×10<sup>-4</sup></td><td>1.27251×10<sup>-6</sup></td><td>3.91885×10<sup>-8</sup></td></tr><tr><td>E-OBS</td><td>3.69610×10<sup>-2</sup></td><td>1.92087×10<sup>-4</sup></td><td>6.93394×10<sup>-7</sup></td><td>2.54014×10<sup>-8</sup></td></tr></tbody></table>

<table><caption>L-BFGS iterations</caption><thead><tr><th>Basis</th><th>N<sub>b</sub> = 1</th><th>N<sub>b</sub> = 2</th><th>N<sub>b</sub> = 3</th><th>N<sub>b</sub> = 4</th></tr></thead><tbody><tr><td>L<sup>2</sup>-OBS</td><td>4</td><td>13</td><td>48</td><td>219</td></tr><tr><td>H<sup>1</sup>-OBS</td><td>7</td><td>17</td><td>235</td><td>not converged after 500 it</td></tr><tr><td>E-OBS</td><td>6</td><td>19</td><td>52</td><td>134</td></tr></tbody></table>

<caption>TABLE 1. (Top & Middle) Values of the different criteria for the HBS and optimal bases, for increasing values of $N_b$. (Bottom) Number of iterations of L-BFGS required for each criterion to achieve convergence up to requested tolerance ($10^{-7}$ on the $l^2$-norm of the gradient).</caption>
<figure><img src="image_4.png" /><figcaption>FIGURE 3. Energies for the optimal bases obtained with the different criteria. (Top) Dissocia-
tion curve. (Bottom) Errors on the energy on the range of configuration $\mathcal{I} = [1.5, 5]$.</figcaption></figure>
<figure><img src="image_5.png" /><figcaption>FIGURE 4. (Top) Densities for the optimal bases obtained with the different criteria. (Middle) Errors on the density for different norms with $N_b = 1$. (Bottom) Error on the density for different norms with $N_b = 3$.</figcaption></figure>
<figure><img src="image_6.png" /><figcaption>$N_b = 2$</figcaption></figure>

<figure><img src="image_7.png" /><figcaption>$N_b = 3$</figcaption></figure>

<figure><img src="image_8.png" /><figcaption>$N_b = 4$</figcaption></figure>

<figure><img src="image_9.png" /><figcaption>FIGURE 5. Optimal basis functions for different criteria, each of them being optimized for different values of $N_b$.</figcaption></figure>
## 4.2. Influence of numerical parameters

### 4.2.1. Random starting points

In Section 4.1, we used the first $N_b$ Hermite functions as a starting point for the optimization procedures. We obtain the same solutions if we start from a random matrix $R$ on the Stiefel manifold, in the sense that the optimal values reached for each criterion are the same, as well as the error plots. However, the L-BFGS algorithm requires more iterations to converge. The basis functions obtained from the optimization algorithms are different from those observed in Figure 5, but still span the same space as the variational solutions are equal.

### 4.2.2. Extrapolating the parameter space $\mathcal{I}$

In Section 4.1, we chose a probability measure $\mathbb{P}$ supported in the interval $[1.5, 5]$ in order to avoid conditioning issues. Indeed, taking smaller values of $a$ results in the L-BFGS algorithm having convergence problems when $N_b$ increases. This phenomenon was observed already for $N_b = 3$ or $N_b = 4$ when including $a = 1$ in the support of $\mathbb{P}$. In practice, this problem can be solved by using preconditioning or getting rid of overcompleteness by pre-processing the basis $\chi_a$ (e.g. selecting a smaller basis by filtering out the very small singular values of the original overlap matrix), but for brevity we will not elaborate further in this direction.

However, once we have computed optimal bases for a reasonable interval $\mathcal{I}$, it is possible to use these bases to solve the variational problem (13) and extrapolate the energy and the density to smaller values of $a$ that are not in the set $\mathcal{I}$. The results are plotted in Figure 6. We notice that the quantities of interest are better approximated on $\mathcal{I} = [1.5, 5]$, but for smaller $a$'s, there is no more gain in accuracy with respect to the standard HBS.

<figure><img src="image_10.png" /><figcaption>FIGURE 6. Energy and densities error with extrapolation up to $a = 0.5$, with basis functions optimized on $\mathcal{I} = [1.5, 5]$</figcaption></figure>

### 4.2.3. Choice of the probability $\mathbb{P}$

The major drawback of our AO basis optimization lies in the necessity to compute very accurate reference solutions for all configurations in the support of $\mathbb{P}$. This is not an issue for our 1D toy model but it can be very time consuming for real systems if the support of $\mathbb{P}$ is too large. It is therefore crucial to reduce as much as possible the support of $\mathbb{P}$.

In this section, we study the influence of the probability measure $\mathbb{P}$ on the quality of the optimized bases. For simplicity, we restrict ourselves to uniform samplings of the interval $\mathcal{I} = [1.5, 5]$. Numerical tests show
that increasing the sample size above the reference sampling with $N_c = 10$ points used in Section 4.1 (see Eq. (29)) brings no significant accuracy improvement. Therefore we chose to investigate in the following the performance of the optimal AO basis sets obtained with very sparse sampling. Figure 7 pictures the error of approximation of the dissociation curve and densities for three samplings: first, the two extreme points of the interval $\mathcal{I} = [1.5, 5]$; second, two points around the equilibrium distance $a_0 \approx 1.925$; third, a single point near the equilibrium distance. All curves are plotted for a fixed number of basis functions $N_b = 3$.

It appears that the latter sampling already provides satisfactory accuracy. The criteria $J_{L^2}$ and $J_{H^1}$ are equal to $-5 \times 10^{-6}$ for optimized basis to be compared with $-1.8 \times 10^{-3}$ for standard HBS. Hence they provide a gain of accuracy in energy of three orders of magnitude over the whole dissociation curve.

<figure><img src="image_11.png" /><figcaption>FIGURE 7. Error plots for probability measures $P$ corresponding to very sparse samplings of the interval $\mathcal{I} = [1.5, 5]: i) the two endpoints of $\mathcal{I}$ ii) two points near the equilibrium distance and iii) one point near the equilibrium distance. (Top line) Error on energy. (Bottom line) Error on density in $L^1$ norm. (Left) OBS for $J_{H^1}$. (Middle) OBS for $J_{L^2}$. (Right) OBS for $J_E$. The “a” in legends are the sampled configurations a.</figcaption></figure>

### 4.2.4. Number of Hilbert basis functions

We now take the same setting as in Section 4.1, except that we set $\mathcal{N} = 5$ instead of $\mathcal{N} = 10$. This provides similar results as those collected in Table 1, see Table 2. However, the values of the criteria $J_A$ and $J_E$ are higher than for $\mathcal{N} = 10$, in particular for $N_b = 4$, where criterion $J_A$ cannot be optimized further than $-10^{-5}$, which makes sense as the space over which the optimization algorithms are performed is smaller. Calculations with $\mathcal{N} = 15$ were also performed: for $N_b = 1, 2, 3$, the criteria are slightly improved but for $N_b = 4$, convergence issues were noticed, due to ill conditioning of the overlap matrices for $a = 1.5$ as the number $\mathcal{N}$ of functions used to describe the optimal bases is larger.
<table><caption>Value of $J_{L^2}$ for the different basis sets</caption><thead><tr><th>Basis</th><th>$N_b = 1$</th><th>$N_b = 2$</th><th>$N_b = 3$</th><th>$N_b = 4$</th></tr></thead><tbody><tr><td>HBS</td><td>-7.40829</td><td>-7.70051</td><td>-7.74312</td><td>-7.77138</td></tr><tr><td>$L^2$-OBS</td><td>-7.43933</td><td>-7.76304</td><td>-7.77554</td><td>-7.77618</td></tr><tr><td>$H^1$-OBS</td><td>-7.43923</td><td>-7.76258</td><td>-7.77525</td><td>-7.77612</td></tr><tr><td>E-OBS</td><td>-7.39401</td><td>-7.76259</td><td>-7.77545</td><td>-7.77615</td></tr></tbody></table>

<table><caption>Value of $J_{H^1}$ for the different basis sets</caption><thead><tr><th>Basis</th><th>$N_b = 1$</th><th>$N_b = 2$</th><th>$N_b = 3$</th><th>$N_b = 4$</th></tr></thead><tbody><tr><td>HBS</td><td>-10.5613</td><td>-11.0566</td><td>-11.1451</td><td>-11.2402</td></tr><tr><td>$L^2$-OBS</td><td>-10.6237</td><td>-11.2225</td><td>-11.2541</td><td>-11.2577</td></tr><tr><td>$H^1$-OBS</td><td>-10.6240</td><td>-11.2244</td><td>-11.2555</td><td>-11.2581</td></tr><tr><td>E-OBS</td><td>-10.5328</td><td>-11.2234</td><td>-11.2547</td><td>-11.2580</td></tr></tbody></table>

<table><caption>Value of $J_E$ for the different basis sets</caption><thead><tr><th>Basis</th><th>$N_b = 1$</th><th>$N_b = 2$</th><th>$N_b = 3$</th><th>$N_b = 4$</th></tr></thead><tbody><tr><td>HBS</td><td>3.77956×10<sup>-2</sup></td><td>3.98301×10<sup>-3</sup></td><td>1.86537×10<sup>-3</sup></td><td>1.35309×10<sup>-4</sup></td></tr><tr><td>$L^2$-OBS</td><td>6.43832×10<sup>-2</sup></td><td>2.46466×10<sup>-4</sup></td><td>1.58667×10<sup>-5</sup></td><td>1.01128×10<sup>-5</sup></td></tr><tr><td>$H^1$-OBS</td><td>6.13025×10<sup>-2</sup></td><td>2.45930×10<sup>-4</sup></td><td>1.62235×10<sup>-5</sup></td><td>1.00611×10<sup>-5</sup></td></tr><tr><td>E-OBS</td><td>3.69681×10<sup>-2</sup></td><td>1.30365×10<sup>-4</sup></td><td>1.41935×10<sup>-5</sup></td><td>9.74560×10<sup>-6</sup></td></tr></tbody></table>

<caption>TABLE 2. Value of the different criteria for the different local (optimized and Hermite) bases, with $N = 5$ and increasing values of $N_b$.</caption>

ACKNOWLEDGEMENTS

The authors thank Susi Lehtola and Etienne Polack for fruitful discussions. This project has received funding from the European Research Council (ERC) under the European Union's Horizon 2020 research and innovation programme (grant agreement EMC2 No 810367). This work was supported by the French 'Investissements d'Avenir' program, project Agence Nationale de la Recherche (ISITE-BFC) (contract ANR-15-IDEX-0003). The work of the second author was also supported by the Ecole des Ponts-ParisTech.
APPENDIX

In this appendix, we will use extensively the two symmetries of the trace: for any matrices M and N such that MN and NM are defined,

$$ \mathrm{Tr}(MN) = \mathrm{Tr}(NM) \quad \text{and} \quad \mathrm{Tr}(M^T) = \mathrm{Tr}(M). $$

Computation of the gradient of $J_A$

Let $R, H \in \mathbb{R}^{N \times N_b}$ and define $I_H = \begin{pmatrix} H & 0 \\ 0 & H \end{pmatrix}$. One has

$$ J_A(R+H) - J_A(R) = - \int_{\Omega} \mathrm{Tr} (M_A^{\mathrm{offline}}(a) (2I_R[S^A(B_a I_R)]^{-1} I_H^T + I_R [d[S^A]^{-1}(B_a I_R) \cdot (B_a I_H)] I_R^T)) \, d\mathbb{P}(a) \\ + O(\|H\|^2) \qquad (30) $$

Considering that

$$ (M+H)^{-1} - M^{-1} = -M^{-1}HM^{-1} + O(\|H\|^2) \quad \text{and} \quad S^A(BI_{R+H}) - S^A(B_a I_R) = I_H^T S^A(B)I_R + I_R^T S^A(B)I_H + O(\|H\|^2), $$

it follows from the chain rule that

$$ d[S^A]^{-1}(B_a I_R) \cdot (B_a I_H) = -[S^A(B_a I_R)]^{-1} (I_H^T S^A(B_a) I_R + I_R^T S^A(B_a) I_H) [S^A(B_a I_R)]^{-1}. $$

From this computation, we obtain that the integrand in expression (30) writes for all a

$$ \begin{aligned} & 2\mathrm{Tr}(M_A^{\mathrm{offline}}(a) [I_R[S^A(B_a I_R)]^{-1} I_H^T - I_R[S^A(B_a I_R)]^{-1} I_H^T S^A(B_a) I_R [S^A(B_a I_R)]^{-1} I_R^T]) \\ & = 2\mathrm{Tr}(M_A^{\mathrm{offline}}(a) I_R [S^A(B_a I_R)]^{-1} I_H^T - I_H^T S^A(B_a) I_R [S^A(B_a I_R)]^{-1} I_R^T M_A^{\mathrm{offline}}(a) I_R [S^A(B_a I_R)]^{-1}). \end{aligned} \qquad (31) $$

The idea is now to write the expression (31) as the inner product of H with a given matrix of $\mathbb{R}^{N \times N_b}$, which we will identify as the integrand of the gradient of $J_A$. Changing from $I_H$ to H imposes to decompose each matrix by block and to write the trace in (31) as the sum of traces over the diagonal blocks. To this end we introduce the superscripts "+", "-", "-+" and "--" associated with one of the four identically shaped blocks of a generic matrix

$$ M = \begin{pmatrix} M^{++} & M^{+-} \\ M^{-+} & M^{--} \end{pmatrix}. \qquad (32) $$

Expression (31) therefore immediately reads

$$ \begin{aligned} & 2\mathrm{Tr}\left(I_H^T \underbrace{[M_A^{\mathrm{offline}}(a) I_R [S^A(B_a I_R)]^{-1} - S^A(B_a) I_R [S^A(B_a I_R)]^{-1} I_R^T M_A^{\mathrm{offline}}(a) I_R [S^A(B_a I_R)]^{-1}]}_{M_A(a, R)}\right) \\ & = 2\mathrm{Tr}(H^T (M_A(a, R)^{++} + M_A(a, R)^{--})). \end{aligned} \qquad (33) $$

One can verify that $M_A(a, R)^{++} + M_A(a, R)^{--}$ is in $\mathbb{R}^{N \times N_b}$ and we conclude by identification that

$$ \nabla J_A(R) = -2 \int_{\Omega} (M_A(a, R)^{++} + M_A(a, R)^{--}) \, d\mathbb{P}(a). \qquad (34) $$
Computation of the gradient of $J_E$

Let $R, H \in \mathbb{R}^{N \times N_b}$ and define $I_H = \begin{pmatrix} H & 0 \\ 0 & H \end{pmatrix}$. We immediately have that

$$ \nabla J_E(R) = -2 \int_{\Omega} \nabla E_a(R) (E_a - E_a(R)) d\mathbb{P}(a), \qquad (35) $$

where

$$ E_a(R) = \mathrm{Tr}(C_a(R)(C_a(R))^T \times \mathcal{H}_a(R)), \qquad (36) $$

with $C_a(R)$ defined in Section 3.2 and $\mathcal{H}_a(R) := I_R^T M_E^{\mathrm{offline}}(a) I_R$. Therefore, if we define $\mathcal{E}_a(R, C) = \mathrm{Tr}(CC^T \mathcal{H}_a(R))$, then $E_a(R) = \mathcal{E}_a(R, C_a(R))$ and we have, by the chain rule,

$$ \nabla E_a(R) \cdot H = \nabla_R \mathcal{E}_a(R, C_a(R)) \cdot H + \nabla_C \mathcal{E}_a(R, C_a(R)) \cdot (dC_a(R) \cdot H). $$

We now detail the computations of the two gradients of $\mathcal{E}_a$, namely $\nabla_R \mathcal{E}_a$ and $\nabla_C \mathcal{E}_a$. Computation of the first gradient $\nabla_R \mathcal{E}_a$. Using notation (32), we introduce

$$ M_a := M_E^{\mathrm{offline}}(a) \text{ and } \Sigma(H) := I_H^T M_a I_R = \begin{pmatrix} H^T M_a^{++} R & H^T M_a^{+-} R \\ H^T M_a^{-+} R & H^T M_a^{--} R \end{pmatrix} \in \mathbb{R}^{(2N_b) \times (2N_b)}, $$

so that, with $P = CC^T$,

$$ \begin{aligned} \mathrm{Tr}(P[\mathcal{H}_a(R) \cdot H]) &= \mathrm{Tr}(P[\Sigma(H) + \Sigma(H)^T]) = 2\mathrm{Tr}(P\Sigma(H)) \\ &= 2\mathrm{Tr}(H^T (M_a^{++} RP^{++} + M_a^{-+} RP^{+-} + M_a^{+-} RP^{-+} + M_a^{--} RP^{--})) . \end{aligned} $$

In the end,

$$ \nabla_R \mathcal{E}_a(R, C) = 2 (M_a^{++} R (CC^T)^{++} + M_a^{+-} R (CC^T)^{-+} + M_a^{-+} R (CC^T)^{+-} + M_a^{--} R (CC^T)^{--}) \in \mathbb{R}^{N \times N_b}. $$

Computation of the second gradient $\nabla_C \mathcal{E}_a$. The Euler-Lagrange equation of the minimization problem (25) yields that there exist a symmetric matrix $\Lambda_a(R) \in \mathbb{R}^{2 \times 2}$ such that

$$ \nabla_C \mathcal{E}_a(R, C_a(R)) = 2 \mathcal{H}_a(R) = 2S(B_a I_R)C_a(R)\Lambda_a(R), $$

where $\Lambda_a(R)$ is actually a diagonal matrix whose diagonal is composed of the two lowest eigenvalues of $\mathcal{H}_a(R)$. Moreover, if we differentiate the constraint $C_a(R)^T S(B_a I_R) C_a(R) = \mathrm{Id}_2$, we get

$$ C_a(R)^T S(B_a I_R) (dC_a(R) \cdot H) + (dC_a(R) \cdot H)^T S(B_a I_R) C_a(R) = -C_a(R)^T (dS(B_a I_R) \cdot H) C_a(R), $$

so that

$$ \begin{aligned} \nabla_C \mathcal{E}_a(R, C_a(R)) \cdot (dC_a(R) \cdot H) &= 2\mathrm{Tr}((S(B_a I_R)C_a(R)\Lambda_a(R))^T (dC_a(R) \cdot H)) \\ &= -\mathrm{Tr}((dS(B_a I_R) \cdot H) C_a(R)\Lambda_a(R)C_a(R)^T). \end{aligned} $$

Now, let us recall that

$$ dS(B_a I_R) \cdot H = I_H^T S(B_a) I_R + I_R^T S(B_a) I_H. $$

Thus, by denoting $Q_a(R) = C_a(R)\Lambda_a(R)C_a(R)^T$, we get that

$$ \nabla_C \mathcal{E}_a(R, C_a(R)) \cdot (dC_a(R) \cdot H) = -2 \mathrm{Tr} \left( H^T (S(B_a)^{++} R Q_a(R)^{++} + S(B_a)^{+-} R Q_a(R)^{-+} + S(B_a)^{-+} R Q_a(R)^{+-}) \right) $$

which ends the computations of the second gradient.
Final gradient. Compiling the computations of the two previous paragraphs, we obtain

$$ \nabla_R E_a(R) = 2 (M_a^{++} R P_a(R)^{++} + M_a^{+-} R P_a(R)^{-+} + M_a^{-+} R P_a(R)^{+-} + M_a^{--} R P_a(R)^{--}) \\ - 2 (S_a^{++} R Q_a(R)^{++} + S_a^{+-} R Q_a(R)^{-+} + S_a^{-+} R Q_a(R)^{+-} + S_a^{--} R Q_a(R)^{--}) \quad (37) $$

where $P_a(R) = C_a(R)C_a(R)^T$, $M_a = M_A^{\text{offline}}(a)$, $S_a = S(B_a)$ and $Q_a(R) = C_a(R)\Lambda_a(R)C_a(R)^T$, and the gradient of $J_E$ is computed with (35).

## REFERENCES

[1] P.-A. ABSIL, R. MAHONY, AND R. SEPULCHRE, *Optimization algorithms on matrix manifolds*, Princeton University Press, 2009.

[2] J. ALMLÖF AND P. R. TAYLOR, *General contraction of Gaussian basis sets. I. Atomic natural orbitals for first- and second-row atoms*, The Journal of Chemical Physics, 86 (1987), pp. 4070–4077.

[3] M. BACHMAYR, H. CHEN, AND R. SCHNEIDER, *Error estimates for hermite and even-tempered gaussian approximations in quantum chemistry*, Numer. Math., 128 (2014), pp. 137–165.

[4] J. BEZANSON, A. EDELMAN, S. KARPINSKI, AND V. B. SHAH, *Julia: A fresh approach to numerical computing*, SIAM Review, 59 (2017), pp. 65–98.

[5] J. S. BINKLEY, J. A. POPLE, AND W. J. HEHRE, *Self-consistent molecular orbital methods. 21. small split-valence basis sets for first-row elements*, Journal of the American, (1980).

[6] S. BLINDER, *Eigenvalues for a pure quartic oscillator*, arXiv preprint arXiv:1903.07471, (2019).

[7] S. F. BOYS AND A. C. EGERTON, *Electronic wave functions - i. a general method of calculation for the stationary states of any molecular system*, Proc. R. Soc. Lond. A Math. Phys. Sci., 200 (1950), pp. 542–554.

[8] E. CANCÈS, G. DUSSON, G. KEMLIN, AND A. LEVITT, *Practical error bounds for properties in plane-wave electronic structure calculations*, 2021.

[9] C. CANUTO, M. Y. HUSSAINI, A. QUARTERONI, AND T. A. ZANG, *Spectral methods: fundamentals in single domains*, Springer Science & Business Media, 2007.

[10] I. CHARPENTIER, F. DE VUYST, AND Y. MADAY, *A component mode synthesis method of infinite order of accuracy using subdomain overlapping: numerical analysis and experiments*, Publication du laboratoire d’Analyse Numerique, 96002 (1996), pp. 55–65.

[11] I. CHARPENTIER, F. DE VUYST, AND Y. MADAY, *Méthode de synthèse modale avec une décomposition de domaine par recouvrement*, Comptes rendus de l’Académie des sciences. Série 1, Mathématique, 322 (1996), pp. 881–888.

[12] L. E. DAGA, B. CIVALLERI, AND L. MASCHIO, *Gaussian basis sets for crystalline solids: All-purpose basis set libraries vs system-specific optimizations*, Journal of chemical theory and computation, 16 (2020), pp. 2192–2201.

[13] T. H. DUNNING, *Gaussian basis sets for use in correlated molecular calculations. i. the atoms boron through neon and hydrogen*, J. Chem. Phys., 90 (1989), pp. 1007–1023.

[14] W. J. HEHRE, R. F. STEWART, AND J. A. POPLE, *Self-Consistent Molecular-Orbital methods. i. use of gaussian expansions of Slater-Type atomic orbitals*, J. Chem. Phys., 51 (1969), pp. 2657–2664.

[15] T. HELGAKER, P. JORGENSEN, AND J. OLSEN, *Molecular Electronic-Structure Theory*, John Wiley & Sons, Aug. 2014.

[16] W. KUTZELNIGG, *Theory of the expansion of wave functions in a gaussian basis*, Int. J. Quantum Chem., (1994).

[17] S. LEHTOLA, *Communication: Curing basis set overcompleteness with pivoted cholesky decompositions*, arXiv preprint arXiv:1911.10372, (2019).

[18] P.-O. LÖWDIN, *On the nonorthogonality problem*, in Advances in Quantum Chemistry, P.-O. Löwdin, ed., vol. 5, Academic Press, Jan. 1970, pp. 185–199.

[19] ——, *On the nonorthogonality problem*, in Advances in quantum chemistry, vol. 5, Elsevier, 1970, pp. 185–199.

[20] R. McWEENY, *Gaussian Approximations to Wave Functions*, Nature, 166 (1950), pp. 21–22.

[21] Q. MÉRIGOT, F. SANTAMBROGIO, AND C. SARRAZIN, *Non-asymptotic convergence bounds for wasserstein approximation using point clouds*, Advances in Neural Information Processing Systems, 34 (2021), pp. 12810–12821.

[22] P. K. MOGENSEN AND A. N. RISETH, *Optim: A mathematical optimization package for Julia*, Journal of Open Source Software, 3 (2018), p. 615.

[23] J. OLSEN, *An introduction and overview of basis sets for molecular and Solid-State calculations*, in Basis Sets in Computational Chemistry, E. Perlt, ed., Springer International Publishing, Cham, 2021, pp. 1–16.

[24] G. PAGÈS, *Introduction to vector quantization and its applications for numerics*, ESAIM: Proceedings and Surveys, 48 (2015), pp. 29–79.

[25] D. H. PHAM, *Galerkin method using optimized wavelet-Gaussian mixed bases for electronic structure calculations in quantum chemistry*, PhD thesis, Université Grenoble Alpes, June 2017.
[26] D. SÁNCHEZ-PORTAL, E. ARTACHO, AND J. M. SOLER, *Analysis of atomic orbital basis sets from the projection of plane-wave results*, Journal of Physics: Condensed Matter, 8 (1996), p. 3859.

[27] S. SCHOLZ AND H. YSERENTANT, *On the approximation of electronic wavefunctions by anisotropic gauss and Gauss-Hermite functions*, Numer. Math., 136 (2017), pp. 841–874.

[28] R. A. SHAW, *The completeness properties of gaussian-type orbitals in quantum chemistry*, Int. J. Quantum Chem., 120 (2020), p. 93.

[29] J. C. SLATER, *Atomic shielding constants*, Physical Review, (1930).

[30] R. F. STEWART, *Small gaussian expansions of atomic orbitals*, J. Chem. Phys., 50 (1969), pp. 2485–2495.