Large-Scale Simulation of Quantum Computational Chemistry on a New Sunway Supercomputer

VQE approximates the exact eigenstate with a unitary transformation acting on a reference state. The reference state is often chosen to be some state that can be easily prepared on a quantum computer and the unitary transformation is represented by a parametric quantum circuit, with the parameters updated iteratively based on the output energy of each minimization step. A classical simulator for VQE simulates this process on a classical computer. The brute-force approach to simulate a quantum circuit is the state vector simulator, which exactly stores the underlying $n$-qubit quantum state. The memory cost of the state vector simulator scales as $2^{n}$ for pure quantum states (implemented in ProjectQ , [22] Qiskit, [23] QuEST [24] and Qulacs [25]) and $2^{2n}$ for mixed quantum states (implemented in DM-Sim , Qiskit, [23] QuEST [24] and Qulacs [25]). The largest quantum circuit simulation with this approach has reached $45$ qubits, but it would be extremely difficult to go further beyond due to the exponential memory cost [26].

Tensor network based algorithms are important alternatives that reduce the memory requirement by representing the quantum state as a tensor network of low-rank tensors, which have been widely applied to simulate RQCs [12, 13, 14, 15, 16, 17, 18]. The tensor network methods used for RQCs are mostly based on high-dimensional tensor network contractions since the RQCs under study usually have well defined two-dimensional geometries. However, the parametric quantum circuits used in VQE would not have such clear geometrical features in general, especially for those physical motivated ansatz such as unitary coupled-cluster as used in this work. In such instances the MPS (a very special type of one-dimensional tensor network that can efficiently simulate weakly entangled quantum states) based simulator is often preferable, especially for quantum chemistry problems. We also note that MPS has already been widely applied to solve quantum chemistry problems by directly computing the ground states in the context of Density Matrix Renormalization Group (DMRG) algorithm [27, 28]. MPS has also been used to directly simulate Shor’s algorithm up to $60$ qubits [29]. MPS-based simulator has already been implemented, for example in QISKIT. However, large-scale implementation of it for VQE on a top supercomputer has not been demonstrated.

The electronic Hamiltonian of a chemical system can be written in a second-quantized formulation:

$\displaystyle\hat{H}$

$\displaystyle=\sum_{p,q}{h_{q}^{p}a_{p}^{\dagger}a_{q}}+\frac{1}{2}{\sum_{\begin{subarray}{c}p,q\\ r,s\end{subarray}}}{g_{rs}^{pq}a_{p}^{\dagger}a_{q}^{\dagger}a_{r}a_{s}}$

(1)

where $h_{q}^{p}$ and $g_{rs}^{pq}$ are one- and two-electron integrals in molecular orbital basis. In a VQE procedure, the creation and annihilation operators $\{{a_{p}^{\dagger}},a_{q}\}$ in the Hamiltonian should be mapped to weighted Pauli strings using the Jordan-Wigner or Bravyi-Kitaev transformation, therefore the energy can be obtained by measuring and summing the expectation values of Pauli strings $\{P\}\in\{I,\sigma_{x},\sigma_{y},\sigma_{z}\}^{\otimes n}$ as

	$\displaystyle E$	$\displaystyle=\langle{\hat{H}}\rangle=\langle\Psi(\theta)|\hat{H}|\Psi(\theta)\rangle$		(2)
		$\displaystyle=\langle\Psi(\theta)|\sum_{p,q}{\tilde{h}_{q}^{p}P_{q}^{p}}|\Psi(\theta)\rangle$
		$\displaystyle+\langle\Psi(\theta)|\sum_{p,q,r,s}{\tilde{g}_{rs}^{pq}P_{rs}^{pq}}|\Psi(\theta)\rangle.$

The accuracy of VQE heavily depends on the wave function ansatz (or the structure of the parametric quantum circuit) used to approximate the unknown ground state of the problem Hamiltonian. The mostly used form of the wave function ansatz is the physically motivated unitary coupled-cluster (UCC)[19], where the wave function takes the form

$$|\Psi(\theta)\rangle=e^{\hat{T}(\theta)-\hat{T}^{\dagger}(\theta)}|\Phi_{0}\rangle.$$

(3)

Here $|\Phi_{0}\rangle$ is an initial state in the form of a single Slater determinant. When truncated to the single and double excitations (UCCSD), the cluster operators are

