Optimizing Electronic Structure Simulations on a Trapped-ion Quantum Computer using Problem Decomposition

Quantum computers have the potential to advance material design and drug discovery by performing costly electronic structure calculations. A critical aspect of this application requires optimizing the limited resources of the quantum hardware. Here, we experimentally demonstrate an end-to-end pipeline that focuses on minimizing quantum resources while maintaining accuracy. Using density matrix embedding theory as a problem decomposition technique, and an ion-trap quantum computer, we simulate a ring of 10 hydrogen atoms without freezing any electrons. The originally 20-qubit system is decomposed into 10 two-qubit problems, making it amenable to currently available hardware. Combining this decomposition with a qubit coupled cluster circuit ansatz, circuit optimization, and density matrix purification, we accurately reproduce the potential energy curve in agreement with the full configuration interaction energy in the minimal basis set. Our experimental results are an early demonstration of the potential for problem decomposition to accurately simulate large molecules on quantum hardware.


INTRODUCTION
Electronic structure simulation is an essential tool for understanding chemical properties of molecules. It is a basis for contemporary materials design and drug discovery. Performing accurate electronic structure simulations on classical computers requires a great amount of computational resources. In particular, they grow exponentially with the system size when employing the full configuration interaction (full CI) method, which calculates the exact solution of the electronic Schrödinger equation in a given basis set. Quantum computing, a computing paradigm that leverages the laws of quantum physics, has the potential to deliver scalable and accurate electronic structure calculations [1,2] beyond the reach of classical computers.
Quantum computing technologies are rapidly advancing, and there has been major progress in simulating molecular systems in the last two decades [3][4][5][6][7][8][9][10][11][12][13]. However, simulating the electronic structure of industrially relevant molecular systems on today's noisy, * snowbirth@gmail.com † erika.lloyd@1qbit.com ‡ MPCoons@dow.com § nam@ionq.co ¶ takeshi.yamazaki@1qbit.com intermediate-scale quantum (NISQ) devices [14] will require systematically scalable, robust methods that allow a given problem to be represented by a small number of qubits and shallow quantum circuits. The treatment of electron correlation is also necessary for applying electronic structure calculations to make accurate predictions about the process of a chemical reaction. Calculating the molecular energy remains a challenge, even for small systems, without limiting the number of configurations or the number of electrons active in performing the electronic structure calculations. Without being able to treat correlations in a scalable way, the stringent constraints of NISQ devices will inhibit their ability to perform highaccuracy simulations of larger molecular systems, particularly those systems where the electron correlation is strong.
Progress is being made on these issues as hardware develops. The largest calculation of the total energy including electron correlation to date is the simulation of BeH 2 using six qubits for the six-electron problem [7]. Note further that some of the authors of the present manuscript simulated a water molecule on a trapped-ion quantum computer. Using a small number of electron configurations that are known to contribute significantly to the total energy of the molecule, the energy estimates obtained were within the widely used measure for chemical accuracy (1.5936×10 −3 hartrees) when compared to classically simulated results [9]. Other recent research [13] simulated a chain of 12 hydrogen atoms using 12 qubits arXiv:2102.07045v3 [quant-ph] 28 Sep 2021 and the Hartree-Fock method. This is a system with 12 uncorrelated electrons and set the record in terms of the largest number of qubits used for a chemistry simulation.
To approach the simulation of larger molecules, problem decomposition techniques can be used to decompose a given molecular system into small subsystems, without sacrificing the accuracy of the electronic structure calculation for a wide class of chemical systems [15][16][17][18]. These techniques admit a more compact representation of a molecule, enabling the explicit inclusion of more electrons in calculating correlation energies. Although the amount of reduction in the computational cost and resulting accuracy is dependent on the problem decomposition algorithm and the system being studied, these techniques have the potential to substantially reduce the qubit count requirements in electronic structure simulations [19][20][21][22][23][24][25][26][27].
For our experiment, we look at a ring of 10 hydrogen atoms taking all electrons into account and without using the frozen core approximation. We choose density matrix embedding theory [28,29] as the appropriate problem decomposition method based on previous studies using classical simulations for this system [29], and for the method's success in quantum simulations of the Hubbard model [19]. DMET has been studied for molecular systems, ranging from model systems such as a ring or a lattice of hydrogen atoms [30] to more-realistic organic molecules [29,31]. Despite the strength and wide applicability of DMET-based decomposition, more complex molecular systems might require a different choice of problem decomposition technique.
In this work we demonstrate experimentally how, given a limited fragment size, we can extend the applicability of quantum hardware to larger systems with the aid of classical computations that implement problem decomposition. We do this by applying the DMET method to generate fragments, and solve the electronic structure problem of the fragments using a quantum algorithm. The most appropriate algorithm for current hardware is the variational quantum eigensolver (VQE) [4] whose advantage is enabling shallow circuits that can realistically be executed, at the expense of taking more measurements. We use the established qubit coupled cluster method [32] to generate the parametric ansatz required for running the VQE algorithm. To run our circuits we use a trapped-ion quantum computer [6,8,9], a platform that allows for an efficient implementation of quantum circuits via complete qubit connectivity [33][34][35]. We further employ a density matrix purification algorithm to post-process the experimentally determined results and mitigate residual error. Our results demonstrate the success of an end-to-end pipeline using problem decomposition to accurately solve a molecular problem.

