Quantum Solvers for Plane-Wave Hamiltonians: Abridging Virtual Spaces Through the Optimization of Pairwise Correlations

For many-body methods such as MCSCF and CASSCF, in which the number of one-electron orbitals is optimized and independent of the basis set used, there are no problems with using plane-wave basis sets. However, for methods currently used in quantum computing such as select configuration interaction (CI) and coupled cluster (CC) methods, it is necessary to have a virtual space that is able to capture a significant amount of electron-electron correlation in the system. The virtual orbitals in a pseudopotential plane-wave Hartree–Fock calculation, because of Coulomb repulsion, are often scattering states that interact very weakly with the filled orbitals. As a result, very little correlation energy is captured from them. The use of virtual spaces derived from the one-electron operators has also been tried, and while some correlations are captured, the amount is quite low. To overcome these limitations, we have been developing new classes of algorithms to define virtual spaces by optimizing orbitals from small pairwise CI Hamiltonians, which we term as correlation optimized virtual orbitals with the abbreviation COVOs. With these procedures, we have been able to derive virtual spaces, containing only a few orbitals, which are able to capture a significant amount of correlation. The focus in this manuscript is on using these derived basis sets to target full CI (FCI) quality results for H2 on near-term quantum computers. However, the initial results for this approach were promising. We were able to obtain good agreement with FCI/cc-pVTZ results for this system with just 4 virtual orbitals, using both FCI and quantum simulations. The quality of the results using COVOs suggests that it may be possible to use them in other many-body approaches, including coupled cluster and Møller–Plesset perturbation theories, and open up the door to many-body calculations for pseudopotential plane-wave basis set methods.