$$\hat{T}(\theta)=\sum_{a,i}\theta_{i}^{a}\hat{a}^{\dagger}_{a}\hat{a}_{i}+\frac{1}{4}{\sum_{a,b,i,j}}\theta_{ij}^{ab}\hat{a}_{a}^{\dagger}\hat{a}_{b}^{\dagger}\hat{a}_{i}\hat{a}_{j},$$

(4)

where $\{i,j,k,\dots\}$, $\{a,b,c,\dots\}$ and $\{p,q,r,\dots\}$ indicate the occupied, virtual, and general orbitals, respectively. The UCCSD method has no exact finite truncation for the Baker-Campbell-Hausdorff expansion, which indicates that it can not be effectively implemented on a classical computer using polynomial expansion methods or the explicit matrix-vector method. UCCSD can be implemented as a parametric quantum circuit by Suzuki-Trotter decomposition of the unitary evolution operator into one- and two-qubit gates [30].

Many contributions have been devoted to reducing the complexity of the quantum circuit ansatz for quantum chemistry problems. Recently, the excited states of Ir(ppy)₃ have been computed with up to $56$ qubits [31]. This is the largest classical simulation of VQE for quantum computational chemistry in terms of the number of qubits up to date. However, to achieve this simulation a very shallow parametric quantum circuit ansatz (with around $60$ two-qubit gates) was employed to reduce the computational complexity.

Following the philosophy that ”divide each of the difficulties under examination into as many parts as possible, and as might be necessary for its adequate solution”, a general VQE framework combined with the divide-and-conquer method (DMET) [32] for large-scale simulation of quantum chemistry has been proposed [33]. Using this approach, the C₁₈ molecule with $144$ qubits is studied classically, where the VQE in each fragment uses at most $16$ qubits [34]. A similar idea has been adopted in a hybrid tensor network approach, in which a bottom-layer tensor network is introduced to describe the strongly correlated intrasubsystem interactions and a top-layer tensor network is used to describe the intersubsystem interactions [35]. This embedding method can also be applied to excited state calculation [36, 37].

Overall, the current state-of-the-art quantum computational chemistry simulators for large-scale systems are limited with respect to computational efficiency and scalability. Thus, the classical simulation of large-scale quantum chemistry systems are severely limited. This represents the key bottleneck in simulating the physical properties of real systems.

The goal we want to achieve is to faithfully simulate the large-scale VQE experiments for practical Chemical systems running on quantum computers, using a world-leading supercomputer. Since we are aiming for a parallelization scale of up to tens of millions of cores, the algorithm-wise challenge is to identify a suitable algorithm which could map the computation onto all the available cores of the underlying supercomputer, thus extending the simulation range of VQE to quantum chemistry systems which are of practical relevance.

The state vector simulator is out of choice due to the huge memory requirement. We thus take the tensor network based simulator which reduces the exponential memory growth. In practice we find that MPS based simulator gives the best performance when used in combination with the UCC ansatz. To spread the computations of the MPS-VQE simulator over the new Sunway supercomputer, there are three implementation-wise challenges:

• •

  To parallelize over the $20$ million cores, we need a scheme which can first efficiently decompose the VQE simulation into sub-tasks (independent quantum circuit evolution and measurements) with little communication overhead and of similar complexity, and then distribute the tensor operations involved in a single circuit evolution in each process.

• •

  To achieve maximal efficiency for the major tensor operations, which are SVD and tensor contraction, we need a well tuned map of these tensor operations to the underlying CPE mesh.

• •

  Given the highly complex nature of the program including a hybridization of DMET, MPS and VQE as well as the application for practical problems, a variety of scientific libraries in high-level programming languages such as Python needs to be used to avoid reinventing the wheels. To build such a complex application with the simplest coding and also maintain a computational efficiency close to the low level C language, we pay a huge effort to develop the Julia compilation framework with SACA runtime support on heterogeneous processors to achieve functional implementation and performance improvement.

The new-generation Sunway supercomputer is used for performance assessment in this work, which is the successor of the Sunway TaihuLight supercomputer. Similar to the Sunway TaihuLight system, the new Sunway supercomputer adopts a new generation of domestic high-performance heterogeneous many-core processors (SW26010Pro) and interconnection network chips in China.

