Periodic Plane-Wave Electronic Structure Calculations on Quantum Computers

A procedure for defining virtual spaces, and the periodic one-electron and two-electron integrals, for plane-wave second quantized Hamiltonians has been developed and demonstrated using full configuration interaction (FCI) simulations and variational quantum eigensolver (VQE) circuits on Quantinuum's ion trap quantum computers accessed through Microsoft's Azure Quantum service. This work is an extension to periodic systems of a new class of algorithms in which the virtual spaces were generated by optimizing orbitals from small pairwise CI Hamiltonians, which we term as correlation optimized virtual orbitals with the abbreviation COVOs. In this extension, the integration of the first Brillouin zone is automatically incorporated into the two-electron integrals. With these procedures we have been able to derive virtual spaces, containing only a few orbitals, that were able to capture a significant amount of correlation. The focus in this manuscript is on comparing the simulations of small molecules calculated with plane-wave basis sets with large periodic unit cells at the $\Gamma$-point, including images, to results for plane-wave basis sets with aperiodic unit cells. The results for this approach were promising as we were able to obtain good agreement between periodic and aperiodic results for an LiH molecule. Simulations performed on the Quantinuum H1-1 quantum computer were able to produce surprisingly good energies, reproducing the FCI values for the 1 COVO Hamiltonian to within 11 milliHartree (6.9 kcal/mol), when corrected for noise.