INTRODUCTION
Quantum chemistry is one of the first and most successful scientific applications of digital computers Boys, 1950;Parr and Mulliken, 1950;Hall, 1951;Roothaan, 1951;Boys et al., 1956;Nesbet, 1960;Allen and Karo, 1962;Nesbet, 1963;Pople et al., 1965;Kohn and Sham, 1965;Reeves, 1966;Pulay, 1969). This success has led to a large number of research, open-source (Dupuis et al., 1989;Stanton et al., 1992;Schmidt et al., 1993;Briggs et al., 1996;Challacombe, 2000;Gygi, 2008;Giannozzi et al., 2009;Deslippe et al., 2012;Hutter et al., 2014;Gonze et al., 2016;Harrison et al., 2016;Apra et al., 2020), and commercial codes (Kresse and Furthmüller, 1996;te Velde et al., 2001;Betteridge et al., 2003;Clark et al., 2005;Werner et al., 2012;Shao et al., 2015;Frisch et al., 2016;Neese, 2018) (for a larger list of quantum chemistry software, see (Wikipedia, The Free Encyclopedia, 2020)), which are used on a regular basis by tens of thousands of scientists, engineers, and students from a variety of scientific and engineering domains. With Moore's law as a backdrop (Moore, 2006), the cycle of new machines leading to new algorithms stimulated the field for many decades, and as a consequence, a large number of quantum chemistry methods were developed along with a variety of numerical methods to solve them. However, in recent decades, the maturity and success of these codes coupled with the imminent death of Moore's law (Dubash, 2005;Rotman, 2020) that made numerical software development much more difficult and less accessible to the average scientist have resulted in the field having priorities other than just new science, such as porting and optimizing these codes to the next generation of computers (Bylaska et al., 2017a, Bylaska et al., 2017bRichard et al., 2018;van Dam et al., 2020), standardization of methods (Crawford et al., 2017;Wilkins-Diehr and Crawford, 2018), and marketing (Goldbeck, 2017;Hocquet and Wieber, 2017).
With the advent of quantum computing, there is excitement again, and quantum chemists are beginning to rethink how they carry out quantum chemistry calculations, in particular very accurate and very expensive instances of systems containing strong electronelectron correlations. This is because it is anticipated that quantum computers with 50-100 qubits will be able to surpass classical digital computers for these types of calculations (Preskill, 2018). Quantum computing has thus emerged as an alternative avenue to the continuity of quantum chemistry in the long run (Wasielewski et al., 2020) but poses several challenges that demand careful consideration in order to eventually mature into a viable replacement for classical computers and large, highly parallelizable high-performance computing clusters.
Present quantum devices are plagued by short coherence times and vulnerability to environment interference, i.e., noise. Albeit quantum algorithms have been developed with proved exactness, such as quantum phase estimation, these are not a viable option in the present/near-term time frame. Therefore, it is desirable to limit the operation of quantum processors to a complementary concerted execution with classical counterparts, whereby each of these components is only in charge of those tasks for which it is more suitable. This has materialized into the variational quantum eigensolver (VQE) (Peruzzo et al., 2014) and other hybrid algorithms. Briefly, this class of algorithms strives to find the lowest eigenvalue of a given observable by assuming that the associated quantum state can be accurately represented by a trial wave function and whose parameters are varied according to the Rayleigh-Ritz method (variational principle), with these parameters being updated by the classical computer. The burden on the quantum processor can be further alleviated with strategies such as Trotterization, which in turn introduce other challenges (Evangelista et al., 2019;Grimsley et al., 2020) but can be successfully exploited in the construction of favorable ansatz, as long as a predefined form for the trial wave function is imposed. This is at the heart of the ADAPT-VQE (Grimsley et al., 2019).
Most high-level methods for strongly correlated systems in use today (e.g., full configuration interaction (CI), coupled cluster (CC) and Green's function (GF) approaches) are based on second-quantized Hamiltonians, which are written in terms of creation and annihilation operators for fermion orbitals. These methods are amenable to quantum computers because fermionic creation and annihilation operators can be readily mapped to qubits through the use of some established transformation, among which Jordan-Wigner (Jordan and Wigner, 1928), Bravyi-Kitaev (Bravyi and Kitaev, 2002), and binary codes (Steudtner and Wehner, 2018) stand out, where the number of qubits scales with the number of orbitals in the second-quantized Hamiltonian. In principle, converting the full many-body electronic Hamiltonian to a second-quantized form is exact and popular CC and GF approximations based on this form are very accurate. However, this conversion has a drawback in that it requires the introduction of a basis set, which, for computational cost reasons, needs to be small. Typically, these basis sets are composed of atomic-like orbitals generated with heuristics based on an atom calculation for each kind of atom in the system. An example of this type of basis set is the popular Dunning correlation consistent basis set (Dunning and Hay, 1977;Dunning, 1989) in which the atomic orbitals are optimized at the CISD (configuration interaction method with single and double excitations) level of theory. While the size of this basis set is small compared to other basis sets used in quantum chemistry, such as plane waves, it still needs to contain a large number of atomic orbitals to produce a truly accurate result.
Solving relevant chemistry problems analogously to what is classically done with MCSCF or FCI on near-term quantum computers that contain 10s to 100s of noisy qubits (Reiher et al., 2017), in which only limited numbers of operations can be performed, is a monumental challenge. One way to reduce the cost of these calculations is to develop new procedures for optimizing basis sets. In this manuscript, a new method is presented for generating a plane-wave derived correlation optimized orbital basis sets. These derived basis sets can also be used in other many-body approaches, including CC theory, and can easily be generalized to work with recently developed Filon's Integration Strategy for two-electron integrals in periodic systems . This method is different from other planewave derived optimized orbital basis sets (Shirley, 1996;Prendergast and Louie, 2009;Chen et al., 2011); in that, it is based on optimizing small select CI problems rather than fitting one-electron eigenvalue spectra and band structures.
The paper is organized as follows. In Section 2, a brief description of the second-quantized Hamiltonian and the double unitary CC downfolding method that can be used with the pseudopotential plane-wave method is given, followed by comparisons between restricted Hartree-Fock (RHF) calculations using plane-wave and Gaussian basis sets. Using this framework, CI calculations up to 20 virtual orbitals, generated from plane-wave Hartree-Fock and one-electron Hamiltonians, are shown for the H 2 molecule. The VQE quantum computing algorithms used in this work are described in Section 3. Section 4 presents a new class of algorithms for generating a virtual space in which the orbitals are generated by minimizing small pairwise CI Hamiltonians, and a complete set of equations for implementing these optimizations is given in Subsections 4.1-4.4. Using this new type of virtual space, CI calculations up to 18 virtual orbitals for the ground state energy curve of the H 2 molecule are presented in Section 5 followed by results using quantum computing simulations in Section 6, and lastly, the conclusions are given in Section 7.

PSEUDOPOTENTIAL PLANE-WAVE MANY-BODY HAMILTONIAN
The nonrelativistic electronic Schrödinger eigenvalue equation of quantum chemistry can be written as where H is the electronic structure Hamiltonian under the Born-Oppenheimer approximation and Ψ(x 1 , x 2 , . . . , x N e )〉 is the quantum mechanical wave function that is a function of the spatial and spin coordinates of the N e electrons, x i (r i , σ i ). When solving this equation, the Pauli exclusion principle constraint of particle exchange must be enforced, in which the wave function changes sign when the coordinates of two particles, x i and x j , are interchanged; i.e., For the Born-Oppenheimer Hamiltonian, the interaction between the electrons and nuclei is described by the proper potentials Ze |ri−RA| , which for plane-wave solvers can cause trouble with convergence because of the singular behavior at |r − R A |. A standard way to remove this issue in plane-wave calculations is to replace these singular potentials with pseudopotentials. By making this replacement, the Hamiltonian, H, in Eq. 1 can be written as where the first term is the kinetic energy operator, the second term contains the local and nonlocal pseudopotentials, V (A) local and V (A),lm NL represent the electron-ion interactions, and the last term is the electron-electron repulsion.
Instead of writing the many-electronic Hamiltonian in the traditional Schrödinger form, as in the equations above, it is more common today to write it in an alternative representation, known as the second-quantization form, defined using the creation, a † p , and annihilation, a p , operators. The secondquantized Hamiltonian is written as where ϕ p (x) represent the one-electron spin-orbital basis. A nice feature about this form of the Hamiltonian is that the antisymmetry of wavefunction requirement as given in Eq. 2 is automatically enforced through the standard fermionic anticommutation relations {a p , a † q } δ pq and {a p , a q } {a † p , a † q } 0. In this formulation, the choice of the one-electron spinorbital basis is nebulous and requires some care in its choosing in order to obtain accurate results with this type of Hamiltonian. Typically, in quantum chemistry, one uses the filled and virtual orbitals from a Hartree-Fock calculation. For methods that utilize linear combinations of atomic orbitals (LCAO) as the basis, the size of the basis set and subsequently generated Hartree-Fock orbitals is fairly small. However, for plane-wave solvers and other grid-based solvers, the size of the basis set is very large and the number of the one-and two-electron integrals in Eq. 4 will become prohibitive if all possible Hartree-Fock orbitals are used.
One approach to this problem is to only include virtual orbital up to a certain energy threshold, and another related approach is to use the plane-wave derived optimized orbital basis set, e.g., the Shirley approach. While the number of these orbitals needed to accurately describe eigenvalue spectra over a range of ∼100 eV is significantly smaller than the number of plane waves, it is still significantly larger than the number of orbitals generated by an LCAO method. The reason for this is that the virtual orbitals in a plane-wave Hartree-Fock calculation, because of Coulomb repulsion, are often unbound scattering states that interact very weakly with the filled orbitals. As a result, very little correlation energy is captured from them. In contrast, LCAO basis methods can only describe bound states, and hence, Hartree-Fock calculations on this basis do not generate these types of scattering states.

Many-Body Downfolding Techniques
One technique for reducing the dimensionality of large planewave calculations is to construct effective Hamiltonians that capture correlation effects of the full calculation in manageable active spaces. That way, all many-body effects are retained as opposed to simply truncating the orbital space. In Bauman et al.
where σ int and σ ext are the general type anti-Hermitian operators defined by amplitudes defining action within and outside of the predefined active space, respectively; i.e., the amplitudes defining the σ ext operator must carry at least one inactive spin-orbital index whereas all amplitudes defining the σ int operator carry active spin-orbital indices only. In Eq. 5, |Φ〉 designates properly chosen reference function (usually chosen as a Hartree-Fock (HF) Slater determinant). The exactness of the expansion 5 has been recently discussed in Kowalski and Bauman (2020), where it was also shown that the standard UCC expansions can provide a basic approximation of the exact σ int and σ ext operators, i.e., where T int and T ext are single-reference-type internal and external cluster amplitudes (in the sense defined above). Using ansatz in Eq. 5 we have shown (Kowalski and Bauman, 2020) that the exact energy of the systems can be reproduced by the diagonalization of the effective (or downfolded) Hamiltonian in the corresponding active space: where In Eq. 11, P and Q int are the projection operators onto the reference function and all active-space excited Slater determinants (with respect to |Φ〉).
We will discuss the utility of the downfolding techniques in the next section for the ground state calculations of H 2 . This is just one of the two approaches presented in this paper for reducing the dimensionality of the quantum problem ( Figure 1).

Results for the 1 Σ + g Ground State of H 2 Using Virtual Space From Hartree-Fock and One-Electron Hamiltonians
The NWChem program package (Kendall et al., 2000;Valiev et al., 2010;Bylaska et al., 2011;Bylaska, 2017;Apra et al., 2020) was used for all calculations in this study, except for the FCI calculations, which used the TINYMRCC suite by Jiří FIGURE 1 | Schematic representation of the dimensionality reduction algorithms considered in this paper: (1) discretization of the many-body problem by employing efficient single-particle basis sets (in this paper, we consider correlation optimized virtual orbitals (COVOs)) and (2) downfolding techniques based on the double unitary coupled cluster (DUCC) formalism Kowalski and Bauman, 2020); in this step, the many-body problem is rerepresented in a subspace of entire Hilbert space.
FIGURE 2 | The ground state energy curves for H 2 with RHF, CCSD, and DUCC/QDK methods using plane-wave and LCAO Gaussian basis sets. It should be noted that for the two-electron H 2 molecule CCSD gives the same answer as FCI.
Frontiers in Chemistry | www.frontiersin.org March 2021 | Volume 9 | Article 603019 Pittner. The plane-wave calculations used a simple cubic box with L 26a 0 and cutoff energy of 100 Ry. The valence electron interactions with the atomic core are approximated with generalized norm-conserving Hamann (Hamann, 1989) pseudopotentials modified to the separable form suggested by Kleinman and Bylander (Kleinman and Bylander, 1982). The pseudopotentials used in this study were constructed using the following core radii: H: rcs 0.8 a.u. and rcp 0.8 a.u.; Be: rcs 1.219 a.u., and rcp 1.219 a.u. (vide infra). The RHF and coupled cluster singles and doubles (CCSD) LCAO calculations used the Dunning cc-pVTZ Gaussian basis set.
As can be seen in Figure 2, the RHF ground state energy curve of the H 2 molecule using plane-wave and LCAO Gaussian basis sets gives nearly identical results. However, when we performed plane-wave FCI calculations (not shown) for this system using up to 20 RHF virtual orbitals, the amount of correlation energy calculated was nearly zero (<1.0 e-4 Hartree). This result was not surprising since most of the virtual states were scattering states as shown in Figure 3.
Instead of using virtual states of the RHF Hamiltonian, virtual states were also generated using the 1-electron part of the RHF Hamiltonian, H 1 (i.e., just the kinetic energy and pseudopotential terms). As shown in Figure 3, the H 1 Hamiltonian generated virtual orbitals that were bound and looked like the virtual orbitals generated in an LCAO RHF calculation. Using these H 1 generated orbitals, we performed plane-wave CI calculations using 19 of these virtual orbitals. As seen in Figure 2 and Table 1, a significant improvement was seen using these orbitals as they were able to capture a nonzero amount of the correlation energy; however, it was still significantly less than that found in LCAO calculations. In addition, results using the quantum phase estimation (QPE) algorithm in the Microsoft QDK package (Svore et al., 2018;Low et al., 2019) in which the number of orbitals was reduced to 4 and 6 orbitals using the DUCC method are shown. These results showed that the DUCC QDK QPE method produces total energies that are within a few milli-Hartrees of the 20 orbital FCI result with only 4 or 6 orbitals.