The new SW processor is designed for massive thread and data parallelism and to deliver high performance on parallel workloads. The architecture of the SW26010Pro processor is shown in Fig. 1. Each processor contains 6 core-groups (CGs), with 65 cores in each CG, and in total 390 cores. Each CG contains one management processing element (MPE), one cluster of computing processing elements (CPEs) and one memory controller (MC). The MPE within each CG is used for computations, management and communication. The CPEs are organized as an $8\times 8$ mesh (64 cores) and are designed to maximize the aggregated computing power and to minimize the complexity of the micro-architecture. The CPEs are organized with a mesh network to achieve high-bandwidth data communication (P2P and collective communications) among the CPEs in one CG, which is referred to as the remote scratchpad memory access (RMA).

Each SW26010Pro processor contains 96 GB memory, with 16 GB memory in each CG. The MPE and the CPEs within the same CG share the same memory which is controlled by the MC. Each CPE has a 32 KB L1 instruction cache, and a 256 KB scratch pad memory (SPM, also called the Local Data Memory (LDM)), which serves the same functionality as the L1 cache. The data storage space can also be configured as a local data cache (LDCache) which is automatically managed by the hardware. Data transfer between LDM and main memory can be realized by direct memory access (DMA), and data transfer between LDCache and main memory can also be realized by conventional load/storage instructions.

In quantum chemistry, the FCI generates a correlated wave function by considering all the excitations above the reference Hartree-Fock state. The FCI wave function $\left|\Psi\right\rangle$ can be written as:

$$\left|\Psi\right\rangle=\sum_{i_{1}...i_{N}}c_{i_{1}i_{2}i_{3}\ldots i_{N}}|i_{1}i_{2}i_{3}\ldots i_{N}\rangle$$

(5)

where $|i_{1}i_{2}i_{3}\ldots i_{N}\rangle$ is the computation basis used, the coefficient $c_{i_{1}i_{2}i_{3}\ldots i_{N}}$ is a rank-$N$ tensor of $2^{N}$ complex numbers. Finding the tensor $c$ corresponding to the ground state is in generally exponentially hard. Various simplified ansatz has been used instead of the exact form in Eq.(5) to reduce the exponential scaling, among which the tensor network states (TNSs) constitute an outstanding class [20, 38]. TNS represents the underlying high-rank tensor as the product of a tensor network made of low-rank tensors, and the number of parameters in a TNS typically scales linearly with the system size (number of qubits). MPS is a special class of one-dimensional TNS which not only allows efficient representation, but also enables efficiently manipulations, such as the local gate operations and evaluating local observables, on it. In the context of random circuit sampling, it has been shown that the two-dimensional variant of TNS, namely the Projected Entangled Pair States (PEPS), would have better performance than MPS since current random circuits are mostly two-dimensional [12]. Here we will use the MPS ansatz for the quantum state due to the complex geometry of our parametric circuit ansatz, as shown in Fig. 2(a), which can be written as

$$c_{i_{1}i_{2}\ldots i_{N}}=\sum_{a_{1},\dots,a_{N-1}}B^{i_{1}}_{a_{1}}B^{i_{2}}_{a_{1},a_{2}}\dots B^{i_{N}}_{a_{N-1}},$$

(6)

with $i_{n}$ the ‘physical’ index and $a_{n}$ the ‘auxiliary’ index. The largest size of the auxiliary indices is referred to as the bond dimension of MPS, denoted as $D=\dim_{1\leq n<N}(a_{n})$. Eq.(6) can represent arbitrary quantum states if $D$ is allowed to grow exponentially. A quantum state is said to be efficiently representable as an MPS if $D$ remains almost constant when $N$ grows. An MPS based algorithm in general has a complexity of $O(D^{3})$. In practice, it is often advantageous to keep MPS in a canonical form, for example, the right-canonical form by ensuring that $\sum_{i_{n},a_{n}}(B^{i_{n}}_{a_{n-1},a_{n}})^{\ast}B^{i_{n}}_{a_{n-1}^{\prime},a_{n}}=\delta_{a_{n-1},a_{n-1}^{\prime}}$ (for one thing, the algorithm would be numerically more stable since the tensors on each site are well normalized).

