QuaSiMo: A Composable Library to Program Hybrid Workflows for Quantum Simulation

We present a composable design scheme for the development of hybrid quantum/classical algorithms and workflows for applications of quantum simulation. Our object-oriented approach is based on constructing an expressive set of common data structures and methods that enable programming of a broad variety of complex hybrid quantum simulation applications. The abstract core of our scheme is distilled from the analysis of the current quantum simulation algorithms. Subsequently, it allows a synthesis of new hybrid algorithms and workflows via the extension, specialization, and dynamic customization of the abstract core classes defined by our design. We implement our design scheme using the hardware-agnostic programming language QCOR into the QuaSiMo library. To validate our implementation, we test and show its utility on commercial quantum processors from IBM and Rigetti, running some prototypical quantum simulations.


INTRODUCTION
Quantum simulation is an important use case of quantum computing for scientific computing applications. Whereas numerical calculations of quantum dynamics and structure are staples of modern scientific computing, quantum simulation represents the analogous computation based on the principles of quantum physics. Specific applications are wide-ranging and include calculations of electronic structure [1][2][3][4], scattering [5], dissociation [6], thermal rate constants [7], materials dynamics [8], and response functions [9]. This manuscript has been authored by UT-Battelle, LLC under Contract No. DE-AC05-00OR22725 and Lawrence Berkeley National Laboratory under Contract No. DE-AC02-05CH11231 with the U.S. Department of Energy. The United States Government retains and the publisher, by accepting the article for publication, acknowledges that the United States Government retains a non-exclusive, paid-up, irrevocable, world-wide license to publish or reproduce the published form of this manuscript, or allow others to do so, for United States Government purposes. The Department of Energy will provide public access to these results of federally sponsored research in accordance with the DOE Public Access Plan. (http://energy.gov/downloads/doe-publicaccess-plan).
Presently, this diversity of quantum simulation applications is being explored with quantum computing despite the limitations on the fidelity and capacity of quantum hardware [10][11][12]. These applications are tailored to such limitations by designing algorithms that can be tuned and optimized in the presence of noise or model representations that can be reduced in dimensionality. Examples include variational methods such as the variational quantum eigensolver (VQE) [13][14][15][16], quantum approximate optimization algorithm (QAOA) [17], quantum imaginary time evolution (QITE) [18], and quantum machine learning (QML) among others.
The varied use of quantum simulation raises concerns for efficient and effective programming of these applications. The current diversity in quantum computing hardware and low-level, hardware-specific languages imposes a significant burden on the application user. For instance, IBM provides Aqua [19], which is part of the Qiskit framework, targeting high-level quantum applications such as chemistry and finance. The emphasis of Aqua is on providing robust implementations of quantum algorithms, yet the concept of reusable and extensible workflows, especially for quantum simulations, is not formally supported. The user usually has to implement custom workflows from lower-level constructs, such as circuits and operators, available in Qiskit Terra [20]. Similarly, Tequila [21] is another Python library that provides commonly-used functionalities for quick prototyping of variational-type quantum algorithms. Orquestra [22], a commercially-available solution from Zapata, on the other hand, orchestrates quantum application workflows as blackbox phases, thus requires users to provide implementations for each phase.
The lack of a common workflow for applications of quantum simulation hinders broader progress in testing and evaluation of such hardware. A common, reusable and extensible programming workflow for quantum simulation would enable broader adoption of these applications and support more robust testing by the quantum computing community.
In this contribution, we address development of common workflows to unify applications of quantum simulation. Our approach constructs common data structures and methods to program varying quantum simulation applications, and we leverage the hardware-agnostic language QCOR and programming framework XACC to implement these ideas. We demonstrate these methods with example applications from materials science and chemistry, and we discuss how to extend these workflows to experimental validation of quantum computation advantage, in which numerical simulations can benchmark programs for small-sized models [12,[23][24][25][26].
Upon release of the QuaSiMo library, we became aware of the additional application modules available to the Qiskit framework which share a number of features with our proposed composable workflow design. In particular, the Qiskit Nature module [27] introduces new concepts to model and solve quantum simulation problems, and reiterates the need for modular, domain-specific, and workflow-driven quantum application builders, which the QuaSiMo work addresses.
The key differentiators for our work lie in the performance and extensibility of the implementation. Since QuaSiMo is developed on top of the C++ QCOR compiler infrastructure, it can take advantage of optimal performance for classical computing components of the workflow, such as circuit construction [28] or post-processing. Moreover, the pluginbased extension model of QuaSiMo is distinct from that of Qiskit Nature. Rather than requiring new extensions being imported, QuaSiMo puts forward a common interface for all of its extension points [29], and therefore enables the development of portable user applications w.r.t. the underlying library implementation.

SOFTWARE ARCHITECTURE
Cloud-based access to quantum computing naturally differentiates programming into conventional and quantum tasks [31,32]. The resulting hybrid execution model yields a loosely integrated computing system by which common methods have emerged for programming and data flow. We emphasize this concept of workflow to organize programming applications for quantum simulation. Figure 1 shows the blueprint of our Quantum Simulation Modeling (QuaSiMo) library. The programming workflow is defined by a QuantumSimulationWorkflow concept which encapsulates the hybrid quantum-classical procedures pertinent to a quantum simulation, e.g., VQE, QAOA, or dynamical quantum simulation. A quantum simulation workflow exposes an execute method taking as input a Quantum-SimulationModel object representing the quantum model that needs to be simulated. This model captures quantum mechanical observables, such as energy, spin magnetization, etc., that we want the workflow to solve or simulate for. In addition, information about the system Hamiltonian, if different from the observable operator of interest, and customized initial quantum state preparation can also be specified in the QuantumSimulationModel.
By separating the quantum simulation model from the simulation workflow, our object-oriented design allows the concrete simulation workflow to simulate rather generic quantum models. This design leverages the ModelFactory utility, implementing the object-oriented factory method pattern. A broad variety of input mechanisms, such as those provided by the QCOR infrastructure or based on custom interoperability wrappers for quantum-chemistry software, can thus be covered by a single customizable polymorphic model. For additional flexibility, the last createModel factory method overload accepts a polymorphic builder interface ModelBuilder the implementations of which can build arbitrarily composed QuantumSimulationModel objects.
QuantumSimulationWorkflow is the main extension point of our QuaSiMo library. Built upon the CppMicroServices framework conforming to the Open Services Gateway Initiative (OSGi) standard [33], QuaSiMo allows implementation of a new quantum workflow as a plugin loadable at runtime. At the time of this writing, we have developed the QuantumSimulationWorkflow plugins for the VQE, QAOA, QITE, and time-dependent simulation algorithms, as depicted in Fig. 1. All these plugins are implemented in the QCOR language [28,34] using the externally-provided library routines.
At its core, a hybrid quantum-classical workflow is a procedural description of the quantum circuit composition, pre-processing, execution (on hardware or simulators), and post-processing. To facilitate modularity and reusability in workflow development, we put forward two concepts, AnsatzGenerator and CostFunctionEvaluator. AnsatzGenerator is a helper utility used to generate quantum circuits based on a predefined method such as the Trotter decomposition [35,36] or the unitary coupled-cluster (UCC) ansatz [37]. CostFunctionEvaluator automates the process of calculating the expectation value of an observable operator. For example, a common approach is to use the partial state tomography method of adding change-ofbasis gates to compute the operator expectation value in the Z basis. Given the CostFunctionEvaluator interface, quantum workflow instances can abstract away the quantum backend execution and the corresponding post-processing of the results. This functional decomposition is particularly advantageous in the NISQ regime since one can easily integrate the noise-mitigation techniques, e.g., the verified quantum phase estimation protocol [38], into the QuaSiMo library, which can then be used interchangeably by all existing workflows.
Finally, our abstract QuantumSimulationWorkflow class also exposes a public validate method accepting a variety of concrete implementations of the abstract Quantum-ValidationModel class via a polymorphic interface. Given the quantum simulation results produced by the execute Figure 1: The class UML diagram of the quantum simulation application. The fully typed version is provided separately (see [30]). method of QuantumSimulationWorkflow, the concrete implementations of QuantumValidationModel must implement its accept _ results method based on different validation protocols and acceptance criteria. For example, the acceptance criteria can consist of distance measures of the results from previously validated values, or from the results of validated simulators. The measure may also be taken relative to experimentally obtained data, which, with sufficient error analysis to bound confidence in its accuracy, can serve as a ground truth for validation. A more concrete example in a NISQ workflow includes the use of the Quan-tumSimulationWorkflow class to instantiate a variational quantum eigensolver simulator, followed by the use of validate to instantiate a state vector simulator. Results from both simulators can be passed to the QuantumValidation-Model accept _ results method which evaluates a distance measure method and optionally calls a decision method which returns a binary answer. Other acceptance criteria include evaluation of formulae with input data, application of curve fits, and user-defined criteria provided in the concrete implementation of the abstract QuantumValidationModel class. The validation workflow relies on the modular architecture of our approach, which effectively means that writing custom validation methods and constructing userdefined validation workflows is achieved by extending the abstract QuantumValidationModel class.
In our opinion, the proposed object-oriented design is well-suited to serve as a pattern for implementing diverse hybrid quantum-classical simulation algorithms and workflows which can then be aggregated inside a library under a unified object-oriented interface. Importantly, our standardized polymorphic design with a clear separation of concerns and multiple extension points provides a high level of composability to developers interested in implementing rather complex quantum simulation workflows.

TESTING AND EVALUATION
Our implementation of the programming workflow for applications of quantum simulation is available online [39]. We have tested this implementation against several of the original use cases to validate the correctness of the implementation and to evaluate performance considerations.

Dynamical Simulation
As a first sample use case, we consider a non-equilibrium dynamics simulation of the Heisenberg model in the form of a quantum quench. A quench of a quantum system is generally carried out by initializing the system in the ground state of some initial Hamiltonian, , and then evolving the system through time under a final Hamiltonian, . Here, we demonstrate a simulation of a quantum quench of a onedimensional (1D) antiferromagnetic (AF) Heisenberg model using the QCOR library to design and execute the quantum circuits. Our AF Heisenberg Hamiltonian of interest is given by where > 0 gives the strength of the exchange couplings between nearest neighbor spins pairs ⟨ , ⟩, > 0 defines the anisotropy in the system, and is the -th Pauli operator acting on qubit . We choose our initial Hamiltonian to be the Hamiltonian in equation 1 in the limit of → ∞. Thus, setting = 1, = +1 , where is an arbitrarily large constant. The ground state of is the Néel state, given by which is simple to prepare on the quantum computer. We choose our final Hamiltonian to have a finite, positive value of , so observable of interest is the staggered magnetization [40], which is related to the AF order parameter and is defined as using namespace QuaSiMo; // AF Heisenberg model auto problemModel = ModelFactory::createModel( "Heisenberg", {{"Jx", 1.0}, {"Jy", 1.0}, // Jz == g parameter {"Jz", g}, // No external field {"h _ ext", 0.0}, {"num _ spins", n _ spins}, {"initial _ spins", initial _ spins}, {"observable", "staggered _ magnetization"}}); // Time-dependent simulation workflow auto workflow = std::make _ shared<TimeDependentWorkflow>(); workflow->initialize({{"dt", dt}, {"steps", n _ steps}}); // Execute the workflow auto result = workflow->execute(problemModel); // Get the observable values // (staggered magnetization) auto obsVals = result.get<std::vector<double>>("exp-vals"); where is the number of spins in the system. Fig. 2 shows sample results for = 9 spins for a three different values for in . The qualitatively different behaviours of the staggered magnetization after the quench for < 1 and > 1 are apparent, and agree with previous studies [40]. We present a listing of the code expressing this implementation in Fig. 3.
We develop QuaSiMo on top of the QCOR infrastructure, as shown in Fig. 1; thus, any quantum simulation workflows constructed in QuaSiMo are retargetable to a broad range of quantum backends. The results that we have demonstrated in Fig. 2 are from a simulator backend. The same code as shown in Fig. 3 can also be recompiled with a -qpu flag to target a cloud-based quantum processor, such as those available in the IBMQ network.
Currently available quantum processors, known as noisy intermediate-scale quantum (NISQ) computers [41], have relatively high gate-error rates and small qubit decoherence times, which limit the depth of quantum circuits that can be executed with high-fidelity. As a result, long-time dynamic simulations are challenging for NISQ devices as current algorithms produce quantum circuits that increase in depth with increasing numbers of time-steps [42]. To limit the circuit size, we simulated a small AF Heisenberg model, eq. 1, with only three spins on the IBM's Yorktown (ibmqx2) and Casablanca (ibmq _ casablanca) devices.
The simulation results from real quantum hardware for values of 0.0 and 4.0 are shown in Fig. 4, where we can see the effects of gate-errors and qubit decoherence leading to a significant impairment of the measured staggered magnetization (circles) compared to the theoretical values (solid lines). In particular, Fig. 4 demonstrates how the quality of the quantum hardware can affect simulation performance. The Yorktown backend has considerably worse performance metrics than the Casablanca backend 1 . Specifically, compared to the Casablanca backend, Yorktown has a slightly higher two-qubit gate-error rate, nearly double the read-out error rate, and substantially lower qubit decoherence times. While identical quantum circuits were run on the two machines, we see much better distinguishability between the results for the two values of in the results from Casablanca than those from Yorktown.
The staggered magnetization response to a quench for a simple three-qubit AF Heisenberg model in Fig. 4, albeit noisy, illustrates non-trivial dynamics beyond that of decoherence (decaying to zero). Improvements in circuit construction (Trotter decomposition) and optimization, noise mitigation, and, most importantly, hardware performance (gate fidelity and qubit coherence) are required to scale up this time-domain simulation workflow for large quantum systems.

Variational Quantum Eigensolver
As a second use case demonstration, we apply the Variational Quantum Eigensolver (VQE) algorithm to find the ground state energy of H 2 . The VQE is a quantum-classic hybrid algorithm used to find a Hamiltonian's eigenvalues, where the quantum process side is represented by a parametrized quantum circuit whose parameters are updated by a classical optimization process [43]. The algorithm updates the quantum circuit parameters to minimize the Hamiltonian's expectation value until it converges. The performance of the VQE algorithm, as any other quantum-classical variational algorithms [44], depends in 1   part on the selection of the classical optimizer and the circuit ansatz. The design scheme implemented in this work allows us to tune the VQE components to pursue better performance. We present a listing of the code expressing this implementation in Fig. 5, in which we define the different parameters of the VQE algorithm in a customtailored way. In the code snippet, @qjit is a directive to activate the QCOR just-in-time compiler, which compiles the kernel body into the intermediate representation, and  QuaSiMo.getWorkflow is a utility function of the library to construct and initialize workflow objects of various types. In Fig. 6(a), we present simulations considering the Simultaneous Perturbation Stochastic Approximation (SPSA) for ansatz updating. There, we show the energy as a function of the quantum circuits used in the learning process with a budget of 200 function evaluations with 5000 shots (experimental repetitions) per evaluation by using the QCOR's VQE module, QuaSiMo.getWorkflow('vqe'), and by using the Qiskit class aqua.algorithms.VQE [19]. In both cases, we consider (0., 0., 0.) as the initial set of parameters, therefore the energy values in the first iterations are similar, due to the SPSA hyperparameter's calibration. As it is expected, since we are using the same optimizer in both quantum simulations, there is not a relevant difference between the learning paths. However, the execution time in QCOR is shorter. This feature is presented in Fig. 6(b), where we show the execution time rate between QCOR and Qiskit as a function of the classical optimization algorithm. To evaluate the typical execution time of each optimizer, we consider 51 runs of VQE towards the generation of 2 quantum ground state.

Symmetry reduction.
The presence of symmetries, such as rotations, reflections, number of particles, etc., in the Hamiltonian model allows us to map the model to a model with fewer qubits [45]. In 2 , we apply the QCOR's function operatorTransform('qubit-tapering', H) to reduce the four-qubit Hamiltonian, see Fig. 5, to a onequbit Hamiltonian model. We present this implementation in Fig. 7, in which we transform the Hamiltonian introduced in Fig. 5 into a one-qubit model, and we redefine the anzats. In the output section of Fig. 7, we present the one-qubit model and the VQE's output that converge to a similar value of the four-qubit model.

Fermion-qubit map.
An important feature included in the QCOR compiler is the fermion-to-qubit mapping that facilitates the quantum state searching in VQE. In Fig. 8, we present an example of how to use OpenFermion operators [47] in the VQE workflow. In that implementation, we define the ansatz by using Scipy and OpenFermion, the QCOR compiler decomposes the ansatz into quantum gates; we follow the same structure presented in Fig. 5 for the VQE workflow.

Quantum Approximate Optimization Algorithm
To further demonstrate the utility of QCOR we present an implementation of the quantum approximate optimization algorithm (QAOA) [17]. QAOA translates a classical cost function into a quantum operator then uses a variational quantum-classical optimization loop to find quantum states that minimize the expectation value ⟨ ⟩. The optimized quantum states are then prepared and measured to obtain bitstrings that correspond to classical solutions to the optimization problem. Figure 9 shows example python code that uses QCOR simulations to find optimized quantum states with QAOA. Opening lines define the number of qubits and construct the problem Hamiltonian for MaxCut on a star graph . The Hamiltonian is then used to create the QuaSiMo from qcor import * @qjit def ansatz(q : qreg, x : List[float]): X(q[0]) with decompose(q, kak) as u: from scipy.sparse.linalg import expm from openfermion.ops import QubitOperator from openfermion.transforms import \ get _ sparse _ operator qop = QubitOperator('X0 Y1') \ -QubitOperator('Y0 X1') qubit _ sparse = get _ sparse _ operator(qop) u = expm(0.5j * Here we depict how to use OpenFermion operators to construct the state-preparation kernel (ansatz) for the VQE workflow. We use SciPy [46] and OpenFermion [47] to construct the exponential of ( 0 1 − 0 1 ) operator as a matrix. The matrix will be decomposed into quantum gates by the QCOR compiler.
model and workflow to simulate -step QAOA and return the expectation ⟨ ⟩, similar to the previous VQE example of Fig. 5. Figure 10 shows optimized energies ⟨ ⟩ computed in QCOR for star graphs with numbers of qubits = 2, ..., 9 at various QAOA depth parameters . Each time the program is run it begins with random initial parameters to be optimized, so we used several random starts for some of the graphs, keeping the smallest result as the ⟨ ⟩ shown in the figure. The QCOR results are in perfect agreement with the optimized standards [48] of Lotshaw et al. [49]. The results have interesting features, for example when is odd perfect ground state energies ⟨ ⟩ = 0 are obtained at = 2, while when is even the results need = 3 to approach close to 0 . The simplicity of the QCOR code and simulations make this an attractive route to studying interesting features like these in future research on QAOA.    Δ , we approximate the imaginary time evolution based on the results of the previous steps.

Quantum Imaginary Time Evolution Algorithm
As shown in Fig. 1, the QITE workflow is integrated into the QuaSiMo library. In this example, we demonstrate the use of QITE to find the ground state energy of a three-qubit transverse field Ising model (TFIM), where = ℎ = −1.0 and = 3 is the number of spins (qubits). Fig. 11 is the code snippet to set up the QuaSiMo problem description and workflow for this problem. The problem model is captured by a single Hamiltonian operator constructed by direct Pauli operator algebra. We also want to note that the Pauli operator algebra of QuaSiMo allows for hierarchical construction of the Hamiltonian, e.g., via for loops, suitable for generic Hamiltonians similar to the one described in eq. (3).
For QITE workflow configurations, we set step-size to 0.45 and steps to 20, for a total imaginary time of = 9.0 (ℏ = 1). The system begins in different initial states by using a state-preparation circuit, as shown in Fig. 11. Similar to the Python API, users can also use the utility function getWorkflow to retrieve an instance of the QITE workflow from the QCOR service registry with the name key "qite" as shown in Fig. 11. The main drawback of QITE algorithm is that the propagating circuit size increases during the imaginary timestepping procedure. To alleviate this constraint, especially for execution on NISQ hardware, QuaSiMo's QITE workflow implementation support custom, externally-provided circuit optimizers that will be invoked during the algorithm execution to minimize the circuit depth. In this demonstration, we take advantage of the QSearch [51] optimizer from BQSKit library, which is capable of synthesizing constant-depth circuits for a variety of common use cases including the TFIM model in (3). For instance, in this example, the final QITE circuit (step = 20) has approximately 8,000 gates (3,000 CNOT gates), which is clearly beyond the capability of current NISQ devices. Thanks to QSearch, we can always re-synthesize a constant depth circuit with only 14 CNOT gates (87 gates in total) for any of the time steps, which is the theoretical lower bound [52], 1 4 (4 − 3 − 1) = 13.5, for CNOT gate count in three-qubit circuits. The results of the QITE workflow execution are shown in Fig. 12, where we can see the energy value exponentially decays to the analytically computed ground state energy of -3.49396 for all initial states.

CONCLUSIONS
We have presented and demonstrated a programming workflow for applications of quantum simulation that promotes common, reusable methods and data structures for scientific applications.
We note that while the framework presented here is readily applicable to use cases in the NISQ era -with the Quan-tumSimulationWorkflow in particular being extendable to NISQ simulation algorithms such as VQE -the workflow is also extendable to universal algorithms. Code portability is ensured through intermediate representation that is then implemented for each backend according to specific APIs, which are accessible in the QuaSiMo library. This enables rapid porting of a simulation algorithm to multiple machines (write once, run everywhere paradigm). This is an ideal environment in which to construct both benchmarks and validation protocols, along with short depth quantum simulations that can quickly be run on multiple processors and directly compared. While we maintain a focus on rapid prototyping of quantum simulation on today's NISQ devices, we note that the nature of the intermediate representation and the modular backend structure enable targeting fault tolerant devices as well. A fault tolerant (FT) architecture may be represented as an additional backend with encoding transpiler preprocessing. While this configuration is ideal for QuaSiMo's unencoded qubit targets, we also note that the parent language QCOR, is capable of expressing fully quantum-errorcorrection-encoded algorithms as well. Therefore, we expect the framework to be extendable and to find use in workflows involving universal or FT applications in the future.