VARIATIONAL QUANTUM EIGENSOLVER METHODS
VQE is a method to find the quantum state that minimizes a cost function defined in operator form (Peruzzo et al., 2014;O'Malley et al., 2016). This is a hybrid computational approach in which the preparation of the quantum circuit is tuned using feedback from classical evaluations of the cost function. Reduction of a given problem to minimization, such as solving for the ground state energy (lowest energy eigenvalue) of a molecular Hamiltonian, may then rely on the variational principle to affirm that only the true ground state could satisfy the minimum energy (Kandala et al., 2017;McCaskey et al., 2019).
Formally, we may consider the problem of solving for the ground state energy, E g , as Notice that the isovalues of RHF LUMO 1 are very close to zero, which indicates that it is a scattering state. The orientation of H 2 is rotated by 90°in the bound LUMO 3 relative to the bound LUMO 2.
The second equality is pertinent to the current context as it makes explicit the fact that 1) H P , the Hamiltonian in terms of strings of Pauli operators, relates to H through some transformation that maps fermionic creation and annihilation operators to qubits operators, and 2) the trial wave function emerges from the action of the parameterized unitary operator U( θ → ) which builds entanglement, usually starting from an unentangled wave function, such as Hartree-Fock. For practical considerations, H is transformed into H P with the Jordan-Wigner transformation (Jordan and Wigner, 1928), but alternatives have been reported in the literature (Bravyi and Kitaev, 2002;Setia and Whitfield, 2018;Steudtner and Wehner, 2018). The classical search for the quantum state that minimizes the energy represents the conventional computing task, while evaluation of the expectation value is performed using the quantum computer. In particular, the quantum state is prepared by executing a quantum circuit, which is expressed formally as a series of unitary operators acting on a well-defined initial state. The action of a specific sequence of unitaries is to prepare a given state that is subject to the measurements necessary to recover the desired expectation value.
In practice, the quantum state that minimizes the energy is unknown, and, therefore, a search over possible unitaries is necessary to find the form that minimizes the energy. This variational approach to circuit synthesis underlies the VQE method and an essential choice is the selection of a quantum circuit ansatz which defines the range of unitaries that may be formed to prepare and evaluate a quantum state. For electronic structure calculations, seemingly randomized unitaries may offer advantages for efficient circuit construction, but they lack much of the intuition available from theoretical chemistry (Kandala et al., 2017). Rather, ansatz circuits derived from unitary coupled cluster theory offer a convenient connection to the expected unitary forms of the minimal quantum state (Romero et al., 2017).
VQE has been applied previously to recover the electronic energy from the Hamiltonian presented in Eq. 4. The literature provides several examples of usage of VQE for problems of chemical interest, in terms of both simulation and implementation on actual quantum hardware. Given the current limitations faced by present quantum computers, these instances are usually accompanied by strategies that reduce the effective Hilbert space, thus leading to a decrease in the computational expense, such as the use of active spaces and natural orbitals (Verma et al., 2020), as well as downfolding techniques introduced earlier. Another route is to modify the form of the ansatz; an example of this alternative would be the so-called Trotterization, which can be used in conjunction with Hilbert space-reducing techniques.
Recently, the principle of VQE was extended to use ansatz circuits that are tailored to computational chemistry applications and specifically the unitary coupled cluster singles and doubles (UCCSD) ansatz state. Adaptive ansatz construction is attractive because it obeys the underlying complexity of the electronic structure in question, whereas a predefined form for the trial wave function in Eq. 13 may fall short of the flexibility necessary for intricate problems. The prime example of this class of algorithms is the ADAPT-VQE, which iteratively assembles a circuit according to the expected energy gain signaled by the gradient with respect to the variational parameters.
An important consideration in the performance of both VQE and ADAPT-VQE is the depth of the ansatz circuit and the time required to construct the optimal variational circuit. For electronic structures dominated by weak correlation, ADAPT-VQE tends to be very economical, adding only operators that make a meaningful contribution toward the lowest eigenvalue in the spectrum of the Hamiltonian in Eq. 14. On the other hand, the usual UCCSD, by virtue of being defined ahead of time, may contain operators with little impact on the energy, but the classical optimizer will still need to perform a number of calls to the cost function in order to find their best values. Also, the gates originating from these operators, even if they are deemed unimportant because of a small associated parameter, will nevertheless be present in the circuit, adding to its depth. If high accuracy is sought, then ADAPT-VQE may require an ansatz comprised of a large number of operators, which in turn adds to the depth of the underlying circuit. More operators also mean more variational parameters, leading to an onerous optimization process. A more detailed analysis of this trade-off can be found in Grimsley et al. (2019).