VQE requires to apply single-qubit and two-qubit gate operations sequentially onto the quantum state and then compute the expectation value of a (local) observable on the quantum state, which is similar to the time-evolving block decimation algorithm used for time evolution of one-dimensional Hamiltonians [39]. Assuming that the initial MPS is prepared in the right-canonical form, a nearest-neighbour two-qubit gate operation, denoted as $O^{i_{n},i_{n+1}}_{i_{n}^{\prime},i_{n+1}^{\prime}}$, can be applied onto a nearest-pair of tensors on sites $n$ and $n+1$ (also referred to as the $n$-th bond) as follows [40] (The implementation of single-qubit gate operation is straightforward, but it is not necessary since single-qubits gate can be absorbed into two-qubits gate using gate fusion). As shown in Fig. 2(b), first, one obtains the rank-$4$ tensor $M$ by the contraction of the two-qubit gate with the tensors $B^{i_{n}}_{a_{n-1},a_{n}}$ and $B^{i_{n+1}}_{a_{n},a_{n+1}}$

$$M_{a_{n-1},i_{n},i_{n+1},a_{n+1}}=\sum_{a_{n}}O^{i_{n},i_{n+1}}_{i_{n}^{\prime},_{n+1}^{\prime}}B^{i_{n}^{\prime}}_{a_{n-1},a_{n}}B^{i_{n+1}^{\prime}}_{a_{n},a_{n+1}}.$$

(7)

Now to restore the MPS form, one could perform a singular value decomposition (SVD) on $M$ to obtain two rank-$3$ site tensors, however one of them will no longer be right-canonical. To preserve the right canonical form, one could perform SVD on another rank-$4$ tensor

$$M^{\prime}_{a_{n-1},i_{n},i_{n+1},a_{n+1}}=\sum_{a_{n-1}}\lambda_{a_{n-1}}M_{a_{n-1},i_{n},i_{n+1},a_{n+1}}$$

(8)

instead, where $\lambda_{a_{n-1}}$ are the Schmidt numbers on the $n-1$-th bond saved from the previous steps, that is,

$${\rm SVD}(M^{\prime}_{a_{n-1},i_{n},i_{n+1},a_{n+1}})=\sum_{a_{n}}U^{i_{n}}_{a_{n-1},a_{n}}\lambda_{a_{n}}^{\prime}V^{i_{n+1}}_{a_{n},a_{n+1}},$$

(9)

where $\lambda_{a_{n}}^{\prime}$ are the new Schmidt numbers on the $n$-th bond which will be used to replace the previous ones. In case the number of nonzero Schmidt numbers is larger than $D$, we will also truncate them back to $D$ by only reserving the $D$ largest Schmidt numbers. This procedure will produce truncation errors if the maximum $D$ we choose is too small or the circuit is too deep (Nevertheless this error can be monitored as an indication of the imprecision of the computing). $V^{i_{n+1}}_{a_{n},a_{n+1}}$ is already right-canonical and will be used to replace $B^{i_{n+1}}_{a_{n},a_{n+1}}$, the new site tensor on site $n$ can be updated as

$$B^{i_{n}}_{a_{n-1},a_{n}}=\sum_{i_{n+1},a_{n+1}}M_{a_{n-1},i_{n},i_{n+1},a_{n+1}}(V^{i_{n+1}}_{a_{n},a_{n+1}})^{\ast}$$

(10)

which can be verified to also be right-canonical. With a right-canonical form of MPS, the expectation value of a local observable $\hat{O}_{n}$ (with its matrix representation $O_{i_{n}^{\prime}}^{i_{n}}$ in the computational basis) on the MPS can be evaluated as

$$\langle\Psi|\hat{O}|\Psi\rangle=\sum_{a_{n-1},a_{n},i_{n},i_{n}^{\prime}}\lambda_{a_{n-1}}^{2}O_{i_{n}^{\prime}}^{i_{n}}B^{i_{n}^{\prime}}_{a_{n-1},a_{n}}(B^{i_{n}}_{a_{n-1},a_{n}})^{\ast}.$$

(11)

The expectation value of two-qubit or higher observables can be efficiently computed accordingly. In this work, the expectation values of Pauli strings are calculated using Hadamard test, which is equivalent to evaluating the expectation value of the local projector $|0\rangle\langle 0|$ or $|1\rangle\langle 1|$ on the ancillary qubit.

To this end, we note that since our VQE is based on the MPS representation, the expressiveness of the VQE ansatz used in our classical simulator will not exceed the expressiveness of the underlying MPS. As a result one may well substitute the VQE simulator by another MPS based optimization algorithm such as DMRG and a similar or even higher precision would be expected if the same $D$ is used. However, the MPS based VQE simulator, as will be shown later, allows straightforward parallelization to millions of cores with almost no data communications, in comparison in current parallel algorithms for DMRG data communication can not be avoided and the parallelizability is limited [41].