RESULTS AND DISCUSSION
Numerical results of the density matrix embedding theory and the qubit coupled cluster methods We begin by describing our DMET-based methodology for simulating the electronic structure of a molecular system (see Methods section for further details). Consider a system described by a second-quantized Hamiltonian H = H 1-e + H 2-e , where H 1-e is the one-electron interaction and H 2-e is the two-electron interaction of the entire molecule. In contrast to simulating the entire system usingĤ, in DMET, the system to be simulated is divided into small fragments. Each fragment is treated as an open quantum system, entangled with its surrounding environment, or bath. Here we use the mean-field, Hartree-Fock (HF) solution of the entire molecular system as a pre-processing step to find local orbitals that we use to fragment the molecule. Then, the following iterative process, which we call the DMET cycle, initiates. Once the cycle terminates, the electronic structure calculation of the entire molecule is complete. The DMET algorithm is shown in Fig. 1 The DMET cycle begins by constructing the bath orbitals for each fragment. The bath orbitals describe the environment that is active for the electronic structure calculation of the fragment, by virtue of Schmidt decomposition [36]. Note that, if the bath is large, the description of the bath can be greatly simplified. With the simple description of the bath, the Hamiltonian H A for each fragment A, along with its specific bath B, is constructed according to the equation where H 1-e,AB denotes the one-electron interaction within and across the fragment and the bath, H 2-e,A denotes the two-electron interaction within the fragment, µ is the chemical potential, and N A is the number of electrons in the fragment. See Methods section for more details. We use a quantum computer in the DMET calculation to accurately evaluate the minimal expectation value of H A , as well as the number of electrons N A in fragment A. Once both the energy expectation value and the number of electrons for each fragment are computed, we combine them to compute the total system energy, and check for self-consistency. In particular, we choose to compute the sum of the number of electrons in each fragment and check whether the sum is equal to the total number of electrons in the system. If the sum is within a pre-specified range with respect to the total number of electrons, the DMET cycle terminates. If the sum is not within the specified range, we run the DMET cycle again, with the chemical potential µ updated as the difference between the sum and the total number of electrons.
Our explicit example of calculating the electronic structure of a molecule using DMET, is to simulate a ring of 10 hydrogen atoms, H 10 . For this molecule we take advantage of its symmetry to create identical subproblems.   Tables 2 and 3 in Methods section for the numerical values of the gate parameters. (ii.) Post-optimization circuits for YY, required for the classical post-processing of our simulation data. See Table 4 in Methods section for the numerical values of the gate parameters.
This allows us to use a single fragment to solve for the entire system. If this symmetry was not present, we would need to solve each subproblem individually which could be done in parallel. To minimize quantum resources, we choose to divide the molecule into 10 one atom fragments, allocating two spin-orbitals for the fragment and for the bath. This may be compared to a total of 20 spin-orbitals in the simulation of the entire molecule, which shows a large reduction in the problem size: a 20-qubit problem is reduced to 10 two-qubit problems. Classical simulations show that this decomposition reproduces full CI energy within chemical accuracy in all regions of the dissociation curve, except for the repulsive wall. This is consistent with the studies in [29], which also show how correlations can be added to improve the approximation by increasing the fragment size, or by adding an additional self-consistency loop to optimize the correlation potential in the DMET approach. For the points along the dissociation curve that we will explore experimentally, we calculate the total energy per atom of the DMET fragments with a full CI solver, and denote the results under DMET-FCI in Table 1.
To estimate the expectation value of the fragment energy and the number of particles per fragment on a quantum device, we use VQE [4] with the qubit coupled-cluster (QCC) ansatz [32]. For all calculations we use the symmetry-conserving Bravyi-Kitaev transformation [37,38] to transform from a fermion to a qubit basis. The QCC ansatz operatorÛ (τ ) is specified according to the equation where τ k is a variational parameter, n g is the number of multi-qubit Pauli operatorsP k , defined aŝ where n q is the number of qubits and X, Y, Z, and I are the Pauli matrices and a single-qubit identity operator, respectively. While the depth of the circuit rapidly increases as the size of the molecule increases, the QCC ansatz admits a low-depth quantum circuit compared to a widely used unitary coupled-cluster single and double ansatz. The details of the QCC circuit simulations can be found in Supplementary Note 1.
We first consider an ansatz state that is a product state of n q arbitrary single-qubit states, following the method described in the work of Ryabinkin et al. [39]. We evaluate the expectation value of the mean-field Hamiltonian with respect to the ansatz state and use VQE, simulated on a classical computer, to minimize the value. The variational parameters that result from the optimization correspond to the optimal wavefunction for the mean field. We next consider a QCC ansatz operatorÛ (τ ) applied to the previously determined, mean-field optimized state. Note that we aim to minimize the expectation value of the fragment Hamiltonian with respect to the ansatz operator parameters τ . We find whichP k to include in our ansatz by computing the derivative of the expectation value with respect to τ k , which can be computed in a straightforward way when τ = 0. We removeP k terms that have small derivative values. Applied to our example molecule H 10 , we find thatP k of XY and YX have large, identical derivative values. Note that we carry out the computation of the derivatives on a classical computer.
We now investigate the performance of DMET with the QCC ansatz, which we denote as DMET-QCC, on a classical computer. Specifically, we consider the potential energy curve of the symmetric expansion (i.e., increasing the bond length R while maintaining it for each pair of neighbouring atoms) of H 10 . We choose 10 points R = 0.7, 0.85, 1.0, 1.1, 1.3, 1.6, 1.8, 2.0, 2.5, and 5.0 A along the potential energy curve and compare the total energies resulting from the two-qubit DMET-QCC ansatz with DMET-FCI, all calculated classically. The total energies per atom of the H 10 molecule are listed in Table 1 and we include the results obtained from other known methods, such as HF and full CI. The DMET-QCC results almost exactly reproduce the DMET-FCI results indicating that the VQE and QCC methods accurately represent the fragments.
We next describe the simulation of our DMET method on a trapped-ion quantum computer. Instead of variationally optimizing the energies, then running through the DMET cycles, here we focus on the evaluation of the total energy from the quantum simulation for the classically pre-computed optimal parameters. Note that for the DMET cycles, we require XX , Y Y , ZZ , XZ , ZX , Z 0 , Z 1 , X 0 , and X 1 to be simulated, where the subscripts denote the qubit index. See Supplementary Note 2 for details on the Hamiltonian and the DMET energy expressions used in the experiments.
Compiling the circuits for the trapped ion hardware The circuits are executed on IonQ's 11-qubit trappedion quantum computer, which is described in detail elsewhere [40]. In the quantum computer, 15 171 Yb + ions, aligned to form a linear crystal with spacing of about 4 µm, are suspended in a chip trap with a radial pseudopotential frequency of ∼3.1 MHz. We cool the crystal to its motional ground state and use the 11 ions in the middle as qubits, and initialize them to the |0 state. Two of these qubits are used in our experiment. Singleand two-qubit gate fidelities are nominally calibrated to greater than 99.5% and 96.5%, respectively. Counterpropagating laser beams capable of illuminating individual ions are used to implement quantum gates, leveraging the ion-ion coupling mediated by the collective radial motional modes. We read out the quantum state by fluorescing the ions using a detection laser.
We extract the expectation values of the Pauli terms using the three compiled, optimized circuits shown in Figure 1(c)(2)(i.). We use the XX and ZZ circuits for XX and ZZ , and the XZ circuit for XZ and ZX , as the two pre-optimization circuits for XZ and ZX reduce to the same circuit upon optimizing compilation. The single-qubit Pauli terms X 0 , X 1 , Z 0 , and Z 1 , where B i denotes the i-th qubit's expectation value in the B basis, are computed using all of the statistics available from the three circuit executions. For instance, for X 0 , we use the results from both XX and XZ . Likewise approaches are used for X 1 , Z 0 , and Z 1 .