ALGORITHM FOR DEFINING A VIRTUAL SPACE WITH A SMALL CI HAMILTONIAN
In this section, we present a downfolding method to define virtual orbitals for expanding the second-quantized Hamiltonian given in Eq. 4. These new types of orbitals are able to capture significantly more correlation energy than the virtual orbitals coming from Hartree-Fock and one-electron Hamiltonians tested in Section 2.2. The basis of this method is Frontiers in Chemistry | www.frontiersin.org March 2021 | Volume 9 | Article 603019 6 to define a set of virtual orbitals, {ψ (n) e (r)} with n 1 . . . N virtual , which we call correlation optimized virtual orbitals or COVOs for short, by optimizing a small select CI Hamiltonian with respect to a single virtual orbital and then the next virtual orbitals in sequence, subject to them being orthonormal to the filled and previously computed virtual orbitals. The algorithm to calculate these new types of orbitals can be formulated as follows: 1. Set n 1. 2. Using the ground state one-electron orbital, ψ g (r) (or ground state orbitals for many-electron systems), and the virtual orbital to be optimized, ψ (n) e (r), generate a CI matrix. 3. Calculate the select CI expansion coefficients by diagonalizing the CI matrix. 4. Using the CI coefficients associated with the lowest eigenvalue, calculate the gradient with respect to the ψ (n) e (r) and then update with a conjugate-gradient or similar method while making sure that ψ (n) e (r) is normalized and orthogonal to ψ g (r) and ψ (m) e (r) for m 1, . . . , n − 1. 5. If the gradient is small, then n n + 1. 6. If n ≤ N virtual , go to step 2; otherwise, finish.
In the case of the H 2 molecule, a small CI wave function for the 2 electron system composed of 2 one-electron orbitals, ψ g (r) and ψ (n) e (r), can be written as a linear combination of 6 determinant wave functions, or just 3 determinant wave functions for just singlet (or triplet) states, , ψ e (r) + . . . .
Using this small CI ansatz, the energies of the system can be obtained by diagonalizing the following eigenvalue equation: Note that the overlap matrix, S, is the identity matrix for orthonormal ψ g and ψ e . The variation with respect to ψ e (r) can be simply obtained using the following formula: It should be noted that the above formulas can be generalized to work beyond two-electron systems by using corresponding orbitals techniques (King et al., 1967;Bylaska and Rosso, 2018). The next two Subsections 4.1-4.4 provide formulas that can be used to generate the matrix elements in Eq. 14 and the gradients with respect to ψ p e (r) in Eq. 15. We also note that the COVO approach proposed in this work is similar in spirit to the optimized virtual orbital space (OVOS) approach developed over 30 years ago by Adamowicz and Bartlett (Adamowicz and Bartlett, 1987;Adamowicz et al., 1988). The differences in our approach compared to this previous work is that the variational space used by COVOs is significantly bigger because plane-wave basis sets are used instead of LCAO Gaussian basis sets and that a second-order Hylleraas functional (Hylleraas, 1928;Hylleraas, 1929;Hylleraas, 1930;Hylleraas, 1964;Koga, 1992) was used to describe the correlation in the OVOS procedure rather than a small CI Hamiltonian. Other differences with the COVOs approach are that the orbitals are optimized one at a time and the cost to generate them is similar to generating regular RHF virtual orbitals (just 4 to 9 times more expensive relative to RHF). Moreover, the resulting electronic gradient is non-Hermitian, which in addition to requiring more involved optimizers can result in extended energy plateaus that occur during the initial stages of the geodesic line searches in a conjugate-gradient or quasi-Newton optimization method.

