Q²Chemistry: A quantum computation platform for quantum chemistry

1.   Introduction

As the application of accurate classical methods is severely limited by the fast growing computational cost, quantum computation provides a promising pathway to solve quantum chemistry problems^[1-3]. By encoding wave functions into the Hilbert space of qubits, the Schrödinger equation for molecular systems can be solved on a quantum computer. In recent years, various quantum algorithms, such as quantum phase estimation (QPE) and variational quantum eigensolver (VQE)^[2-19], are developed for quantum chemistry. There are some important issues that should be studied for these algorithms. For example, the basis set used to construct the electronic Hamiltonian and the fermion-to-qubit encoding technique may lead to significant gate complexity and measurement overhead in QPE^[20]. The optimization of a typical wave function ansatzes using VQE can suffer from “barren plateaus”^[21,22] or local minimum traps^[23]. A quantum computation platform that can provide extensive functionalities and step-by-step cross verification is therefore desirable for designing novel quantum algorithms for chemistry applications.

Noisy intermediate quantum (NISQ) devices have been used to demonstrate the possibility of studying the ground and excited states of molecular systems using currently available algorithms^{[10-15,24,25]}. However, such experimental demonstrations are limited to tiny molecules with an artificially small basis set. This is because NISQ experiments are limited by the available quantum resources and the error associated with each quantum gate. Therefore, low gate fidelity, short coherence time, and insufficient qubit resource prohibit a systematic study of quantum chemistry oriented quantum algorithms on NISQ devices. The largest quantum computation experiment for chemistry to date uses 16 qubits with 160 two-qubit gates^[25], while a simple VQE circuit of the commonly used unitary coupled-cluster (UCC)^[26,27] ansatz for a small molecule has millions of CNOT gates (Table 1), which is far beyond the capability of NISQ devices. Therefore, it is important to have the capability to simulate the quantum circuit on a classical computer at this stage. Even beyond the NISQ era, such a capability can help us to develop quantum-inspired classical algorithms.

Table  1.  Computational resources required to perform VQE simulations for a number of molecules using the unitary coupled-cluster ansatz truncated up to double excitations (UCCSD) and the minimum basis set STO-3G. UCCGSD means that generalized excitation operators (not distinguishing occupied and virtual orbitals) are used. The way to count CNOT gates is described in the supporting information.

Molecule	Qubits	Parameters		CNOT gates
		UCCSD	UCCGSD	UCCSD	UCCGSD
H$_2$	4	2	5	64	144
LiH	8	14	72	1632	10720
H$_2$O	14	90	630	26272	$2.1\times 10^5$
NH$_3$	16	135	1064	46480	$4.4\times 10^5$
CH$_4$	18	230	1692	95200	$8.2\times 10^5$
C$_2$H$_4$	28	1224	9737	$8.6\times 10^5$	$7.9\times10^6$
C$_3$H$_6$	42	5994	48930	$6.6\times 10^6$	$6.3\times 10^7$

 | Show Table

DownLoad: CSV

Several quantum computation packages have been reported, for example, C++ based ProjectQ^[28] and Qiskit^[29], GPU-enabled Qulacs^[30], and the differentiable simulator Yao^[31] implemented in Julia. Most of these codes are developed as standalone quantum circuit simulators or compilers, which are not dedicated to quantum chemistry applications. At the same time, the performance of simulating the quantum circuit is not very satisfactory in these packages. Most of these packages only implement the brute-force simulating method, which leads to an exponential computational cost. Tensor-based methods have been implemented in packages such as Qiskit and PennyLane^[32], however, without an efficient distributed parallelization algorithm. As a reference, the largest simulated quantum circuit to date contains 28 qubits, which is used to study the ground state of ethylene molecules using VQE^[33].

Based on the above considerations, we develop a versatile and extensible quantum computation platform for quantum chemistry, named Q²Chemistry (pronounced as “Q squared chemistry”) to highlight the two quantum dimensions (the systems to be studied and the tools used to study them). Q²Chemistry adopt a modular design and a mixed-language programming model to achieve versatility together with high performance, where Python is used as its application program interface (API) while C++ and Julia^[34] are used for computationally intensive tasks such as quantum circuit simulators. Q²Chemistry provides interfaces to quantum chemistry packages such as PySCF^[35] to generate the required parameters for the qubit Hamiltonian. The quantum algorithms for solving the eigenstates of the electronic Hamiltonian are programmable at a high level of abstraction by adopting internal templates and/or defining custom circuit generation procedures. Some popular VQE and post-VQE algorithms have already been implemented for ground and excited states of molecular and periodic systems. Q²Chemistry provides native modules to run the circuit generated by a specific quantum algorithm, either via various high performance classical circuit simulation algorithms for simulating the circuit on a classical computer or via extensible interfaces reserved for the upcoming actual quantum processors. An efficient matrix product state (MPS)^[36,37] based circuit simulation engine is implemented to perform quantum simulations for systems with 50–100 qubits.

