Adaptive Basis Sets for Practical Quantum Computing

Electronic structure calculations on small systems such as H$_2$, H$_2$O, LiH, and BeH$_2$ with chemical accuracy are still a challenge for the current generation of the noisy intermediate-scale quantum (NISQ) devices. One of the reasons is that due to the device limitations, only minimal basis sets are commonly applied in quantum chemical calculations, which allow one to keep the number of qubits employed in the calculations at minimum. However, the use of minimal basis sets leads to very large errors in the computed molecular energies as well as potential energy surface shapes. One way to increase the accuracy of electronic structure calculations is through the development of small basis sets better suited for quantum computing. In this work, we show that the use of adaptive basis sets, in which exponents and contraction coefficients depend on molecular structure, provide an easy way to dramatically improve the accuracy of quantum chemical calculations without the need to increase the basis set size and thus the number of qubits utilized in quantum circuits. As a proof of principle, we optimize an adaptive minimal basis set for quantum computing calculations on an H$_2$ molecule, in which exponents and contraction coefficients depend on the H-H distance, and apply it to the generation of H$_2$ potential energy surface on IBM-Q quantum devices. The adaptive minimal basis set reaches the accuracy of the double-zeta basis sets, thus allowing one to perform double-zeta quality calculations on quantum devices without the need to utilize twice as many qubits in simulations. This approach can be extended to other molecular systems and larger basis sets in a straightforward manner.


