\section{Introduction}
Electronic structure~\cite{szabo2012modern, arrazola2021differentiable} has been the intersection among multiple disciplines, including quantum physics, quantum chemistry, and material science. The goal of the field is to study the states and dynamics of many-electron systems, or more specifically, the electronic Hamiltonian~\cite{szabo2012modern} determined by a nuclear field. Consider a quantum system consisting of $N$ electrons and $N'$ nuclei, in which the $a'$th nucleus with charge $Z_{a'}$ is located on position ${\bm{R'_{a'}}}$. The information of all the nuclei forms a set that specifies the nuclear field of the system:
\begin{equation}
    \mathbb{M}\equiv\{({\bm{R'_{a'}}}, Z_{a'})\,|\,a'\!=\!1,\dots ,N'\}.
\end{equation}
The electronic Hamiltonian of this system can be separated from its total Hamiltonian $H_{\rm tot}$ when the Born-Oppenheimer approximation is applied:
\begin{equation}
    H_{\rm elec}(\mathbb{M}; N^\uparrow, N^\downarrow) = H_{\rm tot}(\mathbb{M}; N) - H_{\rm nuc}(\mathbb{M})
\end{equation}
where $N^\uparrow$ and $N^\downarrow$ are the number of spin-up electrons and spin-down electrons, respectively, such that $N\!=\!N^\uparrow\!+\!N^\downarrow$, and $H_{\rm nuc}$ is the sum of the nuclear kinetic energy and nuclear-nuclear repulsion energy. Instead of being the direct sum of single-electron Hamiltonians, $H_{\rm elec}$ (in the Hartree atomic units) has the form:
\begin{equation}
    \begin{aligned}
        &H_{\rm{elec}}(\mathbb{M}; N^\uparrow, N^\downarrow)\\
        =& \sum_{a=1}^N(-\frac{1}{2}\nabla_{\bm{A}}^2 - \sum_{a'=1}^{N'}\frac{Z_{a'}}{|\bm{r_a} - \bm{R'_{a'}}|}) + \hspace{-0.8em}\sum_{a=1, b>a}^{N,N}\frac{1}{|\bm{r_a} - \bm{r_b}|} \\

        \equiv& \sum_ah_1(a;\, \mathbb{M}) + \sum_{b>a}h_2(a,b)
        \label{eq:Helec}
    \end{aligned}
\end{equation}
where $h_1(a;\, \mathbb{M})$ is the one-electron operator that represents the kinetic energy of the $a$th electron and the Coulomb attraction between it and the nuclei in the system; $h_2(a, b)$ is the two-electron operator that represents the interaction between the $a$th and $b$th electrons.

One of the basic goals of electronic structure study is to solve the eigenvalue problem of the electronic problem:
\begin{equation}
    H_{\rm elec}\Phi_n = \mathcal{E}_n\Phi_n
\end{equation}
where $\Phi_n$ is the eigenfunction corresponding to the $n$th eigenenergy $\mathcal{E}_n$. $\Phi_n$ is a many-electron wavefunction that obeys the anti-symmetric constraint for identical fermions.

In an effort to alleviate the exponentially growing computational cost with respect to the total number of electrons, various ansatzes and numerical optimizations for the many-electron eigenfunctions have been proposed. Within the domain of classical numerical simulation, mean-field methods such as the Hartree--Fock method~\cite{szabo2012modern, slater1951simplification, valatin1961generalized} were proposed to approximate the ground state with a non-interacting fermionic state (i.e., a single Slater determinant state), while post-Hartree--Fock methods~\cite{szabo2012modern, bartlett1989alternative,  szalay2012multiconfiguration, mcardle2020quantum} try to further represent the ground state with the superposition of Slater determinant states based on one reference state (Hartree--Fock state) or multiple reference states.

With the successful demonstration on a quantum device to solve for the ground state of He--$\rm{H^+}$~\cite{peruzzo_variational_2014}, variational quantum eigensolver (VQE) has established its potential in utilizing noisy intermediate-scale quantum (NISQ) devices~\cite{preskill_quantum_2018, leymann2020bitter} for calculating electronic structure~\cite{mcclean2016theory, kandala2017hardware, cerezo_variational_2021, tilly2022variational}. However, due to the limited scale of NISQ devices, the applicable basis sets to perform qubit encoding for VQE have been restricted to small-size basis sets~\cite{tilly2022variational}. On the other hand, the interest around basis set optimization for electronic structure has grown recently~\cite{lehtola2015automatic, tamayo2018automatic, kasim2022dqc} with the development of computer algebra techniques such as automatic differentiation (AD)~\cite{griewank1989automatic, baydin2018automatic}.

In this paper, we propose a framework for customizing and optimizing composite basis sets based on Gaussian-type orbital (GTO). This is realized by the open-source software toolkit we developed in the Julia programming language~\cite{bezanson2012julia, bezanson2018julia}: Quiqbox.jl v0.5.2 (Quiqbox)~\cite{Quiqbox2022}.

We believe that a unified control and comparable design of basis sets play a crucial role in developing and benchmarking electronic structure algorithms in the present NISQ era when the exploration of the boundary between classical and quantum computational power is important for exploring practical quantum primacy~\cite{preskill_quantum_2018, daley2022practical, whitfield2022quantum}. We show that with the help of Quiqbox, one can control, modify, and improve the existing Gaussian basis sets composed of atomic orbitals (AOs), which are used to form a linear combination of atomic orbitals (LCAO) as a localized molecular orbital. Moreover, one can also construct and optimize more general basis sets beyond AOs, which outperform typical Gaussian basis sets or achieve similar performance while meeting specific requirements such as minimal basis set size.

The paper is organized as follows. In Sec.~\ref{sec2}, we elaborate on the proposal to carry out a more systematic study of basis set design in the NISQ era. In Sec.~\ref{sec3}, we describe the additional freedom our composite basis set generation framework brings for general electronic structure algorithms. We then derive the general expression of the basis set parameter gradient under the said framework and explain Quiqbox's specialized optimization architecture in Sec.~\ref{sec4}, followed by comparison and benchmarking against existing packages in Sec.~\ref{sec:performance}. In Sec.~\ref{sec:DO}, we demonstrate an example of using Quiqbox to construct basis sets with exotic properties. We conclude the paper and discuss the outlook of Quiqbox's development in Sec.~\ref{sec:conclusion}. Additionally, we include tables to summarize all the essential acronyms and notation conventions we use throughout this paper in the Appendix (Sec.~\ref{sec:app1}).

\section{Basis set design for the NISQ era}\label{sec2}

Let $E_0$ be the approximate ground-state energy of $H_{\rm elec}$ associated with a quantum many-electron system. Given an non-orthogonal basis set $\{\phi_n \,|\,  n=1,\dots,W\}$ that contains $W$ spatial orbitals with an overlap matrix $\bm{S}$ defined by $S_{ij}=\braket{\phi_i|\phi_j}$, we can construct an orthonormal basis set
\begin{equation}
    \{\varphi_m \,|\, \braket{\varphi_i|\varphi_j}=\delta_{ij}, m=1,\dots,W\},
\end{equation}
by performing a symmetric orthogonalization
\begin{equation}
    \varphi_m=\sum_nX_{nm}\phi_n, \;{\rm{with}}\; \bm{X} \equiv \bm{S}^{1/2}.
    \label{eq:no2o}
\end{equation}
Thus, a trial state $\ket{\Psi_0(\bm{c})}$ parameterized by $\bm{c}$ can be constructed from $\{\varphi_m\}$ following an ansatz of the ground state. Depending on the selected ansatz, $\bm{c}$ may not only characterizes the single-electron modes but also contains the amplitudes of multiple Slater determinant states formed by the modes. According to the variational principle, $E_0$ would be the lower bound of the energy expectation $\tilde{E_0}(\bm{c})$ with respect to $\ket{\Psi_0(\bm{c})}$:
\begin{equation}
    E_0 \leq \tilde{E_0}(\bm{c}) = \frac{\braket{\Psi_0(\bm{c})|H_{\rm elec}|\Psi_0(\bm{c})}}{\braket{\Psi_0(\bm{c})|\Psi_0(\bm{c})}}.
    \label{eq:E0v}
\end{equation}
Regardless of the choice of ansatz (e.g., the Hartree--Fock approximation) used to construct $\ket{\Psi_0(\bm{c})}$, $\tilde{E_0}(\bm{c})$ in Eq.~(\ref{eq:E0v}) can always be decomposed into three parts:
\begin{equation}
    \begin{aligned}
            \quad \bm{x} &= \bm{f_1}({\bm{A}}[\{\varphi_m\}]),\\
            \quad \bm{y} &= \bm{f_2}({\bm{B}}[\{\varphi_m\}]),\\
            \tilde{E_0}(\bm{c}) &= f_3(\bm{c};\, \bm{x},\, \bm{y}),
        \label{eq:E0general}
    \end{aligned}
\end{equation}
where vector functions $\bm{f_1}$: $\mathbb{R}^{W^2}\!\to\!\mathbb{R}^{W'}$ and $\bm{f_2}$: $\mathbb{R}^{W^4}\!\to\!\mathbb{R}^{W''}$ generate intermediate results $\bm{x}$ and $\bm{y}$ mapped from tensor ${\bm{A}}$ and tensor ${\bm{B}}$ which discretize the one-electron and two-electron operators respectively onto $\{\varphi_m\}$. The tensor elements of ${\bm{A}}$ and ${\bm{B}}$ are called electronic integrals (or molecular integrals in the context of molecular systems):
\begin{align}
    \begin{split}
        A_{mn} &= \hspace{-0.5em}\int \hspace{-0.4em}d\bm{r_a}\varphi_m^*(\bm{r_a})h_1(a;\, \mathbb{M})\varphi_n(\bm{r_a})\\
    &\equiv (\varphi_m|\varphi_n), \label{eq:A}
    \end{split}\\
    \begin{split}
    B_{mnpq} &= \hspace{-0.5em}\int \hspace{-0.4em}d\bm{r_a}d\bm{r_b}\varphi_m^*\hspace{-0.1em}(\hspace{-0.1em}\bm{r_a}\hspace{-0.1em})\varphi_p^*\hspace{-0.1em}(\hspace{-0.1em}\bm{r_b}\hspace{-0.1em})h_2(a,b)\varphi_q\hspace{-0.1em}(\hspace{-0.1em}\bm{r_b}\hspace{-0.1em})\varphi_n\hspace{-0.1em}(\hspace{-0.1em}\bm{r_a}\hspace{-0.1em})\\
    &\equiv (\varphi_m\varphi_n|\varphi_p\varphi_q).\label{eq:B}
    \end{split}
\end{align}
The third mapping function $f_3$ in Eq.~(\ref{eq:E0general}) takes t$\bm{x}$ and $\bm{y}$ as parameters, and $\bm{c}$ as the input variable. In this way, the calculation of $E_0$ is transformed into an optimization problem of finding the optimal value of $\bm{c}$:
\begin{equation}
    \bm{c^0} = \underset{\bm{c}}{{\rm arg}\,{\rm min}}\,\tilde{E_0}(\bm{c}).
    \label{eq:minimizeE}
\end{equation}
This can be done using either a self-consistent field (SCF) method~\cite{szabo2012modern, pulay1982improved, kudin2002black, hu2010accelerating}, or a gradient-based method realized by AD~\cite{yoshikawa2022automatic}, as the optimizer of $\bm{c}$.

One can see from Eq.~(\ref{eq:E0general}) that the expression of $\bm{f_1}$, $\bm{f_2}$, and $f_3$ depend on the choice of the ground state ansatz. More specifically, elements of $\bm{x}$ and $\bm{y}$ can be considered as the portions of ground-state energy contributed by single-Slater determinant states. Thus, for a multi-configuration ansatz, $W'$ and $W''$ would both be larger than one, as opposed to Hartree--Fock methods. Consequently, the vector space of $\bm{c}$ will be larger to presumably contain a better approximation of the true ground-state energy $\mathcal{E}_0$, which generally also makes it harder for the optimizer to find $\bm{c^0}$. An example of $\bm{f_1}$, $\bm{f_2}$ and $f_3$'s expression in the case of restricted Hartree--Fock as the applied electronic structure algorithm is shown in Eq.~(\ref{eq:dHF}) in Sec.~\ref{sec:generaldE} where we show the analytical differentiation process implemented in Quiqbox.

Eqs.~(\ref{eq:A}) and (\ref{eq:B}) indicate that the discretization of $H_{\rm elec}$ is affected by the choice of $\{\varphi_m\}$, which can not only introduce error due to the incompleteness of finite basis set but also affect the hardness of the electronic structure problem~\cite{ogorman_electronic_2021}. Following Eqs.~(\ref{eq:E0general}--\ref{eq:minimizeE}), the process of solving for the electronic structure properties, such as the ground-state energy, can be generalized into a standard procedure shown in FIG.~\ref{fig:ElectronicStructure1}. The electronic Hamiltonian of the system of interest is first discretized by a finite set of basis functions. After that, based on a model ansatz, a further approximation is applied (usually about the many-electron wavefunction) to construct the algebraic expression for the desired observable as an objective function (i.e., the ground-state energy) with respect to the discretized Hamiltonian. Thus, the objective function is implicitly a functional of the basis set. Finally, a numerical eigensolver (typically a variational optimizer) finds the approximate solution $\bm{c^0}$. Accordingly, the procedure of VQE for molecular electronic structure also falls into four similar steps as shown in FIG.~\ref{fig:ElectronicStructure2}.

\begin{figure}[htbp]
    \centering
    \begin{subfigure}[b]{0.48\textwidth}
         \centering
         \includegraphics[width=\textwidth]{images/ElectronicStructureWorkflow3p1.drawio.pdf}
         \caption{The general steps of solving a many-electron system.}
         \label{fig:ElectronicStructure1}
     \end{subfigure}
     \hfill
     \begin{subfigure}[b]{0.48\textwidth}
         \centering
         \includegraphics[width=\textwidth]{images/ElectronicStructureWorkflow3p2.drawio.pdf}
         \caption{Applying VQE to the electronic structure problem of molecular systems.}
         \label{fig:ElectronicStructure2}
     \end{subfigure}
     \hfill
     \begin{subfigure}[b]{0.48\textwidth}
         \centering
         \includegraphics[width=\textwidth]{images/ElectronicStructureWorkflow3p3.drawio.pdf}
         \caption{Using basis set design and molecular ensemble to construct well-defined computational tasks for NISQ or classical algorithms.}
         \label{fig:ElectronicStructure3}
     \end{subfigure}
    \caption{The overall procedure of solving the electronic Hamiltonian can be split into four steps for both classical and quantum algorithms as shown in FIG.~\ref{fig:ElectronicStructure1} and FIG.~\ref{fig:ElectronicStructure2}. A more systematic design and optimization of basis sets can aid the development of these algorithms for electronic structure problems, as shown in FIG.~\ref{fig:ElectronicStructure3}.}
    \label{fig:ElectronicStructure}
\end{figure}

From the generalized procedure of solving an electronic Hamiltonian, we explained the effect that a basis set has on both computational accuracy and resources. In the context of using VQE for electronic structure, this effect becomes more significant. Due to the limited number of qubits in NISQ devices for fermionic encoding~\cite{bravyi2002fermionic, havlivcek2017operator, setia2019superfast}, there is often a trade-off between the size of the basis set and the size of the studied many-electron systems~\cite{mcardle2020quantum, aspuru2005simulated}. Moreover, the applied Gaussian-based AOs (e.g., STO-3G~\cite{stewart1970small}, cc-pVDZ~\cite{dunning1989gaussian}) are, in fact, not fully optimized to minimize electronic structure properties such as the ground state energy given a fixed number of single-electron modes. This can be supported by the literature about the optimization of AOs~\cite{tachikawa1999simultaneous, tamayo2018automatic}. For the purpose of improving the practicality of using NISQ algorithms to solve electronic structure problems, we believe the systematic study of customizable basis set generation is of high importance. For instance, including basis set optimization and design as part of VQE's architecture can introduce the physical characteristics of many-electron systems in real space for VQE. This can potentially improve various aspects of the algorithm, such as trial qubit state preparation or optimizer configuration for the landscape of the corresponding parameter space.

Additionally, customizable basis set generation enables a better formulation of studying the electronic Hamiltonian as a computational task. It has been shown that approximating the ground state of the electronic Hamiltonian with a fixed basis set in real space (and fixed particle number) is QMA-complete~\cite{ogorman_electronic_2021, bookatz2012qma}. This suggests that the way a basis set discretizes the real space can affect the hardness of the electronic structure problem. We can further incorporate this idea into basis set designs. Given an ensemble of electronic Hamiltonians with the same electron number, by designing customized basis sets, we shall generate instances of electronic structure problems classified by different computational complexities (see FIG.~\ref{fig:ElectronicStructure3}). Not only does this help study the complexity of electronic Hamiltonians, but it also provides another metric for benchmarking the capability of quantum simulation algorithms (or hybrid algorithms that rely on quantum resources too)~\cite{manin1980computable, feynman2018simulating, aspuru2005simulated, lanyon2010towards, cerezo2021variational, daley2022practical} against pure classical approaches.

\section{More freedom in Gaussian basis set design} \label{sec3}
The basis sets used in electronic systems typically can be classified into three main categories: delocalized basis sets that rely on periodic functions like plane waves as basis functions~\cite{kresse1996efficient, kresse1996efficiency}, real-space basis sets that use localized functions such as wavelets distributed on grid points~\cite{pask1999real, genovese2008daubechies, ratcliff2020flexibilities}, and Gaussian basis sets that utilize GTO to form contracted functions to approximate Slater-type orbitals~\cite{dunning1977gaussian, hill2013gaussian}.

We also utilize GTO as the building block of more complicated basis functions for several reasons. First and foremost, GTO currently is one of the common choices for constructing basis set for molecular systems, the scale of which is applicable to VQE in the current state of NISQ devices~\cite{peruzzo_variational_2014, mcardle2020quantum, tilly2022variational}. This allows us to optimize and build new basis sets on the existing progress of GTO research. Secondly, the orbital integrals in Eqs.~(\ref{eq:A}) and (\ref{eq:B}) in the case of GTO can be calculated analytically for the most part~\cite{head1988method, petersson2009detailed}, making the numerical error of the tensor elements of the discretized electronic Hamiltonian easily controlled. Last but not least, with proper manipulations of Gaussian functions, various Gaussian-based basis sets with different features from AOs, such as floating Gaussian basis sets~\cite{hurley1988computation, frost1967floating, huber1979geometry}, Gaussian-based real-space basis sets~\cite{white2017hybrid, white2019multisliced}, and even-tempered basis sets~\cite{cherkes2009spanning}, can also be constructed.

By careful design, the variety of basis functions built from GTO can be greatly improved compared to existing Gaussian basis set generation methods in the standard software packages such as PySCF~\cite{sun2020recent} and Psi4~\cite{smith2020psi4}. In this section, we will formally introduce GTO (Gaussian-type orbital). Then, we focus on describing the proposed framework that provides extra freedom in using GTO to build composite basis functions that improve the customizability of Gaussian-based basis sets.

In general, the form of a GTO in Cartesian coordinates can be expressed as
\begin{equation}
    \begin{alignedat}{2}
        & &&\phi^{\rm{GTO}}(\bm{r};\bm{R},\mathpzc{i},\mathpzc{j},\mathpzc{k},\alpha)\\
        &\,\equiv\,&& N(x\!-\!X)^{\mathpzc{i}}(y\!-\!Y)^{\mathpzc{j}}(z\!-\!Z)^{\mathpzc{k}}{\rm{exp}}(-\alpha|\bm{r\!-\!R}|^2),\\
        \mathpzc{l} &\,\equiv\,&& \mathpzc{i} + \mathpzc{j} + \mathpzc{k},\quad

        \bm{r} \,\equiv\, [x,\, y,\, z],\quad
        \bm{R} \,\equiv\, [X,\, Y,\, Z],
    \end{alignedat}
    \label{eq:GTO}
\end{equation}
where $\alpha$ is the exponent coefficient that determines how diffused the GTO is, $\mathpzc{l}$ represents the orbital angular momentum, and $\bm{R}$ is the center position of the orbital. The most common use of GTOs is to form a contracted Gaussian-type orbital (CGTO)
\begin{equation}
    \begin{aligned}
        &\phi_m^{\rm{CGTO}}(\bm{r}; \bm{R_m}, \mathpzc{i}_m, \mathpzc{j}_m, \mathpzc{k}_m)\\
        \equiv& \sum_n d_n\,\phi^{\rm{GTO}}_n(\bm{r}; \bm{R_m}, \alpha_n, \mathpzc{i}_m, \mathpzc{j}_m,\mathpzc{k}_m).
    \end{aligned}
    \label{eq:CGTO}
\end{equation}
$d_n$ are called ``contraction coefficients'' as multiple concentric GTOs with equal angular momentum and different weights are ``contracted'' to approximate a Slater-type orbital (STO). STO is a type of orbital that provides an accurate representation of AOs~\cite{szabo2012modern}. The function curve of STO has a cusp near the orbital center due to its non-vanishing gradient, which GTO is incapable of reconstructing. However, forming a CGTO with appropriate contraction coefficients instead of a single GTO can alleviate this problem while retaining the computational efficiency of GTO~\cite{stewart1970small}. In our basis set generation framework, a basis function need not be a CGTO composed of concentric GTOs with the same angular momentum but can also be the linear combination of multi-center GTOs that may or may not have the same angular momentum. We define this type of generalized Gaussian-based orbital as \textbf{mixed-contracted} Gaussian-type orbital (MCGTO):
\begin{equation}
    \phi_{m'}^{\rm{MCGTO}}(\bm{r}) \equiv \sum_m \phi_m^{\rm{CGTO}}(\bm{r}; \bm{R_m}, \mathpzc{i}_m, \mathpzc{j}_m, \mathpzc{k}_m).
    \label{eq:MCGTO}
\end{equation}
MCGTO provides the first additional freedom in our framework, which makes it possible to construct more delocalized basis functions than a GTO to incorporate the interatomic bonding information before applying the SCF methods. For instance, one can use the molecular orbitals of small systems as the basis functions for larger combined systems.

Primitive variables that correlate multiple CGTOs through customized mapping functions are the second additional freedom to basis set generation. These correlations can have direct physical interpretations. For example, homonuclear diatomic molecules have spatial symmetry that can be described by the point group $D_{\infty h}$, i.e., they have rotational symmetry about their internuclear axis and reflection symmetry about a plane perpendicular to the axis~\cite{bishop1993group}. Such molecular symmetry affects the formation of their molecular orbitals. Therefore, it is reasonable to correlate the parameters of AOs during the AO-based basis set generation so that such symmetry is imposed to reduce the computational cost for later parameter optimization. Taking a step further, by using the primitive variables that correlate basis functions, grid-based Gaussian basis sets become available. In such a scenario, the center position of each basis function would also be the location of a grid point controlled by the spacing of the grid box. A mapping function from the single spacing parameter (i.e., the sparsity of gird points) to the CGTO center positions greatly lower the total number of basis set parameters while maintaining the topology of the grid box. Moreover, the exponent and contraction coefficients of the CGTOs can be correlated to reduce parameters further. One example of such type of basis set is shown in FIG.~\ref{fig:gb}.

\begin{figure}[htbp]
    \centering
    \includegraphics[width=0.45\textwidth]{images/gridBsetL.pdf}
    \caption{A $1\!\times\!2\!\times\!1$ grid-based Gaussian basis set for $\rm{H_2}$.

    Based on the geometric symmetry of $\rm{H_2}$, constraints on CGTO parameters can be imposed. For example, the CGTOs located on the grid points with the same color can have identical exponent and contraction coefficients. Correlated grid points can be directly generated using the helper function \texttt{GridBox} in Quiqbox.}
    \label{fig:gb}
\end{figure}

Combining the design of a generalized Gaussian-based orbital, MCGTO, and the additional correlation configuration among CGTO parameters, we complete our design of a more customizable generation of Gaussian-based basis sets. To provide a straightforward comparison between it and the standard basis set generation method, the computational graphs of both methods in the context of electronic ground-state energy computation are shown in FIG.~\ref{fig:EANN}.

\begin{figure*}[htpb]
    \centering
    \begin{subfigure}[b]{0.65\textwidth}
         \centering
         \includegraphics[width=0.99\textwidth]{images/ADchainN8_0.drawio.pdf}
    \caption{The complete computational graph of an algorithm solving for the electronic ground-state energy in standard electronic structure packages.}
         \label{fig:NN1}
     \end{subfigure}
     \hfill
     \begin{subfigure}[b]{0.65\textwidth}
         \centering
         \includegraphics[width=0.99\textwidth]{images/ADchainN8.drawio.pdf}
         \caption{The computational graph with extra connectivity and layers after replacing the standard basis set generation method in FIG. \ref{fig:NN1} with our framework. In the light blue box, the leftmost added layer represents the  primitive parameters $\{\theta_i\}$ that control CGTO parameters $\{p_i\}$ through mapping functions, and the rightmost layer represents MCGTOs as the output basis functions. Additionally, since the CGTOs can share parameters in Quiqbox, more connectivity (blue lines) between the CGTO parameters and CGTOs becomes available.}
         \label{fig:NN2}
     \end{subfigure}
    \caption{The comparison between the standard basis set generation method and our framework for finding the electronic ground-state energy. The light blue box represents the basis set construction. The light purple box represents the electronic algorithm to compute the approximate ground-state energy $E_0$. The yellow dashed box represents the chosen eigensolver that produces $\bm{c^0}$ as the final input to obtain $E_0$.}
    \label{fig:EANN}
\end{figure*}

Aside from the common input parameters indicated in Eq.~(\ref{eq:Helec}), the parameters from the generated basis set $\{\phi_n\}$ can be considered the primary parameters that affect $E_0$. Specifically, for the standard Gaussian basis set generation method (FIG.~\ref{fig:NN1}), only the CGTO parameters $\{p_i\,|\,i=1\dots N_p\}$ can be the basis set parameters. However, in our basis set generation framework (FIG.~\ref{fig:NN2}), such restrictions are lifted by allowing more primitive parameters $\{\theta_i\,|\,i=1\dots N_\theta\}$ to control the orbital parameters. On top of this, correlations between CGTOs can be implemented with shared orbital parameters among them, and MCGTOs serve as the finally generated basis functions. As we fixed the coefficients generated by the eigensolver as shown in the yellow dashed box (along with other parameters in the grey circled boxes in FIG.~\ref{fig:EANN}), $E_0$ is mapped from $\{\theta_i\}$ through layers of differentiable nodes connected by various mapping functions. Thus, $\{\theta_i\}$ as the basis set parameters can be optimized using variational methods based on the analytical gradient of $E_0$ with respect to them, which forms the process of basis set optimization.

\section{Variational optimization of basis set parameters} \label{sec4}

\subsection{Generalized analytical electronic structure gradient with respect to basis set parameters}\label{sec:generaldE}

Let $\{\theta_i\,|\,\theta_i\in \mathbb{R}, i=1,\dots, N_\theta\}$ be a set of unique real parameters (meaning they are distinguishable from each other by their symbols and indices regardless of whether they have the same value) that $E_0$ explicitly or implicitly depends on. To minimize $E_0$ with respect to $\{\theta_i\}$, we need to calculate the gradient of $E_0$ with respect to the vectorized parameters $\bm{\theta} = [\theta_1, \dots, \theta_i, \dots, \theta_{N_\theta}]^T$, the analytical form of which can be expressed as~\cite{jo1983ab, hurley1988computation}:
\begin{equation}
    \frac{dE_0}{d\bm{\theta}} = \sum_m\frac{\partial E_0}{\partial \varphi_m}\frac{\partial \varphi_m}{\partial \bm{\theta}} + \braket{\Psi_0|\frac{\partial H_{\rm elec}}{\partial \bm{\theta}}|\Psi_0}.
    \label{eq:dE0}
\end{equation}
The first term in Eq.~(\ref{eq:dE0}) represents the explicit dependency of the parameters based on the ground state approximation, and the second term represents the effective force from the Hellmann--Feynman theorem~\cite{hurley1988computation, politzer2018hellmann}. For the purpose of a systematic basis set optimization, $\{\theta_i\}$ represents a set of unique and independent basis set parameters to be optimized. Assuming the electronic Hamiltonian is independent of these parameters (e.g., the basis function centers can be detached from the nuclei locations, forming floating basis functions~\cite{hurley1988computation, frost1967floating, huber1979geometry}), the second term in Eq.~(\ref{eq:dE0}) is eliminated. Thus, the partial derivative of $E_0$ with respect to a single basis set parameter $\theta_i$ reduces to
\begin{equation}
    \frac{\partial E_0}{\partial \theta_i} = \sum_m\frac{\partial E_0}{\partial \varphi_m}\frac{\partial \varphi_m}{\partial \theta_i}.
    \label{eq:E0}
\end{equation}

Combining Eq.~(\ref{eq:E0}) with Eqs.~(\ref{eq:E0general}--\ref{eq:B}), we can derive a more detailed formula of the parameter partial derivative for a general electronic structure algorithm:
\begin{equation}
    \begin{aligned}
        \frac{\partial E_0}{\partial \theta_i} = &\sum_j\frac{\partial E_0}{\partial x_j}\frac{\partial x_j}{\partial \bm{A}}\cdot\frac{\partial \bm{A}}{\partial \theta_i} + \sum_k\frac{\partial E_0}{\partial y_k}\frac{\partial y_k}{\partial \bm{B}}\cdot\frac{\partial \bm{B}}{\partial \theta_i}\\
        = &\sum_{mn} (\nabla_{\hspace{-0.3em}\bm{A}} E_0 \hspace{-0.2em}\odot\hspace{-0.2em} \partial_{\theta_i} \bm{A})_{mn} \hspace{-0.3em}+\hspace{-0.6em} \sum_{mnpq}\hspace{-0.2em}(\nabla_{\hspace{-0.3em}\bm{B}} E_0 \hspace{-0.2em}\odot\hspace{-0.2em} \partial_{\theta_i} \bm{B})_{mnpq} \label{eq:dE1}
    \end{aligned}
\end{equation}
where $\odot$ denotes element-wise multiplication, and
\begin{align}
    \begin{split}
        \hspace{-1.2em}(\nabla_{\bm{A}}E_0)_{mn} \hspace{-0.2em}=& (\sum_j\frac{\partial E_0}{\partial x_j}\frac{\partial x_j}{\partial \bm{A}})_{mn},\label{eq:dEdA}
    \end{split}\\
    \begin{split}
        \hspace{-1.2em}(\nabla_{\bm{B}}E_0)_{mnpq} \hspace{-0.2em}=& (\sum_k\frac{\partial E_0}{\partial y_k}\frac{\partial y_k}{\partial\bm{B}})_{mnpq},\label{eq:dEdB}
    \end{split}\\
    \begin{split}
        (\partial_{\theta_i} \bm{A})_{mn} \hspace{-0.2em}=& (\partial_{\theta_i} \varphi_m|\varphi_n) + (\varphi_m|\partial_{\theta_i} \varphi_n),\label{eq:dA}
    \end{split}\\
    \begin{split}
        (\partial_{\theta_i} \bm{B})_{mnpq} \hspace{-0.2em}=&
        (\partial_{\theta_i} \varphi_m\varphi_n|\varphi_p\varphi_q)\hspace{-0.2em}+\hspace{-0.2em}
        (\varphi_m\partial_{\theta_i} \varphi_n|\varphi_p\varphi_q) + \\&
        (\varphi_m\varphi_n|\partial_{\theta_i} \varphi_p\varphi_q)\hspace{-0.2em}+\hspace{-0.2em} (\varphi_m\varphi_n|\varphi_p\partial_{\theta_i} \varphi_q).
        \label{eq:dB}
    \end{split}
\end{align}
The term $\sum_l\frac{\partial E_0}{\partial c_l}\frac{\partial c_l}{\partial \theta_i}$ is not included in Eq.~(\ref{eq:dE1}) because the value of $\nabla_{\bm{c}}E_0$ is calculated at $\bm{c}=\bm{c^0}$, which is always a null vector as
$E_0$ is the global minimum of $f_3(\bm{c})$ at $\bm{c^0}$ according to Eqs. (\ref{eq:E0v}) and (\ref{eq:minimizeE}).

To further expand Eqs.~(\ref{eq:dA}) and (\ref{eq:dB}), based on Eq.~(\ref{eq:no2o}), the partial derivatives of the orthonormal $\{\varphi_m\}$ can be expressed as~\cite{jo1983ab}:
\begin{equation}
    \begin{aligned}
        \partial_{\theta_i}\varphi_m &=  \partial_{\theta_i}\phi_i X + \phi_i \partial_{\theta_i} X,\\
        (\partial_{\theta_i} X)_{mn} &= -\sum_{pq} V_{mp}\frac{(V^\dagger\partial_{\theta_i} SV)_{pq}}{(\lambda_p\lambda_q)^{1/2}(\lambda_p^{1/2}+\lambda_q^{1/2})}V^*_{nq}, \\
        (\partial_{\theta_i} S)_{ij} &= \braket{\partial_{\theta_i}\phi_i|\phi_j} + \braket{\phi_i|\partial_{\theta_i}\phi_j},
    \end{aligned}
    \label{eq:dE3}
\end{equation}
where $\bm{V_p} = [V_{1p},\dots, V_{Wp}]^T$ is the $p$th eigenvector of $\bm{S}$ with corresponding eigenvalue $\lambda_p$.

We can even further expand the partial derivative of each $\phi_n$ as the product of two terms within the summation:
\begin{equation}
    \partial_{\theta_i}\phi_n =\!\! \sum_{p_j\in\{p_j^n\}}\!\!\frac{\partial\phi_n}{\partial p_j}\frac{\partial p_j}{\partial\theta_i}
    \label{eq:dphi}
\end{equation}
where $\{p_j^n\}$ are all the differentiable CGTO parameters (i.e., the exponent coefficients, the contraction coefficients, and the GTO center) $\phi_n$ explicitly depends on
\begin{equation}
    \begin{aligned}
        \{p_j^n\} = \{&\alpha_1^n,\mydots,\alpha^n_{N^n_{\rm{GTO}}}, d_1^n,\mydots,d^n_{N^n_{\rm{GTO}}}, \\
        &X^n_1, Y^n_1, Z^n_1, \mydots, X^n_{N^n_{\rm{cen}}}, Y^n_{N^n_{\rm{cen}}}, Z^n_{N^n_{\rm{cen}}}\}.
    \end{aligned}
\end{equation}
The superscript $n$ of the included CGTO parameters (e.g., $\alpha_1^n$) indicates that they are corresponding to $\phi_n$. Specifically, $N^n_{\rm{GTO}}$ is the total number of GTOs in $\phi_n$, and $N^n_{\rm{cen}}$ is the total number of GTO centers that group the GTOs in $\phi_n$ into concentric GTOs. Thus, all the explicit CGTO parameters from $\{\phi_n\}$ would be
\begin{equation}
    \{p_j\}  = \{p_j^1\} \cup \{p_j^2\} \cup \dots \cup \{p_j^W\}.
    \label{eq:ebsp}
\end{equation}
In the recent Gaussian basis set optimization method~\cite{tamayo2018automatic}, all basis functions in
$\{\phi_n\}$ are CGTOs, so $N^n_{\rm{cen}}=1$, and the basis set parameters are just a subset of all the explicit parameters as shown in Eq.~(\ref{eq:ebsp}). This effectively reduces Eq.~(\ref{eq:dphi}) to only one term, a partial derivative with respect to $\theta_i$, as $\partial p_j/\partial \theta_i = \delta_{ij}$.

By calculating the analytical parameter gradient, one can perform variational optimization for the generated basis set. In fact, optimization of minimal AO-based Gaussian basis sets with respect to the Hartree--Fock energy has been implemented in a Python software package named ``DiffiQult''~\cite{tamayo2018automatic}. Optimizing the basis set parameters according to the Hartree--Fock energy is a good starting point for many more advanced electronic structure algorithms, such as multi-configuration and multi-reference computations that are based on the Hartree--Fock state. Thus, when we try to implement the generalized basis set optimization in Quiqbox, we also set the objective function of basis set parameter optimization to the Hartree--Fock energy. Specifically, in the case of the restricted closed-shell Hartree--Fock method (RHF), Eq.~(\ref{eq:E0general}) would be specialized as~\cite{szabo2012modern}
\begin{equation}
    \begin{aligned}
        \bm{T^x} =& f_1(\bm{A}) := \bm{A},\\
        \bm{T^y} =& f_2(\bm{B}) := \bm{B},\\
         \tilde{E_0}(\bm{D}) =&  f_3(\bm{D};\,\bm{T^x},\,\bm{T^y})\\ :=&\sum_{\mu,\nu}\hspace{-0.2em}D_{\mu\nu}({T^{\bm{x}}}_{\hspace{-0.8em}\mu\nu}\!+\! \sum_{\lambda,\sigma}\hspace{-0.2em}D_{\lambda\sigma}[2{T^{\bm{y}}}_{\hspace{-0.7em}\mu\nu\sigma\lambda}\!-\!{T^{\bm{y}}}_{\hspace{-0.7em}\mu\lambda\sigma\nu}]),
    \end{aligned}
    \label{eq:dHF}
\end{equation}
where $\bm{D}$ is the density matrix of the electrons with the same spin (assuming $N^\uparrow = N^\downarrow$ and each occupied molecular spatial orbital is doubly occupied), which is a matrix format of $\bm{c}$. Similarly, $\bm{T^x}$ and $\bm{T^y}$ are, respectively, $\bm{x}$ and $\bm{y}$ reformatted as tensors.

\subsection{Variational parameter optimization based on hybrid differentiation design}

If one can find all the elementary differentiation rules that will be used to compute Eq.~(\ref{eq:dE1}), along with additional necessary AD rules to implement the calculation as a programming function in an existing AD software library, they will be able to complete the code for differentiation-based parameter optimization with minimal effort. Specifically, reverse-mode AD would be the preferable AD variant to implement as the number of basis set parameters is normally larger than the output dimension of the objective function (e.g., $E_0$ for the Hartree--Fock energy)~\cite{baydin2018automatic}. However, reverse-mode AD is not well defined when the diagonalization process of matrices with degenerate eigenvalues emerges in the programming code. This led to the compromise of using the rather inefficient forward-mode AD to perform the gradient calculation in DiffiQult~\cite{tamayo2018automatic}.

Instead of fully relying on forward-mode AD or adapting reverse-mode AD with customized AD rules as a workaround, we designed a hybrid differentiation engine in Quiqbox. This differentiation engine combines both AD and symbolic differentiation (SD) so that it can provide both efficiency and extensibility to Quiqbox's parameter optimization procedure. The fundamental design of it is to divide the parameter differentiation, i.e., Eq.~(\ref{eq:dE1}), into three parts:
\begin{enumerate}
    \item compute $\nabla_{\bm{A}}E_0$ and $\nabla_{\bm{B}}E_0$;
    \item compute $\partial_{\theta_i} \bm{A}$ and $\partial_{\theta_i} \bm{B}$;
    \item compute $\partial_{\theta_i} E_0 = f_4(\nabla_{\bm{A}}E_0,\:\nabla_{\bm{B}}E_0,\: \partial_{\theta_i} \bm{A},\:\partial_{\theta_i} \bm{B})$
\end{enumerate}
where
\begin{equation}
    \begin{aligned}
        &f_4(\bm{T^1},\:\bm{T^2},\:\bm{T^3},\:\bm{T^4})\\
        :=& \sum_{mn} (\bm{T^1} \odot \bm{T^3})_{mn} + \sum_{mnpq}\hspace{-0.2em}(\bm{T^2} \odot \bm{T^4})_{mnpq}.
    \end{aligned}
\end{equation}
The first part requires the computation of Eqs.~(\ref{eq:dEdA}) and (\ref{eq:dEdB}), and hence only depends on the choice of electronic structure algorithm, whereas the second part requires the computation of Eqs.~(\ref{eq:dA}) and (\ref{eq:dB}), and thus only depends on the generated basis set. The third part depends only on tensor contractions.

By dividing the chain of analytical differentiation into computational modules, we lose the capability of utilizing a single full AD computation. However, since each module can be independently maintained and updated, doing so can reduce the risk of code breaking due to bugs coming from external AD libraries. This is crucial from the perspective of developing sustainable scientific programs. More importantly, we now gain the flexibility to design tailored optimization for each module.

Particularly, for part one, one can choose between AD and SD depending on the complexity of Eqs.~(\ref{eq:dEdA}) and (\ref{eq:dEdB}) determined by the applied electronic structure algorithm. As for part two, the differentiation of each $\phi_n$ can be optimized by transforming it into the efficient Gaussian-based electronic integral computation. Due to the generality of MCGTO, as it can represent any linear combination based on Eqs.~(\ref{eq:GTO}--\ref{eq:MCGTO}), any partial derivative of MCGTO, including a CGTO, can be fully expressed by a new MCGTO. This means the expression of MCGTO is complete for arbitrary order of gradients of MCGTO. Thus, we can define a new basis function $\psi_n^{\theta_i}$ with respect to $\theta_i$ and $\phi_n$ such that
\begin{equation}
    \psi_n^{\theta_i} \equiv \partial_{\theta_i}\phi_n.
\end{equation}
This effectively transforms the calculation of Eqs.~(\ref{eq:dA}) and (\ref{eq:dB}) into the calculation of electronic integrals of MCGTOs that compose a new \textbf{derivative basis set}:
\begin{equation}
    \mathbb{D} \equiv \{\phi_n\} \cup \{\psi_n^{\theta_1}\} \cup \dots \cup \{\psi_n^{\theta_{N_\theta}}\}.
\end{equation}
Quiqbox utilizes AD to compute all the derivatives that directly contribute to contraction coefficients but uses SD to generate new MCGTOs as the partial differentiation of basis functions. Due to the powerful multiple dispatch feature in the Julia programming language~\cite{bezanson2012julia, bezanson2018julia}, $\mathbb{D}$ can be generated and stored efficiently using the code that supports ahead-of-time (AOT) code compilation.

The overall architecture of Quiqbox's parameter optimization procedure applied with the aforementioned hybrid-differentiation design is shown in FIG.~\ref{fig:po}. The computation of the objective function is also separated from the computation of the gradient, as it provides a stationary point of $\bm{c}$ (i.e., $\bm{c^0}$) before each differentiation computation. As a result, the choice of electronic structure algorithm can also be expanded separately by programming new computational modules. This architecture, by design, makes it easier for Quiqbox to adapt more advanced post-Hartree--Fock algorithms for basis set optimization. As for its performance, we shall discuss it in the next section.

\begin{figure}[htbp]
    \centering
    \includegraphics[width=0.99\linewidth]{images/QuiqboxVO5.drawio.pdf}
    \caption{The diagram of the basis set optimization procedure. Within each optimization iteration cycle, both AD and SD are applied to parameter differentiation. Particularly, partial derivatives of basis functions with respect to basis set parameters are then transformed into a new derivative basis set $\mathbb{M}$ such that the cost of gradient computation can mostly rely on the direct electronic integral calculation in the teal box.}
    \label{fig:po}
\end{figure}

\section{Performance} \label{sec:performance}
To demonstrate the performance of Quiqbox (version 0.5.2 running on Julia 1.8.2), we compare it in three scenarios with existing packages that have limited similar features. The tests were run on a cluster node equipped with an Intel Xeon E7-8891 v4 CPU.

\subsection{The Hartree--Fock method} \label{sec:HF}

As the Hartree--Fock method is currently the main objective function for basis set optimization in Quiqbox, it is crucial to verify the accuracy and efficiency of Quiqbox in performing such computations. We used the equal-spacing hydrogen chain, with the distance between adjacent nuclei being $1.0$ a.u., as a test system for how the computation accuracy and efficientcy scale with respect to the growth of system size. As for the package to benchmark against, we chose the widely-used ab initio electronic structure package PySCF (version 2.1.1 running on Python 3.9.12). The Hartree--Fock method, specifically RHF, was applied using each package with all the default settings except that the convergence threshold of RHF SCF iteration was manually set to be $1\!\times\!10^{-9}\,\rm{Ha}$. The multi-threading optimization was enabled for both packages, and the total number of available threads was set to 16. The benchmarking results are shown in FIG.~\ref{fig:HFc}. In FIG.~\ref{fig:HFc_dE}, we show that with the hydrogen chain lengthening, Quiqbox can maintain better accuracy than PySCF, especially for smaller basis sets. Additionally, as shown in FIG.~\ref{fig:HFc_dE}, Quiqbox, on average, needs less time to complete the SCF process than PySCF but has worse scaling. This is mainly due to the increasing relative inefficiency of Quiqbox's electronic integral computation compared to PySCF.

\begin{figure}[htbp]
    \centering
    \begin{subfigure}[b]{0.48\textwidth}
         \centering
         \includegraphics[width=0.99\textwidth]{images/HF_HCdE_Thread-16.0_Time-2022-11-09T13-07-11-307.pdf}
    \caption{The Hartree--Fock energies computed by PySCF minus the same quantities computed by Quiqbox.}
         \label{fig:HFc_dE}
     \end{subfigure}
     \hfill
     \begin{subfigure}[b]{0.48\textwidth}
         \centering
         \includegraphics[width=0.99\textwidth]{images/HF_HCt_Thread-16.0_Time-2022-11-09T13-07-11-779.pdf}
         \caption{The relative mean run time for Quiqbox's Hartree--Fock SCF iteration with respect to PySCF. The upper error bar indicates the mean time plus one standard deviation, and the lower error bar records the minimal run time in each testing case.}
         \label{fig:HFc_dt}
     \end{subfigure}
    \caption{The benchmark for Quiqbox's RHF function against PySCF.}
    \label{fig:HFc}
\end{figure}

\subsection{Electronic integrals} \label{sec:Int}

Instead of directly relying on an external Gaussian-type electronic integral library, we have written our own integral engine (also in Julia) as part of the core functions of Quiqbox, based on reference~\cite{petersson2009detailed}. Additionally, we have added specialized optimization methods to avoid repetitive integral computation from the resulting tensor elements and the orbital expressions with shared parameters. On top of that, the integral engine is fully compatible with all kinds of MCGTOs. By contrast, PySCF relies on the C-based electronic integral library Libcint~\cite{libcint2015}.

\begin{figure}[htp]
    \centering
    \begin{subfigure}[b]{0.48\textwidth}
         \centering
         \includegraphics[width=0.99\textwidth]{images/Int_HC_t_Hcore_Thread-16.0_Time-2022-11-13T16-41-25-675.pdf}
    \caption{The mean run time comparison of core Hamiltonian (one-electron tensor) computation.}
         \label{fig:Intc_dE}
     \end{subfigure}
     \hfill
     \begin{subfigure}[b]{0.48\textwidth}
         \centering
         \includegraphics[width=0.99\textwidth]{images/Int_HC_t_eeI_Thread-16.0_Time-2022-11-13T16-41-25-698.pdf}
         \caption{The mean run time comparison of electron-electron interaction (two-electron tensor) computation.}
         \label{fig:Intc_dt}
     \end{subfigure}
    \caption{The relative run time of electronic integral computation of Quiqbox's native electronic integral engine with respect to Libcint used in PySCF. Same as in FIG.~\ref{fig:HFc_dt}, the upper error bar indicates the mean time plus one standard deviation, and the lower error bar records the minimal run time in each testing case.}
    \label{fig:Intc}
\end{figure}

We compared the electronic integral computations between Quiqbox and Libcint for the same basis sets and molecules used in the Hartree--Fock tests. Quiqbox and Libcint agreed on all the computation results within a difference of $5\!\times\!10^{-7}$ a.u. In respect of computation efficiency, as shown in FIG.~\ref{fig:Intc}, for smaller basis sets such as STO-3G~\cite{stewart1970small} and 6-31G~\cite{hehre1972self}, the minimal run time difference between Quiqbox's integral engine and Libcint is within the same magnitude in best cases. For the cc-pVDZ~\cite{dunning1989gaussian} basis set (with optimized general contraction), the differences increase significantly, particularly for electron-electron interaction integrals. One possible reason is that Quiqbox's integral engine does not take the spatial distribution of identical floating atomic orbitals (only $\bm{R}$ changes) into consideration to optimize the integral computation. For example, in an $\rm{H_8}$ chain, the integral computation of orbitals located at the 5th, 6th, 7th, and 8th hydrogen nuclei, should produce the same result as the integral computation of orbitals located at the 1st, 2nd, 3rd, 4th hydrogen nuclei after appropriate relabeling and reordering; thus they need not be recomputed. This type of issue can be avoided more easily by implementing more powerful integral computation reduction functions. Moreover, the author of Libcint has mentioned in their technical literature that they implemented hardware-level optimization such as CPU caching and addressing optimization in Libcint, as well as sparse matrix computation optimization on the software level~\cite{libcint2015}, both of which are currently lacking in Quiqbox's integral engine.

We believe that with more optimization in the future, we can improve Quiqbox's integral engine to shorten the efficiency gap to Libcint. It is worth mentioning that implementing our native integral engine allows us to realize proper normalization for Cartesian-coordinate-based GTO integrals, which has been a legacy problem for Libcint~\cite{libcintNorm}, is crucial for realizing MCGTOs in Quiqbox.

\subsection{Basis set parameter optimizations} \label{sec:po}

To benchmark the performance of Quiqbox's basis set optimization, we compared it with DiffiQult~\cite{tamayo2018automatic}. Since DiffiQult does not support basis sets beyond LACO, the selected basis set was the minimal basis set STO-3G~\cite{stewart1970small} for both packages. Moreover, DiffiQult does not support multi-threading computation, so in the case of $\rm{H_2}$, Quiqbox is constrained to only use one thread both when running native Julia functions and when calling the default back-end linear algebra library OpenBLAS~\cite{wang2013augem}.

Nevertheless, Quiqbox's optimization function shows a significant advantage against DiffiQult as presented in TABLE~\ref{tb:poc}. Using the same numerical optimizer as DiffiQult, L-BFGS~\cite{liu1989limited} with a line search algorithm using the strong Wolfe (SW) conditions, Quiqbox provided a 300 times speed-up. By replacing the line search algorithm with the one used in the nonlinear conjugate gradient method proposed by Hager and Zhang (HZ)~\cite{hager2005new}, Quiqbox was able to achieve even faster convergence both in terms of the number of steps and overall computing time. This allowed a total of about 450 times speed-up relative to DiffiQult. When only using a first-order stochastic gradient descent optimizer like Adam~\cite{kingma2014adam}, Quiqbox still outperformed DiffiQult with over 50 times speed-up. To be noted, we included two testing cases for $\rm{H_2}$: one was for optimizing only the exponent and contraction coefficients ($\alpha$ and $d$), and another one also added the GTO centers ($\bm{R}$). The reason for doing so is that DiffiQult does not support optimizing GTO centers while simultaneously optimizing the other two, which is not a limitation for Quiqbox.

The drastic improvement of basis set optimization's efficiency pushes forward the practicality of variational basis set optimization, especially for larger systems. For example, in the case of optimizing the STO-3G basis set for LiH, Quiqbox was able to finish the optimization within two hours. In contrast, DiffiQult had not achieved convergence even after running for more than 11 hours. This result is included in TABLE~\ref{tb:poc} as well.

\begin{center}
    \begin{table*}[htbp]
    \setlength{\tabcolsep}{4.5pt}
    \renewcommand{\arraystretch}{1.3}
    \begin{tabular}{ccccccrrrrr}
    \toprule[1.5pt]
    \multirow{2}{*}{Molecule} & \multicolumn{2}{c}{Parameters} & \multirow{2}{*}{Package} & \multirow{2}{*}{Thread} & \multirow{2}{*}{Optimizer} & \multicolumn{1}{c}{\multirow{2}{*}{$E_0$ (Ha)}} & \multicolumn{1}{c}{\multirow{2}{*}{Step}} & \multicolumn{1}{c}{\multirow{2}{*}{Time (sec)}} & \multicolumn{1}{c}{\multirow{2}{*}{Time/Step}} & \multicolumn{1}{c}{\multirow{2}{*}{Speed-up}} \\
     & Types & Quantity &  &  &  & \multicolumn{1}{c}{} & \multicolumn{1}{c}{} & \multicolumn{1}{c}{} & \multicolumn{1}{c}{} & \multicolumn{1}{c}{} \\ \midrule[1.2pt]\midrule[1.2pt]
    \multirow{8}{*}{$\rm{H_2}$} & \multirow{4}{*}{$\alpha$, $d$} & 12 & DiffiQult & 1 & L-BFGS (SW) & $-1.837307$ & 30 & 901.5 & 30.05 & 1.00 \\ \cline{3-11}
     &  & \multirow{3}{*}{12 (6)} & \multirow{3}{*}{Quiqbox} & \multirow{3}{*}{1} & L-BFGS (SW) & $-1.837307$ & 34 & 3.0 & 0.09 & 300.50 \\
     &  &  &  &  & L-BFGS (HZ) & $-1.837307$ & 18 & 2.0 & 0.11 & 450.75 \\
     &  &  &  &  & Adam & $-1.837306$ & 489 & 12.4 & 0.03 & 72.70 \\ \cmidrule[1.2pt]{2-11}
     & \multirow{4}{*}{$\alpha$, $d$, $\bm{R}$} & 18 & DiffiQult & 1 & L-BFGS (SW) & $-1.840628$ & 37 & 1122.9 & 30.35 & 1.00 \\ \cline{3-11}
     &  & \multirow{3}{*}{18 (9)} & \multirow{3}{*}{Quiqbox} & \multirow{3}{*}{1} & L-BFGS (SW) & $-1.840820$ & 35 & 3.4 & 0.10 & 330.26 \\
     &  &  &  &  & L-BFGS (HZ) & $-1.840820$ & 21 & 3.1 & 0.15 & 362.22 \\
     &  &  &  &  & Adam & $-1.840782$ & 502 & 16.8 & 0.03 & 66.84 \\ \midrule[1.2pt]
    \multirow{3}{*}{LiH} & \multirow{3}{*}{$\alpha$, $d$, $\bm{R}$} & 30 & DiffiQult & 1 & L-BFGS (SW) & N/A & $> \!19$ & $> \!4\!\times\!10^4$ & N/A & 1.00 \\ \cline{3-11}
     &  & \multirow{2}{*}{30 (28)} & Quiqbox & 1 & \multirow{2}{*}{L-BFGS (HZ)} & \multirow{2}{*}{$-8.963143$} & \multirow{2}{*}{433} & 6330.1 & 14.62 & $> \!6.32$ \\
     &  &  & Quiqbox & 8 &  &  &  & 5239.0 & 12.10 & $> \!7.64$ \\
     \bottomrule[1.5pt]
    \end{tabular}
    \caption{The benchmark for Quiqbox's parameter optimization of STO-3G basis set concerning the RHF energy against DiffiQult. In the column ``Parameters'', $\alpha$, $d$, and $\bm{R}$ respectively represent all exponent coefficients, all contraction coefficients, and all orbital centers from the basis functions in the basis set. Each row of the sub-column ``Quantity'' shows the total number of parameters, with the extra number in parentheses (when they appear) being the reduced number of unique parameters to be optimized after imposing parameter correlations using Quiqbox's basis set generation framework.}
    \label{tb:poc}
    \end{table*}
\end{center}

\section{Example: building completely delocalized orbitals} \label{sec:DO}

Now we showcase the extra customizability Quiqbox provides. Consider a basis set $\{\phi'_n\}$, where the Coulomb interaction between every two different spatial orbitals is equal to a constant $\mathcal{J}$, and the exchange-correlation between every two different spatial orbitals is equal to a constant $\mathcal{K}$, i.e.,
\begin{equation}
    \begin{aligned}
        \forall \:\phi'_i,\: \phi'_j &\in \{\phi'_n\}\: {\rm{s.t.}} \:\phi'_i \neq \phi'_j,\\
        \mathcal{J} &\equiv (\phi'_i\phi'_i|\phi'_j\phi'_j),\\
        \mathcal{K} &\equiv (\phi'_i\phi'_j|\phi'_i\phi'_j).
    \end{aligned}
\end{equation}
Such basis sets can be defined as basis sets with ``completely delocalized orbitals (CDOs)'', resulting in permutation invariance of electronic modes with the same spin configuration.

One example of such a basis set is three pairs of identical floating CGTOs that all intersect at the same point such that every line segment that connects two CGTOs (at their centers) in the same pair is a perpendicular bisector of other pairs. The two CGTOs of each pair are combined as a new orbital in the basis set. If each CGTO contains $n$ GTOs, then in total, only $n$ floating GTOs with unique exponents and contraction coefficients are used to construct all three basis functions of the basis set.

The geometry of this basis set is shown in FIG.~\ref{fig:CDO3}, and the code for constructing an instance of such type of basis set with $n\!=\!4$ (i.e., CDO3-4G) in Quiqbox is shown in line 10-36 of Listing \ref{lst:DObasis}.

\begin{figure}[htbp]
    \centering
    \includegraphics[width=0.45\textwidth]{images/DOBsetL.pdf}
    \caption{The geometry of a delocalized basis set of three orbitals (CDO3) concerning $\rm H_2$.}
    \label{fig:CDO3}
\end{figure}

{\large \jlinputlisting[caption={The Julia code of generating and optimizing a CDO3-4G basis set for $\rm{H_2}$ using Quiqbox.}, label={lst:DObasis}]{DObasis_code.jl}}

According to the fermionic de Finetti theorem~\cite{krumnow2017fermionic}, a fermionic quantum state of a finite-size system that is invariant under permutations of modes can be well described by a non-interacting fermionic state produced by mean-field methods. It would be of interest to study whether a CDO-based basis set that enforces the orbital permutation invariance would improve the performance of Hartree--Fock methods without losing much information from discretizing the molecular electronic Hamiltonian. Hence, we can use Quiqbox to perform parameter optimization of CDO3-4G (see lines 39--46 of Listing \ref{lst:DObasis}). The result of this optimization is compared to typical AOs in TABLE~\ref{tbCDOHF}. CDO3-4G is able to produce lower Hartree--Fock energy for $\rm{H_2}$ with one-sixth as many orbitals as cc-pVDZ. This suggests the potential of using CDO-based basis sets to compress the information of the true ground state to a single Slater determinant state, as opposed to traditional AOs.

\begin{center}
    \begin{table}[ht]
    \setlength{\tabcolsep}{5pt}
    \renewcommand{\arraystretch}{1.25}
    \centering
    \begin{tabular}{lccrr}
    \toprule[1.5pt]
    Basis set   & $W$ & $N_{\rm{GTO}}$ & \multicolumn{1}{c}{$E_0\:({\rm{Ha}})$} & \begin{tabular}[c]{@{}r@{}}$E_0$ w.r.t.\\ cc-pVDZ\end{tabular} \\ \midrule[1.2pt]\midrule[1.2pt]

    STO-3G      & $2$   & $6$               & $-1.83100$                            & $+0.651\%$                                                                 \\
    6-31G       & $4$   & $6$               & $-1.84103$                            & $+0.107\%$                                                                 \\
    cc-pVDZ     & $18$  & $24$              & $-1.84300$                            & $0.000\%$                                                                     \\
    CDO3-4G      & $3$   & $24$              & $-1.84617$                            & $-0.172\%$                                                                 \\
    CDO3-nG      & $3$   & $\geq\!30$          & $-1.84700$                            & $-0.217\%$                                                                 \\
    \bottomrule[1.5pt]
    \end{tabular}
    \caption{The Hartree--Fock ground state energy approximations of ${\rm{H_2}}$ (with the bond length set to $1.4\;{\rm{a.u.}}$) using different basis sets. $N_{\rm{GTO}}$ represents the total number of GTO used for each basis set regardless of whether they share the same parameters. The CDO3 basis sets are optimized by Quiqbox to provide the lowest possible $E_0$. CDO3-nG is the asymptotic limit of the optimized CDO3 basis sets as the number of unique floating GTOs keeps increasing until $E_0$ converges.}
    \label{tbCDOHF}
    \end{table}
\end{center}

As we see through this example, with Quiqbox, one can quickly construct and optimize a customized basis set to meet specific research needs. In other words, Quiqbox aims to be a ``quick toolbox'' for basis set study.

\section{Conclusion and perspective} \label{sec:conclusion}

On the one hand, the ideas about Gaussian basis set design and optimization~\cite{tachikawa1999simultaneous, kresse1996efficiency, truhlar1998basis, hurley1988computation, frost1967floating, huber1979geometry} have been developing along the advance of classical ab initio calculation for electronic structure problems. On the other hand, the existing basis sets do not provide a smooth transition for NISQ devices to solve electronic structure problems practically. Nor do they permit proper benchmarks of VQE against existing classical algorithms. Tailoring new basis sets for specific scenarios becomes more in demand.

Aside from the need for a more customizable basis set design, we explained the significance of a more systematic study of basis sets beyond providing optimal numerical results for specific many-electron systems in Sec.~\ref{sec2}. Customizing and modifying basis sets to a higher degree across different systems will provide finer control of the electronic structure problem as a general computational problem and help study the effect of basis set discretization of electronic Hamiltonians.

To meet the abovementioned requirements, we proposed a highly customizable Gaussian-based basis set generation framework. It contains two main parts: basis set generation and parameter optimization. We introduced in Sec.~\ref{sec3} a type of more generalized Gaussian-based orbital named ``mixed-contracted Gaussian-type orbital (MCGTO)'' such that constructing basis sets beyond AOs using GTOs becomes possible. In Sec.~\ref{sec4}, we then introduced a variational parameter optimization procedure for the basis sets composed of MCGTOs, which can be implemented with general electronic structure algorithms.

Moreover, this framework is already realized by the open-source software package, Quiqbox.jl, which we have developed using the high-performance dynamic programming language Julia~\cite{bezanson2012julia, bezanson2018julia}. In Sec.~\ref{sec:performance}, we compared the performance of Quiqbox against other existing electronic structure packages. The results showcase Quiqbox's capability and efficiency of variational basis set optimization and provide helpful feedback regarding Quiqbox's potential improvement. Last but not least, we demonstrated using Quiqbox to construct basis sets of completely delocalized orbitals (CDOs). A CDO3-4G basis set consisting of three orbitals built from four unique floating GTOs can provide $0.2\%$ improvement for computing the Hartree--Fock energy of $\rm{H_2}$ compared to cc-pVDZ basis set with 18 orbitals (see TABLE~\ref{tbCDOHF}).

In the future, we would like to add more post-Hartree--Fock methods as the objective functions for the basis set optimization. Doing so will have multiple benefits. First, it will help optimize and tailor the basis set for more electronic structure algorithms. Secondly, we can start studying the correlation between basis set geometry (through basis function designs and parameters) and the performance of corresponding algorithms. For instance, given a ground state ansatz, can one design a basis set that provides low-error discretization of the electronic Hamiltonian and reduces the computational cost of finding the ground state? If such a basis set exists, then the parameter space characterized by the basis set will contain a global minimum closer to the true ground-state wavefunction and allows the optimizer to find it easily. Additionally, we can incorporate the basis set completeness optimization technique into the basis set construction and parameter optimization~\cite{truhlar1998basis, manninen2006systematic, lehtola_2014}, as it can help develop better automatic basis set generation procedures.

Finally, we would like to optimize Quiqbox's performance further. Notably, we want to improve Quiqbox's integral engine's efficiency. In this way, besides being a basis set generation and optimization packages, it can also be a strong option for doing ab initio electronic structure computation in quantum chemistry and quantum physics.

\section{Acknowledgements} \label{sec:ac}
This paper was supported by the ``Quantum Chemistry for Quantum Computers'' project sponsored by the U.S. Department of Energy under Award No. DESC0019374. The authors are thankful for Brent Harrison's proofreading and comments. JDW was supported by the US NSF grant PHYS-1820747 and NSF EPSCoR-1921199 and by the Office of Science, Office of Advanced Scientific Computing Research under programs Fundamental Algorithmic Research for Quantum Computing (FAR-QC) and the Optimization, Verification, and Engineered Reliability of Quantum Computers (OVER-QC) projects. JDW holds concurrent appointments at Dartmouth College and as an Amazon Visiting Academic. This paper describes work performed at Dartmouth College and is not associated with Amazon. WW would like to express his gratitude to his mother, father, and older sister. This paper would not have been possible without their utmost support. He would also like to thank his friends for their constant encouragement.

\bibliographystyle{unsrt}