Many-body techniques based on the double unitary coupled cluster (DUCC) ansatz can be used to downfold electronic Hamiltonians into low-dimensional active spaces. It can be shown that the resulting dimensionality reduced Hamiltonians are amenable for quantum computing. Recent studies performed for several benchmark systems using phase estimation (PE) algorithms for quantum computers demonstrated that these formulations can recover a significant portion of ground-state dynamical correlation effects that stem from the electron excitations outside of the active space. These results have also been confirmed in studies of ground-state potential energy surfaces using quantum simulators. In this letter, we study the effectiveness of the DUCC formalism in describing excited states. We also emphasize the role of the PE formalism and its stochastic nature in discovering/identifying excited states or excited-state processes in situations when the knowledge about the true configurational structure of a sought after excited state is limited or postulated (due to the specific physics driving excited-state processes of interest). In this context, we can view PE algorithms as an engine for verifying various hypotheses for excited-state processes and providing statistically meaningful results that correspond to the electronic state(s) with the largest overlap with a postulated configurational structure. We illustrate these ideas on examples of strongly correlated molecular systems, characterized by small energy gaps and high density of quasidegenerate states around the Fermi level.

## I. INTRODUCTION

There is significant interest in applying quantum computing techniques to describe and simulate chemical systems and processes.^{1–17} This approach brings hope to address the exponential barriers limiting the applicability of exact diagonalization procedures [or full configuration interaction (CI) methods] and also to provide access to complicated multiconfigurational electronic states, which often cannot be identified or described by vast classes of approximate methods used in routine simulations. Of special interest is the application to strongly correlated molecular systems characterized by small energy gaps between occupied and unoccupied orbitals where multiple electronic states that are defined by complex collective phenomena exist, usually involving higher than single excitations in the corresponding wavefunction expansions.

Even though impressive progress has been achieved in the development of wavefunction-based excited-state approaches such as complete-active-space perturbation theory (CASPT),^{18,19} multireference NEVPT,^{20,21} configuration interaction,^{22} equation-of-motion coupled cluster (EOMCC),^{23–29} multireference coupled cluster (MRCC),^{30–47} and the density matrix renormalization group (DMRG),^{48–51} problems with the description of complicated states dominated by high-rank excitations still exist. For example, in order to capture these states with EOMCC formalisms, one needs to include higher-than-double excitations.^{52} For doubly excited states, the “minimum” level of theory to tackle these states is EOMCC with singles, doubles, and triples (EOMCCSDT),^{53,54} although in several cases, it may not provide a quantitative level of accuracy.^{55} These problems commonly occur even for small molecular systems, and it is reasonable to expect that they may intensify for larger systems and strongly correlated systems (transition metal oxides, metal clusters, and actinides) with a high density of states located around the Fermi level. Recently, significant progress in addressing these problems has been achieved by integrating stochastic configuration interaction (CI) Quantum Monte Carlo (QMC) framework^{56,57} with deterministic EOMCC formulations.^{58} Another interesting aspect of modeling these complex excited states is the attainability of these states in situations where their initial configurational structure cannot be easily obtained to commence convergent iterative procedures.

Progress in the development of quantum algorithms may provide alternative solutions to these problems. The recent application of Variational Quantum Eigensolvers (VQEs)^{5,8,11,12,59–63} and Phase Estimation (PE)^{6,64–71} to excited states^{61,72–75} demonstrates that excited states can be effectively simulated on the quantum computers while at the same time bypassing the problems of conventional computing and approximate formulations.

In this paper, we present an excited state extension of recently developed techniques for active-space downfolding of the electronic Hamiltonian based on the double unitary coupled cluster (DUCC) transformation. This is combined with PE simulations available in the Microsoft Quantum Development Kit (QDK)^{74,76} to illustrate the excited-state version of DUCC formalism on the examples of H_{2} at equilibrium and stretched bond lengths and two H_{4} models: (1) trapezoidal H_{4}, which is a popular benchmark system for studying quasidegenerate states,^{77} and (2) a linear H_{4} molecule, which was intensively studied in the context of singlet fission processes.^{78} For each of these systems, we perform PE simulations to characterize the structure of excited states using DUCC effective Hamiltonians. Moreover, we carry out a number of PE simulations to investigate the role and effect of different initial guesses, which are based on limited knowledge about the excited state of interest.

## II. THEORY OF DUCC DOWNFOLDED HAMILTONIANS

In Ref. 79, we introduced the unitary extension of the subsystem embedding subalgebra CC approach (SES-CC),^{80} which utilizes the double unitary CC expansion,

The characteristics of the expansion (1) is similar to the expansion discussed in the single-reference formulation of the active-space coupled cluster formalism^{81,82} (see also Refs. 83 and 84), which also utilizes the decomposition of the cluster operator into internal and external parts.

In analogy to Ref. 79, *σ*_{int} and *σ*_{ext} are the anti-Hermitian operators ($\sigma int\u2020=\u2212\sigma int$ and $\sigma ext\u2020=\u2212\sigma ext$) defined by excitations/de-excitations within and outside of active space, respectively. To be more precise, the amplitudes defining the *σ*_{ext} operator must carry at least one inactive spin-orbital index. Using (1) in Schrödinger’s equation, one obtains equations for cluster amplitudes and the corresponding energy,

where *Q* is a projection operator on the space spanned by determinants orthogonal to the reference function $\Phi $. In Eqs. (1)–(17) of this section, we consider the case of the exact limit (*σ*_{int} and *σ*_{ext} include all possible excitations). In Ref. 79, we showed that when *σ*_{int} contains all possible excitations/de-excitations within the complete active space, the energy of system (2) can be obtained by diagonalizing the DUCC effective Hamiltonian,

where

and

In the above eigenvalue problem, the $e\sigma int\Phi $ expansion defines a corresponding eigenvector and *P* and *Q*_{int} are projection operators onto the reference function |Φ⟩ and excited determinants in the active space orthogonal to $\Phi $, respectively.

To show this property, it is sufficient to introduce the resolution of identity $e\sigma inte\u2212\sigma int$ to the left of the $H\xafextDUCC$ operator in

where we explicitly used the property of the $e\sigma int|\Phi $ expansion

and to note that $e\u2212\sigma intH\xafextDUCCe\sigma int=e\u2212\sigma inte\u2212\sigma extHe\sigma exte\sigma int$. Next, using matrix representation of the *σ*_{int} operator in the CAS space denoted as *σ*_{int}, this equation can be rewritten as