The remainder of this article is organized as follows. After introducing the framework of Q²Chemistry in Section 2, we present functionalities for handling ab-initio chemistry quantities in Section 3. In Section 4 we explain the implementation of the circuit simulation on a classical computer or running on the upcoming quantum hardware. Section 5 provides a general introduction to the natively implemented VQE-based algorithms. Finally, in Section 6, we provide some examples to demonstrate the power of Q²Chemistry in chemistry applications, and provide the road map for further extensions including the support for more quantum algorithms, classical circuit simulation backends, and circuit optimization algorithms.

2.   The framework of Q²Chemistry

As shown in Fig. 1a, Q²Chemistry contains three modules. (ⅰ) $\mathtt{q2chem.qchem}$: the quantum chemistry module that defines quantum chemistry problems in the qubit space typically with the help of external classical quantum chemistry packages; (ⅱ) $\mathtt{q2chem.qcirc}$: the quantum circuit module that provides the quantum circuit related functionalities, including circuit construction, visualization, optimization, and execution on a virtual or real quantum computer; (ⅲ) $\mathtt{q2chem.qalgo}$: the quantum algorithm module, which includes native quantum algorithms and provides tools to implement new algorithms to solve chemistry problems.

Figure  1.  (a) The framework of Q²Chemistry. (b) A typical workflow of solving a chemical problem using a quantum algorithm.

DownLoad: Full-Size Img PowerPoint

High-level modules in Q²Chemistry are implemented using the scripting language Python and provide base implementations for classes in submodules. Benefiting from the modular design and Python’s extensibility, extending existing submodules are easily achieved by constructing derived classes and implementing only a small number of virtual functions, which requires no modifications to higher level modules. Core functions in the backends are implemented using programming languages including C++ and Julia. These functions are integrated into the Python interfaces using just-in-time (JIT) technology, which compiles the code at runtime to provide machine-specific optimizations and deliver outstanding performance. Most of the low-level data structures accessible from Python are stored using NumPy’s^[38] $\mathtt{ndarray}$. Therefore, auxiliary operations, such as exact diagonalization for the qubit Hamiltonian, can be realized by using NumPy-compatible packages such as SciPy^[39] or PyTorch^[40].

Using such a modular framework, a general workflow to solve a quantum chemistry problem on a quantum computer is briefly illustrated in Fig. 1b. (ⅰ) Collect classical data such as electron integrals and mean-field orbital coefficients to generate Hamiltonian and an initial quantum state. (ⅱ) Choose or develop a suitable quantum algorithm and generate corresponding quantum circuits, which may depend on the Hamiltonian, such as in the QPE algorithm, or be system independent, such as in some VQE algorithms. (ⅲ) Perform lower-level circuit optimizations, including eliminating redundant quantum gates to reduce circuit depth or reconstructing the circuit to fit a specific quantum processor. (ⅳ) Execute the quantum circuits and perform some measurements to extract necessary information, using either a virtual or real quantum computer.

Currently, Q²Chemistry can be routinely used for performing VQE simulations on a classical computer, which is powered by Hamiltonian generation for molecular and periodic systems, hardware-efficient and UCC-based ansatzes for parametric circuit construction, and a scalable noise-free tensor network backend for classical simulations of large quantum circuits.

3.   The quantum chemistry module

The $\mathtt{q2chem.qchem}$ module mainly handles the system Hamiltonian and wave function mapping. The second-quantized Hamiltonian is constructed from classically calculated quantities such as molecular orbital coefficients. Wave function mapping determines how the orbitals of the simulated wave function are mapped onto qubits therefore influences the measurement strategy.

3.1   Hamiltonian

For a Hamiltonian that is expressed as the linear combination of Pauli strings

where $\hat{P}_{i}$ is the product of Pauli operators $\{\hat{\sigma}_{x}, \hat{\sigma}_{y}, \hat{\sigma}_{z}\}^{\otimes}$, and the expectation value $E$ can be obtained through quantum measurement techniques such as the Hadamard test as

Given the Hartree-Fock orbitals, the second-quantized electronic Hamiltonian can be written as

with

$h_{pq}$ and $g^{pq}_{rs}$ are one- and two-electron integrals. To convert the second-quantized Hamiltonian Eq. (3) into the qubit form as given in Eq. (1), a fermion-to-qubit mapping such as the Jordan-Wigner or Bravyi-Kitaev mapping is required^[41-43].

In Q²Chemistry, an interface with the PySCF package is provided to calculate one- and two-electron integrals for molecular and periodic systems^[44,45]. Orbital optimization is also supported by linking to the $\mathtt{pyscf.lo}$ module or using a custom cost function. A unitary matrix is then obtained to transform the integrals as follows:

where $U^{*}$ is the elementwise complex conjugate of $U$. This step is carried out efficiently by calling the optimized tensor contraction package $\mathtt{opt}$_$\mathtt{einsum}$^[46]. For fermion-to-qubit mapping, Q²Chemistry implements an efficient Jordan-Wigner transformation written in pure Julia. Other methods, such as Bravyi-Kitaev, are currently provided through the interface of OpenFermion^[47]. The qubit Hamiltonian can then be used to construct quantum gates for measurement in quantum algorithms such as VQE or QPE.

3.2   Wave function

To represent a many-electron quantum state on a quantum computer, the most commonly used strategy is to map its molecular orbitals onto qubits, which can span the Fock space. Such a straightforward orbital-to-qubit approach can be written as

where $i_{j} \in \{0, 1\}$ represents both the occupation of orbitals and the quantum state $|0\rangle$ or $|1\rangle$ of the corresponding qubit. Eq. (6) actually describes the quantum state corresponding to the eigenstates of a qubit Hamiltonian obtained from the Jordan-Wigner transformation. If other fermion-to-qubit mapping algorithms such as Bravyi-Kitaev are used for the Hamiltonian, such a correspondence does not necessarily exist. With a quantum state mapped to qubits, Q²Chemistry can provide a Hadamard test to evaluate the expectation value of the quantum state with respect to an operator.

In simple orbital-to-qubit mapping, the number of qubits required to simulate the wave function has a linear dependence on the number of basis functions, which prohibits the use of a large basis set on NISQ devices. Q²Chemistry provides another strategy that maps a classical tensor network (TN) state onto quantum circuits^[48,49]. In this way, the qubits determine the classical bond dimension of the tensor network. For chemical systems that contain weak electron correlations or have a special symmetry, such a TN-based strategy provides a possible solution to effectively reduce the number of qubits at the expense of performing more measurements. Generally, different wave function mapping strategy leads to distinctive structures of the quantum circuits, therefore, may bring special restrictions to subsequent quantum algorithms and measurement protocols. The $\mathtt{q2chem.qchem}$ module passes the mapping strategy to quantum algorithms in $\mathtt{q2chem.qalgo}$ to generate an abstract circuit class for initialization. The abstract circuits are then extended according to the adopted quantum algorithm and instantiated to a common readable quantum circuit that can be used in $\mathtt{q2chem.qcirc}$.

4.   The quantum circuit module

The $\mathtt{q2chem.qcirc}$ module provides quantum circuit execution functionalities. It mainly contains two parts: an interface reserved for quantum computer manufacturers and a classical simulator that simulates state evolution determined by quantum circuits on a classical computer.

4.1   Interfaces with quantum devices

Functionalities in $\mathtt{q2chem.qcirc}$ are inherited from the _$\mathtt{Base}$ class. High-level Python APIs enable efficient extensions of Q²Chemistry, which can be used to interface with different quantum device operation systems. Currently, there are multiple competing technical routes of quantum computers, e.g., photonic qubits have long coherence time, while superconducting platforms have good scalability. To efficiently and accurately study quantum chemical problems, various quantum hardwares may be required depending on the characteristics of the quantum algorithm. Therefore, we provide a flexible interface for connecting to different quantum hardware platforms as illustrated in Code Example 1. In a general _$\mathtt{BaseHardware}$, the circuit generated by Q²Chemistry is first represented by quantum assembly language such as OpenQASM via _$\mathtt{compile}$_$\mathtt{circuit()}$. This step generally should include circuit optimization and reconstruction to fit the architecture of specific hardware. Then, the sequence of gates is converted into signals or pulses and sent to the quantum devices through _$\mathtt{quantum}$_$\mathtt{hardware}$_$\mathtt{evolution()}$ using a compiler provided by the vendor of the hardware system. Finally, the results such as measurement statistics or bit strings are collected and postprocessed in Q²Chemistry to obtain the required quantities.

4.2   Classical quantum circuit simulator

Q²Chemistry implements several simulators to run quantum circuits on a classical computer. The quantum state of multiple qubits can be expanded using a certain basis set

where $N$ is the number of qubits, and $|i_1 i_2 \ldots i_N\rangle$ is the basis state. The coefficients $c_{i_1 i_2 \ldots i_N}$ form an $N$-dimensional tensor that contains $2^N$ amplitudes. Eq. (7) is similar to the correlated wave function in quantum chemistry represented by a configuration interaction (CI) expansion.