Experimental results
The total energies obtained from the quantum computer (points plotted in blue), along with the classically computed potential energies (lines), are shown in Fig. 2. We calculate the values of the error bars using the boot-strapping method described in the Methods section. The experimentally determined energies nearly coincide with the full CI energies (black line) and DMET-QCC reference values, and agree within the margin of error.
To improve our experimental results, we perform classical post-processing using McWeeney's density matrix purification technique [41] applied to the reduced density matrices (RDM). Its effectiveness in improving energy estimation has been shown in [13] for uncorrelated oneparticle reduced density matrices (1-RDM), and in [10] for two-particle reduced density matrices (2-RDM) of correlated two-electron systems. As in the latter case, ours is a two-electron system whose 2-RDM is the full density matrix P pqrs . Only because of this aspect, we are able to purify our experiment iteratively using the following procedure: P new pqrs = 3(P old pqrs ) 2 − 2(P old pqrs ) 3 until the convergence criterion of Tr(P 2 pqrs −P pqrs ) < 1.0×10 −2 has been met, see Methods section for details.
We emphasize that this purification only works when applied to the 2-RDM for our system because our fragment only has two electrons -thus the 2-RDM happens to be the full density matrix, allowing idempotency to be imposed. To extend this method to larger correlated systems, we would need to purify our 2-RDMs using the more general n-representability conditions [42][43][44]. This requires further study and is outside the scope of this work. Obtaining P pqrs requires us to run an additional circuit, YY, on the trapped-ion quantum computer. The corresponding circuit is shown in Fig. 1(c)(2)(ii.) comprising three single-qubit and one two-qubit native gate operations; see Methods section for details on the optimized compilation. Table 1 and the points plotted in red in Fig. 2 show the post-processed total energies for all points, along with their respective errors. All points are in agreement with the simulated DMET-QCC values, and, with the exception of the repulsive wall and the point R = 1.8 A, all values are within chemical accuracy. Again excluding the repulsive wall, the purification method improves the accuracy of the energy estimate for all points along the potential energy curve, and brings the point R = 2.5 A within chemical accuracy of DMET-QCC and full CI. This post-processing is enabled by the high quality of the two-qubit gates of the hardware; see Supplementary Note 3 for further details.