where the first component of the [** y**] vector is equivalent to $\u27e8\Phi |e\u2212\sigma inte\u2212\sigma extHe\sigma exte\sigma int|\Phi \u27e9\u2212E$, while the remaining components correspond to projections of $e\u2212\sigma inte\u2212\sigma extHe\sigma exte\sigma int|\Phi $ onto excited configurations belonging to

*Q*

_{int}. The $[e\sigma int]$ matrix is also nonsingular, which is a consequence of the formula

and the anti-Hermitian characteristics of the *σ*_{int} matrix, i.e., Tr(*σ*_{int}) = 0 (where the real characteristics of *σ*_{int} cluster amplitudes is assumed). Given the nonsingular characteristics of the $[e\sigma int]$ matrix (see also Ref. 79), this proves the equivalence of these two representations.

The proof of the above property is not limited to the ground state and can be extended to any electronic state described by (1)–(3). Assuming that this ansatz can describe excited states, the DUCC effective Hamiltonian formalism can be used in the context of excited-state simulations. We will denote the general DUCC solution corresponding to the *K*th state as

where the *K*th state energy can be obtained from diagonalizing the state-specific effective Hamiltonian,

where

and

Similar to the ground-state effective/downfolded Hamiltonians, the operators $H\xafexteff(DUCC)(K)$ are Hermitian and therefore amenable to the real-time simulation on a quantum computer. In analogy to the ground-state representation of DUCC, we will assume that

and

where *S*_{int}(*K*) and *S*_{ext}(*K*) are CC-type cluster operators producing excitations within and outside the active space, respectively, when acting on the reference function $\Phi $.

If the exact form of the operator *σ*_{ext}(*K*) [or *S*_{ext}(*K*)] is known, the effective Hamiltonian $H\xafexteff(DUCC)(K)$ can be diagonalized to find corresponding exact energy *E*_{K}. In practice, in likeness to the ground-state formulation, we obtain an approximate model *σ*_{ext}(*K*) operator by way of calculations with excited-state CC models. Additionally, we previously employed in ground-state DUCC^{79} the many-body form of the effective Hamiltonian $H\xafexteff(DUCC)$ driven by the perturbative analysis of the ground-state energy expansion. However, the same arguments cannot be invoked in the context of the excited-state variant of $H\xafextDUCC(K)$. Instead, in the analysis of the excited states, we use a basic expansion for $H\xafextDUCC(K)$ given by an expression involving a single commutator,

In the future studies, we will explore the impact of the higher-order commutators on the quality of DUCC energies.

In this paper, we explore a simple strategy based on the utilization of the excited-state wavefunction in the equation-of-motion CC parameterization $\Psi KEOMCC$,

as a reference to extract the relevant information about *σ*_{ext}(*K*). In Eq. (19), the cluster operator *T* satisfies the CC equations and the excitation operator *R*_{K} (corresponding to the *K*th excited state) is obtained through diagonalization of the similarity transformed Hamiltonian $H\xaf=e\u2212THeT$. Since (11) represents a normalized state, in order to compare with the corresponding EOMCC expansion in the exact wavefunction limit, one needs to use a normalized form $\Psi \u0303KEOMCC$ of (19),

where

In the simplest approximate variants, we use the EOMCCSD approximations (EOMCC with singles and doubles^{24}) to extract singly excited [*S*_{ext, 1}(*K*)] and doubly excited [*S*_{ext, 2}(*K*)] components of the *S*_{ext}(*K*) operator,

where *Q*_{ext, 1} and *Q*_{ext, 2} are projection operators on subspaces of singly and doubly excited external excitations, respectively. In (22) and (23), we approximate $\Psi \u0303KEOMCC\u2248\Psi KEOMCCSD(A)$, which is defined in the following way:

where *Q*_{1} and *Q*_{2} operators are projection operators on spaces of singly and doubly excited configurations, *R*_{K,i} (*i* = 0, 1, 2) represent EOMCC excitation operators producing *i*-tuply excited configurations when acting onto reference functions, and *T*_{1} and *T*_{2} are the singly and doubly excited cluster operators, respectively. In our approximation, referred to as the DUCC-ex(*K*) approach, we use *T*_{1} and *T*_{2} from standard CCSD calculations and *R*_{K,i} (*i* = 0, 1, 2) from the standard EOMCCSD diagonalization procedure. Consequently, *S*_{ext, 1}(*K*) and *S*_{ext, 2}(*K*) take the following form:

where

As in the ground-state case, the DUCC-ex(*K*) approximation is defined by the length of the commutator expansion and the source of the *σ*_{ext}(*K*) amplitudes. Since the EOMCCSD approximation is used to define *σ*_{ext}(*K*), one should expect that this scheme works in cases where the EOMCCSD approach delivers reliable results and active-space used in the state-selective DUCC-ex(*K*) formalism capable of capturing main configurations corresponding to the *K*th state.

## III. ALGORITHM

Our simulation of the DUCC Hamiltonian consists of two broad components: (1) prepare a suitable classical representation of the DUCC Hamiltonian and the CCSD and EOMCCSD state to be simulated and (2) perform a quantum simulation of the system using quantum algorithms. The proposed workflow is analogous to the one introduced by Wang *et al.*,^{72} up to a difference in initial state(s), the form of the diagonalized Hamiltonian, and the numerical simulation of the quantum algorithm.

The CCSD and EOMCCSD calculations are performed for the ground and excited states using Tensor Contraction Engine (TCE)^{85} CCSD and EOMCCSD implementations available in NWChem.^{86} The quantum simulations are performed using the QDK-NWChem interface and targets a classical emulation of a quantum computer. Here, NWChem generates a yaml file containing the orbital approximation of the Hamiltonian and quantum state. Subsequently, eigenenergy estimates are obtained using QDK implementations of QPE and simulation based on either Trotterization^{66} or qubitization.^{76}

We now describe these components and estimate the quantum resources required. Note that the number of quantum gates executed is polynomial in the number of orbitals. This is efficient on a quantum computer but of course requires exponential time on a classical simulation.

### A. Representation of Hamiltonian and state

The following steps 1–4 are performed using classical resources.