The MPS based simulator takes advantage of the low entanglement nature of the underlying quantum state to dramatically reduce the required computational resources, which is often true for NISQ devices for which decoherence limits the amount of entanglement that can be generated (except for RQC which quickly entangles the quantum states but usually has no clear physical meaning). An illustrative comparison between the performances of state vector (SV), density matrix (DM) and MPS simulators is shown in Fig.2(c), with a quantum circuit that generates a quantum state which can be written as an MPS with $D=8$. This circuit implements unitaries each of which entangles 4 consecutive qubits, corresponding to the correlations between neighbouring orbitals in a chemical system. It is also used as the MPS-inspired ansatz[42, 43] for slightly correlated molecules. For the SV or DM method, the exponentially growing memory storage and computation complexity strictly limited the simulation scale, while the MPS method approximates the low-entangled quantum states with smaller tensors, that only spans a very tiny fraction of the overall Hilbert space, so the required computational resources is dramatically reduced. However, faithfully simulating of the whole VQE calculation is still numerically demanding since a huge number of quantum circuits ($\mathcal{O}(N_{q}^{4})$ for $N_{q}$) has to be evaluated. Nevertheless, this is an ideal situation for massive parallelization, which will be discussed in Sec. III-C.

To facilitate the quantum simulation of the large systems, we use density matrix embedding theory (DMET) to reduce large systems into smaller fragments and use MPS-VQE as a fragment solver, as shown in Fig. 3. The total energy of the system is the sum of the energies of the fragments, each of which is embedded in an effective environment (quantum bath) compressed from the rest of the whole system. DMET provides an economic and highly accurate way to describe the interaction between the fragment and the effective environment. In this way, the electronic structure of the system can be computed by iteratively solving a bunch of smaller problems, which can be done using the state vector or MPS simulators (or ultimately using a quantum computer). In practice, the standard procedures of DMET-MPS-VQE method are listed as follows:

1. 1.

   Perform low-level calculation for the whole system.

2. 2.

   Divide the whole system into fragments with reasonable number of orbitals for each fragment.

3. 3.

   For every fragment, construct bath orbitals and the reduced Hamiltonian by Schmidt decomposition and projection.

4. 4.

   Get the fragment energy and the wave function with VQE, and then compute the number of electrons with 1-RDM (reduced density matrix).

5. 5.

   Check if the sum of the number of electrons in the fragments agrees with the number of electrons for the entire system. If not, modify the global chemical potential and go back the step 3.

In DMET, a high-level calculation for each fragment is carried out individually until the self-consistency criterion has been met: the sum of the number of electrons of all of the fragments (the number of fragments are around 10 to 200) agrees with the number of electrons for the entire system.

Three levels of parallelization are used for the DMET-MPS-VQE simulator. As shown in Fig. 4, the calculation of different fragments can be performed in an embarrassingly parallel manner, that is, the first level of parallelization. As no communication between different fragment calculations is required, we split the whole CPU pool into different sub-groups and sub-communicators. Within each sub-group, the total energy is calculated with the VQE method. Within each VQE calculation, another two levels of parallelization have been used, that is, the second level over circuits and the third level over linear algebra operations. In the VQE calculation, the electronic molecular Hamiltonian $H$ can be written as a sum of the Pauli strings: $H=\sum_{i}{a_{i}P_{i}}$, where each Pauli string $P_{i}\in\{I,\sigma_{x},\sigma_{y},\sigma_{z}\}^{\otimes n}$ is a tensor product of $n$-qubit Pauli operators. We can estimate each $\langle\psi|P_{i}|\psi\rangle$ independently and then sum over all those results with corresponding coefficients ($a_{i}$) to get $\langle\psi|H|\psi\rangle$. This enables high parallel scalability with adapted dynamical load balancing algorithm. The parallel simulation algorithm based on distributed memory over the circuits, just “mimic” the actual quantum computers (for quantum computers each term $\langle\psi|P_{i}|\psi\rangle$ also has to be evaluated individually), so our method can offer a good reference for VQE running on the quantum computers. For tensor operations including tensor contraction and singular value decomposition (SVD), a low-level multi-threaded parallelism on the CPE is used to further boost the performance.