Introduction
With the arrival of quantum computers, researchers are actively developing new algorithms to carry out quantum chemistry calculations on these platforms, in particular for calculations containing strong electron-electron correlations (aka high-level quantum chemistry methods). This is because it is anticipated that quantum computers with 50-100 qubits will eventually surpass classical digital computers for these types of calculations (Preskill 2018). However, in order for quantum computing to reach its full potential, there are hardware and software challenges that need to be addressed before it can become a viable replacement (Wasielewski et al. 2020) for existing high-performance classical computers and the associated cutting-edge parallel software that have been developed in the last two decades.
Most high-level quantum chemistry methods in use today (e.g., full configuration interaction (FCI) (Szabo and Ostlund 2012;Ross 1952;Gan et al. 2006;McArdle et al. 2020;Tubman et al. 2020;Sugisaki et al. 2021;Kawashima et al. 2021), coupled cluster (CC) (Coester 1958;Coester and Kummel 1960;Čížek 1966;Paldus et al. 1972;Mukherjee et al. 1975;Purvis and Bartlett 1982;Bishop and Kümmel 1987;Paldus and Li 1999;Crawford and Schaefer 2000;Bartlett and Musiał 2007;Peng et al. 2021), and Green's function (GF) (Green 1854(Green , 2014Feynman 1949Feynman , 1948Martin and Schwinger 1959;Baym and Kadanoff 1961;Peng et al. 2021) approaches) are based on second-quantized Hamiltonians, which are written in terms of the creation and annihilation operators of the Fermion orbitals along with the one-electron and two-electron integrals for the system. In principle, this formulation is exact, however, conventional computing methods are restricted in their accuracy due to the prohibitive computational cost for exact modeling of the exponentially growing wavefunction from the basis set that is introduced. As a result, these basis sets are typically highly engineered. One of the first, and still popular, class of basis sets used in quantum chemistry methods are atomic-like orbitals or the linear combination of atomic orbitals (LCAO) basis set. Pioneered by J. Lennard-Jones (1929), L. Pauling (1931), and J.C. Slater (1930), the atomic orbitals are generated by carrying out an atom calculation for each kind of atom in the system; guided by the heuristic that says the electronic states of a molecule or solid can be thought of as a superposition of atomic orbitals. For high-level methods, a popular basis set is the Dunning correlation consistent basis set (Dunning and Hay 1977;Dunning 1989;Prascher et al. 2011), in which the atomic orbitals are optimized at the configuration interaction singles and doubles (CISD) level of theory (Handy 1980). While the size of this intuitive, optimistically a priori, class of basis set is small compared to modern style basis sets that are more complete, e.g., plane-wave basis sets, it still needs to contain enough atomic orbitals to produce a truly accurate result.
Another challenge is calculating the two-electron integrals for condensed phase systems, since one typically wants to use periodic boundary conditions to carry out the simulation. While this is natural for plane-wave DFT methods (Bylaska et al. 2011a, b;Bylaska 2017;Pickett 1989;Ihm et al. 1979;Car and Parrinello 1985;Payne et al. 1992;Remler and Madden 1990;Kresse and Furthmüller 1996;Hutter 2000, 2009;Martin 2004;Bylaska et al. 2009;Chen et al. 2016;Gygi 2008;Peng et al. 2021) with low levels of theory, it is significantly more complicated to calculate exact exchange (Gygi and Baldereschi 1985, 1986Bylaska et al. 2011;Bylaska 2017;Chawla and Voth 1998;Sorouri et al. 2006;Marsman et al. 2008;Görling 1996;DiStasio Jr et al. 2014;Peng et al. 2021) and the other two-electron integrals (Bylaska and Rosso 2018) with periodic boundary conditions, as it requires special integration strategies to handle the integration of the Brillouin zone. At first glance, periodic many-body calculations would appear to be intractable because the expansion of one-electron orbitals in terms of Bloch states leads to a large number of orbitals describing the first Brillouin zone  where ψ σ ,nk (G) are the expansion coefficients, is the volume of the primitive cell ( � = [a 1 , a 2 , a 3 ] = a 1 · (a 2 × a 3 ) ), r is the position in real space, G are the reciprocal lattice vectors, σ and n are the spin and orbitals indexes, and k is a vector in the first Brillouin zone (Kittel 2005;Ashcroft and Mermin 1976). Simple approximations to the integration over the Brillouin zone in the exact exchange and other two-electron integrals lead to very inaccurate results, e.g. a straightforward Ŵ-point approximated calculation results in the two-electron integrals being infinite (Bylaska and Rosso 2018;. To overcome these limitations, we have recently developed new methods for generating optimized orbital basis sets, called COVOs . 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. In this work, the COVOs method is extended to periodic systems at the Ŵ-point using the recently developed Filon integration strategy ) for calculating exact exchange energies and two-electron periodic integrals in electron transfer calculations (Bylaska and Rosso 2018;Simonnin et al. 2021), in which the integration of the first Brillouin zone is automatically incorporated.
In addition, present quantum devices are plagued by short coherence times and vulnerability to environment interference, i.e., noise. Although quantum algorithms such as quantum phase estimation can calculate molecular energies with proved exactness, these are not yet viable to run on near-term intermediate scale (NISQ) devices (Preskill 2018;Reiher et al. 2017). 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 development of Variational Quantum Algorithms (VQA) (Peruzzo et al. 2014;Cerezo et al. 2021). Briefly, this class of algorithms strives to find the lowest eigenvalue of a given observable by assuming 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 paper is organized as follows. In "Pseudopotential Plane-Wave Second-Quantized Hamiltonian" section, a brief description of the second-quantized Hamiltonian and one-electron and two-electron integrals with periodic boundary conditions is given, followed in "Algorithm for defining a virtual space with a small CI Hamiltonian" section in which a new class of algorithm for generating a virtual space in which the orbitals are generated by minimizing small pairwise CI Hamiltonians. A complete set of equations for implementing these optimizations is given in "One-electron Virtual and N-Filled Orbitals"-"Matrix elements from the two-electron operators" subsections. Using this new type of virtual space, CI calculations up to 18 virtual orbitals for the ground state energy curve of the LiH molecule in a periodic box are presented in "Results for the Ground State of the LiH Molecule Using Periodic Boundary Conditions" section. LiH is a commonly used test case in quantum computing (Kandala et al. 2017;Low et al. 2019). In "Quantum Computer Calculations for the Ground State of the LiH Molecule Using Periodic Boundary Conditions" section, results from quantum computing simulations using variational quantum eigensolver (VQE) quantum computing algorithms are presented, and lastly, the conclusions are given in "Conclusion" section.

Pseudopotential Plane-Wave Second-Quantized Hamiltonian
The non-relativistic 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 wavefunction 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 wavefunction 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 are described by the proper potentials Ze |r i −R A | , 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 by 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 non-local pseudpotentials, V (A) local and V (A),lm NL , that represent the electron-ion interactions, and the last term is the electron-electron repulsion. (1) Instead of writing the many-body 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. In this form, single particle (electron) creation a † p |0� = |1� and annihilation a p |1� = |0� operators are introduced, where the occupation of a specified state p is defined as |1� and |0� for the occupied and unoccupied orbitals respectively. The second-quantized Hamiltonian is written as  where ψ p (x) represent 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 anti-commutation 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 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. 3 will become prohibitive if all possible Hartree-Fock orbitals are used.
We note the formulae for the one-electron and two-electron integrals in "Periodic One-Electron Integrals using the Pseudopotential Plane-Wave Method", "Periodic Two-Electron Integrals using the Pseudopotential Plane-Wave Method" and "Periodic Ion-Ion Electrostatic Energy using the Pseudopotential Plane-Wave Method" subsections are given in terms of the spatial orbitals rather than spin orbitals. The spin functions α and β are integrated out in the standard way, to involve only spatial functions and integrals (Szabo and Ostlund 2012). Many of the periodic forms presented in the following sections are written in terms of Fourier space using periodic plane-wave basis sets, rather than real space. Descriptions of the plane-wave methods used in this work can be found in the following references (Bylaska et al. 2011(Bylaska et al. , 2011Bylaska 2017;Bylaska and Rosso 2018;Pickett 1989;Ihm et al. 1979;Car and Parrinello 1985;Payne et al. 1992;Remler and Madden 1990;Kresse and Furthmüller 1996;Marx and Hutter 2000;Martin 2004;Chen et al. 2016

Periodic One-Electron Integrals using the Pseudopotential Plane-Wave Method
The one-electron integrals in the pseudopotential plane-wave method can be written as a sum of the kinetic, local pseudopotential and non-local pseudopotential energies (Bylaska 2017).
The kinetic energy can be written as where ψ p (G) and ψ q (G) molecular orbitals in Fourier space. The local pseudopotential energy can be evaluated as where the valence overlap electron density in reciprocal space ρ pq (G) is obtained from taking the fast Fourier transform of its real-space representation, ρ pq (r) = ψ * p (r)ψ q (r) . The local potential is defined to be periodic and is represented as a sum of piecewise functions on the Bravais lattice by where R I is the location of the atom, I, in the unit cell, L is a Bravais lattice vector, and V I local (r) is a radial local pseudopotential for the atom obtained from a Kleinman-Bylander expansion of a norm-conserving pseudopotential (Kleinman and Bylander 1982;Hamann 1989). The local pseudopotential in reciprocal space can be generated by using an ( l = 0 ) spherical Bessel transform.
where j 0 (x) = sin(x) x is the l = 0 spherical Bessel function of the first kind. The non-local pseudopotential energy can be evaluated as where P I nlm (G) is the reciprocal space representation of the nonlocal projector obtained from the Kleinman-Bylander (or generalized Kleinman-Bylander Vanderbilt (1985); Blöchl (1990)) expansion of the pseudopotential, which can be obtained from spherical Bessel transforms.
where T l,m is a real space spherical harmonic or Tesseral harmonic , and j l (x) are the spherical Bessel functions of the first kind of degree l.

Periodic Two-Electron Integrals using the Pseudopotential Plane-Wave Method
The two-electron integrals are written as where the periodic Coulomb and screened Coulomb energies are and where the Fourier representation of the densities are The periodic exchange term in Eq. 7 is approximated by The filter potential is approximated using the cutoff Coulomb kernel from our prior exchange paper based on the Wannier orbitals (Bylaska et al. 2011), written in real-space as where R = |r i − r j | , and N and R cut are adjustable parameters. The design of this cutoff kernel is chosen to remove the interactions between redundant periodic images of Wannier orbitals, because of the long-range nature of the Coulomb potential. Recently, we developed a Filon integration strategy , which showed that filter potential for periodic exchange can be formulated as where V BZ = 2π 2 � is the volume of the first Brillouin zone, and moreover this potential can be approximated by the cutoff Coulomb kernel.
To derive the form of the Eqs. 4 and 7, we compared the results from the "corresponding orbital transformation" developed by King et al. King et al. (1967) (and generalized to periodic boundary conditions, see "Algorithm for defining a virtual space with a small CI Hamiltonian" section and reference (Bylaska and Rosso 2018)) to the results using the one-electron and two-electron integrals for the electronic structure Hamiltonian integrals, H AB = � A |H | B � between two determinants | A � and | B �.

Periodic Ion-Ion Electrostatic Energy using the Pseudopotential Plane-Wave Method
The ion-ion electrostatic energy for a periodic system can be calculated using the Ewald decomposition (Ewald 1921).
where Z I are atom charges, ǫ is a constant (typically on the order of 1) and L is a lattice vector.

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. 3. As previously shown, these new types of orbitals are able to capture significantly more correlation energy than with virtual orbitals coming from Hartee-Fock . The basis of this method is 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 configuration interaction (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 type of orbitals can be formulated as follows: 1 Set n = 1 2 Using the ground state one-electron orbitals for many electron systems, ψ f 1 (r) , ψ f 2 (r) , · · · , ψ f N (r) , 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) then update with a conjugate gradient or similar method while making sure that ψ (n) e (r) is normalized and orthogonal to ψ f 1 (r) , ψ f 2 (r) , · · · , ψ f N (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 finished.
A small CI wavefunction is constructed by varying the top orbitals to produce 3 determinant wavefunctions for the 2N-electron system composed of (N+1) one-electron orbitals, ψ f 1 (r) and ψ (n) e (r) , can be written as a linear combination of 6 determinant wavefunctions, or just 3 determinant wavefunctions for just singlet (or triplet) states.
Using this small CI ansatz, the energies of the system can be obtained by diagonalizing the following eigenvalue equation.
where Note the overlap matrix, S, is the identity matrix for orthonormal ψ g and ψ e . The variation with respect to ψ e (r) can simply be 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, "One-electron Virtual and N-Filled Orbitals"-"Matrix elements from the two-electron operators", provide formulas that can be used to generate the matrix elements in Eq. 9 and the gradients with respect to ψ * e (r) in Eq. 10. We also note the COVOs approach 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 main difference is that the variational space used by COVOs is significantly larger, because plane-wave basis sets are used instead of LCAO Gaussian basis sets used by OVOS. Another difference between the approaches is that the correlation is described by a small CI Hamiltonian for COVOs, and a second-order Hylleraas functional (Hylleraas 1928(Hylleraas , 1929(Hylleraas , 1930(Hylleraas , 1964Koga 1992) for OVOS. The cost to generate COVOs is similar to the cost to generate regular RHF virtual orbitals (just 4 to 9 times more expensive than RHF virtual orbitals). However, because the orbitals are generated one at a time, the resulting electronic gradient is non-Hermitian, which requires more advanced optimizers.

One-electron Virtual and N-Filled Orbitals
The one-electron spin orbitals of (N+1)-state Hamiltonian are where the spatial orbitals and spin functions are orthonormalized. et al. Materials Theory (2023)

Fillings
For the (N+1)-state system, there are six 2N-electron wavefunctions, two of which are singlet, two of which are triplet, and two of which contain a mixture of singlet and triplet character. These wavefunctions can be written as Note that a and b cannot be written as a product of a spatial wavefunction 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 wavefunctions that are separable in spatial and spin functions, and at the same time being eigenfunctions of S 2 and S z .

Incorporating Brillouin Zone Integration
For systems with periodic boundary conditions, the matrix elements for calculating H AB are used with the Bloch states, i.e.
where k 1 , k 2 , . . . are points in the first Brillouin zone, and a ik j (x) and b ik j (x) are the one-electron Bloch orbitals of A and B , where the orbitals in each determinant are taken from the same orthonormal set. The corresponding orbital transformation (King et al. 1967) can be used to generalize for different orthonormal sets. Since the overlaps between orbitals with different k-points vanishes, the one-electron operators can be carried out per k-point (i.e., block by block). The matrix elements, however, for the two-electron operators are in general not block diagonal with respect to the k-points. In cases, where the two-electron matrix elements of the spin-orbitals have a double noncoincidence (King et al. 1967) the matrix elements are again block diagonal, otherwise the matrix elements can be represented as a sum of periodic Coulomb and exact exchange energies, where the Filon integration strategy  can be used to fold in the first Brillouin zone integration present in the exact exchange energies.
To compare the energy states, E i , between calculations with periodic and free space boundary conditions, it is convenient when calculating the � A |H | A � , � B |H | B � , � A |H | B � , and � B |H | A � matrix elements to shift the Hamiltonian by a constant equal to the Ewald ion-ion energy, Eq. 8, plus the charge correction ( Q 2 M 2r s ) for systems with periodic conditions, and a constant equal to the free space ion-ion energy for free-space boundary conditions. It should be noted that the constant shift does not affect energy differences.

Matrix elements from the one-electron operators
The H 1 operator for a periodic system written in reciprocal-space containing N-electrons per unit cell is where the h(G) function/operator is The periodic form of local pseudopotential is given in Eq. 5, and based on Eq. 6 the nonlocal pseudopotential kernel is defined as The one-electron matrix elements between | g � , | e � , and | m � states of the periodic 3 × 3 select CI Hamiltonian can be written using the corresponding orbitals formulas (Bylaska and Rosso 2018) as the following. Song et al. Materials Theory (2023) 7:2 The variation of these integrals with respect to ψ * e (G) are then Song et al. Materials Theory (2023) 7:2 where V I local (G ′ ) and V I NL (G, G ′ ) are given in Eqs. 5 and 11.

Matrix elements from the two-electron operators
The two-electron matrix elements between | g � , | e � , and | m � states of the periodic 3 × 3 select CI Hamiltonian can be written using the Slater-Condon rules or the corresponding orbitals formulas (Bylaska and Rosso 2018;) as the following.
The variation of the two-electron integrals with respect to ψ * e (G) are then et al. Materials Theory (2023)

Results for the Ground State of the LiH Molecule Using Periodic Boundary Conditions
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 planewave calculations used a simple cubic box with L = 15 Å, and a cutoff energy of 100 Ry. The web application EMSL Arrows ) was used to set up and perform the plane-wave calculations. The valence electron interactions with the atomic core are approximated with generalized norm-conserving Hamann (1989) pseudopotentials modified to the separable form suggested by Kleinman and Bylander (1982). The pseudopotentials used in this study were constructed using the following core radii: H: r cs =0.800 a.u and r cp =0.800 a.u.; Li: r cs =1.869 a.u, and r cp =1.551 a.u.. The results for PW FCI calculations of LiH with 1, 4, 8, 12, and 18 COVOs are shown in Fig. 1 and Table 2 in Appendix A. Our code produced the whole energy curves that show inversion symmetry about the central point at R = 7.5 Å, i.e., the energy at the distance (15 Å −R ) produced the same energy found at R with the simple-cubic supercell size of 15 Å due to the periodicity. The average difference error for the 1, 4, 8, and 12 COVOs calculations from the 18 COVOs calculation is 12.9 kcal/mol, 2.7 kcal/mol, 1.0 kcal/mol, and 0.4 kcal/mol respectively. While the error is significant for 1 virtual orbital, the difference is quite small by 4 virtual orbitals, and the error steadily decreases as the number of virtual orbitals is increased. Another measure of the error is the extensivity error. The energy for large R should be the same as the combined energy of the isolated H and Li atoms. The aperiodic PW FCI energies for the dissociated atoms (at R = 7 Å) were found to be -0.66372, -0.68739, -0.68945, -0.68946, and -0.69011 Hartrees for 1, 4, 8, 12, and 18 optimized virtual orbitals, respectively. The sequence of numbers shows the convergence to the combined Hartree-Fock energy of the isolated H and Li atoms which is -0.691388 Hartrees ( E(H) = −0.498883 Hartrees and E(Li) = −0.192505 Hartrees) calculated by the pseudopotential plane-wave method.
In Fig. 2 we compared the total energies from aperiodic (see Table 3 in Appendix A) and periodic plane-wave FCI calculations for the LiH molecule with 1 and 18 correlation optimized virtual orbitals. The energies from periodic plane-wave FCI calculations are lower than the energies from aperiodic calculations from R = 1.3 Å to R = 3.5 Å, while the former are higher than the latter from R = 4.0 Å to R = 7.0 Å. The average difference error for the 1 and 18 COVOs calculations between the aperiodic and periodic energies is 1.2 kcal/mol and 1.3 kcal/mol respectively, which suggests that periodic results agree with the aperiodic ones. However, at large R a significant difference between aperiodic and periodic calculations can be observed. The comparison between the total energies from aperiodic and periodic plane-wave FCI calculations for the H 2 molecule with 8 correlation optimized virtual orbitals is shown in Fig. 3. The difference in the agreement between periodic and aperiodic energies at large R for LiH and H 2 molecules is due to the dipoles in molecules. Since H 2 is a non-polar molecule, there are no dipoles that affect the total energies in both aperiodic and periodic systems while for the polar LiH molecule at large R, the dipoles between Li and H atoms and their images in periodic systems cancel each other in the periodic systems and thus the energy becomes higher than the energies in the aperiodic system.

Quantum Computer Calculations for the Ground State of the LiH Molecule Using Periodic Boundary Conditions
The COVOs optimized orbital basis sets can reduce the circuit depth and complexity for quantum algorithms, opening up applications in chemistry and physics. This is particularly meaningful as current and near-future quantum computers are noisy intermediate-scale quantum (NISQ) devices. These devices are restricted in the number of qubits, qubit connectivity, and fidelity of single-and multi-qubit entangling gates. To effectively utilize such quantum hardware, one must employ algorithms that minimize gate count and withstand noise, which is any undesired internal or external factor that changes the quantum system. Thanks to a grant from Microsoft, we were provided credits for jobs submitted on quantum computers and their corresponding emulators to demonstrate the applicability of the COVOs method in quantum computing for the periodic LiH system.
In this project, we had a few goals in mind in our endeavor to illustrate the applicability of a generated COVOs basis in quantum computing. First, we wanted to compute ground-state energy on an actual quantum computer. We also wanted ground-state energies at various Li-H internuclear distances to measure the performance of a quantum algorithm in different regimes or electronic correlation. Through Microsoft's Azure Quantum cloud computing service (Azure quantum 2022), we had access to the H1-1 quantum computer provided by Quantinuum (a company formed from the combination of Honeywell Quantum Solutions and Cambridge Quantum) (Quantinuum 2022). The H1-1 quantum computer is the latest hardware in Quantinuum's H1 generation of quantum computers with high-fidelity, fully connected qubits. The high fidelity of qubits corresponds to lower errors brought on by noise. In addition to computing ground-state energies at different bond lengths on a quantum computer, we also wanted to measure the effects of noise on the corresponding circuit evaluations. With these goals in mind, we set out on this venture while being mindful of the available computational funds granted to us, which was achieved through a threestage process.
We probed the potential energy surface of LiH at 1.7, 3.0, and 7.0 Å inter-nuclear distances employing VQE, one of the most widely used near-term applications for quantum computing that has successfully been deployed to various kinds of quantum hardware (Peruzzo et al. 2014;Kandala et al. 2017). The VQE method is a hybrid quantum-classical approach in which energies are evaluated on quantum hardware or simulators, and classical computers perform the algorithm to optimize the variational parameters. The repeated evaluation of the quantum circuits can be costly, especially if there are slow convergences of the variational parameters. In the first stage of our calculations, we carried out noise-free VQE simulations to obtain optimal variational parameters. These simulations employed Qiskit's Aer simulator with the simultaneous perturbation stochastic approximation optimizer and EfficientSU2 two-level circuit as the ansatz (Treinish et al. 2022). At this point, we limited ourselves to the 1 COVO basis and three internuclear distances to evaluate for the sake of computational cost. A batch of four circuits, shown in Fig. 4, were evaluated to compute the energy. These circuits all require two qubits and consist of 8 R y and 8 R y rotation gates for the 16 variational parameters in the ansatz, 3 Controlled-X (or Controlled-NOT) gates, and 2 measurements on the two qubits, along with either 0, 1, or 2 Hadamard gates (McMahon 2007). These results reproduced FCI energies to less than a milliHartree for the three geometries. Currently, proposals for robust quantum error correction require qubit numbers and performance that are not yet available via Cloud-based NISQ devices today (Shor 1995;Cory et al. 1998;Reed et al. 2012), so before executing the circuits on the H1-1 quantum computer, we wanted to ensure that noise played a manageable role in computing the ground-state energies. So, for the second stage of the calculations, we used the optimal variational parameters obtained from the noise-free simulation with the noisy Quantinuum emulator that mimics the actual behavior of the Quantinuum H1-1 quantum computer. The error from noise was corrected using the post-processing mitigation technique called the full calibration measurement correction fitter, which measures a circuit with an expected result several times to construct a calibration matrix. The corresponding circuits can be seen in Fig. 5. There was a limit to the number of times a circuit could be executed, which was 500 times, so we performed 500 shots in the simulation. Between the circuits for the energy evaluation and the error mitigation, each complete run consumed nearly 80 credit units, a significant portion of our allocation. The number of credit units required is computed using the following formula: Fig. 4 The two-level circuits consisting of 8 R y and 8 R y rotation gates, 3 Controlled-X (or Controlled-NOT) gates, and 2 measurements on the two qubits, along with either 0 (top), 1 (middle two), or 2 (bottom) Hadamard gates. The values in parentheses represent the 16 variational parameters in the circuit where C is the shot count, N 1q is the number of one-qubit operations, N 2q is the number of two-qubit operations, and N m is the number of state preparation and measurement operations per circuit. After convincing evidence that the error from the noise for this circuit can be well tempered, we performed the last stage of the calculations, where the same energy evaluation and error mitigation technique was performed for 500 runs on the Quantinuum H1-1 quantum computer.
For the three points, energies obtained on the H1-1 quantum computer and simulations reproduced the FCI values to less than 11 milliHartree (6.9 kcal/mol) when corrected for noise (see Table 1 and Fig. 6). These errors are expected to improve with a larger number of circuit runs, which was not feasible then. Given the advertised highfidelity of qubits for the H1-1 quantum computer, error mitigation played a small role in hardware calculations and simulations, reducing all energies by 1-3 milliHartree. Overall, our results are promising, especially considering that another study of LiH showed error mitigated results ranging from ∼10-60 milliHartree along the potential energy surface ).

Conclusion
In summary, we have extended the COVOs method to periodic systems at the Ŵ point using the recently developed Filon integration strategy  for two-electron periodic integrals, in which the integration of the first Brillouin zone is automatically incorporated. We would also like to note that Fig. 1 in reference  illustrates how these integrations can be generalized to include explicit integrations over the first Brillouin zone by breaking it up into patches (see https:// mater ialst heory. sprin gerop en. com/ artic les/ 10. 1186/ s41313-020-00019-9/ figur es/1). For an (N+1)-state Hamiltonian, the method is based on optimizing the virtual orbitals to minimize a small select CI Hamiltonian (i.e., COVOs)

Fig. 5
The four circuits used to construct the calibration matrix for the error mitigation that contains configurations containing all N 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 filled orbitals and the previously optimized virtual orbitals. The method was applied to the simple, but non-trivial, LiH molecule in a periodic system, and we were able to obtain good agreement between the total energies from aperioidic and periodic plane-wave FCI calculations. Also, as shown in Fig. 7, the shapes of the periodic COVOs are basically the same as what is found for the COVOs from aperiodic calculations, which indicates that this extended periodic COVOs method can reproduce the results by the aperiodic Fig. 6 VQE simulations and calculations on the Quantinuum H1-1 quantum computer. Red triangles correspond to energies obtained with the H1-1 quantum computer simulator, while the blue squares correspond to the energies obtained using the H1-1 quantum computer. Open triangles and square are used to represent energies before error mitigation, while the filled shapes are the error corrected values. We note both the quantum computer and simulator results are very good, and error mitigation has very little effect on the overall results. The energies are plotted with the FCI potential energy curve, given by the solid black line. A simple-cubic supercell (L=15.0 Å) was used COVOs method in our previous work . 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 1 kcal/mol.
To test the validity of the periodic COVOs method on a NISQ device, we carried out VQE simulations on the H1-1 quantum computer and its simulator. It was found that the energies obtained using the H1-1 quantum computer were able to reproduce the FCI values to less than 11 milliHartree (6.9 kcal/mol) with a modest number of 500 shots performed; slightly less when corrected for noise. These errors are expected to improve with a larger number of circuit runs. For both simulation and hardware calculations, it was found that error mitigation played a small role, only reducing the energies by 1-3 milliHartree. These results were promising, and open the door to running larger molecular and crystalline systems on NISQ devices in the near future. Fig. 7 The 1 filled RHF orbital and 18 COVOs for the LiH molecule from periodic and aperiodic plane-wave FCI calculations are shown on 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 1.6 Å for LiH. The positive and negative isosurfaces are colored in blue and orange respectively