A Hybrid Classical/Quantum Approach for LargeScale Studies of Quantum Systems with Density Matrix Embedding Theory
Abstract
Abstract
Determining ground state energies of quantum systems by hybrid classical/quantum methods has emerged as a promising candidate application for nearterm quantum computational resources. Short of largescale faulttolerant quantum computers, smallscale devices can be leveraged with current computational techniques to identify important subspaces of relatively large Hamiltonians. Inspired by the work that described the merging of dynamical meanfield theory (DMFT) with a smallscale quantum computational resource as an impurity solver [Bauer et al., arXiv:1510.03859v2], we describe an alternative embedding scheme, density matrix embedding theory (DMET), that naturally fits with the output from the variational quantum eigensolver and other hybrid approaches. This approach is validated using a quantum abstract machine simulator [Smith et al., arXiv:1608.03355] that reproduces the ground state energy of the Hubbard model converged to the infinite limit.
[inlinelist,1]label=(), \pdfstringdefDisableCommands
I Introduction
Many current computational studies of strongly correlated chemical or material systems involve a high accuracy simulation of their electronic structure. A ubiquitous theme in classical computational methods is to identify a relevant subspace of the full Hamiltonian such that a high fidelity simulation in this subspace results in a wavefunction with a high degree of overlap with the true full system wavefunction. Examples of this methodology are active space methods, popular in quantum chemistry, and embedding or impurity problems, prevalent in studies of strongly correlated materials. The embedded problem commonly requires a fullconfiguration interactionlike solution (FCI, an exact diagonalization of the embedded Hamiltonian) to significantly improve the description of the electronic structure. Unfortunately, FCI scales exponentially in the number of orbitals, restricting its application to small problems.
The integration of active space or impurity methods with quantum computation offers a route to improving wavefunction approximations by increasing the size of the active space to include relevant lowenergy states. Bauer et al. Bauer et al. (2015) recently proposed the integration of dynamical meanfield theory (DMFT) and the quantum phaseestimation algorithm which was demonstrated in a proofofconcept study by Kreula et al. Kreula et al. (2016). In a similar vein, Bravyi and Gosset studied the complexity of impurity methods integrated with quantum computation Bravyi and Gosset (2016). An integration of complete active space selfconsistent field (CASSCF) and the phaseestimation algorithm is an example of one such integration that has already been suggested Reiher et al. (2016).
The variational quantum eigensolver (VQE) is a hybrid classical/quantum algorithm to approximate the ground state eigenvalues and density matrices of a Hamiltonian Peruzzo et al. (2014); McClean et al. (2016a); McClean (2015); McClean et al. (2016b); Wang et al. (2015); Wecker et al. (2015). The VQE is a classical optimization loop invoking a quantum abstract machine (QAM) Smith et al. (2016) used for state preparation and operator measurement. Recent work investigating this approach yielded quantum circuits for ansatz state preparation that are short and demonstrated that the algorithm is potentially robust to standard noise models McClean et al. (2016b); Sawaya et al. (2016).
The VQE using a wavefunction ansatz that is exponentially expensive on a classical computer, such as unitary coupledcluster McClean et al. (2016a) or approximate Hamiltonian evolution Wecker et al. (2015), offers a potentially advantageous alternative to other high fidelity solvers (FCI, density matrix renormalization group, etc.) for active space calculations Reiher et al. (2016) and embedding schemes Bauer et al. (2015); DallaireDemers and Wilhelm (2016). In general, the VQE can be integrated into any computational chemistry methodology that requires a highfidelity solver with an interface requiring the oneparticle and twoparticle reduced density matrices ( and RDM, respectively) Helgaker et al. (2014); Werner and Knowles (1985); Gidofalvi and Mazziotti (2007).
Density matrix embedding theory (DMET) is a powerful alternative to DMFT that replaces the problem of finding the impurity’s local Green’s function with determining the embedded problem’s RDM Knizia and Chan (2012, 2013); Wouters et al. (2016); Bulik et al. (2014a, b); Bulik (2015). DMET maps the problem of finding the ground state of a large particle system onto many interacting impurity problems. This is achieved by selfconsistently matching the RDM from a lowlevel wavefunction for the entire system with the RDMs of the impurities calculated with a highlevel technique. Though it is not necessarily clear that such an iterative scheme would result in accurate local RDMs, DMET has been demonstrated to be highly accurate for one and twodimensional lattice models Knizia and Chan (2012); Bulik et al. (2014a), extended chemical systems Zheng and Chan (2016); Bulik et al. (2014b); Knizia and Chan (2013), and for ground state and transition state molecules Wouters et al. (2016). Furthermore, spectral functions (on the realfrequency axis) are accessible from DMET without any bath discretization error allowing for the construction of arbitrary dynamic correlation functions Booth and Chan (2015). These studies demonstrate the amazing quality of DMET with modestly sized impurities.
A key component of DMET is the embedded problem solver. Ideally, this solver produces a highaccuracy RDM for the embedded problem such that the fragment piece of the meanfield RDM can be matched to the correlated RDM for the fragment. In this work we propose the use of VQE for the embedded problem solver as an alternative to FCI or DMRG routines. This change potentially allows for {enumerate*}
the study of a larger embedding Hamiltonian, allowing access to longer range correlation functions, and
generally accelerates the DMET algorithm by quickly finding the solution to the embedded problem on a quantum computer. As a proof of principle, the ground state energy of the repulsive Hubbard model is determined for small rings and the thermodynamic limit.
Ii Density Matrix Embedding
In this section we review the salient features of density matrix embedding and discuss how the method naturally incorporates the use of the variational quantum eigensolver as an embedded Hamiltonian solver. References Wouters et al. (2016); Bulik (2015) provide a more detailed review on density matrix embedding while Bulik et al. (2014a, b) present alternatives to the selfconsistent procedure that appears in DMET.
ii.1 Preliminaries
For a large system , its wavefunction can be arbitrarily bipartitioned into a fragment (or impurity) and bath (or environment). Some examples of a fragment piece include a single or multiple sites of a lattice model or localized atomic orbitals corresponding to a piece of a larger molecular basis set. The total wavefunction is naturally expressed in the tensor product of basis states of the fragment and the bath, , which has a linear dimension of where (resp. ) is the dimension of the Hilbert space of the fragment (resp. bath). Utilizing the Schmidt decomposition, an eigenstate of can be written as a sum over the tensor products of the Schmidt basis states^{1}^{1}1This Schmidt basis is constructed by the singular value decomposition of the coefficient tensor. Wouters et al. (2016); Peschel (2012); Nielsen and Chuang (2010).
(1) 
This reformulation demonstrates that there is a local fragment basis and environment basis of the same size that produces an equivalent representation of . Without loss of generality, we can assume the Hilbert space of the fragment is smaller than the environment. In the Schmidt basis, no matter how large the bath, a fragment can only entangle with bath states. The Schmidt states can be used to project the Hamiltonian into a combined impurity/bath basis that necessarily has the same ground state as the original Hamiltonian but is significantly smaller in size:
(2) 
For a general largescale quantum system the wavefunction can only be computed approximately. This constraint naturally leads to the fundamental approximation in DMET: the embedded Hamiltonian is approximated by constructing bath states from the Schmidt decomposition of the ground state of an approximate quadratic Hamiltonian for the total system. As a consequence, the embedded Hamiltonian now contains an interacting fragment embedded in a noninteracting bath. The exact embedding Hamiltonian is approximated by matching the fragment’s RDM with the lowlevel meanfield density matrix of the system by varying the embedding potential that appears in the quadratic Hamiltonian Knizia and Chan (2012).
ii.2 The LowLevel Quadratic Hamiltonian
For an arbitrary chemical system, the Hamiltonian contains terms that are at most four creation/annihilation operators corresponding to one and twoparticle interactions:
(3) 
In order to vary the form of the bath states used when constructing the embedded Hamiltonian, a potential is added to the full chemical Hamiltonian. The augmented total system Hamiltonian is then approximated using the Hartree–Fock method to find an optimal singleparticle basis:
(4) 
ii.3 Bath Oribtals from a Slater Determinant
An alternative way to construct the fragment and bath embedding states that is similar to the Schmidt decomposition in (II.1) is a contraction over the bath basis:
(5)  
(6)  
(7) 
Orthogonalizing into a set of states gives:
(8) 
Note that while each , there are no more than of them. With this form for the embedding basis, one simply needs to find the form of the bath states constructed from the meanfield wavefunction. This is achieved by projecting the occupied states with overlap on the fragment onto the environment orbitals and normalizing.
There are no formal restrictions on the lowlevel wavefunction representing the total system , though it is desirable that one can quickly construct an approximate eigenfunction of . Alternative approximate wavefunctions that have been used are the antisymmeterized geminal power wavefunctions Tsuchimochi et al. (2015) that introduce some correlation into the bath, products of Bogoliubov quasiparticle states Zheng and Chan (2016), and a Slater determinant Knizia and Chan (2012). The following is a derivation of the embedding basis from a Slater determinant wavefunction. This lowlevel wavefunction for the full system can be expressed as the product of fermionic operators in the occupied subspace acting on the true vacuum:
(9) 
The single particle orbitals can be constructed from linear combinations of site orbitals, , in the lattice:
(10) 
The matrix corresponds with the Hartree–Fock transform from local lattice spinorbitals to the basis that minimizes the Hartree–Fock Hamiltonian of the system. A rotation to the fragment and bath basis can be constructed in a similar manner to the canonical orthogonalization of atomic orbitals that appears in the Hartree–Fock procedure applied to molecules Szabo and Ostlund (1989); Helgaker et al. (2014).
The overlaps of the occupied orbitals projected onto the fragment space are calculated using the fragment projection operator acting on the occupied set of orbitals indexed by and :
(11) 
The eigenvectors of this overlap matrix correspond to an orbital rotation that forms the fragment and bath basis. Diagonalizing with a unitary gives
(12) 
Naturally, the eigenvalues of the overlap matrix are between and . Zero eigenvalues correspond to occupied states that have zero overlap on the fragment. The bath states can be constructed by projecting the occupied orbitals with overlap on the fragment into the bath:
(13) 
Finally, the rotation to the Schmidt basis is defined by the direct sum of the bath transformation and an identity matrix of dimension equal to the number of orbitals on the fragment, :
(14) 
The core states with zero overlap on the fragment are eliminated from the embedded Hamiltonian by including their interaction at the meanfield level similar to generating an active space interacting with a frozen set of core states Hosteny et al. (1975).
ii.4 Embedded Hamiltonian
The embedded Hamiltonian is constructed by projecting into the fragment and bath states
(15) 
which, for a Slater determinant wavefunction, corresponds to a single and doubleparticle integral transform Knizia and Chan (2013):
(16)  
where  
(18)  
(19) 
In the above transformation is the interaction with the nonentangled nonoverlapping core states, is a matrix from (14) whose columns are the transformation vectors to the embedding basis, and and are the one and twoparticle integral tensors defined in the lattice basis.
ii.5 Determining the Optimal Embedding Potential
The key step in the DMET procedure requires improving the bath states so they approximate the bath states of the true embedded Hamiltonian. Schmidt bath state tuning is achieved by varying the potential that was added to the system Hamiltonian in Sec. II.2. Variation in this potential is the connection between the embedded Hamiltonian and the quadratic Hamiltonian representing the entire system. The total DMET procedure can be summarized as follows:

Find the ground state Slater determinant wavefunction for the entire system with the HartreeFock Hamiltonian augmented with the embedding potential.

For all fragments, use to construct the embedding basis and project the system Hamiltonian into each embedded Hamiltonian. Solve for the RDM and RDM with a highlevel method.

Adjust the embedding potential in the meanfield Hamiltonian such that the RDM of the impurity calculated with the highlevel wavefunction matches the RDM calculated from the quadratic Hamiltonian.

Repeat until particle conservation and the embedding potential do not change between iterations.
Step 3 is carried out by minimizing the squared norm between the fragment piece of the meanfield RDM and the fragment RDM computed from the ground state wavefunction for the embedded Hamiltonian:
(20) 
In Wouters et al. (2016) this cost function is known as the fragmentonly density matrix fitting. Alternatives such as full matching of the RDMs in the embedding basis or matching the electron densities on the fragment have also been studied Bulik et al. (2014a, b).
The exponentially scaling FCI has prompted studying alternative methods for finding the correlated ground state RDM Wouters et al. (2016); Zheng et al. (2016). Though this may not lead to an optimal approximation to the embedded Hamiltonian fragment and bath states it does allow for the study of significantly larger fragments–important for assessing longer range correlation functions and more accurate local expectation values Zheng and Chan (2016). The ability to study larger fragments at higher fidelity has many implications for accurate simulations of materials and chemistry. For example, if one were studying a reaction center in a metalorganic framework a fragment would likely need to contain the metal site and ligand orbitals for an accurate description of the energy and properties. The VQE enters in step when solving for the RDM and RDM of the embedded Hamiltonian. DMET’s need for the RDM of the embedded Hamiltonian makes VQE an ideal quantum algorithm for this use case.
Iii The Unitary Cluster Ansatz and Trotterization in VQE
The VQE is a functional minimization scheme that leverages fast construction of the wavefunction to find expectation values of operators O’Malley et al. (2016); McClean et al. (2016a). The energy of the system is minimized by varying over parameters for the wavefunction ansatz (21):
(21) 
The expectation value is determined by summing the expectation of each term in the Hamiltonian:
(22) 
For a general quantum chemical Hamiltonian, scales quadratically with respect to basis set size. In order to evaluate the Hamiltonian expectation value, the wavefunction is prepared many times with a subsequent measurement of the operator. In this work the unitary coupled cluster (UCC) state is used as the function for the energy functional though others could be considered such as the variational adiabtic ansatz.
UCC is a manybody expansion wavefunction ansatz parameterized by the cluster coefficients associated with each generator. We use the UCCsinglesdoubles ansatz corresponding to antiHermitian generators performing single and doubleparticle excitations and deexcitations:
(23)  
where  
(25)  
(26) 
In quantum computing, the antisymmetric fermionic creation/annihilation operators are represented by distinguishable qubits. We utilize the Jordan–Wigner transformation (JW) for mapping fermionic creation/annihilation operators to Pauli spinoperators that preserve the anticommutation relations and parity of the second quantized operators (27) Whitfield et al. (2011). It should be noted that the representation of second quantized operators in Pauli matrices by the JW transform scales linearly in the number of Pauli terms^{2}^{2}2There are alternatives such as the Bravyi–Kitaev transformation that have logarithmic scaling when representing second quantized operators Tranter et al. (2015); Seeley et al. (2012). as the basis index is increased. The maps are defined as
(27) 
where
(28) 
For each shot of the VQE algorithm the wavefunction of the system—according to a UCC ansatz—is constructed by first preparing the Hartree–Fock reference, in the number of gates required to prepare the state where is the number of orbitals, then evolving the initial wavefunction according to (23), using in the number of gates for the exponentiation of the unitary generators.
The  and RDM are measured after an optimal set of parameters are determined by measuring expectation values corresponding to the RDM:
(29) 
The RDM is obtained by contraction from the RDM:
(30) 
The fragment RDM of the total system is matched to the fragment piece of the RDM obtained from the VQE algorithm by varying the embedding potential in (3).
iii.1 Commutation Relations for Generators in Unitary CoupledCluster
Unlike the classical cluster operators that commute at all orders Crawford and Schaefer (2000) the unitary coupledcluster operators (25)–(26), henceforth denoted as , do not commute with each other at all orders. Therefore approximating the exponential by a particular order of the Suzuki–Trotter decomposition Suzuki (1994); Trotter (1959) may be important for constructing accurate distributions. The antiHermitian generators do not commute because the deexcitation operators do not commute with the excitation operators. This can be seen for the first order operators. Here we use to denote the particle space and to denote the hole space.
(31)  
(32) 
The first and last commutators vanish while the middle two terms survive when there is an overlapping index in and :
(33)  
(34) 
As a result, Trotter error must be considered when generating circuits for UCC state preparation. In Sec. IV we discuss the performance of UCCSD with Trotter orders and with Trotter steps beyond .^{3}^{3}3The term Trotter order refers to the order of the series approximation to the exponential of two noncommuting operators. The term Trotter steps refers to the number of slices that each order is broken into in order to minimize the the correction term.
Iv Results and Discussion
Exact  