INTRODUCTION
Quantum computing has seen a recent uptick in investment and interest, largely due to the promise that quantum computational tasks could be completed exponentially faster than their classical counterparts. 1 One of the most direct applications for quantum computing is in the area of ab initio quantum chemical calculations 2,3 Particularly, the full configuration interaction (full CI) method, which is a nondeterministic polynomial-time complete (NP-complete) problem, is a difficult problem for classical computers as the dimension of the CI space increases factorially with the size of the basis set. For example, a full CI calculation for an N2 molecule with ANO [4s3p1d] basis set needs 10 10 Slater determinants requiring 80 GB of disk memory. 4 With the development of quantum circuit algorithms, full CI may become solvable even for large molecules within a reasonable time and resources available. [5][6][7] However, the modern noisy intermediate-scale quantum (NISQ) devices are currently far too limited to accurately simulate large chemical systems. Many current quantum devices fall under this NISQ label, meaning they are incapable of large-scale error correction. Full error correction will only be achieved with at least several thousand quantum bits (qubits). 8 Given the current number of qubits available in a NISQ computer, quantum chemistry algorithms have been limited to simulating molecules with only a handful of atoms. Accurate modeling of small systems such as H2, H2O, LiH, and BeH2 is an important step towards the use quantum devices for electronic structure calculations of larger systems, and develops a better understanding of the NISQ systems capabilities. 9,10 While the structure optimizations on quantum computers are still challenging, quantum devices can be relatively easily employed to perform a series of single point energy calculations that allow one to construct a potential energy surface (PES) for a molecule of interest. 11,12 The PES can then be utilized to obtain useful information about the chemical system of interest, such as the molecular structure at a stable conformation (i.e., structures at PES minima), or a dissociation energy. The PES can be an important tool for modeling reaction dynamics or light-induced dynamics of chemical systems. [13][14][15][16][17][18] It can also be employed to perform vibrational analysis, allowing one to construct molecular partition functions and thermodynamic corrections to enthalpy and entropy of the system of interest. 19,20,21 Therefore, even with the current NISQ limitations, one can start collecting useful chemical information with relevance to the molecular structure and reactivity from electronic structure calculations on quantum devices.
Due to NISQ limitations, minimal basis sets, such as STO-3G, are most commonly employed for quantum chemical calculations as the minimal basis set requires the smallest number of qubits. 9,10,22 The STO-nG basis sets possess a single basis function for each core and valence shell of a particular atom, composed of n primitive Gaussian basis functions, that were contracted to reproduce the shape of a Slater-type orbital. 23,24 In the case of the STO-3G basis set, the s-orbital is described by a linear combination of three primitive Gaussian functions (assuming the atom is centered at the origin of the Cartesian system): (1) where di are the contraction coefficients, ai represent spatial spread of the wavefunction, and r is the distance from nucleus. The STO-3G basis set is available for atoms from hydrogen to xenon. 25 Traditionally, most basis sets are atom-centered and do not change with the molecular conformation. However, molecular orbital shapes and energies depend on relative atom positions. This is not a problem for large basis sets, but minimal basis sets do not have enough flexibility to describe the shapes of the molecular orbitals with the same accuracy at each molecular conformation. Thus, while the minimal basis sets are computationally efficient thanks to their small size, their accuracy can vary depending on the atom positions. On the other hand, the accuracy of complete or near-complete basis sets does not depend on molecular conformation, but they are computationally more expensive due to their larger size. 26,27 For example, the calculated CISD/STO-3G energy of H2 at 0.7 Å is 22.8 kcal/mol higher in energy than the energy calculated using a more complete basis set at CISD/aug-cc-pVQZ level, while the energy difference is 46.0 kcal/mol at 2.0 Å, 28 suggesting that the errors due to the basis set size are not uniform with respect to the molecular conformation (see Figure 1). Typically, one would increase the accuracy of the calculation by increasing the basis set size. For example, Figure 1 shows the significant improvement in both the potential energy surface (PES) energy and shape going from a minimal STO-3G and STO-6G basis sets to double-zeta basis sets (3-21G and 6-31G), and double-zeta polarized basis set (6-31G**). Another way to increase the accuracy of small basis sets, without increasing the number of basis functions, is to change (or adapt) the basis set based on the molecular conformation, creating a so-called "adaptive" basis set. Adaptive basis sets employ basis functions with exponents and contraction coefficients that depend on the molecular structure. Decades ago, adaptive basis sets were employed to provide a better representation of the molecular orbitals, using the idea of extracting atomic orbitals from molecular orbitals with population analysis, which is useful for calculating electron correlation with quantitative ab initio formulation and the interpretation of reactions. [29][30][31][32][33][34][35] Adaptive basis sets were also utilized to reduce the computational cost of SCF algorithms. [36][37][38][39][40][41] Recently, Schütt et al. employed machine-learning to construct adaptive basis sets for large scale molecular dynamics simulations using polarized atomic orbitals and quasi-atomic minimal basis representation. [42][43][44] Due to the limitations on the number of available qubits in the current NISQ devices, as well as the increase in noise leading to substantial errors when using a larger number of qubits, adaptive basis sets provide an alternative way to increase the computational accuracy. Adaptive basis sets could improve the basis set quality without increasing the number of basis functions and thus the number of qubits utilized in the calculations, avoiding the additional noise and error that comes from using larger quantum circuits with bigger depth.
The aim of this work is to explore the use of adaptive minimal basis sets in quantum chemical calculations on NISQ devices, taking the H2 molecule as our test system due to its small size and simplicity. We present an adaptive STO-3G basis set (labeled as A-STO-3G) that was constructed to produce the lowest energy at each H-H distance along the H2 potential energy surface. Basis set contraction coefficients and exponents were optimized at the Hartree-Fock (HF), configuration interaction with single and double excitations (CISD), and density functional (B3LYP) levels of theory, resulting in three adaptive sets A-STO-3G_HF, A-STO-3G_CISD, and A-STO-3G_B3LYP. 45-47 All A-STO-3G basis sets are more accurate than STO-3G when compared to the benchmark aug-cc-pVQZ basis set. They show similar accuracy as unmodified STO-6G around the H2 PES minimum, and a significant improvement over both STO-3G and STO-6G at long H-H distances. The PES of H2 at the UCCSD/A-STO-3G_CI level of theory was also generated with several IBM-Q quantum devices to check its validity. The ability to produce accurate molecular PESs with quantum computers would help to extract useful chemical properties, such as dissociation energies and vibrational frequencies, and to enable studies of the chemical reactivity.