The second level of parallelism is based on the fact that the electronic structure Hamiltonian can be expressed as the summation of a polynomial number of mutually uncorrelated Pauli strings. Expectation values of each Pauli string can thus be calculated independently. Fig.5 gives an example for the hydrogen molecule. The $4$-qubit Hamiltonian of hydrogen molecule contains $15$ Pauli strings, therefore $15$ circuits are required for each of them. For each circuit, the first $2$ $X$ gates prepare a Hartree-Fock reference state. The next $120$ gates (with variational parameters in the $RZ$ gates) on q₀, q₁, q₂ and q₃ constitute the ansatz circuit, which corresponds to the first-order Trotterized UCC ansatz and is identical for each Pauli string. $q_{4}$ is the ancillary qubit for Hadamard test used to compute the expectation value, and the gates on $q_{4}$ in different circuits corresponding to the measurement part are constructed according to each Pauli string. As shown in Fig.4, during the VQE optimization, circuits corresponding to mutually exclusive subsets of Pauli strings in Eq.2 are mapped to different processes. The energy $E(\theta)={\langle\psi(\theta)|\hat{H}|\psi(\theta)\rangle}$ under current parameters $\{\theta\}$ is obtained by reducing the results on different processes and then provided to the optimizer.

It should be noted that the number of Pauli strings in the electronic Hamiltonian scales as $\mathcal{O}(N_{q}^{4})$, where $N_{q}$ is the number of qubits. Even for a medium-size molecule such as benzene (72 qubits, 330816 Pauli strings), storing all the $330816$ circuits on each process will bring a lot of pressure on the memory space of CGs and synchronizing the circuits after each optimization step will also result in remarkable overhead on communication. Noticing that the ansatz part of different circuits are identical and the measurement part can be constructed on-the-fly in the first energy evaluation and then kept constant during the VQE optimization, a straightforward optimization for memory usage can be achieved by keeping only one replica of the ansatz circuit on each process. The speedup of such a memory-optimized simulation is given in Sec. IV-B.

Julia is designed to take advantage of modern techniques for executing dynamic languages efficiently. It has the performance of a statically compiled language while providing interactive dynamic behavior and productivity[44]. The key ingredients of performance are:

• •

  Rich type information, provided naturally by multiple dispatch;

• •

  Aggressive code specialization against run-time types;

• •

  JIT compilation using the LLVM compiler framework[45];

We build the core components of Q²Chemistry based on the Julia language to take advantage of the above advantages. Julia codes can be highly extensible thanks to its type system and multiple dispatch mechanism, and the meta-programming ability makes developing customized syntax and device-specific programs simple[46]. Generic programming in Julia helps us to reach optimized performance while still keeping the code general and concise. Julia integrates well with other programming languages. It is generally straightforward to use external libraries written in other languages in Q²Chemistry. We have also a python version of Q²Chemistry, which will be reported elsewhere. In this work, the electronic structure code PySCF[47] and openfermion[48] are used to calculate chemical quantities such as the orbital occupations and qubit Hamiltonian for quantum state evolution and measurements.

We unleash the Julia compiler on the new Sunway supercomputer and deploy Q²Chemistry on it. Fig. 6 shows the schematic overview of our Sunway Julia compiler. The MPE compiler is the main part which serves as a JIT compiler for MPEs. The SACA.jl library focuses on providing definitions for CPE operations that are required for writing effective applications based on the SACA programming model[49]. The main interface for calling functions on CPEs resembles a call to an ordinary Julia function: @saca (config...) function(args...), where the config tuple indicates the spawn configuration similar to the arguments in SACA APIs. In addition, Julia features powerful meta-programming and reflection capabilities. For example, the high-level Julia IR is accessible with the code_lowered and code_typede reflections, and can be modified with generated functions. The CPE Julia Compiler implements multiple types of extension based on SW26010Pro architecture features. The SIMD extensions support the definition of non-uniform width short vector data types and the built-in operations in Julia language layer. The multi-level storage extensions make it easy for CPE Julia code to migrate date between multi-level storage, such as LDM, TLS, CROSS, etc. The async-operation extensions provide flexible interfaces for DMA, RMA and other asynchronous components. To facilitate interactions with LLVM, we adopt LLVM.jl[50] package which provides a high-level wrapper to the LLVM C API, using Julia’s powerful Foreign Function Interface (FFI) to interact efficiently with the underlying libraries. To meet the application requirements of Q²Chemistry, we port a number of Julia standard and third-party libraries, such as MPI.jl[51], Distributed.jl, PyCall.jl, etc. MPI.jl and Distributed.jl are two commonly used distribution modes of Julia. Distributed.jl is the method to spawn new processes in separate memory spaces. Like how multi-threading is set up, its master-slave distributed architecture makes it nice and convenient for smallish parallelization. MPI.jl is a basic Julia wrapper for the Message Passing Interface (MPI)[52]. The MPI libraries are versatile and highly optimized for SW26010Pro, so that the MPI.jl could be easy to use it from Julia by preparing functions with almost the same name as C functions and get similar performance. Given the performance and scalability of Q²Chemistry on the new Sunway supercomputer, we finally chose MPI.jl as the main way for distributed computing. We use swBLAS as the backend for matrix-vector and matrix-matrix multiplications. PyCall.jl is used to directly call python packages from Julia.