UCCSD  
DMET(1)ED  
DMET(2)ED  
DMET(1)UCCSD 
We considered the onedimensional Hubbard model (35) with repulsive interactions and antiperiodic boundary conditions as a representative test system for DMET.
(35)  
The DMET selfconsistency loops were run using the QCDMET Wouters (2015) code with an interface to PySCF Su (2016) for exact diagonalization of the embedded Hamiltonian. The VQE algorithm for solving the embedded Hamiltonian used Quil and Rigetti Computing’s pyQuil package to construct the UCCSD ansatz and evolved the wavefunction on a noiseless quantum virtual machine Smith et al. (2016). We describe the programming environment in more detail in Appendix VII.1. All UCCSDVQE runs started with a secondorder Møller–Plesset perturbation theory guess for the cluster amplitudes and required about 20 iterations of BFGS Nocedal and Wright (2006) with the gradient numerically approximated. For singlesite DMET calculations the lowlevel and highlevel fragment density matrices were equivalent and thus the DMET loop only involved setting the chemical potential such that the number of electrons in each fragment summed to the total system’s electron count. For twosite models the fragment RDM was generally not equal to the fragment Slater determinant RDM and thus required a numerical search for the optimal embedding potential.
The noncommuting nature of the cluster operators do not drastically effect the accuracy of the UCCSD ansatz. Figure 1 plots the convergence of the energypersite as the Trotter order and Trotter steps are increased for . Trotter order and Trotter steps indicate the structure of the wave function ansatz for UCCSD. For example, the first order Trotter with Trotter steps produces an ansatz of form:
(36) 
where and correspond to the single and double anitHermitian generators with their coefficients which are described in Eq. (25) and Eq. (26). The persistence of the energy gap as Trotter order and Trotter steps are increased demonstrates that the single reference nature and the truncated set of generators are the main sources of error.
The DMET methodology introduces three sources of error: {enumerate*}
the fragment size and thus the number of noninteracting bath states that can entangle with the fragment,
the embedded Hamiltonian solver, and
the pieces of the embedded Hamiltonian RDM that are chosen to match for updating the embedding potential. The effects of (i) and (ii) are determined by studying small lattices and comparing the ground state energy against the FCI solution and the full UCCSD solution. Table 1 compares the ground state energypersite of a foursite ring at various interaction strengths between the FCI solution, the UCCSD solution, and various fragment sizes for DMET. The DMET solutions are labeled by DMET(), where is the number of sites considered in the fragment, followed by an acronym for the type of embedded solver–i.e. ED corresponds to FCI. The accuracy of the UCCSD solution to the foursitering lattice generally decreases as is increased demonstrating the expected result that singles and doubles generators cannot parametrize the full unitary group. A foursite lattice requires eightspinorbitals and thus eight qubits. The UCCSD ansatz contains 14 cluster amplitudes when restricting the wavefunction to the singlet subspace Bulik et al. (2015); Scuseria et al. (1988). Under a greedy parallelization of commuting instructions for breaking the unitary evolution into time steps, a single state evolution at first order Trotter involves time steps where multiple qubit operations are performed at each step. Each time step has an average of 0.8692 onequbit gates and 0.8101 twoqubit gates. The coupledcluster reference is usually chosen in a basis where the oneparticle piece of the Hamiltonian is diagonal. For the Hubbard model this involves a Fourier transform on the oneparticle and twoparticle integral tensors resulting in onebody terms and twobody terms where is the number of lattice sites. The twobody terms are proportional to .
A single fragment site results in a trivial twoorbital embedded Hamiltonian. For this system, UCCSD + VQE is equivalent to FCI, CISD, and CCSD solvers. Considering a twosite DMET at and produced energiespersite of and , respectively. Though not exactly equal to the DMET(2)ED result, the error is on the order of the difference between the UCCSD solution for the full lattice and the FCI solution for the full lattice.
In Figure 2, we plot the energypersite for a site lattice calculated with DMET and the Bethe ansatz Lieb and Wu (2003); Shiba (1972). The top subfigure demonstrates that DMET(1)UCCSD method converges to the correct solution at high. The bottom subfigure describes the error for DMET(1)ED, DMET(1)UCCSD, DMET(2)ED, and DMET(4)ED relative to the Bethe ansatz. As expected the quality of the solution improves as fragment size is increased. An site fragment requires qubits and thus we are restricted to small fragment sizes when performing a classical simulation.
V Conclusion
The hybrid classical/quantum computation model naturally fits with current stateoftheart methods for studying the electronic structure of molecules and materials. Many of the current methods involve selecting an active space or embedded piece of the Hamiltonian to treat at higher fidelity. These techniques are especially important when describing the electronic structure of correlated materials and molecules such as metaloxides, heterogeneous catalysis, multielectron redox reactions, and catalysts. In contrast to exact diagonalization solvers commonly used in impurity or active space problems, quantum computation offers a technique to treat these Hamiltonians with high accuracy in polynomial time.
The first steps towards integrating a quantum computational resource with embedding schemes was outlined in Bauer et al. (2015). That work integrated quantum phaseestimation with DMFT as the impurity solver. Their results indicate that a modestly sized quantum computational resource of about 100 qubits can be used to simulate a relevant physical problem. Though their work describes DMFT as the impurity methodology they do discuss the use of alternative embedding schemes along with alternative algorithms to use on the quantum computation side.
In this work we have elaborated on one alternative in particular: DMET integrated with the VQE. DMET works by mapping the problem of finding the particle ground state of a large system to many smaller interacting impurity problems. Therefore, DMET is applicable to very large scale problems be it either molecules or materials. Unlike DMFT that requires the frequency dependent twoparticle Green’s function, DMET requires the RDM of the embedded problem with high accuracy. This requirement makes DMET and the VQE a natural pair. The RDM and RDM are computed at each iteration of the VQE optimization.
The VQE algorithm is a general functional minimization routine that allows for the use of any quantum state preparation method. In this work we utilized the UCC ansatz and demonstrated that despite the noncommuting nature of the generators, the error for the ansatz comes from the truncation of the series and not Trotter error. Therefore, at least in the case of simulating local interaction models, loworder Trotter approximation and Trotter number are sufficient for accurate state preparation.
Current classical techniques for describing quantum chemistry or correlated materials are naturally accelerated with a hardware quantum processing unit (QPU). Many of the classical methods that require a RDM or RDM of an active space can be integrated with the VQE algorithm or other hybrid classical/quantum techniques. Further integration can lead to methodologies to study large quantum systems with higher resolution than what is currently possible with state of the art classical methods.
Vi Acknowledgments
The author acknowledges useful discussions of VQE and DMET with Jarrod McClean and members of the Rigetti Computing software team: Will Zeng and Robert Smith. A special thanks to Robert Smith for writing the Appendix on quantum programming and detailed editing of the draft.
Vii Appendix
vii.1 Details of the Quantum Programming Environment
In this work, execution of quantum programs was done using the quantum abstract machine and Quil, its quantum instruction language Smith et al. (2016). In this section, we outline elements of this quantum programming environment, including examples of essential techniques used to compute the  and RDMs.
The VQE is, at its core, a classical optimizer around a quantum execution unit. In particular, a classical optimization loop produces the next parameters of (21) for minimization. Mathematically, parameterizes a wavefunction ansatz, but computationally, they are regarded as parameters to generate a new quantum program. These programs are run on a quantum execution unit modeled as a restricted quantum abstract machine^{4}^{4}4The machine is one with quantum state , classical state , static gates , parametric gates , program , and program counter . with
(37)  
(38) 
where the gates are understood to be the collection of operators acting on each qubit combination of the quantum state . The classical memory is as large as —the number of qubits—to hold the basis measurements for each loop of the VQE, which are collected to approximate . Alternatively, by using a quantum virtual machine, the amplitudes of can be observed in the computational basis and can be computed directly.
The object of interest to simulate is usually specified in terms of sums of fermionic operations, such as the Hamiltonian in (16). After a fermionic transform such as (27), said object is represented as a sum of products of Pauli spin operators. We represent these Pauli sums algebraically as firstclass objects in a canonical form. We wish to compute a program which acts identically to the exponentiation of this Pauli representation—as in (23)—on the quantum abstract machine’s . Exponentiating a Pauli sum requires two steps: factorization via Trotterization, and exponentiation of Pauli terms. Consider the Pauli terms
(39) 
These are represented as objects in Python using the pyQuil library:
> A = 2.0*PauliTerm(’X’, 0)*PauliTerm(’X’, 1) > B = 0.5*PauliTerm(’X’, 0)*PauliTerm(’Z’, 2) > print "A =", A, "\nB =", B A = 2.0*X0*X1 B = 0.5*X0*Z2
Exponentiation of a Pauli term has three parts: a change to the basis, an entanglement of the qubits on which the operator nontrivially acts, and a rotation about . For example, we can see this by looking at a Quil program which allows one to compute the action of .
> print exponentiate(A) H 0 H 1 CNOT 0 1 RZ(4.0) 1 CNOT 0 1 H 0 H 1
This Quil program itself is also a firstclass object. Indeed, the result of a Suzuki–Trotter approximation of, say, must combine programs which compute and . We can see this by looking at the Quil program which computes in a firstorder decomposition:
> print trotterize(A, B) H 0 H 1 CNOT 0 1 RZ(4.0) 1 CNOT 0 1 H 0 H 1 H 0 CNOT 0 2 RZ(1.0) 2 CNOT 0 2 H 0
Quil programs are written as a serial list of instructions, but in fact the semantics are not changed if the order between successive commuting instructions is changed. As such, we can think of such instructions occurring in a single time slice^{5}^{5}5If the ^{th} commuting instruction executes in time, then the time slice itself executes in time and is the worst case for arbitrary with maximal parallelization. If all are known for each time slice, then it is possible to do aggressive parallelization by overlapping commuting instructions from subsequent time slices, despite the time slices as a whole (i.e., its action on the total Hilbert space) being noncommuting.. A straightline Quil program itself is parallelized if it is transformed according to a maximal semanticspreserving parallelization. Consider the following 16instruction sample from a state evolution according to the UCCSD ansatz with firstorder Trotter and one time slice for the foursite Hubbard model:
X 2 X 3 H 3 RX(1.5707963267948966) 5 CNOT 3 4 CNOT 4 5 RZ(0.00001) 5 CNOT 4 5 CNOT 3 4 H 3 RX(1.5707963267948966) 5 RX(1.5707963267948966) 1 H 5 CNOT 1 2 CNOT 2 3 CNOT 3 4
Instruction sequences like X 2
and X 3
commute and may execute in the same time slice. However, X 3
and H 3
do not commute and must execute in different time slices. The parallelized Quil program is as follows:
Time Slice #1: X 2 X 3 RX(1.5707963267948966) 5 RX(1.5707963267948966) 1 Time Slice #2: H 3 CNOT 1 2 Time Slice #3: CNOT 3 4 Time Slice #4: CNOT 4 5 Time Slice #5: RZ(0.00001) 5 Time Slice #6: CNOT 4 5 Time Slice #7: CNOT 3 4 RX(1.5707963267948966) 5 Time Slice #8: H 3 H 5 Time Slice #9: CNOT 2 3 Time Slice #10: CNOT 3 4
The 16instruction straightline program has been reduced to a 10timeslice program averaging onequbit gates and twoqubit gates per time slice. All instructions being equal, this gives a 38% improvement in timing.
When a final Quil program is prepared through the various means described, it is dispatched to a representation of a QAM: either a hardware QPU or a software quantum virtual machine. In this work, a remotely deployed quantum virtual machine was used for all Quil execution.
References
 Bauer et al. (2015) B. Bauer, D. Wecker, A. J. Millis, M. B. Hastings, and M. Troyer, arXiv preprint arXiv:1510.03859 (2015).
 Kreula et al. (2016) J. M. Kreula, L. GarcíaÁlvarez, L. Lamata, S. R. Clark, E. Solano, and D. Jaksch, EPJ Quantum Technology 3, 11 (2016).
 Bravyi and Gosset (2016) S. Bravyi and D. Gosset, arXiv preprint arXiv:1609.00735 (2016).
 Reiher et al. (2016) M. Reiher, N. Wiebe, K. M. Svore, D. Wecker, and M. Troyer, arXiv preprint arXiv:1605.03590 (2016).
 Peruzzo et al. (2014) A. Peruzzo, J. McClean, P. Shadbolt, M.H. Yung, X.Q. Zhou, P. J. Love, A. AspuruGuzik, and J. L. OâBrien, Nature communications 5 (2014).
 McClean et al. (2016a) J. R. McClean, J. Romero, R. Babbush, and A. AspuruGuzik, New Journal of Physics 18, 023023 (2016a).
 McClean (2015) J. R. McClean, Algorithms Bridging Quantum Computation and Chemistry, Ph.D. thesis, Harvard University (2015).
 McClean et al. (2016b) J. R. McClean, M. E. Schwartz, J. Carter, and W. A. de Jong, arXiv preprint arXiv:1603.05681 (2016b).
 Wang et al. (2015) Y. Wang, F. Dolde, J. Biamonte, R. Babbush, V. Bergholm, S. Yang, I. Jakobi, P. Neumann, A. AspuruGuzik, J. D. Whitfield, et al., ACS nano 9, 7769 (2015).
 Wecker et al. (2015) D. Wecker, M. B. Hastings, and M. Troyer, Phys. Rev. A 92, 042303 (2015).
 Smith et al. (2016) R. S. Smith, M. J. Curtis, and W. J. Zeng, arXiv preprint arXiv:1608.03355 (2016).
 Sawaya et al. (2016) N. P. D. Sawaya, M. Smelyanskiy, J. R. McClean, and A. AspuruGuzik, Journal of Chemical Theory and Computation 12, 3097 (2016), pMID: 27254482, http://dx.doi.org/10.1021/acs.jctc.6b00220 .
 DallaireDemers and Wilhelm (2016) P.L. DallaireDemers and F. K. Wilhelm, Phys. Rev. A 93, 032303 (2016).
 Helgaker et al. (2014) T. Helgaker, P. Jorgensen, and J. Olsen, Molecular electronicstructure theory (John Wiley & Sons, 2014).
 Werner and Knowles (1985) H. Werner and P. J. Knowles, The Journal of Chemical Physics 82, 5053 (1985).
 Gidofalvi and Mazziotti (2007) G. Gidofalvi and D. A. Mazziotti, The Journal of Chemical Physics 127, 244105 (2007), http://dx.doi.org/10.1063/1.2817602.
 Knizia and Chan (2012) G. Knizia and G. K.L. Chan, Phys. Rev. Lett. 109, 186404 (2012).
 Knizia and Chan (2013) G. Knizia and G. K.L. Chan, Journal of Chemical Theory and Computation 9, 1428 (2013), pMID: 26587604, http://dx.doi.org/10.1021/ct301044e .
 Wouters et al. (2016) S. Wouters, C. A. JimÃ©nezHoyos, Q. Sun, and G. K.L. Chan, Journal of Chemical Theory and Computation 12, 2706 (2016), pMID: 27159268, http://dx.doi.org/10.1021/acs.jctc.6b00316 .
 Bulik et al. (2014a) I. W. Bulik, G. E. Scuseria, and J. Dukelsky, Phys. Rev. B 89, 035140 (2014a).
 Bulik et al. (2014b) I. W. Bulik, W. Chen, and G. E. Scuseria, The Journal of Chemical Physics 141, 054113 (2014b), http://dx.doi.org/10.1063/1.4891861.
 Bulik (2015) I. W. Bulik, Electron correlation in extended systems via quantum embedding, Ph.D. thesis, Rice University (2015).
 Zheng and Chan (2016) B.X. Zheng and G. K.L. Chan, Phys. Rev. B 93, 035126 (2016).
 Booth and Chan (2015) G. H. Booth and G. K.L. Chan, Physical Review B 91, 155107 (2015).
 Peschel (2012) I. Peschel, Brazilian Journal of Physics 42, 267 (2012).
 Nielsen and Chuang (2010) M. A. Nielsen and I. L. Chuang, Quantum computation and quantum information (Cambridge university press, 2010).
 Tsuchimochi et al. (2015) T. Tsuchimochi, M. Welborn, and T. Van Voorhis, The Journal of Chemical Physics 143, 024107 (2015), http://dx.doi.org/10.1063/1.4926650.
 Szabo and Ostlund (1989) A. Szabo and N. S. Ostlund, Modern quantum chemistry: introduction to advanced electronic structure theory (Courier Corporation, 1989).
 Hosteny et al. (1975) R. P. Hosteny, T. H. Dunning, R. R. Gilman, A. Pipano, and I. Shavitt, The Journal of Chemical Physics 62, 4764 (1975).
 Zheng et al. (2016) B.X. Zheng, J. S. Kretchmer, H. Shi, S. Zhang, and G. K. Chan, arXiv preprint arXiv:1608.03316 (2016).
 O’Malley et al. (2016) P. J. J. O’Malley, R. Babbush, I. D. Kivlichan, J. Romero, J. R. McClean, R. Barends, J. Kelly, P. Roushan, A. Tranter, N. Ding, B. Campbell, Y. Chen, Z. Chen, B. Chiaro, A. Dunsworth, A. G. Fowler, E. Jeffrey, E. Lucero, A. Megrant, J. Y. Mutus, M. Neeley, C. Neill, C. Quintana, D. Sank, A. Vainsencher, J. Wenner, T. C. White, P. V. Coveney, P. J. Love, H. Neven, A. AspuruGuzik, and J. M. Martinis, Phys. Rev. X 6, 031007 (2016).
 Whitfield et al. (2011) J. D. Whitfield, J. Biamonte, and A. AspuruGuzik, Molecular Physics 109, 735 (2011).
 Tranter et al. (2015) A. Tranter, S. Sofia, J. Seeley, M. Kaicher, J. McClean, R. Babbush, P. V. Coveney, F. Mintert, F. Wilhelm, and P. J. Love, International Journal of Quantum Chemistry 115, 1431 (2015).
 Seeley et al. (2012) J. T. Seeley, M. J. Richard, and P. J. Love, The Journal of chemical physics 137, 224109 (2012).
 Crawford and Schaefer (2000) T. D. Crawford and H. F. Schaefer, Reviews in computational chemistry 14, 33 (2000).
 Suzuki (1994) M. Suzuki, Communications in mathematical physics 163, 491 (1994).
 Trotter (1959) H. F. Trotter, Proceedings of the American Mathematical Society 10, 545 (1959).
 Wouters (2015) S. Wouters, “QCDMET,” https://github.com/SebWouters/QCDMET (2015).
 Su (2016) Q. Su, “PySCF,” https://github.com/sunqm/pyscf (2016).
 Nocedal and Wright (2006) J. Nocedal and S. Wright, Numerical optimization (Springer Science & Business Media, 2006).
 Bulik et al. (2015) I. W. Bulik, T. M. Henderson, and G. E. Scuseria, Journal of Chemical Theory and Computation 11, 3171 (2015), pMID: 26575754.
 Scuseria et al. (1988) G. E. Scuseria, C. L. Janssen, and H. F. Schaefer, The Journal of Chemical Physics 89, 7382 (1988).
 Shiba (1972) H. Shiba, Phys. Rev. B 6, 930 (1972).
 Lieb and Wu (2003) E. H. Lieb and F. Wu, Physica A: Statistical Mechanics and its Applications 321, 1 (2003).