Currently, there are two strategies to simulate the evolution of the quantum state on a classical computer. One is brute-force simulation which uses a ($2^{N}\times 1$) state vector (SV) or a ($2^{N} \times 2^{N}$) density matrix (DM) to represent the quantum state. The other is tensor network methods which approximate the quantum state by a set of low-rank tensors. The brute-force methods are supported by most of the existing packages (Table 2). Q²Chemistry also supports the state vector and density matrix methods. These brute-force methods implemented in the $\mathtt{q2chem.qcirc}$ module use the compressed sparse row (CSR) format for sparse matrix storage and OpenMP for multithreaded calculation. OpenACC-enabled C++ code is also implemented and can be selected at compile time to utilize GPU, which can provide $2-10$ times speed-up over the multithreaded CPU version.

Table  2.  Backends for circuit simulations in present open-source quantum computation softwares. External indicates that one or more third-party packages are required to enable the functionality. SA stands for single-amplitude simulation.

Software	Backend			Parallelism		GPU	Automatic differentiation
	SV	DM	TN	Threaded	Distributed
ProjectQ[28]	$\checkmark$	Ongoing	$\times$	$\checkmark$	SV	$\times$	$\times$
HiQ[50]	$\checkmark$	Ongoing	SA	$\checkmark$	SV	Ongoing	$\times$
MindQuantum[51]	$\checkmark$	$\checkmark$	$\times$	$\checkmark$	$\times$	$\checkmark$	$\checkmark$
Qulacs[30]	$\checkmark$	$\checkmark$	$\times$	$\checkmark$	$\times$	$\checkmark$	$\times$
Yao[31]	$\checkmark$	$\checkmark$	Ongoing	$\checkmark$	$\times$	$\checkmark$	$\checkmark$
Qiskit[29]	$\checkmark$	$\checkmark$	MPS	$\checkmark$	SV	$\checkmark$	$\times$
Cirq[52]	$\checkmark$	$\checkmark$	External	$\checkmark$	$\times$	$\checkmark$	$\times$
PaddleQuantum[53]	$\checkmark$	$\checkmark$	$\times$	$\checkmark$	$\times$	$\checkmark$	$\checkmark$
PennyLane[32]	$\checkmark$	$\checkmark$	MPS	$\checkmark$	$\times$	External	$\checkmark$
QuEST[54]	$\checkmark$	$\checkmark$	$\times$	$\checkmark$	SV, DM	$\checkmark$	$\times$
TEQUILA[55]	$\checkmark$	$\checkmark$	$\times$	$\checkmark$	$\times$	$\checkmark$	$\checkmark$
QFORTE[56]	$\checkmark$	$\times$	$\times$	$\checkmark$	$\times$	$\times$	$\times$
XACC[57]	$\checkmark$	$\checkmark$	External	$\checkmark$	$\checkmark$	External	$\times$
Q2Chemistry	$\checkmark$	$\checkmark$	MPS	$\checkmark$	$\checkmark$	$\checkmark$	$\checkmark$

 | Show Table

DownLoad: CSV

A typical drawback of brute-force simulation is that the memory usage and computational complexity both scale as $\mathcal{O}(c^{N})$, where $c$ is a constant between 2 and 4 and $N$ is the number of qubits. This exponential scaling prevents quantum simulations for larger molecules with more than ~30 qubits. Recently, matrix product states (MPS) and projected entangled pair states (PEPS) have been used to simulate large scale random quantum circuits^[58-61]. Since the tensor contraction pattern is mostly fixed and does not have NP-hard path optimization problems^[62], the one-dimensional MPS is preferable as a high-performance tensor network backend for a general quantum circuit simulator.

The MPS ansatz factorizes the rank-$N$ coefficients into lower rank tensors, which can be written as

where ${^k T^{i_k}_{u_{k-1} u_{k}}}$ is a rank-$3$ tensor with $i_k$ called the physical index and $u_k$ the auxiliary index. The maximum size of the auxiliary indices is defined as the bond dimension of the MPS, which is denoted as $D=\max_{0\le k \le N}{\{u_k\}}$. Algorithms based on MPS generally have a complexity of $\mathcal{O}(ND^3)$. If the bond dimension $D$ is allowed to grow exponentially, the MPS can exactly represent any quantum state using Eq. (8).

The $\mathtt{q2chem.qcirc}$ module implements the MPS simulation algorithm on top of the machine learning framework PyTorch^[40] based on the algorithm proposed by Vidal et al.^[63,64]. A demonstrative procedure for simulating quantum circuits using the MPS algorithm is presented in Fig. 2. Different from common classical MPS-based methods such as the density matrix renormalization group (DMRG)^[37], a set of auxiliary matrices are inserted between the rank-3 tensors and stored to maintain a normalized quantum state after truncated SVD.