**Form**$|\Psi \u0303KEOMCC$**.**This state approximates the targeted excited eigenstate, call it $\Psi K$, and is obtained by normalizing the $|\Psi KEOMCCSD(A)$ vector and forming*S*_{ext,1}(*K*) and*S*_{ext,2}(*K*) amplitudes as defined by (25) and (26).**Form the**$H\xafextDUCC(K)$**Hamiltonian in the particle-hole representation.**As in Ref. 79, we calculate the normal product form of $H\xafextDUCC(K)$ (18), which is approximated by the scalar, one-, and two-body interactions, i.e.,

where *N*[⋯]^{ph} denotes the particle-hole (ph) normal product form of a given operator expression and *P*, *Q*, *R*, *S*, … denote general active spin-orbital indices. The algebraic form of the $\gamma (K)QP$ and $\gamma (K)RSPQ$ tensors for the situation when all occupied orbitals are deemed active is given in Tables I and II.

**Transform**$H\xafextDUCC(K)$**from particle-hole representation****(28)****to Fermi vacuum representation.**

using identities from Table IV of Ref. 79. In analogy to Ref. 79 as an input for the QDK, we are using a symmetric orbital approximation of the spin orbital *χ* tensor, where it is assumed that Mulliken orbital-type dressed integrals (**PQ**|**RS**) (as well as (**PQ**|**SR**), (**QP**|**RS**), …, (**SR**|**QP**)) are obtained from $\chi R\alpha S\beta P\alpha Q\beta $ (where **P**, **Q**, **R**, and **S** are the orbitals corresponding to active spin orbitals *P*, *Q*, *R*, and *S*, respectively). In general, the number *M* of unique integrals (accounting for the symmetries (**PQ**|**RS**) = (**PQ**|**SR**), (**QP**|**RS**), …, (**SR**|**QP**)) is much fewer than the worst-case of $M=O(N4)$, where *N* is the number of active-space spin-orbitals.

**Transform Fermi vacuum representation to Jordan-Wigner representation.**The targeted quantum computer applies quantum gates on qubits rather than fermions. Thus, we map fermion operators to Pauli operators {*X*_{P},*Y*_{P},*Z*_{P}} acting on qubits*P*∈ [0,*N*− 1]. It suffices for our purposes to use the Jordan-Wigner representation,

One may verify that the transformed operators satisfy the expected fermion anticommutation relationships. Thus, the DUCC Hamiltonian is of the form

where *P*_{j} are Pauli operators with coefficients *α*_{j}.

Amplitude . | Expression . | Amplitude . | Expression . |
---|---|---|---|

H_{N} term | |||

$\gamma AB=fAB$ | $\gamma IJ=fIJ$ | ||

$\gamma AI=fAI$ | $\gamma IA=fIA$ | ||

$\gamma IABC=vIABC$ | $\gamma IJKA=vIJKA$ | ||

$\gamma ABCI=vABCI$ | $\gamma KAIJ=vKAIJ$ | ||

$\gamma IAJB=vIAJB$ | $\gamma ABCD=vABCD$ | ||

$\gamma IJKL=vIJKL$ | $\gamma ABIJ=vABIJ$ | ||

$\gamma IJAB=vIJAB$ |

Amplitude . | Expression . | Amplitude . | Expression . |
---|---|---|---|

H_{N} term | |||

$\gamma AB=fAB$ | $\gamma IJ=fIJ$ | ||

$\gamma AI=fAI$ | $\gamma IA=fIA$ | ||

$\gamma IABC=vIABC$ | $\gamma IJKA=vIJKA$ | ||

$\gamma ABCI=vABCI$ | $\gamma KAIJ=vKAIJ$ | ||

$\gamma IAJB=vIAJB$ | $\gamma ABCD=vABCD$ | ||

$\gamma IJKL=vIJKL$ | $\gamma ABIJ=vABIJ$ | ||

$\gamma IJAB=vIJAB$ |

Amplitude . | Expression . |
---|---|

$(HNText)C,open$ term | |

$\gamma AB+=veAMBsMe\u221212veAMNsMNeB$ | |

$\gamma IJ+=veIMJsMe+12vefMJsMIef$ | |

$\gamma AI+=veAMIsMe$ | |

$\gamma IA+=feAsIe+veIMAsMe\u221212veIMNsMNeA+12vefMAsMIef$ | |

$\gamma IABC+=veABCsIe\u2212veAMBsMIeC+veAMCsMIeB$ | |

$\gamma IJKA+=veJKAsIe\u2212veIKAsJe+12vefKAsIJef\u2212veJMKsMIeA+veIMKsMJeA$ | |

$\gamma KAIJ+=veAIJsKe$ | |

$\gamma IAJB+=veAJBsIe\u2212veAMJsMIeB$ | |

$\gamma IJKL+=veJKLsIe\u2212veIKLsJe+12vefKLsIJef$ | |

$\gamma IJAB+=veJABsIe\u2212veIABsJe+feAsIJeB\u2212feBsIJeA+12vefABsIJef$ | |

$\u2212veJMAsMIeB+veIMAsMJeB+veJMBsMIeA\u2212veIMBsMJeA$ | |

$(Text\u2020HN)C,open$ term | |

$\gamma BA+=vMBeAseM\u221212vMNeAseBMN$ | |

$\gamma JI+=vMJeIsMe+12vMJefsefMI$ | |

$\gamma IA+=vMIeAsMe$ | |

$\gamma AI+=fAeseI+vMAeIseM\u221212vMNeIseAMN+12vMAefsefMI$ | |

$\gamma BCIA+=vBCeAseI\u2212vMBeAseCMI+vMCeAseBMI$ | |

$\gamma KAIJ+=vKAeJseI\u2212vKAeIseJ+12vKAefsefIJ\u2212vMKeJseAMI+vMKeIseAMJ$ | |

$\gamma IJKA+=vIJeAseK$ | |

$\gamma JBIA+=vJBeAseI\u2212vMJeAseBMI$ | |

$\gamma KLIJ+=vKLeJseI\u2212vKLeIseJ+12vKLefsefIJ$ | |

$\gamma ABIJ+=vABeJseI\u2212vABeIseJ+fAeseBIJ\u2212fBeseAIJ+12vABefsefIJ$ | |

$\u2212vMAeJseBMI+vMAeIseBMJ+vMBeJseAMI\u2212vMBeIseAMJ$ |

Amplitude . | Expression . |
---|---|

$(HNText)C,open$ term | |

$\gamma AB+=veAMBsMe\u221212veAMNsMNeB$ | |

$\gamma IJ+=veIMJsMe+12vefMJsMIef$ | |

$\gamma AI+=veAMIsMe$ | |