Parameter optimization for adaptive basis sets
All calculations for basis set parameter optimization were carried out using the Gaussian 16 software package (Revision A.03). 48 Since HF and CI are variational methods, the basis set parameters can be optimized by minimizing the energy of the molecular system at each conformation. The pairs of ai and di parameters in Eq (1) are arbitrary, so we constrained a1 > a2 > a3 to obtain consistent results. Coefficient d3 was set equal to one to avoid redundant results (sets of different coefficients and exponents that yield the same energy), brought upon by the normalization condition of the optimized basis function: (2) where N is the normalization constant and yA-STO-3G is the optimized basis function.

Three-point optimizer
We have implemented the method of non-opportunistic coordinate search (see ref. 49 for a discussion of the classical coordinate search methods) in a script called three-point optimizer (TPO) for basis set parameter optimization with the Gaussian 16 software. A flowchart for the entire TPO process is found in Figure 2. TPO starts from the default STO-3G parameters at a given initial geometry. The optimizer then alters one coefficient by a user-defined constant, ± c, calculating single point energies at a set of parameters based on each alteration. If the center point is not the minimum of the three, the optimizer moves the center to the minimum among the three and repeats the step. The second or further iteration at this step requires only one more single-point energy calculation. If the center point between the three alterations is the minimum, then the optimizer moves on to the next coefficient until all five parameters (excluding d3 set to one) yield optimization to the lowest energy, which is tested by another cycle through each parameter. Once all the parameters are optimized at the ± c level, the process is repeated with a higher resolution (c/10) until it meets the criteria defined by the desired decimal place. At the next geometry, exponents and coefficients obtained in the previous step are employed as the starting point for optimization.
In addition to TPO, limited-memory Broyden-Fletcher-Goldfarb-Shanno algorithm with boundaries (L-BFGS-B), 50 and sequential least squares programming (SLSQP) optimizers 51 were tested. In the end, TPO optimized energies obtained were lower than energies from the two other optimizers but the computational cost of using TPO is much higher.

Figure 2.
Flowchart of the TPO process.

Sensitivity Analysis
When constructing a basis set, it may not be known a priori whether all basis set parameters are equally important to determine the energy of the system. One or more parameters could be an order of magnitude more or less influential than the rest. Mathematically, the study of parameter influence is known as sensitivity analysis and is an active area of research. [52][53][54][55][56][57][58] Broadly, sensitivity analysis can be divided into two categories: local sensitivity analysis (LSA) and global sensitivity analysis (GSA). 59 LSA approximates sensitivities locally by evaluating the gradient of the quantity of interest (in this case, ) at a particular parameter value. 59,60 Frequently, LSA will utilize many evaluation points scattered across the domain to overcome the limitations of locality.
Full GSA, on the other hand, constructs an expansion of the quantity of interest in terms of its conditional statistical moments. 61,62 For basis-set construction, LSA was performed on each basis set coefficient parameter to identify which parameter is the most influential.
For ease of explanation and implementation, in this work we use LSA with evaluation points of the parameters, . The evaluation points are sampled from uniform distributions, which we chose because we must avoid negative and because a uniform distribution is considered an unbiased prior distribution. 59 We also want to study parameter sensitivity as the H-H distance changes.
At a particular bond length , we compute the sensitivity of parameter as where p j is the j-th realization of the parameters. For the sampling intervals, we used which we determined by collecting the maximum and minimum optimized parameter values and then padding the relevant intervals even further. We need at least one to be relatively small to yield a wide wavefunction, so we chose without loss of generality. Furthermore, since the coefficients and are normalized, we set , again without loss of generality.

H2 potential energy surface construction
Our adaptive basis set was developed over a range of H-H distances between 0.5 to 6.0 Å. The grid consists of equally spaced points at intervals of 0.1 Å within the specified range. The TPO is used at each grid point to find the parameters that minimize the energy at HF, B3LYP, and CISD levels of theory. The optimized basis set at each point can then be used in single point energy calculations at CISD level to construct the potential energy surface.