Figure  2.  (a) Applying a single qubit gate on the MPS quantum state is simulated by simply a local contraction. (b) Applying a two-qubit gate on neighboring qubits generally has 2 steps: (i) reshape the two-qubit gate into a 4-dimensional tensor and contract with the qubits to form a two-qubit tensor; (ii) perform a singular value decomposition to restore the two-qubit tensor back to the MPS formulation. Postprocessing is usually required to maintain normalization or canonicalization of MPS tensors. (c) Auxiliary matrices that contain truncated and normalized singular values are used to normalize the quantum state.

DownLoad: Full-Size Img PowerPoint

Classical quantum simulator backends can also be extended within the framework of Q²Chemistry in a similar style as the interfaces for hardware. It requires only a minimum effort to implement functions _$\mathtt{quantum}$_$\mathtt{gate}$_$\mathtt{evolution()}$ and _$\mathtt{measure}$_$\mathtt{qubit()}$ for a derived class of _$\mathtt{BaseSimulator}$ and, if necessary, a custom data structure to store the quantum state. No modification to the upper level modules such as expectation value evaluation or higher level quantum algorithms is needed. In Section 6.2, we show the application of Q²Chemistry interfaced with an external MPS simulator QuantumSpins^[59,65], which efficiently simulated a 40-qubit molecule with a high accuracy.

4.3   Reversible automatic differentiation

In many quantum algorithms, a parametric quantum circuit is constructed, and optimization of the circuit parameters is carried out iteratively. On a quantum computer, the gradients of energy or another target function with respect to circuit parameters can be calculated through the parameter-shift rule or finite difference steps, which will introduce an additional complexity factor of $\mathcal{O}({N_{p}})$, where ${N_{p}}$ is the number of parameters. Nevertheless, using a classical quantum circuit simulator the gradients can be evaluated efficiently with an approximately $\mathcal{O}(1)$ complexity, which is extremely helpful just as the reverse-mode automatic differentiation (AD) algorithm used in classical machine learning if a large number of parameters are involved. Different from conventional machine learning, the computation graph of a quantum expectation value with respect to variational parameters is reversible:

which indicates that the intermediate quantum states $\{|\Psi\rangle_{n-1}, |\Psi\rangle_{n-2}, \cdots\}$ can be calculated iteratively and thus do not need to be stored explicitly in memory, enabling a significant reduction in memory consumption during the reverse-mode gradient evaluations^[66-68]. This classical algorithm, which is termed reversible AD, is first implemented in the Julia package Yao.jl, and a detailed description is given in Algorithm 1.

Algorithm 1: Reversible AD algorithm to calculate $g_{i}=\frac{\partial E}{\partial \theta_{i}}$.

Data:$\hat{H}$, $\{U_{0}(\theta_{0}), U_{0}(\theta_{0}), \ldots\}$, $|\Psi\rangle$

Result: $g$: gradients of energy w.r.t. parameters $\{\theta_{0}, \theta_{1}, \ldots\}$

$N_{p}\leftarrow$ number of parameters, $g\leftarrow$empty array of length $N_{p}$;

for i=0; i$\leqslant N_{p}-1$; i+=1 do

$\left| \; \; \; \; \; \right.|\Psi\rangle\leftarrow U_{i}(\theta_{i})|\Psi\rangle$;

end

$|\Psi_{l} \rangle \leftarrow \hat{H} | \Psi \rangle$, $|\Psi_{r} \rangle \leftarrow |\Psi\rangle$;

for i=$N_{p}-1$; i$\geqslant 0$; i-=1 do

$\left| \begin{aligned} & \;\;\;\;|{{\Psi }_{r}}\rangle \leftarrow U_{i}^{\dagger }({{\theta }_{i}})|{{\Psi }_{r}}\rangle; \\ & \;\;\;\;g\left[ i \right]=2\times \text{Re}( \langle {{\Psi }_{l}}|\frac{\partial {{U}_{i}}({{\theta }_{i}})}{\partial {{\theta }_{i}}}{{\Psi }_{r}}\rangle); \\ & \;\;\;\;|{{\Psi }_{l}}\rangle \leftarrow U_{i}^{\dagger }({{\theta }_{i}})|{{\Psi }_{l}}\rangle; \\ \end{aligned} \right.$

end

Q²Chemistry implements reversible AD for brute-force backends. It should be noted that the reversible AD algorithm is invalid on real quantum devices even if the original equation is used:

since the derivative gates $\left\{ {\dfrac{\partial U_{i}(\theta_{i})}{\partial \theta_{i}} } \right\}$ are generally nonunitary. At the same time, using reversible AD with the MPS simulator requires careful modifications to the algorithm and specially designed techniques due to SVD truncation during the simulation. A naive implementation of Algorithm 1 probably leads to suboptimal performance and numerical errors.

5.   Quantum algorithms