Discussion
The level of accuracy in our results is due to the ability of the decomposition method to well approximate most of the curve, the success of the VQE algorithm and QCC circuit ansatz to express the fragments, and the small systematic errors of the hardware used in this experiment. We note that the error bars are quite large, which is a result of the limited sampling capacity. The achievable precision will increase as the availability of quantum resources grows. Simulating a molecule of this size was only possible through using problem decomposition to create fragments that would fit the hardware. This demonstrates how we can use classical methods to extend the applicability of quantum hardware that is limited in size, and scale the correlation captured in our fragments as hardware develops.
In order to capture the nature of the entanglement in our simulation, which is the characteristic phenomenon of quantum simulation, we also investigate two types of quantum entanglement. The Hilbert space of the twoqubit system is split into two subsystems: the first and second qubits, and the fragment and the bath orbitals. In both cases, the experimental results are very close to theoretical predictions and clearly show that the states we observe in the present experiment are entangled. See Supplementary Note 4 for more details.

CONCLUSION
In this work, we have experimentally demonstrated an end-to-end pipeline using a PD technique to reduce the size of an electronic structure simulation, making a 20 qubit problem realizable on a trapped-ion NISQ device. Our method involves combining DMET and VQE with a QCC ansatz to compress the quantum simulation circuit for the electronic structure calculation, and circuit optimization techniques that target trapped-ion quantum computers. A density matrix purification method is then applied to mitigate residual errors. All of these steps successfully construct the potential energy curve of the ring of 10 hydrogen atoms, and serves as a proof of concept for the usefulness of problem decomposition to accurately and efficiently represent large molecules. Our results are in agreement with full CI, and almost all points are within chemical accuracy of their DMET-QCC simulated values. This demonstrates that the approach employed herein can describe bond-breaking and bondforming in molecular systems.
Computer-aided molecular design in both materials science and the life sciences is key to our attaining a sustainable future through the accurate prediction of chemical reactions, such as the catalytic reaction of organometallic compounds in advanced materials innovation, enzymatic reactions in the life sciences, and the electrochemical reactions implemented in next-generation batteries. In a previous work by some authors of the present paper [26], it is shown that a quantum simulation of a typical, industrially relevant organometallic compound would require over 2000 qubits. Selecting the appropriate PD technique could result in the reduction of qubit requirements significantly. A reduction by a factor of five for industrially relevant systems is shown in another work [22]. This sort of reduction could place these simulations within the reach of near-future NISQ devices that have hundreds of qubits. The precise reduction in the amount of required resources can be both algorithm dependent and application dependent. For example, if we were to use the DMET algorithm as the PD technique, and if our target system were not symmetric, we would need to run quantum simulations for each fragment that the DMET algorithm would create. If we were to use a larger fragment size for greater accuracy, the quantum resource required to simulate the fragment would also become greater, but would always remain smaller than or equal in size to the case where the entire system were simulated without using any PD.
Our present work experimentally demonstrates the potential of problem decomposition methods as a preprocessing step to reduce the quantum resources required in simulating molecular systems while preserving accuracy. We were able treat all electrons in the system, and include correlations as necessary to produce accurate results using subproblems that fit on hardware. This additional classical component allows us to apply results from quantum devices to a larger family of problems, which we believe will be a valuable tool for using quantum computers to enhance molecular design platforms.