When compiling Q²Chemistry, each computing node has a separate JIT compilation engine, and there is no data interaction during JIT compilation, so it can be easily scaled up to more than 20M cores. For more than 20M cores, the simultaneous loading of the Julia dynamic libraries on computing nodes can become a performance bottleneck. We avoid this problem by using the ROFS (Read Only File System). Compared with the GFS (Global File System), The ROFS has smaller capacity but faster access speed, which significantly reduces the loading time of Julia dynamic libraries.

When running Q²Chemistry, the hotspots are mainly tensor contraction and the SVD functions in the Julia’s basic library LinearAlgebra. In the tensor contraction, the first step is the index permutation of the tensors, followed by the matrix multiplication (ZGEMM) to accomplish the calculation. Here we use the fused permutation and multiplication technique. In the SVD calculation, the matrix is transformed into a bidiagonal matrix using an orthogonal transformation, then the bidiagonal matrix is diagonalized using BDC (bidiagonal divide-and-conquer) or QR decomposition. Thus the core calculations of SVD are mainly the matrix vector multiplication, matrix multiplication, vector multiplication and matrix transpose, which is done with the swBLAS library.

We also note that running Q²Chemistry on non-Sunway architectures needs to address the following issues:

1. 1.

   The compiling and running environment for Julia program on computing nodes: The Julia compilation framework needs to be supported on non-Sunway architectures to meet the compilation and running requirements of the Julia language. Currently the latest Julia version only supports i686 / x86_64 / arm / aarch64 / powerpc64le. Other architectures require porting of the Julia compilation framework;

2. 2.

   MPI.jl distributed environment: The underlying MPI protocol and the MPI.jl library need to be supported on non-Sunway architectures to meet Q²Chemistry’s massively scalable communication requirements. All of the mainstream supercomputer systems support the MPI protocol, and some architectures need to transplant the MPI.jl library;

3. 3.

   Architecture specific optimization: Optimized third-party computing libraries, including BLAS and LAPACK, need to be supported on non-Sunway architectures to meet Q²Chemistry’s computing performance requirements. The current mainstream non-Sunway Architectures basically meet this requirement.

Therefore, in general Q²Chemistry has good portability and can be rapidly and efficiently transplanted to the current mainstream non-Sunway architectures.

We firstly validate our method by comparing with the FCI results as shown in Fig. 7(a). We start with the calculation of the hydrogen ring with 10 atoms to show the effectiveness of the DMET-MPS-VQE method. Restricted Hartree-Fock method is used as the low-level calculation of the whole system, with the STO-3G basis. In the DMET-MPS-VQE calculation, the hydrogen atoms are divided into fragments with two atoms. Here we can see the potential energy curve with the DMET-MPS-VQE method, which shows very good agreement with those from FCI, and the relative errors are within 0.5%. We also compare the simulation results of MPS-VQE with FCI, and the results of the ground state energy are very accurate for the chemical systems of H₂, LiH, and H₂O, with the relative errors at the level of 0.01%.

We further examine the carbon ring molecule (C18) with a large basis set (cc-pVDZ) as shown in Fig. 7(b). The equilibrium geometry of this molecule is the bond-length alternated structure as confirmed in experiment[53]. Here we scan the energy of this molecule with the change of the bond length alternation (BLA), defined as the lengths of the two sets of C-C bonds in the molecule. We show results with cc-pVDZ basis set (for carbon, the 1s orbital is frozen, only the 2s and 2p basis orbitals of carbon atom are considered in the VQE calculation, while the high angular momentum orbitals are treated at the mean-field level). We can see that, similar to the CCSD results, our DMET-VQE results predict that the energy of the bond-length alternating structure is lower, which in good agreement with experimental observations.