Two-Electron Orbitals for a Two-State Hamiltonian
For the two-state system, there are six two-electron wave functions, two of which are singlet, two are triplet, and two contain a mixture of singlet and triplet character. These wave functions can be written as − ψ e (r 1 )β(s 1 )ψ g (r 2 )α(s 2 ) ⎤ ⎥ ⎦ , Note that Ψ a and Ψ b cannot be written as a product of a spatial wave function times a spin function. Moreover, these functions are not eigenfunctions of the spin operators S 2 and S z , and as a result, these determinants contain both singlet and triplet components. However, if we take linear combinations of them, we can get two new wave functions that are separable in spatial and spin functions and at the same time being eigenfunctions of S 2 and S z ,

Matrix Elements From the Two-Electron Operators
The H 2 operator for H 2 molecule is , .

RESULTS FOR 1 Σ + g GROUND STATE OF THE H 2 MOLECULE USING CORRELATION OPTIMIZED VIRTUAL ORBITALS (COVOS)
The results for PW FCI calculations of H 2 with 1, 4, 8, 12, and 18 COVOs are shown in Figure 4 and Table 2. The average difference error for the 1, 4, 8, and 12 COVOs calculations from the 18 COVOs calculation is 11.8 kcal/mol, 1.4 kcal/mol, 0.9 kcal/mol, and 0.3 kcal/ mol, respectively. While the error is significant for 1 virtual, the difference is quite small by 4 virtual orbitals, and the error steadily decreases as the number of virtual orbitals increases. The error seen in the 4 optimized virtual orbitals' calculations is similar to the 1.6 kcal/ mol error seen in the DUCC calculations for the 19 H 1 virtual orbitals' calculations in Section 2.2. Another measure of the error is the extensivity error. The energy for large R should be the same as the energy of twice the energy of an isolated H atom. For the pseudopotential plane-wave method being used, the energy of 2 H atoms is −0.997765 Hartrees (E(1H) −0.498825 Hartrees). This difference at R 7 Å is found to be 11.6, 1.2, 1.1, 0.5, and 0.4 kcal/mol for 1, 4, 8, 12, and 18 optimized virtual orbital calculations, respectively. With only 4 optimized virtual orbitals, the correlation energy at the minimum was found to be −0.035 Hartrees, which is comparable to the −0.039 Hartrees found with CCSD/cc-pVTZ. The correlation energy decreases to −0.037, then −0.038, and finally −0.039 Hartrees as the number of optimized virtual orbitals increases to 8, 12, and 18 orbitals, respectively. These results showed that by 18 COVOs the same amount of correlation energy was recovered as with the cc-pVTZ LCAO basis set calculation. The results also showed that there was systematic convergence toward the benchmark LCAO result as the number of COVOs was increased, and with just 4 COVOs, a significant portion of the correlation energy was recovered.
QUANTUM SIMULATIONS OF THE 1 Σ + g GROUND STATE OF THE H 2 MOLECULE USING COVOS The previous section provides indisputable evidence for the performance of the proposed virtual orbitals for correlation energy recovery. Besides the possible ramifications in quantum chemistry carried out with classical computers, one immediate application is in the realm of quantum simulations. Because the present quantum hardware has not fully matured, hybrid algorithms that leverage classical resources and restrict the workload delegated to quantum computers, namely, state preparation and measurements of highly entangled states, are essential to meaningful quantum computations. The COVOs meet this requirement by decreasing the dimensionality of the problem, i.e., by enabling simulations with fewer qubits.
In order to probe the performance of COVOs in quantum simulations, we use the Hamiltonian with 4 COVOs and simulate the 1 Σ + g ground state of H 2 in the same bond distances shown in Table 2 and Figure 4. The ansatz circuit for the simulations is generated according to the ADAPT-VQE algorithm as implemented in the XACC (McCaskey et al., 2018b, McCaskey et al., 2020 framework for hybrid quantum computing using the tensor network quantum virtual machine (TNQVM) as the noiseless simulator backend (McCaskey et al., 2018a). In the present study, the ADAPT-VQE cycle is repeated until the norm of the gradient vector falls below 1e-2 and we use an operator pool containing all spin-adapted single and double excitation operators (one-and two-body rotations). A detailed account of ADAPT-VQE is exposed elsewhere (Grimsley et al., 2019). Optimization of the parameterized circuit is conducted with the COBYLA (Powell, 1994) optimizer as implemented in the NLOpt package (Johnson, 2007). Results for the simulated potential energy curve are plotted in Figure 5. It is evident from Figure 5 that ADAPT-VQE can generate a circuit capable of reproducing the FCI results in the current active space. These simulations deliver a smooth, continuous potential energy curve that tracks the FCI values strikingly well. The deviations from the corresponding FCI energies are all found below 1e-4 Hartree. This means that not only do these simulations deliver results that are well below the conventional chemical accuracy mark, but also more importantly in the current context is that this error is inconsequential compared to the effect of noise in case of deployment on actual quantum hardware.
It is remarkable that the results degrade little throughout the energy scan, which attests to the aptness and flexibility of ADAPT-VQE in determining an ansatz according to the complexity of the underlying electronic structure. The ansatz in the vicinity of the equilibrium bond length 0.5-1.0 Å is comprised solely by pair excitations as would be expected given a restricted HF reference, which means no determinant obtained via one-body rotations can lower the energy below that of HF. As we approach the Coulson-Fischer point (Coulson and Fischer, 1949), single excitations start to become part of the ansatz, which signals the inadequacy of a restricted reference wave function and that inclusion of these operators enables the ansatz to remain in the 1 Σ + g potential energy curve, which means that this flexibility may come at the expense of deeper circuits. Because one-qubit gates tend to be executed in a short timescale and are fairly insensitive to noise, we can use the number of CNOTs present in the circuit as indicative of the complexity in its implementation, which we provide in Figure 6, showing that the ansatzes generated from ADAPT-VQE are much more affordable than those obtained by ordinary UCCSD VQE simulations.
Along these lines, once the operator composition of the ansatz is defined, by virtue of introducing more parameters, we are likely to experience a more arduous optimization of the corresponding parameterized gates. This has a compound effect with the circuit depth since more measurements are needed, each of which requires the circuit to be implemented and measurements to take place. Figure 7 gives a profile of the optimization performance along the potential energy scan.
It should come as no surprise that the optimization is more difficult in the regime of stronger correlation. This region also demands a more complex ansatz, as the top plot in Figure 7 shows that only in this vicinity (1.5-3.5 Å) we observe ansatzes with more than four operators. Interestingly, the number of objective function calls does not show large deviations for ansatzes with 1-3 operators, regardless of where they are found in the potential energy curve, which is further corroborated by the relatively small error bars in the corresponding columns of the bottom plot. This observation does not hold as more parameters/operators are introduced in the ansatz in order to accommodate a more complex electronic structure. Thus, with four parameters, not only are more calls to the objective function needed, but also there is a more pronounced standard deviation. Ansatzes with five or more operators can only be found in the (1.5-3.5 Å), as we can see that the calls to the objective function coming from them dominate the overall number of optimization cycles. Due to the scarce occurrence of these ansatzes in the current energy scan, the corresponding statistical information that can be derived from these instances is not as reliable. All in all, this plot is valuable in lending additional insight into the resources required to perform these simulations. It is important to mention  that, for every new ansatz, the variational parameters are initialized at zero. Alternatively, the parameters corresponding to the previously optimizer ansatzes could be initialized at their optimal values and the new parameter would be introduced in the ansatz, which would accelerate convergence. Moreover, the convergence profile likely displays pronounced dependence on the chosen optimizer, which is not pursued here.

CONCLUSION
In summary, we have developed a new approach for defining virtual spaces with a pseudopotential plane-wave code for use in many-body methods described by second-quantized Hamiltonians. The method is based on optimizing the virtual orbitals to minimize a small select CI Hamiltonian (i.e., COVOs) that contains configurations containing filled RHF orbitals and the one virtual orbital to be optimized. Subsequent virtual orbitals are optimized in the same way, but with the added constraint of being orthogonal to the previously calculated filled and virtual orbitals. The method was applied to the simple, but nontrivial, H 2 molecule. As summarized in Figure 8, these new types of virtual orbitals were significantly better at capturing correlation in plane-wave calculations than from virtual spaces from Hartree-Fock and one-electron Hamiltonian, and moreover, we were able to obtain good agreement with Gaussian cc-pVTZ basis set results FIGURE 8 | Summary of various plane-wave and Gaussian basis set RHF and FCI calculations for the 1 Σ + g ground state of H 2 molecule. FCI calculations (not shown) using up to 20 RHF virtual orbitals produced only a negligible amount of correlation energy (<1.0 e-4 Hartree, i.e., visually the same as RHF results). It should be noted that for the two-electron H 2 molecule CCSD gives the same answer as FCI.
FIGURE 9 | Potential energy curves in kcal/mol for the Be 2 dimer using plane-wave Hartree-Fock and FCI with 5 COVOs calculations.
with just 4 virtual orbitals for the H 2 molecule. Subsequent calculations showed that the correlation energy converged steadily as more virtual orbitals were included in the calculation. With 18 virtual orbitals, the correlation energies were found to be converged to less than 0.5 kcal/mol. The robustness of the proposed basis sets is corroborated by its ready applicability to quantum simulations, which in the case of ADAPT-VQE show remarkable agreement with the classical, exact diagonalization result (FCI) in the same basis set (4 COVOs).
Because this study is focused on how one might carry out planewave CI calculations on near-term quantum computers in the next few years, we have only shown results for the H 2 dimer. However, we are optimistic that these correlation optimized virtual orbitals open up the door to many-body calculations using pseudopotential planewave calculations, including coupled cluster, Møller-Plesset, and Green's function theories as well as other FCI-approaching methods for quantum computers. We hope in future studies to more thoroughly test the effectiveness of the COVOs procedure on larger and more complicated molecules and materials. To lend credence to this assertion, we show the promising results for Be 2 dimer with a small number of COVOs in Figure 9. Also as shown in Figure 10, the shapes of the COVOs end up being similar to what is found for the virtual orbitals from LCAO calculations. This suggests that new classes of LCAO basis sets might be able to be generated using a simple rotation of the filled orbitals and COVOs. Future work will focus on using this approach on larger molecular and periodic systems. With the validation granted by our quantum simulations, further studies are called for, including the further reduction of dimension by employing active-space DUCC downfolded Hamiltonians, OVOS, and natural orbitals, as well as work in conjunction with VQE methods.

DATA AVAILABILITY STATEMENT
The datasets presented in this study can be found in online repositories. The names of the repository/repositories and accession number(s) can be found below: https://github.com/ ebylaska/PWH2-Data/archive/master.zip.

AUTHOR CONTRIBUTIONS
All authors contributed to the writing and editing of the manuscript. EB and DS developed and implemented in the NWChem the algorithm for defining a virtual space with small CI Hamiltonians. KK wrote a part of the FIGURE 10 | The 1 filled RHF orbital and 18 COVOs for the H 2 molecule and 2 filled RHF orbitals and 5 COVOs for the Be 2 molecule are shown in the left and right panels, respectively. The orbitals are displayed in the order of increasing orbital energy from left to right and bottom to top. The distance between two atoms at which the energy achieves its minimum is 0.7 Å for H 2 and 2.8 Å for Be 2 . The positive and negative isosurfaces are colored in blue and orange, respectively. Notice that, for H 2 , the shapes of COVOs 1, 3, and 4 are the same as those of LUMOs 1, 2, and 3 from the H 1 Hamiltonian, respectively, while COVO 2 has a similar shape to LUMO 1 from the straight HF in Figure 3 although the former is a bound state while the latter is essentially a scattering state.