$\gamma IA+=feAsIe+veIMAsMe\u221212veIMNsMNeA+12vefMAsMIef$ | |

$\gamma IABC+=veABCsIe\u2212veAMBsMIeC+veAMCsMIeB$ | |

$\gamma IJKA+=veJKAsIe\u2212veIKAsJe+12vefKAsIJef\u2212veJMKsMIeA+veIMKsMJeA$ | |

$\gamma KAIJ+=veAIJsKe$ | |

$\gamma IAJB+=veAJBsIe\u2212veAMJsMIeB$ | |

$\gamma IJKL+=veJKLsIe\u2212veIKLsJe+12vefKLsIJef$ | |

$\gamma IJAB+=veJABsIe\u2212veIABsJe+feAsIJeB\u2212feBsIJeA+12vefABsIJef$ | |

$\u2212veJMAsMIeB+veIMAsMJeB+veJMBsMIeA\u2212veIMBsMJeA$ | |

$(Text\u2020HN)C,open$ term | |

$\gamma BA+=vMBeAseM\u221212vMNeAseBMN$ | |

$\gamma JI+=vMJeIsMe+12vMJefsefMI$ | |

$\gamma IA+=vMIeAsMe$ | |

$\gamma AI+=fAeseI+vMAeIseM\u221212vMNeIseAMN+12vMAefsefMI$ | |

$\gamma BCIA+=vBCeAseI\u2212vMBeAseCMI+vMCeAseBMI$ | |

$\gamma KAIJ+=vKAeJseI\u2212vKAeIseJ+12vKAefsefIJ\u2212vMKeJseAMI+vMKeIseAMJ$ | |

$\gamma IJKA+=vIJeAseK$ | |

$\gamma JBIA+=vJBeAseI\u2212vMJeAseBMI$ | |

$\gamma KLIJ+=vKLeJseI\u2212vKLeIseJ+12vKLefsefIJ$ | |

$\gamma ABIJ+=vABeJseI\u2212vABeIseJ+fAeseBIJ\u2212fBeseAIJ+12vABefsefIJ$ | |

$\u2212vMAeJseBMI+vMAeIseBMJ+vMBeJseAMI\u2212vMBeIseAMJ$ |

### B. Simulation with quantum algorithms

We obtain estimates of excited state energies using a phase estimation algorithm on a quantum circuit approximating the unitary time-evolution operator *e*^{−iHt}, acting on a trial wavefunction $\Psi Trial$. As abstractly illustrated in Fig. 1, the trial wavefunction may have overlap with several states ($\Psi 1$ and $\Psi 2$) in Fig. 1 and be orthogonal to others ($\Psi 3$) in Fig. 1. Each run of PE returns an approximate energy estimate for a random eigenstate with probability equal to the amount of overlap |⟨Ψ_{j}|Ψ_{Trial}⟩|^{2} of the trial wavefunction with that state’s corresponding wavefunction. As an example, Fig. 2 illustrates just such a distribution of energies for a strongly correlated system. Importantly, one will always obtain an energy estimate even if this overlap is small. This is sufficient for the purpose of characterizing the excited state spectrum.

We now outline the quantum gate cost of these algorithms.

#### 1. Trial wavefunction preparation

As described in Ref. 74, the distribution of energies for the ground and excited states from PE is determined by the Hamiltonian and a trial wavefunction composed of a superposition of *m* Slater determinants. In the QDK, this state can be prepared nondeterministically with success probability $O(1/N2)$ using $O(mN)$ quantum gates. With quantum amplitude amplification,^{87} the success probability can be boosted close to unity using a total of $O(mN2)$ quantum gates. Note that alternate techniques may efficiently synthesize on a quantum computer trial wavefunction with exponentially many components, such as using the UCC formalism.^{79} However, this is unnecessary for PE methods when the exact wavefunction is dominated by only a few Slater determinants. Unlike VQE techniques, where obtaining energy estimates that are correct to chemical accuracy requires the overlap to also be unity up to at least chemical accuracy, PE methods only need the overlap to be polynomially large, say, |⟨Ψ_{j}|Ψ_{Trial}⟩|^{2} = Ω(1/poly(*N*)).

For this purpose, a wide array of excited-state EOMCC formulations implemented in NWChem can be used to provide necessary information about the configurational structure of excited states in molecular systems and therefore complement a strategy based on the utilization multiconfigurational self-consistent field (MCSCF) wave functions (see Ref. 72). The following part of the discussion in this section will however deal with the hypothetical situation when our knowledge of the initial state, for whatever reason, is not perfect and we cannot precisely determine its structure, instead, we can only hypothesize about the structure of the sought state. This is a situation which may occur in studies of strongly correlated systems. Through the repeated application of PE where a carefully chosen trial wavefunction is reinitialized each time, one can obtain a distribution of energies for several states of interest.

The stochastic nature of PE algorithms in providing information on a distribution of targeted states is unlike VQE algorithms that only provide energy estimates for a single targeted state. As which single state to target may not be obvious *a priori*, PE methods enable the discovery and investigation of exotic and novel states that are unobtainable with conventional computing and current approximate formulations.

#### 2. Phase estimation