RESULTS
In this section, we first discuss the results of the sensitivity analysis for the STO-3G parameters.
Then we present the values of each optimized parameter across the entire H-H distance range.
Next, we analyze the potential energy surfaces produced from these optimized parameters based on electronic structure calculations on classical computers. Finally, we report the results of the PES calculations using the optimized parameters on IBMQ quantum devices and simulators.

Sensitivity Analysis Results
We evaluated the sensitivity at the bond lengths 0.5, 0.6, …, 2.5 (in Å), equally spaced to every 0.1 Å, and the results are shown in Figure 3. The largest sensitivity corresponds to the exponent on the most diffuse Gaussian primitive, which is to be expected. However, the separation between the sensitivity of a3 and the other parameters is not several orders of magnitude, which is the typical heuristic of discarding parameters. 59 Additionally, the sensitivities have low variation as the H-H distance changes. Therefore, we optimize all six parameters, paying the most attention to a3 and the least attention to a1.

Parameters for the A-STO-3G basis sets
TPO optimization was done at the HF, B3LYP, and CISD levels of theory, resulting in the optimized A-STO-3G_HF, A-STO-3G_B3LYP, and A-STO-3G_CISD basis sets, respectively.
While the H2 molecule is small and CISD calculations are relatively inexpensive, such calculations may not be practical for larger systems, where the use of HF or DFT methods may be more accessible and cost-effective.
The TPO optimized A-STO-3G parameters are plotted individually in Figure 4. As seen in

Accuracy of A-STO-3G basis sets
The PES of H2 resulting from CISD/aug-cc-pVQZ calculations acts as the reference for determining if the use of the A-STO-3G basis sets has improved the accuracy. Neither the standard STO-3G, STO-6G, or the A-STO-3G variants match well with the reference around the minimum (H-H distance = 0.773 Å, see Figure 5). However, the energies obtained from the CISD/A-STO-3G calculations are significantly closer to the reference than the default STO-3G at long distances (see Figure 5, H-H distance > 0.9 Å). Importantly, the results obtained from the A-STO-3G_CISD basis set are very close to those obtained from the double-zeta quality basis sets (see Figure 6). results with double-zeta basis sets and the aug-cc-pVQZ reference basis set.
Unmodified STO-6G is the most accurate at calculating the energies around the minima (H-H distance = 0.7 Å), followed by the A-STO-3G_CISD basis set that is ~3 kcal/mol higher in energy.
For H-H distances longer than 1.0 Å, PES obtained using the A-STO-3G_CISD basis set demonstrates a considerable improvement in accuracy compared to the unmodified STO-3G and STO-6G; at H-H distance = 2.5 Å, the energy difference between the A-STO-3G_CISD basis set and aug-cc-PVQZ is only 7.27 kcal/mol, while the energy difference between the STO-6G and aug-cc-pVQZ basis sets is 39.85 kcal/mol. The PESs from HF-and B3LYP-optimized parameters also provide improved results, but are slightly less successful than the CISD-optimized STO-3G, leading to the errors at 2.5 Å of 25.46 kcal/mol and 12.29 kcal/mol, respectively. While it seems that it is preferable to employ CISD for the basis set optimization, parameters optimized at the B3LYP level provide similar accuracy to the CISD-optimized parameters, suggesting that basis   The use of the A-STO-3G basis sets results in shapes closer to the reference, as illustrated by the numerical derivative of each PES calculated displayed in Figure 8 (bottom panel). The STO-nG basis sets perform worse compared to the numerical derivative of the reference until ~1.7 Å. In contrast, the numerical derivative of A-STO-3G_CISD PES is closer to that of aug-cc-pVQZ, meaning its overall shapes are more similar.
PESs can also be used to extract useful quantities such as the molecular structure at minimum (Req) or dissociation energy (De). We have obtained the minimum for each PES by first constructing the sixth-order polynomial fits and then calculating its first derivative and setting it equal to zero. Furthermore, the approximate dissociation energy for each PES curve was obtained as the change in energy from the minimum to the energy from a H-H distance of 6.0 Å. Overall, the use of A-STO-3G basis sets results in PESs that are slightly less accurate in approximating H-H equilibrium bond distance than the STO-3G basis set, but significantly more accurate for determination of the dissociation energies and the overall shape of the H2 PES than the STO-nG basis sets. More importantly, the results obtained with the A-STO-3G basis sets reach the quality of the double-zeta basis sets that use twice as many basis functions. It is also worth noting that A-STO-3G basis set optimized at the CISD and B3LYP levels of theory provide very similar results, while the HF-optimized A-STO-3G is slightly less accurate.