In this section, we compare two other libraries with our proposed algorithm and analyze the test results. The two libraries are qiskit (using the state vector and MPS algorithm) and quimb (using the MPS algorithm). We evaluate the four methods by comparing the simulating time of one circuit for three molecule systems using one process, as shown in Fig. 8. The computing time of quimb (MPS) is on average $3$ times slower than qiskit (state-vector). Our method (Q²Chemistry) is on average around $7$ times faster than quimb and around $2$ times faster than qiskit (both MPS and state vector).

Figure 9 shows the result for the memory-efficient optimisation in the circuit simulation. Here we choose three systems as examples. The (H₂)3, LiH and H₂O molecules have $919$, $630$ and $1085$ circuits respectively. The number of circuits per process is $18$, $19$, $17$. We can see the speedup is around 15x and the memory reduction is around 20x. This indicates that the memory-optimization scheme is effective and efficient.

Figure 10 shows the simulation time with the MPS-VQE simulator for the VQE circuits of the hydrogen chain. The number of the electrons/atoms changes from 6 to 100, while the corresponding number of the qubits changes from 12 to 200. It is clearly shown that the computing time scales linearly with the number of qubits, which demonstrates the power of Q²Chemistry. There are two major bottlenecks for further scaling up to larger systems: 1) the memory requirement, which grows linearly with the number of qubits ; 2) the quantum circuit corresponding to more complex systems could be significantly deeper, for which one needs a larger bond dimension.

In the MPS-VQE simulation, the percentages of time spent in tensor contraction and SVD are around 15% and 82% when the number of qubits changes from 33 to 129. Fig. 11 give the evaluation of the tensor contraction and SVD performance improvements when using different bond dimensions. We observe that compared to the original version that only uses the MPE, the optimized version use both the MPE and the 64 CPEs produces an overall speedup in the tensor contraction calculation ranging from 2.3x to 46.5x and the SVD calculation ranging from 1.04× to 15.5×, and the speedup is more significant when increasing the bond dimension from 256 to 1024.

We also evaluate the performance of Q²Chemistry on the x86 (AMD Ryzen EPYC 7452) CPU, using a quantum circuit similar to that used in Fig. 2(c), with unitaries replaced by 2-qubit gates acting on neighbouring qubits. The initial quantum state is generated randomly according to a bond dimension threshold (which is 512 in this test). We observe that the SW version is approximately $1.1$ times faster than the x86 version based on OpenBLAS and 16.6 times faster than the x86 version based on LAPACK-3.2.

Figure 12 shows the strong scaling performance of the quantum computational chemistry simulation for the hydrogen chain with 1280 atoms (1280 electrons). In all calculations, each MPI sub-group is mapped to 2048 processes. The code shows good strong scaling performance: as the number of Sunway processes increases from $10,240$ to $327,680$ for the calculation (the number of cores increases from $665,600$ to $21,299,200$), the parallelization efficiency exceeds 92% and 30$\times$ speedup with respect to the $10,240$ processes is achieved.

Fig. 13 shows the weak scaling of the quantum computational chemistry simulation for the hydrogen chain. To obtain the weak scaling data, we time the MPS-VQE calculations with an increasing number of fragments ($2048$ processes per fragment). The number of atoms (electrons) are (40, 80, 320, 1280) respectively in the weak scaling calculation. The calculation achieves excellent weak scaling, exhibiting high parallelization efficiency as the number of cores grows. Taking $10,240$ processes as a reference, the code achieves a parallelization efficiency of around 92% with $327,680$ processes ( $21,299,200$ cores).

The reason for the high parallelization efficiency is due to the low data communication overhead in our algorithm. There is almost no data communication on the DMET level parallelization except the final summation over all fragments (a scalar from each process) to get the total DMET energy. On the MPS-VQE level parallelization over 2048 processes, we first need to broadcast the parameters from the root process to all processes using MPI_Bcast, and after the computation, we need to collect all the results to the root process using MPI_Reduce, which are shown in Fig. 4. We use the Julia performance profiler to measure the communication costs, which shows that the amount of data communication per process is around $15.6$ KB and the time spent on data communication is less than $0.001$ second per VQE iteration.