In the current release, Q²Chemistry provides a couple of VQE-based algorithms. The module $\mathtt{q2chem.qalgo}$ adopted several variational wave function ansatzes that can be used individually or collectively to solve for the eigenstates of the given chemical system. New quantum algorithms can also be conveniently implemented by users.

5.1   Variational quantum circuit ansatz

Introducing a parametric wave function $| \Psi(\theta) \rangle$, the lowest eigenvalue $E_{0}$ can be obtained variationally:

Such a protocol is implemented in the $\mathtt{q2chem.qalgo}$ module to obtain the eigenstates of a given Hamiltonian. Properly constructing a parametric quantum circuit, the wave function ansatz $| \Psi(\theta) \rangle$ is encoded into the quantum state of qubits. Combining the measurements for expectation value evaluation on a quantum computer and a numerical optimization algorithm on a classical computer, the variational procedure can then be performed in a hybrid quantum-classical way.

Generally, two broad types of variational quantum circuit ansatzes exist^[8]. They are physically motivated ansatz (PMA), which is inspired by classical wave function methods that systematically approach the exact electronic wave function, and hardware heuristic ansatz (HHA), which considers specific hardware structures and employs entangling blocks. Both types are currently integrated into the $\mathtt{q2chem.qalgo}$ module.

Unitary coupled-cluster is one of the most commonly used PMAs for quantum computing. Generally, the UCC wave function is defined as

In Q²Chemistry, the spin-adapt CCD0 cluster operators^[69] are used to construct the UCCSD and UCCGSD wave functions. A fermion-to-qubit mapping is performed, and first-order Trotter-Suzuki decomposition^[70,71] is implemented to convert Eq. (11) into the product of Pauli strings:

where $\theta_{i}$ is the variational parameter corresponding to the $i$-th anti-Hermitian fermion excitation operator $\hat{T}_i - \hat{T}^{\dagger}_i$, which is transformed into qubit form $\sum\limits_{j} {\text{i} c_{ij} \hat{P}_{ij}}$. Each exponential term is mapped to the quantum circuit following Algorithm 2. An example of the mapped circuit is presented in Fig. 3a.

Figure  3.  (a) The quantum circuit corresponding to the operator $\exp{({\rm{i}}\theta \hat{\sigma^{x}} \hat{\sigma^{y}} \hat{\sigma^{z}} \hat{\sigma^{x}})}$ and (b) the two-layer HEA circuit that entangles all neighboring qubits using the controlled-U gate. Blue squares represent nonparametric gates, while green squares represent parametric quantum gates such as Rz and the three-parameter (controlled-)U gate.

DownLoad: Full-Size Img PowerPoint

The hardware efficient ansatz (HEA) circuit^[14] is implemented as an HHA ansatz. If the type of entanglement gates, the ordering of entanglement qubits and the number of layers are set, Q²Chemistry automatically generates the parametric hardware efficient circuit for simulation. Fig. 3b shows an example of a two-layer HEA circuit with all neighboring qubits entangled by the three-parameter controlled-U gate.

Algorithm 2: Map $\exp{({\rm{i}}\theta \hat{P})}$ to circuit. ${\rm{HY}}$ is the Hadamard-Y gate defined as ${\rm{HY}} = {\sqrt{2}}/{2}\times({\rm{Z}} + {\rm{Y}})$

Data: $\hat{P}$, $\theta$

Result: C: the quantum circuit

$N_{q} \leftarrow$ number of qubits, ${\rm{C}}\leftarrow$ empty circuit;

for i=0; i$\leqslant N_{q}-1$; i+=1 do

$\left| \begin{align} & \ \ \ \ {{p}_{i}}=P\left[ i \right]; \\ & \ \ \ \ \mathbf{if}\ \ {{p}_{i}}=={{{\hat{\sigma }}}_{x}}\ \ \bf{then} \\ & \ \ \ \ \left| \ \ \text{C}+={{\text{H}}_{i}} \right. \\ & \ \ \ \ \bf{else}\ \bf{if}\ {{p}_{i}}=={{{\hat{\sigma }}}_{y}}\ \ \bf{then} \\ & \ \ \ \ \left| \ \ \text{C}+=\text{H}{{\text{Y}}_{i}} \right. \\ & \ \ \ \ \bf{end} \end{align} \right.$

end

for i=$N_{q}-2$; i$\geqslant 0$; i-=1 do

|　　C += ${\rm{CNOT}}_{(i+1, i)}$

end

C += ${\rm{RZ}}(-2\theta)_{N_{q} - 1}$

for i=0; i$\leqslant N_{q}-2$; i+=1 do

|　　C += ${\rm{CNOT}}_{i+1, i}$

end

for i=0; i$\leqslant N_{q}-1$; i+=1 do

