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 H₂ 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 (Mulliken et al., 1941; Mulliken and Rieke, 1941; 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., 2017b; Richard 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 electron-electron 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 (Bylaska et al., 2020). This method is different from other plane-wave 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₂ 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₂ 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

$H|\Psi \left(x_1,x_2,\dots ,x_{N_e}\right)\rangle =E|\Psi \left(x_1,x_2,\dots ,x_{N_e}\right)\rangle ,(1)$

where H is the electronic structure Hamiltonian under the Born–Oppenheimer approximation and $|\Psi \left(x_1,x_2,\dots ,x_{N_e}\right)\rangle$ is the quantum mechanical wave function that is a function of the spatial and spin coordinates of the $N_e$ electrons, $x_i=\left(r_i,{\sigma }_i\right)$. 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.,

$\begin{array}{l}|\Psi \left(x_1,x_2,\dots x_i,\dots x_j,\dots ,x_{N_e}\right)\rangle \\ =-|\Psi \left(x_1,x_2,\dots x_j,\dots x_i,\dots ,x_{N_e}\right)\rangle .\end{array}(2)$

For the Born–Oppenheimer Hamiltonian, the interaction between the electrons and nuclei is described by the proper potentials $\frac{Ze}{\left|r_i-R_A\right|}$, which for plane-wave solvers can cause trouble with convergence because of the singular behavior at $\left|r-R_A\right|$. 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

$\begin{array}{c}H=-\frac{1}{2}{\displaystyle \sum _{i=1}^{N_e}{\nabla }_i^2}\\ +{\displaystyle \sum _{i=1}^{N_e}{\displaystyle \sum _{A=1}^{N_A}}}\left(V_{local}^{\left(A\right)}\left(\left|r_i-R_A\right|\right)+{\displaystyle \sum _{lm}{\widehat{V}}_{NL}^{\left(A\right),lm}}\right)\\ +{\displaystyle \sum _{i=1}^{N_e}{\displaystyle \sum _{j>i}^{N_e}\frac{1}{\left|r_i-r_j\right|}}},\end{array}(3)$

where the first term is the kinetic energy operator, the second term contains the local and nonlocal pseudopotentials, $V_{local}^{\left(A\right)}$ and ${\widehat{V}}_{NL}^{\left(A\right),lm}$ 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 second-quantized Hamiltonian is written as

$\begin{array}{c}H={\displaystyle \sum _{p=1}^{N_{basis}}{\displaystyle \sum _{q=1}^{N_{basis}}{h}_{pq}{a}_p^{†}{a}_q+\frac{1}{2}{\displaystyle \sum _{pqrs}{h}_{pqrs}{a}_p^{†}{a}_r^{†}{a}_s{a}_q}}},\\ h_{pq}={\displaystyle \int dx{\varphi }_p^{\ast }\left(x\right)\left(-\frac{1}{2}{\nabla }^2\right){\varphi }_q\left(x\right)}\\ +{\displaystyle \int dx{\varphi }_p^{\ast }\left(x\right)\left[{\displaystyle \sum }_{A=1}^{N_A}\left(V_{local}^{\left(A\right)}\left(\left|r-R_A\right|\right)+{\displaystyle \sum }_{lm}{\widehat{V}}_{NL}^{\left(A\right),lm}\right)\right]{\varphi }_q\left(x\right)},\\ h_{pqrs}={\displaystyle \int d{x}_{1}d{x}_2{\varphi }_p^{\ast }\left(x_1\right){\varphi }_r^{\ast }\left(x_2\right)\frac{1}{\left|r_1-r_2\right|}{\varphi }_s\left(x_2\right){\varphi }_q\left(x_1\right)},\end{array}(4)$

where ${\varphi }_p\left(x\right)$ 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 $\left\{a_p,a_q^{†}\right\}={\delta }_{pq}$ and $\left\{a_p,a_q\right\}=\left\{a_p^{†},a_q^{†}\right\}=0$.

In this formulation, the choice of the one-electron spin-orbital 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 plane-wave 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. (2019), we introduced downfolding techniques, which utilize the double unitary CC (DUCC) ansatz for exact ground state wave function $|\Psi \rangle$,

$|\Psi \rangle =e^{{\sigma }_{ext}}{e}^{{\sigma }_{\mathrm{int}}}|\Phi \rangle ,(5)$

where ${\sigma }_{\mathrm{int}}$ and ${\sigma }_{ext}$ are the general type anti-Hermitian operators

${\sigma }_{\mathrm{int}}^{†}=-{\sigma }_{\mathrm{int}},(6)$

${\sigma }_{ext}^{†}=-{\sigma }_{ext},(7)$

defined by amplitudes defining action within and outside of the predefined active space, respectively; i.e., the amplitudes defining the ${\sigma }_{ext}$ operator must carry at least one inactive spin-orbital index whereas all amplitudes defining the ${\sigma }_{\text{int}}$ operator carry active spin-orbital indices only. In Eq. 5, $|\Phi \rangle$ 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 ${\sigma }_{\text{int}}$ and ${\sigma }_{ext}$ operators, i.e.,

${\sigma }_{\mathrm{int}}\simeq T_{\mathrm{int}}-T_{\mathrm{int}}^{†},(8)$

${\sigma }_{ext}\simeq T_{ext}-T_{ext}^{†},(9)$

where $T_{\text{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 ${\overline{H}}_{eff}^{\left(DUCC\right)}$ in the corresponding active space:

${\overline{H}}_{eff}^{\left(DUCC\right)}{e}^{{\sigma }_{\mathrm{int}}}|\Phi \rangle =E{e}^{{\sigma }_{\mathrm{int}}}|\Phi \rangle ,(10)$

where

${\overline{H}}_{eff}^{\left(DUCC\right)}=\left(P+Q_{\mathrm{int}}\right)e^{-{\sigma }_{ext}}H{e}^{{\sigma }_{ext}}\left(P+Q_{\mathrm{int}}\right).(11)$

In Eq. 11, P and $Q_{\text{int}}$ are the projection operators onto the reference function and all active-space excited Slater determinants (with respect to $|\Phi \rangle$).

We will discuss the utility of the downfolding techniques in the next section for the ground state calculations of H₂. This is just one of the two approaches presented in this paper for reducing the dimensionality of the quantum problem (Figure 1).

FIGURE 1

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 (Bauman et al., 2019; Kowalski and Bauman, 2020); in this step, the many-body problem is rerepresented in a subspace of entire Hilbert space.

Results for the ¹${\text{Σ}}_g^{+}$ Ground State of H₂ 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ří Pittner. The plane-wave calculations used a simple cubic box with L = 26a₀ 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₂ 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.

FIGURE 2

FIGURE 2. The ground state energy curves for H₂ 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₂ molecule CCSD gives the same answer as FCI.

FIGURE 3

FIGURE 3. The HOMO and first three LUMOs generated from the straight HF calculation and the $H_1$ Hamiltonian are shown in the left and right panels, respectively. The orbitals are displayed in the order of decreasing orbital energy from top to bottom. The isovalues of positive and negative isosurfaces are 0.007226 and −0.007226 for the RHF LUMO 1; 0.01148 and −0.01148 for the RHF LUMO 2; 0.002404 and −0.002404 for the RHF LUMO 3; and 0.03117 and −0.03117 for the others. 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₂ is rotated by 90° in the bound LUMO 3 relative to the bound LUMO 2.

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.

TABLE 1

TABLE 1. Total energies as a function of distance from plane-wave FCI calculations for the H₂ molecule 19 $H_1$ virtual 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

$E_g=\underset{|\Psi \rangle }{\min }\langle \Psi |H|\Psi \rangle ,(12)$

where H represents the second-quantized Hamiltonian of Eq. 4 and $|\Psi \rangle$ labels the electronic configuration. Within the context of quantum computing, the fermionic representations of the Hamiltonian and state are transformed to alternate representations in terms of the spin operators. This transformation recast the molecular Hamiltonian into a representation $H_P$ that is defined in terms of the usual Pauli spin operators. Similarly, the electronic state $|\Psi \rangle$ is represented as a variable unitary operator $\widehat{U}\left(\overrightarrow{\theta }\right)$ acting on a fiducial quantum state (|0>). This yields the equivalent representation of Eq. 13 as

$E_g=\underset{\overrightarrow{\theta }}{\min }\langle 0|{\widehat{U}}^{†}\left(\overrightarrow{\theta }\right)H_P\widehat{U}\left(\overrightarrow{\theta }\right)|0\rangle .(13)$

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 $\widehat{U}\left(\overrightarrow{\theta }\right)$ 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 to define a set of virtual orbitals, $\left\{{\psi }_e^{\left(n\right)}\left(r\right)\right\}$ with $n=1\dots 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, ${\psi }_g\left(r\right)$ (or ground state orbitals for many-electron systems), and the virtual orbital to be optimized, ${\psi }_e^{\left(n\right)}\left(r\right)$, 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 ${\psi }_e^{\left(n\right)}\left(r\right)$ and then update with a conjugate-gradient or similar method while making sure that ${\psi }_e^{\left(n\right)}\left(r\right)$ is normalized and orthogonal to ${\psi }_g\left(r\right)$ and ${\psi }_e^{\left(m\right)}\left(r\right)$ for $m=1,\dots ,n-1$.

5. If the gradient is small, then $n=n+1$.

6. If $n\le N_{virtual}$, go to step 2; otherwise, finish.

In the case of the H₂ molecule, a small CI wave function for the 2 electron system composed of 2 one-electron orbitals, ${\psi }_g\left(r\right)$ and ${\psi }_e^{\left(n\right)}\left(r\right)$, can be written as a linear combination of 6 determinant wave functions, or just 3 determinant wave functions for just singlet (or triplet) states,

$\begin{array}{c}{\Psi }_i\left[{\psi }_g\left(r\right),{\psi }_e\left(r\right)\right]=c_g^{\left(i\right)}{\Psi }_g\left[{\psi }_g\left(r\right)\right]\\ +c_e^{\left(i\right)}{\Psi }_e\left[{\psi }_e\left(r\right)\right]\\ +c_m^{\left(i\right)}{\Psi }_m\left[{\psi }_g\left(r\right),{\psi }_e\left(r\right)\right]+\dots .\end{array}$

Using this small CI ansatz, the energies of the system can be obtained by diagonalizing the following eigenvalue equation:

$H{C}_i=E_{i}S{C}_i,$

where

$\begin{array}{l}H=\left[\begin{array}{ccc}\langle {\Psi }_g|H|{\Psi }_g\rangle & \langle {\Psi }_g|H|{\Psi }_e\rangle & \langle {\Psi }_g|H|{\Psi }_m\rangle \\ \langle {\Psi }_e|H|{\Psi }_g\rangle & \langle {\Psi }_e|H|{\Psi }_e\rangle & \langle {\Psi }_e|H|{\Psi }_m\rangle \\ \langle {\Psi }_m|H|{\Psi }_g\rangle & \langle {\Psi }_m|H|{\Psi }_e\rangle & \langle {\Psi }_m|H|{\Psi }_m\rangle \end{array}\right],\\ S=\left[\begin{array}{ccc}\langle {\Psi }_g|{\Psi }_g\rangle & \langle {\Psi }_g|{\Psi }_e\rangle & \langle {\Psi }_g|{\Psi }_m\rangle \\ \langle {\Psi }_e|{\Psi }_g\rangle & \langle {\Psi }_e|{\Psi }_e\rangle & \langle {\Psi }_e|{\Psi }_m\rangle \\ \langle {\Psi }_m|{\Psi }_g\rangle & \langle {\Psi }_m|{\Psi }_e\rangle & \langle {\Psi }_m|{\Psi }_m\rangle \end{array}\right],\\ C_i=\left[\begin{array}{c}{c}_g^{\left(i\right)}\\ c_e^{\left(i\right)}\\ c_m^{\left(i\right)}\end{array}\right].\end{array}(14)$

Note that the overlap matrix, S, is the identity matrix for orthonormal ${\psi }_g$ and ${\psi }_e$. The variation with respect to ${\psi }_e\left(r\right)$ can be simply obtained using the following formula:

$\begin{array}{c}\frac{\delta E_i}{\delta {\psi }_e^{\ast }\left(r\right)}=c_g^{\left(i\right)}\frac{\delta \langle {\Psi }_g|H|{\Psi }_g\rangle }{\delta {\psi }_e^{\ast }\left(r\right)}{c}_g^{\left(i\right)}+c_g^{\left(i\right)}\frac{\delta \langle {\Psi }_g|H|{\Psi }_e\rangle }{\delta {\psi }_e^{\ast }\left(r\right)}{c}_e^{\left(i\right)}\\ +c_g^{\left(i\right)}\frac{\delta \langle {\Psi }_g|H|{\Psi }_m\rangle }{\delta {\psi }_e^{\ast }\left(r\right)}{c}_m^{\left(i\right)}+c_e^{\left(i\right)}\frac{\delta \langle {\Psi }_e|H|{\Psi }_g\rangle }{\delta {\psi }_e^{\ast }\left(r\right)}{c}_g^{\left(i\right)}\\ +c_e^{\left(i\right)}\frac{\delta \langle {\Psi }_e|H|{\Psi }_e\rangle }{\delta {\psi }_e^{\ast }\left(r\right)}{c}_e^{\left(i\right)}+c_e^{\left(i\right)}\frac{\delta \langle {\Psi }_e|H|{\Psi }_m\rangle }{\delta {\psi }_e^{\ast }\left(r\right)}{c}_m^{\left(i\right)}\\ +c_m^{\left(i\right)}\frac{\delta \langle {\Psi }_m|H|{\Psi }_g\rangle }{\delta {\psi }_e^{\ast }\left(r\right)}{c}_g^{\left(i\right)}+c_m^{\left(i\right)}\frac{\delta \langle {\Psi }_m|H|{\Psi }_e\rangle }{\delta {\psi }_e^{\ast }\left(r\right)}{c}_e^{\left(i\right)}\\ +c_m^{\left(i\right)}\frac{\delta \langle {\Psi }_m|H|{\Psi }_m\rangle }{\delta {\psi }_e^{\ast }\left(r\right)}{c}_m^{\left(i\right)}.\end{array}(15)$

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 ${\psi }_e^{\ast }\left(r\right)$ 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.

One-Electron Orbitals for Two-State Hamiltonian

The four one-electron spin orbitals of two-state Hamiltonian are

$\begin{array}{l}{\chi }_1\left(x\right)={\psi }_g\left(r\right)\alpha \left(s\right),\\ {\chi }_2\left(x\right)={\psi }_g\left(r\right)\beta \left(s\right),\\ {\chi }_3\left(x\right)={\psi }_e\left(r\right)\alpha \left(s\right),\\ {\chi }_4\left(x\right)={\psi }_e\left(r\right)\beta \left(s\right),\end{array}$

where the spatial orbitals and spin functions are orthonormalized,

$\begin{array}{l}{\displaystyle \int {\psi }_g^{\ast }\left(r\right){\psi }_e\left(r\right)dr}={\displaystyle \int {\psi }_e^{\ast }\left(r\right){\psi }_g\left(r\right)dr}=0,\\ {\displaystyle \int {\psi }_g^{\ast }\left(r\right){\psi }_g\left(r\right)dr}={\displaystyle \int {\psi }_e^{\ast }\left(r\right){\psi }_e\left(r\right)dr=1},\\ {\displaystyle \int {\alpha }^{\ast }\left(s\right)\beta \left(s\right)ds}={\displaystyle \int {\beta }^{\ast }\left(s\right)\alpha \left(s\right)ds=0},\\ {\displaystyle \int {\alpha }^{\ast }\left(s\right)\alpha \left(s\right)ds}={\displaystyle \int {\beta }^{\ast }\left(s\right)\beta \left(s\right)ds=1}.\end{array}$

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

$\begin{array}{c}|{\chi }_1{\chi }_2\rangle ={\Psi }_g\left(x_1,x_2\right)\\ =\left[{\psi }_g\left(r_1\right){\psi }_g\left(r_2\right)\right]\frac{1}{\sqrt{2}}\left[\alpha \left(s_1\right)\beta \left(s_2\right)-\alpha \left(s_2\right)\beta \left(s_1\right)\right],\\ |{\chi }_3{\chi }_4\rangle ={\Psi }_e\left(x_1,x_2\right)\\ =\left[{\psi }_e\left(r_1\right){\psi }_e\left(r_2\right)\right]\frac{1}{\sqrt{2}}\left[\alpha \left(s_1\right)\beta \left(s_2\right)-\alpha \left(s_2\right)\beta \left(s_1\right)\right],\\ |{\chi }_1{\chi }_4\rangle ={\Psi }_a\left(x_1,x_2\right)\\ =\frac{1}{\sqrt{2}}\left[\begin{array}{c}{\psi }_g\left(r_1\right)\alpha \left(s_1\right){\psi }_e\left(r_2\right)\beta \left(s_2\right)\\ -{\psi }_e\left(r_1\right)\beta \left(s_1\right){\psi }_g\left(r_2\right)\alpha \left(s_2\right)\end{array}\right],\\ |{\chi }_2{\chi }_3\rangle ={\Psi }_b\left(x_1,x_2\right)\\ =\frac{1}{\sqrt{2}}\left[\begin{array}{c}{\psi }_g\left(r_1\right)\beta \left(s_1\right){\psi }_e\left(r_2\right)\alpha \left(s_2\right)\\ -{\psi }_e\left(r_1\right)\alpha \left(s_1\right){\psi }_g\left(r_2\right)\beta \left(s_2\right)\end{array}\right],\\ |{\chi }_1{\chi }_3\rangle ={\Psi }_u\left(x_1,x_2\right)\\ =\frac{1}{\sqrt{2}}\left[{\psi }_g\left(r_1\right){\psi }_e\left(r_2\right)-{\psi }_e\left(r_1\right){\psi }_g\left(r_2\right)\right]\left[\alpha \left(s_1\right)\alpha \left(s_2\right)\right],\\ |{\chi }_2{\chi }_4\rangle ={\Psi }_d\left(x_1,x_2\right)\\ =\frac{1}{\sqrt{2}}\left[{\psi }_g\left(r_1\right){\psi }_e\left(r_2\right)-{\psi }_e\left(r_1\right){\psi }_g\left(r_2\right)\right]\left[\beta \left(s_1\right)\beta \left(s_2\right)\right].\end{array}$

Note that ${\Psi }_a$ and ${\Psi }_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$,

$\begin{array}{c}{\Psi }_m={\Psi }_{a-b}=\frac{1}{\sqrt{2}}\left(|{\chi }_1{\chi }_4\rangle -|{\chi }_2{\chi }_3\rangle \right)\\ =\frac{1}{\sqrt{2}}\left[{\Psi }_a\left(x_1,x_2\right)-{\Psi }_b\left(x_1,x_2\right)\right]\\ =\frac{1}{\sqrt{2}}\left[{\psi }_g\left(r_1\right){\psi }_e\left(r_2\right)+{\psi }_e\left(r_1\right){\psi }_g\left(r_2\right)\right]\\ \times \frac{1}{\sqrt{2}}\left[\alpha \left(s_1\right)\beta \left(s_2\right)-\beta \left(s_1\right)\alpha \left(s_2\right)\right],\\ {\Psi }_p={\Psi }_{a+b}=\frac{1}{\sqrt{2}}\left(|{\chi }_1{\chi }_4\rangle +|{\chi }_2{\chi }_3\rangle \right)\\ =\frac{1}{\sqrt{2}}\left[{\Psi }_a\left(x_1,x_2\right)+{\Psi }_b\left(x_1,x_2\right)\right]\\ =\frac{1}{\sqrt{2}}\left[{\psi }_g\left(r_1\right){\psi }_e\left(r_2\right)-{\psi }_e\left(r_1\right){\psi }_g\left(r_2\right)\right]\\ \times \frac{1}{\sqrt{2}}\left[\alpha \left(s_1\right)\beta \left(s_2\right)+\beta \left(s_1\right)\alpha \left(s_2\right)\right].\end{array}$

Matrix Elements From the One-Electron Operators

The $H_1$ operator for H₂ molecule is

$H_1=h\left(r_1\right)+h\left(r_2\right),$

where $h\left(r\right)$ is a function/operator of the coordinate r; i.e.,

$\begin{array}{c}h\left(r\right)=-\frac{1}{2}{\nabla }_r^2+{\displaystyle \sum _{A=1}^{N_A}\left(V_{local}^{\left(A\right)}\left(\left|r-R_A\right|\right)+{\displaystyle \sum _{lm}{\widehat{V}}_{NL}^{\left(A\right),lm}}\right)},\\ \langle {\Psi }_g|H_1|{\Psi }_g\rangle =2{\displaystyle \int {\psi }_g^{\ast }\left(r\right)h\left(r\right){\psi }_g\left(r\right)dr},\\ \langle {\Psi }_g|H_1|{\Psi }_e\rangle =0,\\ \langle {\Psi }_e|H_1|{\Psi }_e\rangle =2{\displaystyle \int {\psi }_e^{\ast }\left(r\right)h\left(r\right){\psi }_e\left(r\right)dr},\\ \langle {\Psi }_m|H_1|{\Psi }_m\rangle ={\displaystyle \int {\psi }_g^{\ast }\left(r\right)h\left(r\right){\psi }_g\left(r\right)dr}+{\displaystyle \int {\psi }_e^{\ast }\left(r\right)h\left(r\right){\psi }_e\left(r\right)dr},\\ \langle {\Psi }_g|H_1|{\Psi }_m\rangle =\sqrt{2}{\displaystyle \int {\psi }_g^{\ast }\left(r\right)h\left(r\right){\psi }_e\left(r\right)dr},\\ \langle {\Psi }_m|H_1|{\Psi }_g\rangle =\sqrt{2}{\displaystyle \int {\psi }_e^{\ast }\left(r\right)h\left(r\right){\psi }_g\left(r\right)dr},\\ \langle {\Psi }_e|H_1|{\Psi }_m\rangle =\sqrt{2}{\displaystyle \int {\psi }_e^{\ast }\left(r\right)h\left(r\right){\psi }_g\left(r\right)dr},\\ \langle {\Psi }_m|H_1|{\Psi }_e\rangle =\sqrt{2}{\displaystyle \int {\psi }_g^{\ast }\left(r\right)h\left(r\right){\psi }_e\left(r\right)dr},\\ \frac{\delta \langle {\Psi }_g|H_1|{\Psi }_g\rangle }{\delta {\psi }_e^{\ast }\left(r\right)}=0,\\ \frac{\delta \langle {\Psi }_g|H_1|{\Psi }_e\rangle }{\delta {\psi }_e^{\ast }\left(r\right)}=0,\\ \frac{\delta \langle {\Psi }_e|H_1|{\Psi }_e\rangle }{\delta {\psi }_e^{\ast }\left(r\right)}=2h\left(r\right){\psi }_e\left(r\right),\\ \frac{\delta \langle {\Psi }_m|H_1|{\Psi }_m\rangle }{\delta {\psi }_e^{\ast }\left(r\right)}=h\left(r\right){\psi }_e\left(r\right),\\ \frac{\delta \langle {\Psi }_g|H_1|{\Psi }_m\rangle }{\delta {\psi }_e^{\ast }\left(r\right)}=0,\\ \frac{\delta \langle {\Psi }_m|H_1|{\Psi }_g\rangle }{\delta {\psi }_e^{\ast }\left(r\right)}=\sqrt{2}h\left(r\right){\psi }_g\left(r\right),\\ \frac{\delta \langle {\Psi }_e|H_1|{\Psi }_m\rangle }{\delta {\psi }_e^{\ast }\left(r\right)}=\sqrt{2}h\left(r\right){\psi }_g\left(r\right),\\ \frac{\delta \langle {\Psi }_m|H_1|{\Psi }_e\rangle }{\delta {\psi }_e^{\ast }\left(r\right)}=0.\end{array}$

Matrix Elements From the Two-Electron Operators

The $H_2$ operator for H₂ molecule is

$\begin{array}{l}{H}_2=\frac{1}{r_{12}}=\frac{1}{\left|r_1-r_2\right|},\\ \langle {\Psi }_g|H_2|{\Psi }_g\rangle ={\displaystyle \iint {\psi }_g^{\ast }\left(r\right){\psi }_g\left(r\right)\frac{1}{\left|r-r^{\prime }\right|}{\psi }_g^{\ast }\left(r^{\prime }\right){\psi }_g\left(r^{\prime }\right)drd{r}^{\prime }},\\ \langle {\Psi }_g|H_2|{\Psi }_e\rangle ={\displaystyle \iint {\psi }_g^{\ast }\left(r\right){\psi }_e\left(r\right)\frac{1}{\left|r-r^{\prime }\right|}{\psi }_g^{\ast }\left(r^{\prime }\right){\psi }_e\left(r^{\prime }\right)drd{r}^{\prime }},\\ \langle {\Psi }_e|H_2|{\Psi }_g\rangle ={\displaystyle \iint {\psi }_e^{\ast }\left(r\right){\psi }_g\left(r\right)\frac{1}{\left|r-r^{\prime }\right|}{\psi }_e^{\ast }\left(r^{\prime }\right){\psi }_g\left(r^{\prime }\right)drd{r}^{\prime }},\\ \langle {\Psi }_e|H_2|{\Psi }_e\rangle ={\displaystyle \iint {\psi }_e^{\ast }\left(r\right){\psi }_e\left(r\right)\frac{1}{\left|r-r^{\prime }\right|}{\psi }_e^{\ast }\left(r^{\prime }\right){\psi }_e\left(r^{\prime }\right)drd{r}^{\prime }},\\ \langle {\Psi }_m|H_2|{\Psi }_m\rangle =\left[{\displaystyle \iint {\psi }_e^{\ast }\left(r\right){\psi }_e\left(r\right)\frac{1}{\left|r-r^{\prime }\right|}{\psi }_g^{\ast }\left(r^{\prime }\right){\psi }_g\left(r^{\prime }\right)drd{r}^{\prime }+{\displaystyle \iint {\psi }_e^{\ast }\left(r\right){\psi }_g\left(r\right)\frac{1}{\left|r-r^{\prime }\right|}{\psi }_g^{\ast }\left(r^{\prime }\right){\psi }_e\left(r^{\prime }\right)drd{r}^{\prime }}}\right],\\ \langle {\Psi }_g|H_2|{\Psi }_m\rangle =\sqrt{2}{\displaystyle \iint {\psi }_g^{\ast }\left(r\right){\psi }_g\left(r\right)\frac{1}{\left|r-r^{\prime }\right|}{\psi }_g^{\ast }\left(r^{\prime }\right){\psi }_e\left(r^{\prime }\right)drd{r}^{\prime }},\\ \langle {\Psi }_m|H_2|{\Psi }_g\rangle =\sqrt{2}{\displaystyle \iint {\psi }_e^{\ast }\left(r\right){\psi }_g\left(r\right)\frac{1}{\left|r-r^{\prime }\right|}{\psi }_g^{\ast }\left(r^{\prime }\right){\psi }_g\left(r^{\prime }\right)drd{r}^{\prime }},\\ \langle {\Psi }_e|H_2|{\Psi }_m\rangle =\sqrt{2}{\displaystyle \iint {\psi }_e^{\ast }\left(r\right){\psi }_e\left(r\right)\frac{1}{\left|r-r^{\prime }\right|}{\psi }_e^{\ast }\left(r^{\prime }\right){\psi }_g\left(r^{\prime }\right)drd{r}^{\prime }},\\ \langle {\Psi }_m|H_2|{\Psi }_e\rangle =\sqrt{2}{\displaystyle \iint {\psi }_g^{\ast }\left(r\right){\psi }_e\left(r\right)\frac{1}{\left|r-r^{\prime }\right|}{\psi }_e^{\ast }\left(r^{\prime }\right){\psi }_e\left(r^{\prime }\right)drd{r}^{\prime }},\\ \frac{\delta \langle {\Psi }_g|H_2|{\Psi }_g\rangle }{\delta {\psi }_e^{\ast }\left(r\right)}=0,\\ \frac{\delta \langle {\Psi }_g|H_2|{\Psi }_e\rangle }{\delta {\psi }_e^{\ast }\left(r\right)}=0,\\ \frac{\delta \langle {\Psi }_e|H_2|{\Psi }_g\rangle }{\delta {\psi }_e^{\ast }\left(r\right)}=2\left[{\displaystyle \int \frac{{\psi }_e^{\ast }\left(r^{\prime }\right){\psi }_g\left(r^{\prime }\right)}{\left|r-r^{\prime }\right|}d{r}^{\prime }}\right]{\psi }_g\left(r\right),\\ \frac{\delta \langle {\Psi }_e|H_2|{\Psi }_e\rangle }{\delta {\psi }_e^{\ast }\left(r\right)}=2\left[{\displaystyle \int \frac{{\psi }_e^{\ast }\left(r^{\prime }\right){\psi }_e\left(r^{\prime }\right)}{\left|r-r^{\prime }\right|}d{r}^{\prime }}\right]{\psi }_e\left(r\right),\\ \frac{\delta \langle {\Psi }_m|H_2|{\Psi }_m\rangle }{\delta {\psi }_e^{\ast }\left(r\right)}=\left[{\displaystyle \int \frac{{\psi }_g^{\ast }\left(r^{\prime }\right){\psi }_g\left(r^{\prime }\right)}{\left|r-r^{\prime }\right|}d{r}^{\prime }}\right]{\psi }_e\left(r\right),\\ +\left[{\displaystyle \int \frac{{\psi }_g^{\ast }\left(r^{\prime }\right){\psi }_e\left(r^{\prime }\right)}{\left|r-r^{\prime }\right|}d{r}^{\prime }}\right]{\psi }_g\left(r\right),\\ \frac{\delta \langle {\Psi }_g|H_2|{\Psi }_m\rangle }{\delta {\psi }_e^{\ast }\left(r\right)}=0,\\ \frac{\delta \langle {\Psi }_m|H_2|{\Psi }_g\rangle }{\delta {\psi }_e^{\ast }\left(r\right)}=\sqrt{2}\left[{\displaystyle \int \frac{{\psi }_g^{\ast }\left(r^{\prime }\right){\psi }_g\left(r^{\prime }\right)}{\left|r-r^{\prime }\right|}d{r}^{\prime }}\right]{\psi }_g\left(r\right),\\ \frac{\delta \langle {\Psi }_e|H_2|{\Psi }_m\rangle }{\delta {\psi }_e^{\ast }\left(r\right)}=\sqrt{2}\left\{\left[{\displaystyle \int \frac{{\psi }_e^{\ast }\left(r^{\prime }\right){\psi }_g\left(r^{\prime }\right)}{\left|r-r^{\prime }\right|}d{r}^{\prime }}\right]{\psi }_e\left(r\right)+\left[{\displaystyle \int \frac{{\psi }_e^{\ast }\left(r^{\prime }\right){\psi }_e\left(r^{\prime }\right)}{\left|r-r^{\prime }\right|}d{r}^{\prime }}\right]{\psi }_g\left(r\right)\right\},\\ \frac{\delta \langle {\Psi }_m|H_2|{\Psi }_e\rangle }{\delta {\psi }_e^{\ast }\left(r\right)}=\sqrt{2}\left[{\displaystyle \int \frac{{\psi }_g^{\ast }\left(r^{\prime }\right){\psi }_e\left(r^{\prime }\right)}{\left|r-r^{\prime }\right|}d{r}^{\prime }}\right]{\psi }_e\left(r\right).\end{array}$

Results for ¹${\text{Σ}}_g^{+}$ Ground State of the H₂ Molecule Using Correlation Optimized Virtual Orbitals (COVOs)

The results for PW FCI calculations of H₂ 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.

FIGURE 4

FIGURE 4. Plots of total energies as a function of distance from plane-wave FCI calculations for the H₂ molecule with 1, 4, 8, 12, and 18 correlation optimized virtual orbitals. The top plot shows energy from R = 0.4 Å to R = 7.0 Å, and the bottom plot zooms in near the energy minima.

TABLE 2

TABLE 2. Total energies as a function of distance for the H₂ molecule from plane-wave FCI calculations with 1, 4, 8, 12, and 18 COVOs and ADAPT-VQE simulations with 4 COVOs. Nonparallelity errors (NPE) are evaluated with respect to the calculations with the largest virtual orbital space (PW FCI 18 COVOs) and reported in milli-Hartree.

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 ¹${\text{Σ}}_g^{+}$ Ground State of the H₂ 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 ¹${\text{Σ}}_g^{+}$ ground state of H₂ 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 () 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.

FIGURE 5

FIGURE 5. Potential energy curves for FCI and ADAPT-VQE (top) and the deviations in ADAPT-VQE energies with respect to FCI (bottom).

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 ¹${\text{Σ}}_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.

FIGURE 6

FIGURE 6. CNOT gate count for ADAPT-VQE and ordinary VQE with both singlet-adapted singles and doubles operators.

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.

FIGURE 7

FIGURE 7. Number of objective function calls as a function of the H-H distance for the different ansatz compositions (top) and the average number of objective function calls per ansatz size, with error bars representing one standard deviation (bottom). The bar colors on the top plot represent the ansatz sizes in the bottom plot.

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₂ 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 with just 4 virtual orbitals for the H₂ 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).

FIGURE 8

FIGURE 8. Summary of various plane-wave and Gaussian basis set RHF and FCI calculations for the ¹${\text{Σ}}_g^{+}$ ground state of H₂ 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₂ molecule CCSD gives the same answer as FCI.

Because this study is focused on how one might carry out plane-wave CI calculations on near-term quantum computers in the next few years, we have only shown results for the H₂ dimer. However, we are optimistic that these correlation optimized virtual orbitals open up the door to many-body calculations using pseudopotential plane-wave 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₂ 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.

FIGURE 9

FIGURE 9. Potential energy curves in kcal/mol for the Be₂ dimer using plane-wave Hartree–Fock and FCI with 5 COVOs calculations.

FIGURE 10

FIGURE 10. The 1 filled RHF orbital and 18 COVOs for the H₂ molecule and 2 filled RHF orbitals and 5 COVOs for the Be₂ 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₂ and 2.8 Å for Be₂. The positive and negative isosurfaces are colored in blue and orange, respectively. Notice that, for H₂, 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.

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 NWChem-DUCC interface and codeveloped the original form of the DUCC formalism. NB wrote a part of the NWChem-DUCC interface, codeveloped the original form of the DUCC formalism, carried out FCI and QPE calculations, and prepared formatted integrals for VQE calculations. DC wrote the algorithm for and performed the VQE quantum simulations.

Funding

This material is based upon work supported by the United States Department of Energy (DOE), Office of Science, Office of Basic Energy Sciences, Chemical Sciences, Geosciences, and Biosciences Division through its “Embedding Quantum Computing into Many-body Frameworks for Strongly Correlated Molecular and Materials Systems” project at Pacific Northwest National Laboratory (PNNL). This work was also supported by the Quantum Science Center (QSC), a National Quantum Information Science Research Center of the United States Department of Energy (DOE). We also would like to thank the DOE BES Chem CCS, DOE BES Geochemistry, and DOE Advanced Scientific Computing Research (ASCR) ECP NWChemEx programs for their support of software development for high-performance computers and computer time needed to carry out the work. PNNL is operated for the United States Department of Energy by the Battelle Memorial Institute under Contract DE-AC06-76RLO-1830. This research was also partially supported, thru their support of software development for high-performance computers, by DOE BES Chem CCS and DOE BES Geosciences programs, as well as the Exascale Computing Project (17-SC-20-SC), a collaborative effort of the United States Department of Energy Office of Science and the National Nuclear Security Administration. Calculations have been performed using computational resources at the Pacific Northwest National Laboratory (PNNL). This research used resources of the National Energy Research Scientific Computing Center (NERSC), a User Facility supported by the Office of Science of the United States DOE under Contract No. DE-AC02-05CH11231, and the Argonne ALCF computing center through their early science program. This research also used resources of the Oak Ridge Leadership Computing Facility, which is a DOE Office of Science User Facility supported under Contract DE-AC05-00OR22725. This manuscript has been authored in part by UT-Battelle, LLC under Contract No. DE-AC05-00OR22725 with the United States Department of Energy. The United States Government retains and the publisher, by accepting the article for publication, acknowledges that the United States Government retains a nonexclusive, paid-up, irrevocable, worldwide license to publish or reproduce the published form of this manuscript or allow others to do so, for United States Government purposes. The Department of Energy will provide public access to these results of federally sponsored research in accordance with the DOE Public Access Plan.

Conflict of Interest

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Acknowledgments

We would like to thank the NWChem project team and the people that have helped the progress of the NWChem software over the years. DC would like to thank Alexander McCaskey for discussions and help with the software engineering of the quantum algorithm used in this paper.

References

Adamowicz, L., and Bartlett, R. J. (1987). Optimized virtual orbital space for high‐level correlated calculations. J. Chem. Phys. 86, 6314–6324. doi:10.1063/1.452468

CrossRef Full Text | Google Scholar

Adamowicz, L., Bartlett, R. J., and Sadlej, A. J. (1988). Optimized virtual orbital space for high‐level correlated calculations. II. Electric properties. J. Chem. Phys. 88, 5749–5758. doi:10.1063/1.454721

CrossRef Full Text | Google Scholar

Allen, L. C., and Karo, A. M. (1962). Electronic structure of simple Molecules1. J. Phys. Chem. 66, 2329–2331. doi:10.1021/j100818a011

CrossRef Full Text | Google Scholar

Aprà, E., Bylaska, E. J., De Jong, W. A., Govind, N., Kowalski, K., Straatsma, T. P., et al. (2020). NWChem: past, present, and future. J. Chem. Phys. 152, 184102. doi:10.1063/5.0004997

PubMed Abstract | CrossRef Full Text | Google Scholar

Bauman, N. P., Bylaska, E. J., Krishnamoorthy, S., Low, G. H., Wiebe, N., Granade, C. E., et al. (2019). Downfolding of many-body Hamiltonians using active-space models: extension of the sub-system embedding sub-algebras approach to unitary coupled cluster formalisms. J. Chem. Phys. 151, 014107. doi:10.1063/1.5094643

PubMed Abstract | CrossRef Full Text | Google Scholar

Betteridge, P. W., Carruthers, J. R., Cooper, R. I., Prout, K., and Watkin, D. J. (2003). CRYSTALSversion 12: software for guided crystal structure analysis. J. Appl. Crystallogr. 36, 1487. doi:10.1107/S0021889803021800

CrossRef Full Text | Google Scholar

Boys, S. F., Cook, G. B., Reeves, C. M., and Shavitt, I. (1956). Automatic fundamental calculations of molecular structure. Nature 178, 1207–1209. doi:10.1038/1781207a0

CrossRef Full Text | Google Scholar

Boys, S. F. (1950). Electronic wave functions-i. a general method of calculation for the stationary states of any molecular system. Proc. R. Soc. Lond. A Math. Phys. Sci. 200, 542–554. doi:10.1098/rspa.1950.0036

CrossRef Full Text | Google Scholar

Bravyi, S. B., and Kitaev, A. Y. (2002). Fermionic quantum computation. Ann. Phys. 298, 210–226. doi:10.1006/aphy.2002.6254

CrossRef Full Text | Google Scholar

Briggs, E., Sullivan, D., and Bernholc, J. (1996). Real-space multigrid-based approach to large-scale electronic structure calculations. Phys. Rev. B Condens. Matter 54, 14362. doi:10.1103/PhysRevB.54.14362

PubMed Abstract | CrossRef Full Text | Google Scholar

Bylaska, E. J., and Rosso, K. (2018). Corresponding orbitals derived from periodic bloch states for electron transfer calculations of transition metal oxides. J. Chem. Theor. Comput. 14, 4416–4426. doi:10.1021/acs.jctc.7b01180

CrossRef Full Text | Google Scholar

Bylaska, E. J., Aprà, E., Kowalski, K., Jacquelin, M., De Jong, W. A., Vishnu, A., et al. (2017a). “Transitioning NWChem to the next generation of manycore machines,” in Exascale scientific applications: scalability and performance portability. Editors T. P. Straatsma, K. B. Antypas, and T. J. J. Timothy Williams (Boca Raton, FL: CRC Press), 165. doi:10.1201/b21930-8

CrossRef Full Text | Google Scholar

Bylaska, E. J., Jacquelin, M., de Jong, W. A., Hammond, J. R., and Klemm, M. (2017b). “Performance evaluation of NWChem ab-initio molecular dynamics (AIMD) simulations on the Intel® xeon Phi™ processor,” in High performance computing. Editors J. M. Kunkel, R. Yokota, M. Taufer, and J. Shalf (Cham, Switzerland: Springer International Publishing), 404–418. doi:10.1007/978-3-319-67630-2_30

CrossRef Full Text | Google Scholar

Bylaska, E. J. (2017). “Plane-wave DFT methods for chemistry,” in Annuual Reports in Computational Chemistry. Editor D. A. Dixon (Amsterdam, Netherlands: Elsevier), 13, 185–228. doi:10.1016/bs.arcc.2017.06.006

CrossRef Full Text | Google Scholar

Bylaska, E. J., Waters, K., Hermes, E. D., Zádor, J., and Rosso, K. M. (2020). A filon-like integration strategy for calculating exact exchange in periodic boundary conditions: a plane-wave dft implementation. Mater. Theory. 4, 1–31. doi:10.1186/s41313-020-00019-9

CrossRef Full Text | Google Scholar

Bylaska, E., Tsemekhman, K., Govind, N., and Valiev, M. (2011). “Large-scale plane-wave-based density functional theory: formalism, parallelization, and applications,” in Computational methods for large systems: electronic structure approaches for biotechnology and nanotechnology. Editor J. R. Reimers (Hoboken, NJ: Wiley Online Library), 77–116. doi:10.1002/9780470930779.ch3

CrossRef Full Text | Google Scholar

Challacombe, M. (2000). Linear scaling computation of the fock matrix. v. hierarchical cubature for numerical integration of the exchange-correlation matrix. J. Chem. Phys. 113, 10037–10043. doi:10.1063/1.1316012

CrossRef Full Text | Google Scholar

Chen, M., Guo, G., and He, L. (2011). Electronic structure interpolation via atomic orbitals. J. Phys. Condens. Matter 23, 325501. doi:10.1088/0953-8984/23/32/325501

PubMed Abstract | CrossRef Full Text | Google Scholar

Clark, S. J., Segall, M. D., Pickard, C. J., Hasnip, P. J., Probert, M. I. J., Refson, K., et al. (2005). First principles methods using castep. Z. für Kristallogr. - Cryst. Mater. 220, 567–570. doi:10.1524/zkri.220.5.567.65075

CrossRef Full Text | Google Scholar

Coulson, C. A., and Fischer, I. (1949). Notes on the molecular orbital treatment of the hydrogen molecule. The London, Edinburgh, and Dublin Philosophical Magazine and Journal of Science. 40, 386. doi:10.1080/14786444908521726

CrossRef Full Text | Google Scholar

Crawford, T. D., Clementi, C., Harrison, R., Head-Gordon, T., Jha, S., Krylov, A., et al. (2017). The molecular sciences software institute. Abstr. Pap. Am. Chem. Soc. Vol. 254. WASHINGTON, DC: AMER CHEMICAL SOC 1155 16TH ST, NW.

Google Scholar

Deslippe, J., Samsonidze, G., Strubbe, D. A., Jain, M., Cohen, M. L., and Louie, S. G. (2012). Berkeleygw: a massively parallel computer package for the calculation of the quasiparticle and optical properties of materials and nanostructures. Comput. Phys. Commun. 183, 1269–1289. doi:10.1016/j.cpc.2011.12.006

CrossRef Full Text | Google Scholar

Dubash, M. (2005). Moore’s law is dead, says gordon moore. Techworld. com 13.

Google Scholar

Dunning, T. H. (1989). Gaussian basis sets for use in correlated molecular calculations. i. the atoms boron through neon and hydrogen. J. Chem. Phys. 90, 1007–1023. doi:10.1063/1.456153

CrossRef Full Text | Google Scholar

Dunning, T. H., and Hay, P. J. (1977). Gaussian basis sets for molecular calculations. Methods of electronic structure theory 3, 1–27. doi:10.1007/978-1-4757-0887-5_1

CrossRef Full Text | Google Scholar

Dupuis, M., Watts, J. D., Villar, H. O., and Hurst, G. J. B. (1989). The general atomic and molecular electronic structure system hondo: version 7.0. Comput. Phys. Commun. 52, 415–425. doi:10.1016/0010-4655(89)90116-1

CrossRef Full Text | Google Scholar

Evangelista, F. A., Chan, G. K., and Scuseria, G. E. (2019). Exact parameterization of fermionic wave functions via unitary coupled cluster theory. J. Chem. Phys. 151, 244112. doi:10.1063/1.5133059

PubMed Abstract | CrossRef Full Text | Google Scholar

Frisch, M. J., Trucks, G. W., Schlegel, H. B., Scuseria, G. E., Robb, M. A., Cheeseman, J. R., et al. (2016). Gaussian 16 revision C.01. Wallingford CT: Gaussian Inc.

Google Scholar

Giannozzi, P., Baroni, S., Bonini, N., Calandra, M., Car, R., Cavazzoni, C., et al. (2009). Quantum espresso: a modular and open-source software project for quantum simulations of materials. J. Phys. Condens. Matter 21, 395502. doi:10.1088/0953-8984/21/39/395502

PubMed Abstract | CrossRef Full Text | Google Scholar

Goldbeck, G. (2017). The scientific software industry: a general overview. Goldbeck Consulting, Technical Report.

Google Scholar

Gonze, X., Jollet, F., Araujo, F. A., Adams, D., Amadon, B., Applencourt, T., et al. (2016). Recent developments in the abinit software package. Comput. Phys. Commun. 205, 106–131. doi:10.1016/j.cpc.2016.04.003

CrossRef Full Text | Google Scholar

Grimsley, H. R., Claudino, D., Economou, S. E., Barnes, E., and Mayhall, N. J. (2020). Is the trotterized uccsd ansatz chemically well-defined? J. Chem. Theor. Comput. 16, 1–6. doi:10.1021/acs.jctc.9b01083

CrossRef Full Text | Google Scholar

Grimsley, H. R., Economou, S. E., Barnes, E., and Mayhall, N. J. (2019). An adaptive variational algorithm for exact molecular simulations on a quantum computer. Nat. Commun. 10, 3007. doi:10.1038/s41467-019-10988-2

PubMed Abstract | CrossRef Full Text | Google Scholar

Gygi, F. (2008). Architecture of qbox: a scalable first-principles molecular dynamics code. IBM J. Res. Dev. 52, 137–144. doi:10.1147/rd.521.0137

CrossRef Full Text | Google Scholar

Hall, G. (1951). The molecular orbital theory of chemical valency viii. a method of calculating ionization potentials. Proc. R. Soc. Lond. A Math. Phys. Sci. 205, 541–552. doi:10.1098/rspa.1951.0048

CrossRef Full Text | Google Scholar

Hamann, D. R. (1989). Generalized norm-conserving pseudopotentials. Phys. Rev. B Condens. Matter 40, 2980. doi:10.1103/physrevb.40.2980

PubMed Abstract | CrossRef Full Text | Google Scholar

Harrison, R. J., Beylkin, G., Bischoff, F. A., Calvin, J. A., Fann, G. I., Fosso-Tande, J., et al. (2016). Madness: a multiresolution, adaptive numerical environment for scientific simulation. SIAM J. Sci. Comput. 38, S123–S142. doi:10.1137/15M1026171

CrossRef Full Text | Google Scholar

Hocquet, A., and Wieber, F. (2017). “only the initiates will have the secrets revealed”: computational chemists and the openness of scientific software. IEEE Ann. Hist. Comput. 39, 40–58. doi:10.1109/MAHC.2018.1221048

CrossRef Full Text | Google Scholar

Hutter, J., Iannuzzi, M., Schiffmann, F., and VandeVondele, J. (2014). cp2k: atomistic simulations of condensed matter systems. WIREs Comput. Mol. Sci. 4, 15–25. doi:10.1002/wcms.1159

CrossRef Full Text | Google Scholar

Hylleraas, E. A. (1929). Neue berechnung der energie des heliums im grundzustande, sowie des tiefsten terms von ortho-helium. Z. Phys. 54, 347–366. doi:10.1007/bf01375457

CrossRef Full Text | Google Scholar

Hylleraas, E. A. (1964). The schrödinger two-electron atomic problem. Adv. Quant. Chem. 1, 1–33. doi:10.1016/s0065-3276(08)60373-1

CrossRef Full Text | Google Scholar

Hylleraas, E. A. (1930). Über den Grundterm der Zweielektronenprobleme von H?, He, Li+, Be++ usw. Z. Phys. 65, 209–225. doi:10.1007/bf01397032

CrossRef Full Text | Google Scholar

Hylleraas, E. A. (1928). Über den Grundzustand des Heliumatoms. Z. Phys. 48, 469–494. doi:10.1007/bf01340013

CrossRef Full Text | Google Scholar

Johnson, S. (2007). The NLopt nonlinear-optimization package. Available at: http://github.com/stevengj/nlopt.doi:10.1364/np.2007.ntub1

CrossRef Full Text | Google Scholar

Jordan, P., and Wigner, E. (1928). Über das Paulische äquivalenzverbot. Z. Phys. 47, 631–651. doi:10.1007/BF01331938

CrossRef Full Text | Google Scholar

Kandala, A., Mezzacapo, A., Temme, K., Takita, M., Brink, M., Chow, J. M., et al. (2017). Hardware-efficient variational quantum eigensolver for small molecules and quantum magnets. Nature 549, 242–246. doi:10.1038/nature23879

PubMed Abstract | CrossRef Full Text | Google Scholar

Kendall, R. A., Aprà, E., Bernholdt, D. E., Bylaska, E. J., Dupuis, M., Fann, G. I., et al. (2000). High performance computational chemistry: an overview of nwchem a distributed parallel application. Comput. Phys. Commun. 128, 260–283. doi:10.1016/S0010-4655(00)00065-5

CrossRef Full Text | Google Scholar

King, H. F., Stanton, R. E., Kim, H., Wyatt, R. E., and Parr, R. G. (1967). Corresponding orbitals and the nonorthogonality problem in molecular quantum mechanics. J. Chem. Phys. 47, 1936–1941. doi:10.1063/1.1712221

CrossRef Full Text | Google Scholar

Kleinman, L., and Bylander, D. M. (1982). Efficacious form for model pseudopotentials. Phys. Rev. Lett. 48, 1425. doi:10.1103/physrevlett.48.1425

CrossRef Full Text | Google Scholar

Koga, T. (1992). Hylleraas wave functions revisited. J. Chem. Phys. 96, 1276–1279. doi:10.1063/1.462164

CrossRef Full Text | Google Scholar

Kohn, W., and Sham, L. J. (1965). Self-consistent equations including exchange and correlation effects. Phys. Rev. 140, A1133. doi:10.1103/PhysRev.140.A1133

CrossRef Full Text | Google Scholar

Kowalski, K., and Bauman, N. P. (2020). Publisher's Note: “Sub-system quantum dynamics using coupled cluster downfolding techniques” [J. Chem. Phys. 152, 244127 (2020)]. J. Chem. Phys. 153, 049902. doi:10.1063/5.0008436 doi:10.1063/5.0021137

PubMed Abstract | CrossRef Full Text | Google Scholar

Kresse, G., and Furthmüller, J. (1996). Efficiency of ab-initio total energy calculations for metals and semiconductors using a plane-wave basis set. Comput. Mater. Sci. 6, 15–50. doi:10.1016/0927-0256(96)00008-0

CrossRef Full Text | Google Scholar

Low, G. H., Bauman, N. P., Granade, C. E., Peng, B., Wiebe, N., Bylaska, E. J., et al. (2019). Q# and nwchem: tools for scalable quantum chemistry on quantum computers. Preprint: arXiv:1904.01131.

Google Scholar

McCaskey, A., Dumitrescu, E., Chen, M., Lyakh, D., and Humble, T. (2018a). Validating quantum-classical programming models with tensor network simulations. PloS One 13, e0206704–19. doi:10.1371/journal.pone.0206704

PubMed Abstract | CrossRef Full Text | Google Scholar

McCaskey, A. J., Dumitrescu, E. F., Liakh, D., Chen, M., Feng, W., and Humble, T. S. (2018b). A language and hardware independent approach to quantum-classical computing. Software 7, 245–254. doi:10.1016/j.softx.2018.07.007

CrossRef Full Text | Google Scholar

McCaskey, A. J., Lyakh, D. I., Dumitrescu, E. F., Powers, S. S., and Humble, T. S. (2020). XACC: a system-level software infrastructure for heterogeneous quantum-classical computing. Quantum Sci. Technol. 5, 024002. doi:10.1088/2058-9565/ab6bf6

CrossRef Full Text | Google Scholar

McCaskey, A. J., Parks, Z. P., Jakowski, J., Moore, S. V., Morris, T. D., Humble, T. S., et al. (2019). Quantum chemistry as a benchmark for near-term quantum computers. Npj Quantum Inf. 5, 1–8. doi:10.1038/s41534-019-0209-0

CrossRef Full Text | Google Scholar

Moore, G. E. (2006). Cramming more components onto integrated circuits, Reprinted from Electronics, volume 38, number 8, April 19, 1965, pp.114 ff. IEEE Solid-State Circuits Soc. Newsl. 11, 33–35. doi:10.1109/N-SSC.2006.4785860

CrossRef Full Text | Google Scholar

Mulliken, R. S., Rieke, C. A., and Brown, W. G. (1941). Hyperconjugation. J. Am. Chem. Soc. 63, 41–56. doi:10.1021/ja01846a008

CrossRef Full Text | Google Scholar

Mulliken, R. S., and Rieke, C. A. (1941). Improved computations on conjugation and hyperconjugation. J. Am. Chem. Soc. 63, 1770–1771. doi:10.1021/ja01851a514

CrossRef Full Text | Google Scholar

Neese, F. (2018). Software update: the orca program system, version 4.0. WIREs Comput Mol Sci. 8, e1327. doi:10.1002/wcms.1327

CrossRef Full Text | Google Scholar

Nesbet, R. K. (1963). Computer programs for electronic wave-function calculations. Rev. Mod. Phys. 35, 552. doi:10.1103/RevModPhys.35.552

CrossRef Full Text | Google Scholar

Nesbet, R. K. (1960). Ground state electronic wave function of methane. J. Chem. Phys. 32, 1114–1122. doi:10.1063/1.1730858

CrossRef Full Text | Google Scholar

O’Malley, P. J. J., Babbush, R., Kivlichan, I. D., Romero, J., McClean, J. R., Barends, R., et al. (2016). Scalable quantum simulation of molecular energies. Phys. Rev. X. 6, 031007. doi:10.1103/PhysRevX.6.031007

CrossRef Full Text | Google Scholar

Parr, R. G., and Mulliken, R. S. (1950). LCAOSelf‐consistent field calculation of the π‐electron energy levels of cis‐ andtrans−1,3‐butadiene. J. Chem. Phys. 18, 1338–1346. doi:10.1063/1.1747474

CrossRef Full Text | Google Scholar

Peruzzo, A., McClean, J., Shadbolt, P., Yung, M. H., Zhou, X. Q., Love, P. J., et al. (2014). A variational eigenvalue solver on a photonic quantum processor. Nat. Commun. 5, 4213. doi:10.1038/ncomms5213

PubMed Abstract | CrossRef Full Text | Google Scholar

Pople, J. A., Santry, D. P., and Segal, G. A. (1965). Approximate self‐consistent molecular orbital theory. I. Invariant procedures. J. Chem. Phys. 43, S129–S135. doi:10.1063/1.1701475

CrossRef Full Text | Google Scholar

Powell, M. J. D. (1994). A direct search optimization method that models the objective and constraint functions by linear interpolation. Dordrecht, Netherlands: Springer Netherlands, 51–67. doi:10.1007/978-94-015-8330-5_4

CrossRef Full Text | Google Scholar

Prendergast, D., and Louie, S. G. (2009). Bloch-state-based interpolation: an efficient generalization of the shirley approach to interpolating electronic structure. Phys. Rev. B 80, 235126. doi:10.1103/PhysRevB.80.235126

CrossRef Full Text | Google Scholar

Preskill, J. (2018). Quantum computing in the nisq era and beyond. Quantum. 2, 79. doi:10.22331/q-2018-08-06-79

CrossRef Full Text | Google Scholar

Pulay, P. (1969). Ab initiocalculation of force constants and equilibrium geometries in polyatomic molecules. Mol. Phys. 17, 197–204. doi:10.1080/00268976900100941

CrossRef Full Text | Google Scholar

Reeves, C. M. (1966). An algorithm for generating projective reduction formulas for matrix elements of many-electron wavefunctions. Commun. ACM 9, 276–279. doi:10.1145/365278.365540

CrossRef Full Text | Google Scholar

Reiher, M., Wiebe, N., Svore, K. M., Wecker, D., and Troyer, M. (2017). Elucidating reaction mechanisms on quantum computers. Proc. Natl. Acad. Sci. U.S.A. 114, 7555–7560. doi:10.1073/pnas.1619152114

PubMed Abstract | CrossRef Full Text | Google Scholar

Richard, R. M., Bertoni, C., Boschen, J. S., Keipert, K., Pritchard, B., Valeev, E. F., et al. (2019). Developing a computational chemistry framework for the exascale era. Comput. Sci. Eng. 21, 48–58. doi:10.1109/MCSE.2018.2884921

CrossRef Full Text | Google Scholar

Romero, J., Babbush, R., McClean, J. R., Hempel, C., Love, P. J., and Aspuru-Guzik, A. (2017). Strategies for quantum computing molecular energies using the unitary coupled cluster ansatz. Quantum Sci. Technol. 4, 014008. doi:10.1088/2058-9565/aad3e4

CrossRef Full Text | Google Scholar

Roothaan, C. C. J. (1951). New developments in molecular orbital theory. Rev. Mod. Phys. 23, 69. doi:10.1103/RevModPhys.23.69

CrossRef Full Text | Google Scholar

Rotman, D. (2020). We’re not prepared for the end of moore’s law. MIT Technology Review.

Google Scholar

Schmidt, M. W., Baldridge, K. K., Boatz, J. A., Elbert, S. T., Gordon, M. S., Jensen, J. H., et al. (1993). General atomic and molecular electronic structure system. J. Comput. Chem. 14, 1347–1363. doi:10.1002/jcc.540141112

CrossRef Full Text | Google Scholar

Setia, K., and Whitfield, J. D. (2018). Bravyi-kitaev superfast simulation of electronic structure on a quantum computer. J. Chem. Phys. 148, 164104. doi:10.1063/1.5019371

PubMed Abstract | CrossRef Full Text | Google Scholar

Shao, Y., Gan, Z., Epifanovsky, E., Gilbert, A. T., Wormit, M., Kussmann, J., et al. (2015). Advances in molecular quantum chemistry contained in the q-chem 4 program package. Mol. Phys. 113, 184–215. doi:10.1080/00268976.2014.952696

CrossRef Full Text | Google Scholar

Shirley, E. L. (1996). Optimal basis sets for detailed brillouin-zone integrations. Phys. Rev. B Condens. Matter. 54, 16464. doi:10.1103/PhysRevB.54.16464

PubMed Abstract | CrossRef Full Text | Google Scholar

Stanton, J. F., Gauss, J. r., Watts, J. D., Lauderdale, W. J., and Bartlett, R. J. (1992). The aces ii program system. Int. J. Quant. Chem. 44, 879–894. doi:10.1002/qua.560440876

CrossRef Full Text | Google Scholar

Steudtner, M., and Wehner, S. (2018). Lowering qubit requirements for quantum simulations of fermionic systems. Preprint: arXiv:1712.07067v2.

Google Scholar

Svore, K., Roetteler, M., Geller, A., Troyer, M., Azariah, J., Granade, C., et al. (2018). “Q# enabling scalable quantum computing and development with a high-level dsl,” in Proceedings of the Real World Domain Specific Languages Workshop 2018, Vienna, Austria, February 2018. 1–10. doi:10.1145/3183895.3183901

CrossRef Full Text | Google Scholar

te Velde, G., Bickelhaupt, F. M., Baerends, E. J., Fonseca Guerra, C., van Gisbergen, S. J. A., Snijders, J. G., et al. (2001). Chemistry with adf. J. Comput. Chem. 22, 931–967. doi:10.1002/jcc.1056

CrossRef Full Text | Google Scholar

Valiev, M., Bylaska, E. J., Govind, N., Kowalski, K., Straatsma, T. P., Van Dam, H. J. J., et al. (2010). NWChem: a comprehensive and scalable open-source solution for large scale molecular simulations. Comput. Phys. Commun. 181, 1477–1489. doi:10.1016/j.cpc.2010.04.018

CrossRef Full Text | Google Scholar

van Dam, H., Apra, E., Bair, R., Boschen, J., Bylaska, E., De Jong, W., et al. (2020). NWChemEx–computational chemistry for the exascale era. Bull. Am. Phys. Soc. 65.

Google Scholar

Verma, P., Huntington, L., Coons, M., Kawashima, Y., Yamazaki, T., and Zaribafiyan, A. (2020). Scaling up electronic structure calculations on quantum computers: the frozen natural orbital based method of increments. Preprint: arXiv:2002.07901.

Google Scholar

Wasielewski, M. R., Forbes, M. D. E., Frank, N. L., Kowalski, K., Scholes, G. D., Yuen-Zhou, J., et al. (2020). Exploiting chemistry and molecular systems for quantum information science. Nat. Rev. Chem. 4, 490–15. doi:10.1038/s41570-020-0200-5

CrossRef Full Text | Google Scholar

Werner, H.-J., Knowles, P. J., Knizia, G., Manby, F. R., and Schütz, M. (2012). Molpro: a general-purpose quantum chemistry program package. WIREs Comput. Mol. Sci. 2, 242–253. doi:10.1002/wcms.82

CrossRef Full Text | Google Scholar

Wikipedia, The Free Encyclopedia (2020). List of quantum chemistry and solid-state physics software. San Francisco: Wikimedia Foundation Inc.

Google Scholar

Wilkins-Diehr, N., and Crawford, T. D. (2018). NSF's inaugural software institutes: the science gateways community institute and the molecular sciences software institute. Comput. Sci. Eng. 20, 26–38. doi:10.1109/MCSE.2018.05329813

CrossRef Full Text | Google Scholar