Phase estimation algorithms, such as quantum phase estimation,^{87} return a random eigenphase *θ*_{j} of an arbitrary unitary $U=\u2211jei\theta j\Phi j\Phi j$ with eigenvectors $\Phi j$. The dominant cost of this algorithm is in applying controlled-*U* some number *T* times. The dominant error in the PE energy estimate arises from an analog of the time-energy uncertainty relation. Each application of phase estimation estimates *θ*_{j} with a standard deviation of *σ*_{θ} = Θ(*T*^{−1}), and the robust PE algorithm^{68} we use in the QDK empirically achieves *σ*_{θ} ≈ 2/*T* and uses only one ancillary qubit.

Implementing *U* exactly can be difficult as it must be expressed in terms of primitive quantum gates that may be implemented on a universal quantum computer. Thus, we permit a maximum simulation error *ϵ*_{θ} and instead output an operation *U*′ such that

This implies that the eigenvalues and eigenphases $\theta j\u2032$ of *U*′ differ from *U* by at most *ϵ*_{θ}. The overall error in estimating *θ*_{j} is then a random error of *σ*_{θ} and a systematic shift of *ϵ*.

#### 3. Unitary time-evolution

In order to estimate energy with a standard deviation *σ*_{E}, PE algorithms query the unitary time-evolution operator $e\u2212iH\Delta =\u2211je\u2212iEj\Delta \psi j\psi j$ by the DUCC Hamiltonian. In order to avoid ambiguity in the principle value of the phase, Δ should be chosen so that for the eigenstate of energy, the range of energies |*E*_{max} − *E*_{min}|Δ ≤ 2*π*. Thus, obtaining an energy estimate with standard deviation *σ*_{E} requires *σ*_{E}Δ ≈ 2/*T*. In many cases, the step-size Δ determines the error of approximation. Suppose we approximate *e*^{−iHΔ} by a unitary operator with error,

Then, this contributes a systematic error in energy of *ϵ*_{E} = *ϵ*_{θ}/Δ. At present, the two most popular methods for constructing the operator *U* are Trotter-Suzuki methods^{66} and qubitization,^{76} both of which are available in the QDK.

In the first-order Trotter-Suzuki method, real-time evolution by the DUCC Hamiltonian is approximated by

where $\Vert \alpha \u2192\Vert 1=\u2211j|\alpha j|$. The second-order method with error $O|\Delta |3\Vert \alpha \u2192\Vert 13$ applies a symmetric sequence *U*_{2}(Δ) that is *U*_{1}(Δ/2) followed by its component exponentials in reverse order. In order for *ϵ*_{E} to be commensurate with *σ*_{E}, this implies that

Higher-order *j* methods with error $O5j/2\u22121|\Delta |j+1\Vert \alpha \u2192\Vert 1j+1$ are known but can be impractical due to the exponential prefactor. It is also worth noting that these error bounds represent worst-case behavior. In empirical studies,^{13} the step-size can be much larger than suggested by Eq. (35).

The overall gate cost of obtaining a single energy estimate is then *T* times the gate cost of *U*_{2}(Δ). The overall space requirement is *N* qubits for representing a system with *N* spin-orbitals in the Jordan-Wigner representation and one more for phase estimation. Within fault-tolerant models of quantum computation, the Clifford gates that conjugate to *P*_{j} to a single-qubit Pauli *Z* are considered cheap and dominant cost of *U*_{1}(Δ) is in applying Θ(*M*) arbitrary single-qubit *Z*-rotations. Alternate methods based on quantum signal processing^{88} or qubitization^{76,89,90} achieve even smaller error. We do not consider them further as they require $O(log\u2061N)$ additional ancillary qubits, which is not conducive to a numerical simulation.

To summarize, this is the scheme of phase estimation combined with unitary time-evolution and is efficient as the number of quantum gates is polynomial in the number of spin-orbitals. Furthermore, speedups in gate depth and count could be obtained by using more compressed Hamiltonians. However, progress in this field is contingent on the progress achieved in understanding the many-body structure of downfolded Hamiltonians and quantifying its errors, such as in the representation of the *σ*(*K*)_{ext} operator.

## IV. NUMERICAL TESTS

We performed a series of numerical tests with the QDK simulator^{74} for three systems characterized by strong ground-state correlation effects and low-lying singlet excited states defined by complicated multireference configurational structures: (1) H_{2} systems for H–H separation corresponding to the equilibrium (*R*_{H–H} = 1.4008 a.u.) and a stretched geometry (*R*_{H–H} = 10 a.u.), (2) H_{4} system in the trapezoidal configuration corresponding to geometrical parameter *α* equal to 0.001 [see Fig. 3(a)]—for this geometry, one can observe a strong quasidegeneracy of low-lying electronic states, and (3) linear form of H_{4} used in studies of singlet fission processes [see Fig. 3(b)]. Special attention is paid to singlet doubly excited states that pose a significant challenge to existing many-body methodologies. For systems considered here we employ the cc-pVTZ basis set^{91} and we correlate all electrons in all calculations. Since the QDK cannot exploit spatial symmetry, all PE simulations were performed using *C*_{1} symmetry. For simplicity, we will refer to the *D*_{2h}, *C*_{2v}, and *D*_{2h} irreducible representations of the H_{2}, H_{4} trapezoidal, and H_{4} linear systems, respectively, when reporting full configuration interaction (FCI) numbers obtained with NWChem^{86} (or equivalently, EOMCCSD or EOMCCSDTQ results for the H_{2} and H_{4} models, respectively).

For H_{2} models, we examine the performance of the bare and DUCC-ex transformed Hamiltonians for excited states, which have a significant (A) or partial (PA) amount of leading characteristic excitations in the active space, which contains one occupied and three lowest-lying RHF orbitals. Energies of these states for *R*_{H–H} = 1.4008 a.u. and *R*_{H–H} = 10.0 a.u. are shown in Tables III and IV, respectively. Since the lowest-lying doubly excited states are of *A*_{g} symmetry, in the forthcoming analysis, we will entirely focus on the states corresponding to this symmetry. In Table III, we consider three states 2^{1}*A*_{g}, 3^{1}*A*_{g}, and 5^{1}*A*_{g}, which are either comprehensively described within the active space (2^{1}*A*_{g}) or have a significant overlap with it (3^{1}*A*_{g} and 5^{1}*A*_{g} states). The 2^{1}*A*_{g} is dominated by $\Phi 13$ and $\Phi 1\xaf3\xaf$ single excitations, while the two other states have non-negligible out-of-active-space components. For example, the 3^{1}*A*_{g} state is dominated by $\Phi 11\xaf22\xaf$ excitation but contains important component originating in the $\Phi 17$ and $\Phi 1\xaf7\xaf$ excitations. In a similar way, the 5^{1}*A*_{g} state is dominated by $\Phi 11\xaf33\xaf$ excitations, but it also contains important contributions from $\Phi 11\xaf37\xaf$ and $\Phi 11\xaf73\xaf$ configurations, which do not belong to the active space.

State . | . | DUCC-ex(2^{1}A_{g})
. | DUCC-ex(3^{1}A_{g})
. | DUCC-ex(5^{1}A_{g})
. | Bare Hamiltonian . |
---|---|---|---|---|---|

(characteristics)^{a}
. | FCI^{b}
. | SD: $12(\Phi 13+\Phi 1\xaf3\xaf)$ . | SD: $\Phi 11\xaf22\xaf$ . | SD: $\Phi 11\xaf33\xaf$ . | SD: $\Phi 11\xaf22\xaf$ . |

2^{1}A_{g} | −0.5487 | −0.5478 ± 0.0066 | −0.5197 ± 0.0019 | −0.5034 ± 0.0030 | −0.5306 ± 0.0017 |

(A) | (0.0009) | (0.0290) | (0.0453) | (0.0181) | |

3^{1}A_{g} | −0.1210 | … | 0.0644 ± 0.0017 | −0.0696 ± 0.0016 | −0.0622 ± 0.0018 |

(PA) | (0.0566) | (0.0514) | (0.0588) | ||

5^{1}A_{g} | 0.2310 | 0.2914 ± 0.0021 | 0.2819 ± 0.0065 | 0.2671 ± 0.0031 | 0.2842 ± 0.0014 |

(PA) | (0.0604) | (0.0508) | (0.0361) | (0.0532) |

State . | . | DUCC-ex(2^{1}A_{g})
. | DUCC-ex(3^{1}A_{g})
. | DUCC-ex(5^{1}A_{g})
. | Bare Hamiltonian . |
---|---|---|---|---|---|

(characteristics)^{a}
. | FCI^{b}
. | SD: $12(\Phi 13+\Phi 1\xaf3\xaf)$ . | SD: $\Phi 11\xaf22\xaf$ . | SD: $\Phi 11\xaf33\xaf$ . | SD: $\Phi 11\xaf22\xaf$ . |

2^{1}A_{g} | −0.5487 | −0.5478 ± 0.0066 | −0.5197 ± 0.0019 | −0.5034 ± 0.0030 | −0.5306 ± 0.0017 |

(A) | (0.0009) | (0.0290) | (0.0453) | (0.0181) | |

3^{1}A_{g} | −0.1210 | … | 0.0644 ± 0.0017 | −0.0696 ± 0.0016 | −0.0622 ± 0.0018 |

(PA) | (0.0566) | (0.0514) | (0.0588) | ||

5^{1}A_{g} | 0.2310 | 0.2914 ± 0.0021 | 0.2819 ± 0.0065 | 0.2671 ± 0.0031 | 0.2842 ± 0.0014 |

(PA) | (0.0604) | (0.0508) | (0.0361) | (0.0532) |

^{a}

Characteristics of the electronic states are determined by the active-space contribution: A, dominated by active-space configurations; PA, dominated by configurations not belonging to active space.

^{b}

Full configuration interaction calculations were performed using all 30 molecular orbitals.

State . | . | DUCC-ex(2^{1}A_{g})
. | DUCC-ex(3^{1}A_{g})
. | DUCC-ex(5^{1}A_{g})
. | DUCC-ex(7^{1}A_{g})
. | Bare Hamiltonian . |
---|---|---|---|---|---|---|

(characteristics)^{a}
. | FCI^{b}
. | SD: $\Phi 11\xaf22\xaf$ . | SD: $12(\Phi 13+|\Phi 1\xaf3\xaf)$ . | SD: $\Phi 11\xaf33\xaf$ . | SD: $\Phi 11\xaf44\xaf$ . | SD: $\Phi 11\xaf22\xaf$ . |

2^{1}A_{g} | −0.5981 | −0.6047 ± 0.0018 | −0.5847 ± 0.0018 | −0.5517 ± 0.0018 | −0.5610 ± 0.0013 | −0.5803 ± 0.0019 |

(A) | (−0.0065) | (0.0134) | (0.0464) | (0.0371) | (0.0178) | |

3^{1}A_{g} | −0.4873 | −0.4825 ± 0.0045 | −0.4848 ± 0.0050 | −0.4747 ± 0.0043 | −0.4785 ± 0.0035 | −0.4795 ± 0.0043 |

(A) | (0.0048) | (0.0025) | (0.0126) | (0.0088) | (0.0078) | |

5^{1}A_{g} | −0.1040 | −0.0812 ± 0.0018 | −0.0779 ± 0.0017 | −0.0625 ± 0.0018 | −0.0592 ± 0.0020 | −0.0763 ± 0.0017 |

(A) | (0.0228) | (0.0261) | (0.0415) | (0.0448) | (0.0277) | |

6^{1}A_{g} | 0.0238 | … | … | 0.0405 ± 0.0016 | 0.0365 ± 0.0017 | 0.0368 ± 0.0024 |

(A) | (0.0167) | (0.0127) | (0.0130) | |||

7^{1}A_{g} | 0.2326 | 0.2525 ± 0.0023 | 0.2545 ± 0.0019 | 0.2623 ± 0.0019 | 0.2584 ± 0.0017 | 0.2556 ± 0.0020 |

(A) | (0.0199) | (0.0219) | (0.0297) | (0.0258) | (0.0230) |

State . | . | DUCC-ex(2^{1}A_{g})
. | DUCC-ex(3^{1}A_{g})
. | DUCC-ex(5^{1}A_{g})
. | DUCC-ex(7^{1}A_{g})
. | Bare Hamiltonian . |
---|---|---|---|---|---|---|

(characteristics)^{a}
. | FCI^{b}
. | SD: $\Phi 11\xaf22\xaf$ . | SD: $12(\Phi 13+|\Phi 1\xaf3\xaf)$ . | SD: $\Phi 11\xaf33\xaf$ . | SD: $\Phi 11\xaf44\xaf$ . | SD: $\Phi 11\xaf22\xaf$ . |

2^{1}A_{g} | −0.5981 | −0.6047 ± 0.0018 | −0.5847 ± 0.0018 | −0.5517 ± 0.0018 | −0.5610 ± 0.0013 | −0.5803 ± 0.0019 |

(A) | (−0.0065) | (0.0134) | (0.0464) | (0.0371) | (0.0178) | |

3^{1}A_{g} | −0.4873 | −0.4825 ± 0.0045 | −0.4848 ± 0.0050 | −0.4747 ± 0.0043 | −0.4785 ± 0.0035 | −0.4795 ± 0.0043 |

(A) | (0.0048) | (0.0025) | (0.0126) | (0.0088) | (0.0078) | |

5^{1}A_{g} | −0.1040 | −0.0812 ± 0.0018 | −0.0779 ± 0.0017 | −0.0625 ± 0.0018 | −0.0592 ± 0.0020 | −0.0763 ± 0.0017 |

(A) | (0.0228) | (0.0261) | (0.0415) | (0.0448) | (0.0277) | |

6^{1}A_{g} | 0.0238 | … | … | 0.0405 ± 0.0016 | 0.0365 ± 0.0017 | 0.0368 ± 0.0024 |

(A) | (0.0167) | (0.0127) | (0.0130) | |||

7^{1}A_{g} | 0.2326 | 0.2525 ± 0.0023 | 0.2545 ± 0.0019 | 0.2623 ± 0.0019 | 0.2584 ± 0.0017 | 0.2556 ± 0.0020 |

(A) | (0.0199) | (0.0219) | (0.0297) | (0.0258) | (0.0230) |

^{a}

Characteristics of the electronic states are determined by the active-space contribution: A, dominated by active-space configurations and PA, dominated by configurations not belonging to active space.

^{b}

Full configuration interaction calculations were performed using all 30 molecular orbitals.

When the bare Hamiltonian for H_{2} at equilibrium is diagonalized in the active space, all three excited states can be observed. To improve the distribution and get a better statistical sampling for the excited states (particularly 3^{1}*A*_{g}), the initial wavefunction was chosen to be the $\Phi 11\xaf22\xaf$ configuration, although it is not necessary to track these excited states. Unfortunately, there is a significant difference between the FCI energies and energies obtained by diagonalizing the bare Hamiltonian with errors of over 18 mhartree for 2^{1}*A*_{g} and 50 mhartree for 3^{1}*A*_{g} and 5^{1}*A*_{g}.

The DUCC-ex formalism for the 2^{1}*A*_{g} state (DUCC-ex(2^{1}*A*_{g})), where the similarity transformation is driven by *σ*_{ext}(*K*) corresponding to the 2^{1}*A*_{g} state [see (23) and (24)], provides an excellent agreement with the FCI energies with an error less than 1 mhartree, a substantial improvement over the energy obtained with the bare Hamiltonian. From Table III, one can note that by using the trial wavefunction $12(\Phi 13+\Phi 1\xaf3\xaf)$, there is a chance to describe not only target 2^{1}*A*_{g} (which is the physical solution given the state-specific nature of the DUCC-ex approach) but also other states, which have nonzero overlap with the trial wavefunction (in this case 5^{1}*A*_{g}). For the 3^{1}*A*_{g} and 5^{1}*A*_{g} states, the effect of important out-of-active-space Slater determinants needs to be determined perturbatively and is only approximately captured by the corresponding *σ*_{ext}(*K*) operators. Consequently, the DUCC-ex(3^{1}*A*_{g}) and DUCC-ex(5^{1}*A*_{g}) results are less accurate compared to the DUCC-ex(2^{1}*A*_{g}) case. However, the results obtained with the respective DUCC-ex(3^{1}*A*_{g}) and DUCC-ex(5^{1}*A*_{g}) Hamiltonians are still better compared to the diagonalization of the bare-Hamiltonian. In addition, in the simulations for 3^{1}*A*_{g} and 5^{1}*A*_{g}, all three states are observed when the initial wavefunctions are taken to be the leading excitation for the corresponding state.

Compared to the equilibrium geometry, the H_{2} system corresponding to the *R*_{H–H} = 10.0 a.u. is characterized by a larger number of excited states with a significant overlap with active space. For *R*_{H–H} = 10.0 a.u., one can identify five excited states of *A*_{g} symmetry falling into this category. As shown in Table IV, the low-lying excited states can be efficiently described by the DUCC-ex formalism. For example, the DUCC-ex(2^{1}*A*_{g}) and DUCC-ex(3^{1}*A*_{g}) Hamiltonians provide very accurate estimates for the 2^{1}*A*_{g} and 3^{1}*A*_{g} states, which is quite remarkable given the size of the active space. In both cases, one can observe a significant improvement of the results obtained by the diagonalization of the bare Hamiltonian in the same active space, which provides a good illustration of the efficiency of the DUCC-ex formalism, even in a simple case corresponding to single commutator expansion of the DUCC-ex downfolded Hamiltonian. For higher excited states (5^{1}*A*_{g}, 6^{1}*A*_{g}, and 7^{1}*A*_{g}), the DUCC-ex results are less accurate, but still comparable to the results obtained through the diagonalization of the bare Hamiltonian in the active space. This behavior is justified given the complex excitation manifold describing the higher-lying states and ought to be resolved with improved approximations for *σ*_{ext}(*K*).

Table IV also provides an excellent illustration of the fact that the trial wavefunction can be used to probe various electronic states in quantum simulations (see also Fig. 2). For example, using the trial state $\Phi 11\xaf22\xaf$ for PE simulations of the DUCC-ex(2^{1}*A*_{g}) Hamiltonian, one can obtain a statistically meaningful presence of states different from a challenging doubly excited 2^{1}*A*_{g} state, i.e., 1^{1}*A*_{g} (not shown in Table IV), 3^{1}*A*_{g}, 5^{1}*A*_{g}, and 7^{1}*A*_{g}. This fact is associated with the strong quasidegeneracy of the ground-state and low-lying excited states with contribution from the $\Phi 11\xaf22\xaf$ Slater determinant. It is also worth mentioning that all states considered here are either purely doubly excited states or states of mixed single- and doubly excited characteristics, which usually pose a significant challenge for approximate EOMCC formulations. Analogous situations are naturally occurring in strongly correlated systems, and the fact that PE can deal with these challenging problems provides yet another argument in favor of developing quantum algorithms for excited states. Similar behavior can also be seen in the case of the active-space representation of the bare Hamiltonian. Although the DUCC-ex approach is state-specific in its nature, the appearance of other physically interpretable excited states in the spectra of DUCC-ex Hamiltonians should be explored further.

We exploited the state-specific nature of the DUCC-ex approach in an attempt to describe the lowest-lying fully symmetric 2^{1}*A*_{1} and 2^{1}*A*_{g} states (*vide infra*) of the H_{4} system in trapezoidal (*α* = 0.001) [H_{4}(a) system] and linear configuration [H_{4}(b) system] shown in Fig. 3. For both systems, FCI results were obtained by running EOMCCSDTQ calculations using *C*_{2$v$} and *D*_{2h} symmetries, respectively. The active space used to construct DUCC effective Hamiltonians in both configurations includes the two occupied orbitals and the five lowest-lying virtual orbitals. While the ground state for H_{4}(a) discloses strong quasidegeneracy effects between $\Phi $ and $\Phi 22\xaf33\xaf$ Slater determinants, for the H_{4}(b) system, these effects slightly weaker. In the description of the 2^{1}*A*_{1} state of the H_{4}(a) system, the dominant role is played by the $\Phi 22\xaf33\xaf$ determinant, which almost entirely dictates the corresponding wavefunction expansion. The 2^{1}*A*_{g} state of H_{4}(b) reveals a more multiconfigurational structure where determinants $\Phi 13$, $\Phi 1\xaf3\xaf$, $\Phi 24$, $\Phi 2\xaf4\xaf$, $\Phi 1234$, $\Phi 1\xaf2\xaf3\xaf4\xaf$, $\Phi 11\xaf33\xaf$, $\Phi 22\xaf33\xaf$, and $\Phi 22\xaf44\xaf$ all have a non-negligible contribution to the wavefunction. As shown in Table V, we are able to inspect the lowest-lying doubly excited state for both systems with the bare Hamiltonian. It is important to note that unlike the H_{2} systems, one must use an initial guess that contains the $\Phi 22\xaf33\xaf$ determinant with some significant weight. Otherwise, only the ground state will be observed. Even though these states can be tracked with the bare Hamiltonians, they come with errors of 73 mhartree for the H_{4}(a) system and 38 mhartree for the H_{4}(b) system. As shown in Table V, the DUCC-ex(2^{1}*A*_{1}) approach can very efficiently reproduce the FCI 2^{1}*A*_{1} energy with relative error less than 10 mhartree. As in the case of H_{2} system using $\Phi 22\xaf33\xaf$ as an initial state for quantum simulations, one can observe a statistical presence of states which have a non-negligible overlap with this determinant. It is interesting to note that for the H_{4}(b) system, the DUCC-ex results for the 2^{1}*A*_{g} are less accurate and they provide comparable results to the bare Hamiltonian energy. In our opinion, this is associated with the strong multiconfigurational characteristics of the corresponding wavefunction, where simple approximations (25) and (26) may not provide the desired level of accuracy. This statement can be supported by the analysis of lowest-lying orbitals (shapes of the active orbitals used in our calculations are shown in Fig. 4). As shown in Fig. 5, orbital energies of orbitals not included in the active space (orbitals 8,9, …) for H_{4}(a) system grow faster compared to the H_{4}(b) case, where energies of nonactive orbitals reach a plateau involving orbitals 8–13. This behavior suggests that much larger active spaces are needed for modeling H_{4}(b) with the DUCC-ex formalism.

State . | . | Bare Hamiltonian . | DUCC-ex(2^{1}A_{1})
. |
---|---|---|---|

(characteristics)^{a}
. | FCI^{b}
. | SD: $\Phi 22\xaf33\xaf$ . | SD: $\Phi 22\xaf33\xaf$ . |

α = 0.001 [H_{4}(a) system] | |||

2^{1}A_{1} | −2.0280 | −1.9546 ± 0.0037 | −2.0184 ± 0.0033 |

(A) | (0.0734) | (0.0096) | |

Linear [H_{4}(b) system] | |||

2^{1}A_{g} | −1.9901 | −1.9513 ± 0.0052 | −1.9458 ± 0.002894 |

(A) | (0.0388) | (0.0443) |

State . | . | Bare Hamiltonian . | DUCC-ex(2^{1}A_{1})
. |
---|---|---|---|

(characteristics)^{a}
. | FCI^{b}
. | SD: $\Phi 22\xaf33\xaf$ . | SD: $\Phi 22\xaf33\xaf$ . |

α = 0.001 [H_{4}(a) system] | |||

2^{1}A_{1} | −2.0280 | −1.9546 ± 0.0037 | −2.0184 ± 0.0033 |

(A) | (0.0734) | (0.0096) | |

Linear [H_{4}(b) system] | |||

2^{1}A_{g} | −1.9901 | −1.9513 ± 0.0052 | −1.9458 ± 0.002894 |

(A) | (0.0388) | (0.0443) |

^{a}

Characteristics of the electronic states are determined by the active-space contribution: A, dominated by active-space configurations and PA, dominated by configurations not belonging to active space.

^{b}

Full configuration interaction calculations were performed using all 60 molecular orbitals.

## V. CONCLUSIONS

In this paper, we discussed the excited-state extension of the DUCC formalism and its amenability for quantum computing. Using simple approximation schemes based on the utilization of the EOMCCSD wavefunction representation and lowest-rank contribution stemming from a single commutator electronic Hamiltonian and external *σ*_{ext}(*K*) operator, we were able to demonstrate that the active-space representation of downfolded Hamiltonians can be used to reproduce a large portion of excited-state correlation effects. For H_{2} and H_{4}(a) models, we observed significant improvements in excited-state energies compared to the diagonalization of the bare Hamiltonian in the active space, which was especially true for low-lying states dominated by active-space contributions (usually attributed to singly excited states and low-lying doubly excited states). As expected, for active-space dominated states characterized by higher excitation energies and more complicated configurational structure, the efficiency of simple approximation schemes discussed here deteriorates, which is indicative of the need for the inclusion of higher-order commutators and more efficient estimates of the *σ*_{ext}(*K*) operators. This issue will be explored in future research by coupling DUCC-ex formalism with higher-rank EOMCC formulations. We have also demonstrated that quantum phase estimation can provide an efficient tool for testing various “excited-state” hypotheses for strongly correlated systems, which usually pose a significant challenge for existing many-body formalism due to the need of inclusion of higher rank excitation effects and high density of states in a narrow energy gap, which renders the numerical identification of state of interest numerically unfeasible. Instead, using PE techniques, one can define hypothesis or trial state and obtain in the course of calculation statistical footprint of all states that have non-negligible overlap with the hypothesis state. Although this fact is a well-known foundation of quantum computing, it deserves broader exposure. The stochastic characteristics of PE may be very useful in studies of excited-state processes, especially in the strongly correlated or metallic-like system.

## ACKNOWLEDGMENTS

This work was supported by the “Embedding Quantum Computing into Many-body Frameworks for Strongly Correlated Molecular and Materials Systems” project, which is funded by the U.S. Department of Energy (DOE), Office of Science, Office of Basic Energy Sciences, the Division of Chemical Sciences, Geosciences, and Biosciences. A portion of this research was funded by the Quantum Algorithms, Software, and Architectures (QUASAR) Initiative at the Pacific Northwest National Laboratory (PNNL). It was conducted under the Laboratory Directed Research and Development Program at PNNL. All calculations have been performed using the Molecular Science Computing Facility (MSCF) in the Environmental Molecular Sciences Laboratory (EMSL) at the Pacific Northwest National Laboratory (PNNL). EMSL is funded by the Office of Biological and Environmental Research in the U.S. Department of Energy. PNNL is operated for the U.S. Department of Energy by the Battelle Memorial Institute under Contract No. DE-AC06-76RLO-1830.