The System Hamiltonian
The second-quantized electronic Hamiltonian can be written aŝ where p, q, r, and s are distributed over all spin-orbitals, and a † p and a q are the corresponding creation and annihilation operators. The terms h pq and pq|rs are the one-and two-electron integrals (in chemists' notation), respectively. The evaluation of these integrals, as well as the Hartree-Fock (HF) and full CI calculations, are carried out using PySCF [45]. The minimal basis set MINAO [46] is used in our calculation.
The electronic Hamiltonian is then transformed into the qubitized form by using symmetry-conserving Bravyi-Kitaev transformation (scBK) [37,38]. The qubitized Hamiltonian can be written aŝ where p, q, r,. . . are distributed over all qubits, and σ α p ∈ {σ x , σ y , σ z , I} acts on qubit p. The transformation of the electronic Hamiltonian into a qubit basis is performed using OpenFermion [47].

Density Matrix Embedding Theory
Let us consider a molecular system divided into fragments, in which a small fragment A with N A states is surrounded by a large bath B with N B states. If the wavefunction of the entire system |Ψ is known, Schmidt decomposition [36] can be applied, and we may write where |α i and |β i denote particular many-body bases, and |χ i is the orthogonal set of |χ i . The states |χ i are the states of bath B; however, the number of states coincides with that of fragment A. This shows that, regardless of the size of bath B, only N A states can be entangled with the fragment. This will reduce the size of the problem drastically for large-sized systems. The exact wavefunction of the entire system |Ψ is the eigenfunction of the Hamiltonian for the entire system H. The Hamiltonian H of fragment A embedded in bath B can be defined using a projection operator P , that is, where It is now evident that the electronic structure of the entire system can be described exactly by that of the fragments and their surrounding baths. The electronic structure calculation of the entire system can be solved using this smaller problem. However, the exact wavefunction of the molecular system is usually not known a priori; thus, introduction of an approximation is necessary. The wavefunction of the entire system obtained by low-level, mean-field theory, such as the HF calculation, would be a straightforward approximation of the exact wavefunction of the entire system. Using a low-level wavefunction to construct a bath and solve the reduced problem employing a high-level theory is the principal idea behind the density matrix embedding theory (DMET).
The orbitals of the entire system are transformed into unentangled occupied orbitals, unentangled virtual orbitals, local fragment orbitals, and bath orbitals. The orbital space of the entire system is then greatly reduced, as only the local fragment and bath orbitals are employed for each high-level DMET fragment calculation. Practically, we need to optimize the embedding of a bath. In DMET, a high-level calculation for each fragment is carried out individually until self-consistency has been attained according to a certain criterion: the sum of the 1-RDM of all of the fragments agrees with that of the low-level one for the entire system. The DMET energy is calculated using the 1-RDM and 2-RDM. The DMET algorithm used in this work, the single-shot algorithm [29], can be described as follows: 1. Calculate the wavefunction of the entire molecular system using a low-level method and then localize the orbitals to fragment the molecule.
2. Construct the bath orbitals so as to include the surrounding environmental effect.
3. Construct the Hamiltonian of a fragment (including environmental effect) and calculate the wavefunction using a method based on a high-level theory.
4. Calculate the fragment energy and the number of electrons for each fragment with the 1-and 2-RDMs from the wavefunction obtained in Step 3.

5.
Repeat steps 2-4 for each fragment and obtain the total energy and the high-level 1-RDM of the entire molecular system. 6. Repeat steps 2-5 until the sum of the number of electrons in the fragments agrees with the number of electrons for the entire system.
The initial step of a DMET calculation is to perform a mean-field HF calculation for the entire molecule. The localized orbitals are obtained by localizing the canonical orbitals from the HF calculation to determine how to fragment the molecules. The Meta-Löwdin localization scheme [48] is used in this work.
After the molecule is divided, the DMET cycle (steps 2-6) is initiated. The first step in the cycle is to obtain the bath orbitalsâ † r (the active orbitals in the electronic structure calculation used for fragments to describe environmental effects), and the environment density matrix of fragment A, D env,A , is calculated as follows: where C represents the molecular orbital coefficients obtained from the mean-field calculation of the entire molecule. Now, the embedding Hamiltonian H emb,A can be constructed (step 3). It can be defined as where the h pq are the one-electron integrals, the pq|rs are the two-electron integrals in chemists' notation, L A is the number of orbitals in the fragment, L B is the number of bath orbitals, L is the number of orbitals in the entire molecule, and p, q, r, and s are general orbital indices. We introduce the chemical potential δµ, which is optimized in the DMET cycle. Once the Hamiltonian H emb,A has been obtained, the electronic structure calculation is performed for fragment A. We employ VQE with the QCC ansatz in this work. Following these calculations, the algorithm constructs the one-and two-particle density matrices, D A and P A , respectively, from the QCC wavefunction Ψ A . The fragment energy E A is calculated (step 4) from the RDMs as Note that only the elements with fragment orbital indices (p ∈ A) are used for calculating the fragment energy.
Here, the 1-and 2-RDMs are defined as and respectively. The number of electrons N A in fragment A are calculated (step 4) as The DMET energy is calculated by summing the fragment energy for each fragment (step 5), which is obtained according to the equation where E nuc is the nuclear repulsion energy. The DMET cycle (steps 2-6) iterates until the number of electrons in the DMET calculation given by converges to the total number of electrons in the molecule N tot . Convergence is achieved by updating the chemical potential δµ according to the equation where a is positive number.
For the DMET calculation for the H 10 ring, the high symmetry of the molecular structure allows us to reduce the calculation further. We calculate only a single fragment with the assumption that all of the fragments have both the same energy and number of electrons. The DMET total energy with only one fragment multiplied by 10 (the number of hydrogen atoms) coincides with the energy when using all fragments; thus, we treat only one fragment in the DMET calculation and simply multiply the energy and the number of electrons by 10.

The Qubit Coupled-Cluster Method
An accurate and affordable description of the correlated wave function |Ψ required to evaluate the energy of fragment A according to Eq. (12) is achieved using the qubit coupled-cluster (QCC) method [32]. Within the QCC approach, a mean-field wave function |Ω is determined and subsequently utilized in a heuristic [39] to construct a unitary operator ansatzÛ (τ ). This operator recovers the missing electron correlation for the mean-field state |Ω and results in the QCC wave function according to the equation A parameterized mean-field wave function is defined as a tensor product of n q single-qubit states and can be expressed as where is the set of 2n q mean-field parameters. Each single-qubit state is then represented in the qubit computational basis as and is characterized by the Bloch sphere polar and azimuthal angles θ j and φ j , respectively [32]. The mean-field energy functional is given by Minimization of Eq. (22) with respect to Γ gives the ground-state, mean-field energy E MF and the corresponding optimal set of parameters Γ opt . The QCC unitary operator ansatzÛ (τ ) takes the form whereP k is a multi-qubit Pauli operator defined aŝ τ k is a variational parameter, and n g is the number of Pauli operatorsP k included in the ansatz. The QCC energy functional is defined as where τ = {τ k } ng k=1 is the set of n g variational parameters. The heuristic screening procedure [39] is utilized to construct the set of operators {P k } ng k=1 appearing in Eq. (23). This heuristic approach [32] relies on the gradient of Eq. (25) with respect to τ k evaluated using the optimal mean-field wave function (i.e., Γ = Γ opt and τ = 0), the form of which is where [P k ,Ĥ] =P kĤ −ĤP k . Equation (26) is quantified for a representativeP k from each group of electron correlation generators contained in the direct interaction set (DIS) [39]. The n g representative generators from the DIS with the largest energy gradient magnitudes are utilized in Eq. (23). Minimization of the energy functional given by Eq. (25) with respect to τ gives the QCC ground state energy E QCC and the corresponding optimal set of parameters τ opt . We note that, in the present work, the optimal mean-field parameter set Γ opt is not relaxed during the minimization of Eq. (25). The parameter set T opt of size 2n q + n g that specifies the ground state QCC correlated wave function |Ψ T opt is formed as the union of the separately optimized mean-field and QCC parameter sets: T opt = Γ opt ∪ τ opt .
The DMET-QCC framework is applied to compute the energy of each point along the H 10 potential energy curve (see Table 1 and Fig. 2). The pre-optimized quantum circuits utilized for the trapped ion hardware experiments are constructed with T opt obtained from classical DMET-QCC simulations of the embedding Hamiltonian for fragment A (see Eq. (11)). At each point R along the potential energy curve,Ĥ emb,A (R) is a two-qubit opera-tor (n q = 2) after scBK encoding and takes the form The inter-atom, spacing-dependent expansion coefficients ofĤ emb,A (R) in Eq. (27) are provided in Supplementary Note 1 for each point along the H 10 potential energy curve. ProjectQ [49] is employed to perform simulation on classical hardware. Further details of the DMET-QCC circuits and simulations can also found in Supplementary Note 1.

Density Matrix Purification
We perform density matrix purification based on McWeeny's purification scheme [41]. This iterative method purifies the 2-RDM P pqrs according to The iteration is continued until the convergence criterion is met. In this work, = 1.0×10 −2 is used. Note that the change in the total energy is smaller than a millihartree when we changed the criterion to 1.0 × 10 −7 . The tensor multiplication of P is defined as where the Einstein summation is implied. Note that this purification technique can only be applied to two-electron systems [10,13]. Although we here consider a 10-electron system, the RDM in our current work can be purified, as the fragment calculation for DMET involves a two-electron system. More analysis of the purification method applied to the experimental results can be found in Supplementary Note 4.

The Bootstrap Method
The total energies and their statistical errors are calculated using an empirical bootstrapping method. We follow the procedure described in a previous work [9]. We start from the state preparation and measurement (SPAM)-corrected histograms for each Pauli operator (XX, YY, ZZ, XZ, and ZX), and construct a distribution of total energies.
The mean and the standard deviation (σ) are computed from the distribution of energies. The procedure of the bootstrapping method is as follows.
1. Draw a random bootstrap sample of the same size as the original dataset with replacement of the data.
2. Construct a new histogram based on step 1.

Compute the expectation value of the Pauli terms
using the new histogram.
4. Repeat steps 1-3 for each Pauli term and obtain all expectation values needed to construct RDMs.
5. Construct 1-and 2-RDMs. 6. Calculate the total energy. 7. Repeat steps 1-6 10,000 times and obtain a distribution of total energies. 8. Calculate the mean and the standard deviation from the distribution of energies constructed in step 7.
We follow this procedure to construct a histogram of possible measurements consistent with the empirical data. The calculated value of σ is represented using error bars. The mean of this distribution is the measured energy, and the 1σ error estimate is its standard deviation. The difference from the previous work [9] is the construction of the RDMs (step 5), which is required for DMET energy calculation. Note that when we perform density matrix purification, we purify the 2-RDM between steps 5 and 6.
We combine the circuit compilation and optimization techniques reported in [9,50,51] with the circuit optimization technique for the trapped-ion gate set shown in Fig. 3 to obtain the final, optimized circuits, amenable to implementation on a trapped-ion quantum computer. The pre-and post-optimization circuits are represented in Fig. 1(c). Their gate parameters appear in Tables 2  and 3.
Specifically, to perform the classical post-processing step, we obtain the Y Y expectation value from the trapped-ion quantum computer. The pre-optimization circuit is the same as that used to calculate any other expectation value, except the measurement basis transformation operations M 0 and M 1 in Fig. 1(c)(2)(i.) are s † h. The post-optimization circuit is shown in Fig. 1(c)(2)(ii.). Note that we use a real degree of freedom in the twoqubit gate, that is, φφ (θ) := exp(−iσ φ σ φ θ/2), as has been used elsewhere for efficient and high-fidelity simulations [9,52]. The gate parameter values of θ for each point along the potential energy curve are shown in Table 4.

DATA AVAILABILITY
Atomic coordinates for the experimental points can be found in Supplementary Data 1. All other data that support the findings of this study are available from the corresponding author upon reasonable request.

CODE AVAILABILITY
The code is available from the corresponding author upon reasonable request.

COMPETING INTERESTS
The authors declare no competing interest.   TABLE 3. Post-optimized parameter specification. Parameters of the gates that appear in the post-optimization circuits of ZZ , XZ , and XX in Fig. 1(c)(2)(i.). All parameters are in radians.  TABLE 4. Post-optimized parameter specification. Parameters of the gates that appear in the post-optimization circuits of Y Y in Fig. 1(c)(2)(ii.). All parameters are in radians.

Author Contributions
Supplementary Figure 1: The qubit coupled cluster circuit ansatz. Generic DMET-QCC circuit ansatz with explicit dependence on the mean-field and QCC parameters. Both single-qubit states are initialized as |0i.
The qubit Hamiltonians provided above correspond to each point along the H 10 potential energy curve. These Hamiltonians are used to obtain the pre-optimized circuit using the qubit coupled cluster method described in the main text. Note that the Hamiltonians shown above are not exactly the same as the Hamiltonian shown in Eq.(11) of the main text. The orbitals are transformed into molecular orbitals using the coe cients C optimized in Hartree-Fock calculations of the subsystems. The DMET energy expression is not the same as that of the Hamiltonians shown above. The total energy of the DMET calculation, shown in Eq. (12) of the main text, is calculated from the 1-and 2-RDMs constructed from the wavefunction. The total energy per atom for the five points in the final DMET iteration can be expressed in an analytical form as follows: The above DMET energy expressions are for the points R = 0.7, 0.85, 1.0, 1.1, 1.3, 1.6, 1.8, 2.0, 2.5, and 5.0 A, respectively. Note that in both the qubit Hamiltonian and DMET energy expression for the point R = 5.0 A, the coe cients for ZZ and hZZi appear as zero due of the precision (i.e., the number of digits) to which we are reporting, but are finite in the actual calculation. The total energy per atom can be calculated using these expressions. Note that when density matrix purification is carried out, the 2-RDM is constructed and purified; thus, these equations cannot be used to compute the total energies. The total energies are evaluated using Eq. (12) of the main text with the purified RDMs.
Note that the expectation value hY Y i is not needed to compute the total energies. When we compute the DMET-QCC energy, we first measure the Pauli operators XX, YY, ZZ, XZ, and ZX in the experiment. An additional experiment is needed that measures hY Y i for post-processing. We also note that the coe cients for the terms (X 0 , X 1 , XZ, ZX) in the qubit Hamiltonians are the smallest, whereas in the energy expression they are the largest. We observe that the simulated expectation values for these terms are all very small, well below the sampling accuracy of our experiment. This means that the sampling error has a large e↵ect on the results. A further study on the impact of the experimental energies and post-processing is shown in Supplementary Note 3.

SUPPLEMENTARY NOTE 3: DENSITY MATRIX PURIFICATION
In our preliminary investigation, 2-RDMs with large errors are largely improved by the density matrix purification. However, total energies of all points in this work are not improved, as the errors in them are very small. To understand the performance of the density matrix purification method in this regime, we study the dependency of the purified energies on the accuracy of the individual expectation values. In particular, we explore hXZi, hZXi, hX 0 i, and hX 1 i, whose ideal simulated values are below the sampling resolution of our experiment. We find that the finite values measured via experiment have a significant e↵ect on both the resulting energy and the success of the density matrix purification.
Supplementary Figure 2 shows a simulated analysis of the purified energy estimates and their dependence on hXZi + hZXi and hX 0 i + hX 1 i. In the ideal case, each term along the axes is equal, and contributes equally to the total energy. The experimental values are indicated by the blue dot, and the colour bar represents the error value in the purified energy, with chemical accuracy indicated by the yellow-green-blue region. The area within the dashed red lines shows the region where the purification is e↵ective in improving the experimental results. Here we show two cases. In Supplementary Figure 2(a), for the specific case of the point R = 2.5 A, where the purification improves the experimental results (the blue dot is within the dashed red lines), achieving chemical accuracy. Supplementary Figure 2(b) is an example at the point R = 0.7 A, where the purification process does not improve the energy estimate, and the resulting di↵erence is within the error range of the experimental data. The rest of the experimental points follow this trend. The range where the purification will improve the experimental result (the width between the dashed red lines) varies between (a) and (b), as it is dependent on the relative magnitude of the coe cients in the energy expression for these terms, and the accuracy of the rest of the expectation values.
We investigate how the accuracy of hY Y i, for which a two-qubit circuit is used, a↵ects both the experimental energy and the post-processed energy by manually changing the hY Y i expectation value while keeping the other expectation values fixed. We use the point R = 2. 5 A, which allows us to succeed in improving the total energy value by employing density matrix purification. The analysis of the energy dependence on hY Y i is depicted in Supplementary Figure 3. The experimental values are indicated in blue, the theoretical DMET-QCC values are shown using a black line, and the post-processed (purified) values are plotted in red. The errors for the experimental, post-processed, and relevant experimental hY Y i values are calculated using the bootstrapping method explained in the main text, and are indicated by their respective shaded regions. Finally, the energy values within chemical accuracy of the theoretical DMET-QCC results are shown within the yellow region.
We find that the experimental values that are not purified do not depend on hY Y i. However, we observe that the value of hY Y i a↵ects how the energy is postprocessed: in Supplementary Figure 3(a), we see that the post-processed energy falls outside of the chemical accuracy range when hY Y i is smaller than 0.73 or larger than 0.93. The present analysis shows that density matrix purification with a large error in the hY Y i expectation value will not succeed. In our experiment, however, the trapped-ion hardware provides us with a hY Y i expectation value of 0.859 ± 0.005, which is su ciently accurate for it to bring the post-processed total energy within chemical accuracy of the theoretical DMET-QCC values. The relevant range is shown magnified in Supplementary Figure 3   A for illustration. The ideal total energy per atom in the minimal basis calculated with DMET-QCC is shown using a black line, and the surrounding region shaded in yellow indicates the range of chemical accuracy for these simulated values. The blue line with circles plots the energy as determined by the ion-trap experiment, and the red line with circles plots the post-processed experimental energy. We determine the errors in the calculated values using the bootstrapping method described in the main text, and indicate the error ranges using the shaded regions along the curves. Similarly, we include the experimental value and error for the hY Y i term using the green dotted line and corresponding shaded region. (b) A magnification of the relevant experimental range for the hY Y i term and the resulting energy improvement with post-processing applied for the point R = 2.5 A chosen for this illustration.

SUPPLEMENTARY NOTE 4: ENTANGLEMENT ENTROPY
In this section, we describe the computation of two types of entanglement entropies. The first is the entanglement between di↵erent molecular orbitals (MO) and the second is the entanglement between fragment and bath.
The embedding Hamiltonian given by Eq. (11) of the main text is obtained by projecting the total system onto a fragment described by using the basis of fragment and bath orbitals. From this expression, the orbitals are transformed into the MO basis by using the MO coefficients C optimized in the HF calculation of the subsystem. The Hamiltonians used in the experiment are in this MO basis. Therefore, the entanglement between different MOs corresponds to the entanglement between different qubits. One of the motivations for looking at this quantity is to check whether the ground state of the fragment Hamiltonian is an entangled state, despite the fact that the optimized circuits for evaluating the embedding Hamiltonian contain only single-qubit gates. The other type of entanglement entropy, between the fragment and the bath, will capture the structure of the Schmidt decomposition given by Eq. (6) in the main text and, therefore, the properties of DMET more directly.
We consider the entanglement entropies in the fermion picture. The elements of the density matrices we obtain via experimentation are ha † p a † r a s a q i and ha † p a q i, where the indices p, q, r, and s take values in the spin orbitals {1↵, 1 , 2↵, 2 }. In order to split the Hilbert space based on orbitals, we describe the density matrix ⇢ in the basis of |1↵, 1 , 2↵, 2 ih1↵ 0 , 1 0 , 2↵ 0 , 2 0 |. The RDM in our analysis is obtained by taking the trace over the degrees of freedom for orbital 2: It is straightforward to express the RDM in terms of fermionic operators: Note that o↵-diagonal elements in Eq. (22) are zero due to the conservation of spins and electron numbers. These symmetries are guaranteed once the scBK transformation has been performed. Therefore, experimental noise will not change the structure of the RDM. Note that while there is a systematic way to expand the two-qubit RDM into the four-qubit RDM in the sense that there is no additional experimental information required, the values of the entanglement entropies between these two expressions are di↵erent.
The entanglement entropy is defined by the equation where p i are the eigenvalues of the RDM. First, we choose a basis such that 1↵ and 1 correspond to the first qubit and 2↵ and 2 correspond to the second qubit, both after scBK transformation. The results show that the entanglement entropies between the first qubit and the second qubit increase as the bond lengths increase. We perform a transformation such that 1↵ and 1 represent spin-up and spin-down states of the fragment orbitals, respectively, while 2↵ and 2 represent those of the bath orbitals. To compute the entropies in the fragment-bath basis, we transform the basis into the tensor product of the fragment and the bath Hilbert spaces. In this case, the entanglement entropies decrease as the bond lengths increase. The experimental results are very close to the theoretical predictions for both cases. In the MO basis picture, the ground state in the absence of interactions is a product state (i.e., an HF state). Therefore, the entanglement entropies should be zero. Two-body interactions generate entanglement between di↵erent orbitals. Intuitively speaking, the Hilbert space is divided based on the energy levels (i.e., the Slater determinants) in this picture. The fact that the entanglement entropies in the MO picture are small reflects the fact that the ground states are close to the HF states.

Experiment Theory
On the other hand, in the fragment-bath orbital picture, the the Hilbert space is split depending on the orbitals' position in real space. Therefore, the ground state in the fragment-bath basis is highly entangled even when there is no interaction. More explicitly, the MO 7 coe cients C (the transformation matrix between the fragment-bath basis and the MO basis) for R = 0.7 A is