PESs from quantum devices
The IBM Q Paris, Ourense, Rome, Bogota, and Santiago backends were used to collect data depending on their availability. The calculations were done at the UCCSD level of theory implemented in qiskit. Note that since hydrogen molecule contains only two electrons, UCCSD and CISD methods yield identical results. For the calculations on these quantum devices, parity mapping was used to reduce the noise by reducing the number of operated qubits to two. The SLSQP optimizer was used for VQE. The simultaneous perturbation stochastic approximation (SPSA) optimizer, 73 which is designed to work on noisy devices, showed more accurate results than SLSQP at longer distances. However, SPSA is about tenfold slower than SLSQP and becomes even slower for shorter distances requiring more optimization steps for energy to converge. Thus, we were not able to collect all data to construct the PESs with SPSA. With collected data from SLSQP, the trends in PESs between STO-3G and A-STO-3G basis sets are still evident.
We collected 395 data points for STO-3G basis set and 506 data points for A-STO-3G_CISD basis set from our runs on quantum devices. The PESs generated from the average and minima of the VQE results are presented in Figure 9. The results do not only deviate from the CISD surface, but also show a significant amount of noise even with 8000 shots per run. The major sources of error are likely from hardware such as qubit readout error, qubit instability, and CNOT (two-qubit gate) error. The magnitude of the error is not constant, but seem to increase at short and longer H2 distances (see Figure 10). Based on the PESs fitted from the average values, it seems that the error could be interpolated with a continuous function (Figure 9), suggesting that it could potentially be mitigated by statistical analysis. Some of these errors can also be improved through recently introduced error-mitigation methods such as noise amplification 74 and classical postprocessing. 75 Overall, the data collected on the IBM-Q devices confirm the improvements in the H2 PES at intermediate and long distances obtained through the use of the A-STO-3G basis sets that was seen on a classical computer, suggesting that development of the adaptive basis sets is a viable strategy to improve the accuracy of electronic structure calculations at quantum devices.

Conclusions
This work develops an adaptive basis set based on optimizing the STO-3G basis set for an H2 molecule at multiple geometries. A-STO-3G basis sets are more accurate in calculating the H2 PES than unmodified STO-3G. They are comparable to STO-6G around the PES minimum and lead to vast improvement at longer distances as the interaction between the two atoms diminishes. The A-STO-3G basis sets reach the accuracy of double-zeta basis sets (3-21G, 6-31G) that utilize twice as many basis functions. As a result, the PESs obtained with A-STO-3G basis sets possess a better shape and higher accuracy in reproducing the experimental dissociation energy. CISD/A-STO-3G_CISD calculations provide a PES closest in overall accuracy to the reference PES obtained at the CISD/aug-cc-VQZ level of theory, compared to the calculations with the STO-nG, A-STO-3G_HF or A-STO-3G_B3LYP basis sets.
The VQE calculations at the UCCSD/A-STO-3G_CISD level on quantum devices have a large error (> 0.1 hartree), but the average error may be corrected in post processing. Yet, the results obtained with A-STO-3G basis set still show a significant improvement compared to the STO-3G.
Further statistical analysis on extracting useful chemical parameters such as vibrational frequency from noisy data will be performed in the future. Overall, adaptive basis sets provide a viable strategy to increasing the accuracy of electronic structure calculations on quantum devices without increasing the demands on the number of qubits or operations utilized in the simulations.

Notes
This paper describes work completed while the author ZM was at North Carolina State University.
The authors declare no competing financial interest.