$\left| \begin{align} & \ \ \ \ {{p}_{i}}=P\left[ i \right]; \\ & \ \ \ \ \mathbf{if}\ \ {{p}_{i}}=={{{\hat{\sigma }}}_{x}}\ \ \bf{then} \\ & \ \ \ \ \left| \ \ \text{C}+={{\text{H}}_{i}} \right. \\ & \ \ \ \ \mathbf{else}\ \mathbf{if}\ {{p}_{i}}=={{{\hat{\sigma }}}_{y}}\ \ \mathbf{then} \\ & \ \ \ \ \left| \ \ \text{C}+=\text{H}{{\text{Y}}_{i}} \right. \\ & \ \ \ \ \mathbf{end} \\ \end{align} \right.$

end

Several techniques are implemented in Q²Chemistry to reduce the computational overhead, including qubit tapering for Hamiltonian^[72], qubit excitation based (QEB) operator^[73], Pauli entangler^[74,75], symmetry-based operator selection^[33] and the ADAPT-VQE algorithm^[76]. Q²Chemistry supports combination of above methods, for example, using Pauli entanglers that entangle at most 4 qubits together with ADAPT-VQE leads to an iterative-qubit-coupled-cluster (iQCC)^[75] like algorithm.

In addition to ground state methods, the $\mathtt{q2chem.qalgo}$ module offers a couple of post-VQE algorithms for calculating excited states, including the variational quantum deflation (VQD)^[18], quantum subspace expansion (QSE)^[19] and equation-of-motion (EOM)^[45,77,78] theory. VQD consecutively constructs an effective Hamiltonian, the lowest eigenvalue of which corresponds to the 1st, 2nd, 3rd, $\ldots$ excited state energy:

where $\hat{H}^{0}_{{\rm{eff}}}$ is the original second-quantized electronic Hamiltonian. QSE and EOM use a set of fermion excitation operators or a state-transfer operator to construct and solve a generalized eigenvalue problem ${\boldsymbol{MC}} = {\boldsymbol{SCE}}$ by additional quantum measurements on the top of the ground state circuit. These methods are implemented in Q²Chemistry for calculating electron excitations, ionization potentials and electron affinity energies for both molecules and periodic systems.

5.2   Parallel evaluation of an expectation value

As shown in Eq. (2), the Hamiltonian is expressed as the summation of a polynomial number of mutually uncorrelated Pauli strings. The expectation values of each Pauli string can thus be calculated independently. Fig. 4 gives an example of circuits used for evaluating the expectation value $\sum\limits_{i} { c_{i} \langle \Psi | {\hat{P}_{i}} | \Psi \rangle}$ for a given Hamiltonian containing a number of Pauli strings under some fermion-to-qubit transformation. For each of the circuits, the ansatz parts are the same while the measurement parts are constructed according to the specific form of the Pauli string. During the calculation of the expectation value, Q²Chemistry automatically distributes these circuits to different quantum devices or simulator processes. On each device, a subset of circuits are executed and then measured. The measurement outcomes are then postprocessed to calculate expectation values locally. Finally, the results are reduced across all the devices to obtain the total energy.

Figure  4.  Evaluating the expectation value of a linear combination of Pauli strings using multiple quantum devices or simulator processes. $|0_{N}\rangle$ represents an $N$-qubit quantum register with all qubits initialized to $|0\rangle$. In this example, the measurement parts are the Hadamard test circuits.

DownLoad: Full-Size Img PowerPoint

It should be noted that on real quantum devices, since the quantum states are nonreplicable, all the circuits should be executed. However, on a classical simulator, the memory data of quantum states can be reused. Therefore, using a quantum circuit simulator, the ansatz part only needs to be executed once, and the simulated quantum state can then be copied to each process for later measurements. This strategy is used in the simulations of Section 6.

6.   Applications

We present several simulations to show the power of Q²Chemistry, including a scalability test for the MPS-based quantum circuit simulator and the numerical simulation results for ground- and excited-state calculations.

6.1   Scalability benchmark

The parallel measurement for the expectation value introduced in Section 5.2 is extended to a two-level parallelism for the MPS simulator:

(I) The first level parallelization over Hamiltonian. Subsets of Pauli strings from the Hamiltonian in Eq. (1) are distributed to each process. The expectation values of Pauli strings are calculated independently, and a reduce-sum is performed across all processes to obtain the final energy.

(II) The second low-level parallelization using multithread parallelism on CPU or GPU, to accelerate the calculations of linear algebra routines such as matrix multiplication and singular value decomposition.

The two-level strategy enables good parallel scalability if the adapted dynamical distribution algorithm is used to achieve load balance. To reduce memory usage, in the first level, the quantum circuit is stored and evolved only on the $0$th process. Although for small molecules such as H$_2$ (4 qubits), the circuit evolution may contribute over 80% to the total execution time, the number of Pauli strings in the Hamiltonian will quickly go beyond $10^{5}$ due to the $\mathcal{O}(N^{4})$ scaling, and the time cost of circuit simulation will become negligible if larger systems with more than 12 qubits are involved. The scaling benchmark of Q²Chemistry simulating a Cr$_2$ molecule (STO-3G basis set) using the MPS backend is given in Fig. 5. Q²Chemistry achieves good parallel performance up to 768 CPU cores for this 72-qubit system. Note that the main purpose of this test is for parallelization scaling instead of obtaining reliable results. Therefore, only 51 of the 3131 variational parameters are considered, leading to a quantum circuit with 119884 gates. The upper bound for the bond dimension is also set to a relatively low value of 64.

Figure  5.  Simulating the Cr$_2$ molecule using the MPS backend. The STO-3G basis set and the symmetry-reduced UCCSD ansatz are used, where the qubit Hamiltonian contains 305041 Pauli strings. The time cost refers to one VQE iteration (including evolution of the quantum circuit and calculation of energy) being tested. The distributed parallelization is implemented using OpenMPI and Python’s $\mathtt{mpi4py}$ package.

DownLoad: Full-Size Img PowerPoint

6.2   Numerical results

Fig. 6a shows the potential energy curve of the H$_2$ molecule. The calculation is carried out by extending the simulators in the $\mathtt{q2chem.qcirc}$ module to an external MPS circuit simulator written in Julia^[59,65]. Benefiting from the high-level modular structure introduced in Section 2, this external simulator and the parallelization techniques introduced in Section 5.2 can be easily implemented within the framework of Q²Chemistry by adding a few lines of code. The ground-state energies are variationally optimized using the ccJ-pVDZ basis set^[79] and the symmetry-reduced UCCSD^[33] circuit (leading to 40 qubits and 53 variational parameters). The BOBYQA optimizer is used to perform gradient-free optimizations. For each geometry, 2000 optimization steps are performed within 24 h using 560 CPU cores.

Figure  6.  (a) The VQE optimized potential energy curve of H$_2$ calculated by the MPS backend using the ccJ-pVDZ basis set. The VQE results are calculated by interfacing with the external Julia-implemented MPS simulator. (b) Energy difference between EOM-ADAPT-C and classical EOM-CCSD band structures for Si. Inset gives the EOM-ADAPT-C band structure calculated using Q²Chemistry. The FCI and EOM-CCSD energies are obtained using the PySCF code.

DownLoad: Full-Size Img PowerPoint

Figure  .  Code Example 1. Interface to quantum devices in the qcircuit module.

DownLoad: Full-Size Img PowerPoint

Fig. 6b calculates the quasi-particle band structures for silicon. For such periodic systems, the Hamiltonian and UCC-based wave function ansatz need to include the constraint of crystal momentum conservation^[44], which is automatically handled by the $\mathtt{q2chem.qchem}$ module. The simulation uses the GTH-SVZ basis set with a GTH-PADE pseudopotential. A UCCGSD operator pool with complementary operators ^[45] is employed together with the ADAPT algorithm for the ground state ADAPT-C wave function. With a $1 \times 1 \times 1$ k-point grid, 16 qubits are simulated using the state vector backend. The EOM-ADAPT-C method^[45] is used for band structure calculation, which achieves a mean absolute difference of 0.047 eV from EOM-CCSD and the deviation is as small as $\sim 10^{-3}$ eV at the $\Gamma$ point (Fig. 6b).

7.   Conclusions

In this study, we demonstrate that the Q²Chemistry package is suitable for simulating and developing quantum algorithms for quantum chemistry applications. Q²Chemistry provides versatile functionalities for simulating ground- and excited-state properties. The simulator backend, including the parallelized MPS algorithm, achieves high performance for large-scale simulations using a moderate amount of computational resources. Directions for future development include more classical simulation methods, including high-dimensional tensor network based methods, more integrated quantum algorithms, and high-performance quantum circuit optimization algorithms. With the flexibility to link to different quantum devices, Q²Chemistry can be used as a useful platform in pursuing a practical quantum advantage.

Supporting information

The supporting information for this article can be found online at https://doi.org/10.52396/JUSTC-2022-0118. The supporting information includes one figure. The way to count CNOT gates is described in the supporting information.

Acknowledgements

This work was partially supported by the National Natural Science Foundation of China (21825302), the Fundamental Research Funds for the Central Universities (WK2060000018), the National Supercomputing Center in Jinan, and the USTC Supercomputing Center.

Conflict of interest

The authors declare that they have no conflict of interest.

Preprint statement

Research presented in this article was posted on a preprint server prior to publication in JUSTC. The corresponding preprint article can be found here: arXiv: 2208.10978. https://arxiv.org/abs/2208.10978.