Simulating periodic systems on quantum computer

For the periodic system using atom-centered basis sets, Bloch atomic orbitals (BAO) are defined as

where $\mathbf{R}_{n}$ is the translation vector and $\mathbf{k}$ is a crystal momentum vector sampled in the unit cell. $N$ is the number of unit cells. HF orbitals at $\mathbf{k}$ are expanded as a linear combination of BAOs,

which is also known as the Hartree-Fock-Roothaan approximation in first principles molecular calculations. Given HF orbitals, the corresponding one- and two-electron integrals can be computed VanKraMoh05 . The core electron potential is represented with the norm-conserving HGH pseudopotential KleByl82 ; GoeTetHut96 ; HarGoeHut98 , which removes the Coulomb singularity at the origin.

For each $\mathbf{k}$, the Hartree-Fock eigenvalue equation in the representation of BAOs is expressed as,

The Fock and overlap matrix elements are given by

Here $\mathbf{T}$ is the kinetic energy, $\mathbf{V}^{\mathrm{PP}}$ is the pseudopotential, $\mathbf{J}$ is the Coulomb integral and $\mathbf{K}$ is the exchange integral. In the real space, these matrices are given by

where $\Omega$ indicates that the real space integration is performed in the unit cell. The density matrix can be obtained by averaging over sampling $k$-points in an unit cell

where $N_{k}$ is the number of $k$-points and $N_{o}$ is the number of occupied electrons.

In the VQE algorithm, one key ingredient is the wave function ansatz, which prepares an electronic state with a few parametrized unitary operators McCRomBab16 ,

where $|\psi_{0}\rangle$ is the reference state. The parametrized wave function is optimized through Rayleigh-Ritz variational principle

The Hamiltonian in the second-quantized representation is expressed as

where $h^{p}_{q}$ is the one-electron integral, including kinetic energy and ionic potential (pseudopotential in this work) and $h^{pq}_{rs}$ is two-electron integral

The general one- and two-body excitation operators are defined as Noo00 ; MukKut04

$a_{p}^{\dagger}$ and $a_{p}$ are the second-quantized creation and annihilation operators, satisfying the anticommutation relation.

To extend VQE algorithms for periodic systems, we define the general one- and two-body operator pool with Bloch orbitals PieKowFan03 ; Dav03 ; Nak04 ,

where $\tilde{p}=p\mathbf{k}_{p}$. The Hamiltonian is summed over sampling $k$-points in an unit cell,

Because the Hamiltonian satisfies the translational symmetry based on the periodic boundary condition, general one- and two-body excitation operators must conserve crystal momentum,

where $\mathbf{k}_{p}$ and $\mathbf{k}_{r}$ are crystal momenta of creation operators and annihilation operators, respectively. $\mathbf{G}_{m}$ is a reciprocal lattice vector. The primed summation in Eq. (25) indicates that one of the orbital momenta is fixed according to Eq. (26).

In order to analyze the difference between real and complex HF orbitals, we firstly introduce the contracted Schrödinger equation (CSE), which can be derived by contraction of the Schrödinger equation onto the space of two particles

The anti-Hermitian CSE (ACSE) is expressed as

which is also known as the Brillouin condition Maz06 ; Maz07 . The ACSE can be written as a sum of the real part (ACSE-Re)

and the imaginary part (ACSE-Im)

Therefore, for the real wave function, the variationally optimized ADAPT-VQE wave function is exactly the solution of ACSE. This agrees with the conclusion that ACSE enforce stationary of the energy with respect to a sequence of unitary transformation of the reference state. As the Bloch wave function is introduced, Eq. 29 is still satisfied according to the convergence criteria in Eq. 22 while the residual error of ACSE-Im that results in the deviation of the energy is not minimized in the variational optimization procedure. In order to remove the residual error of ACSE-Im, we transform HF orbitals at sampling $k$-points in an unit cell into orbitals at $\Gamma$-point in the corresponding supercell. This generates a set of real wave function and Hamiltonian for the VQE algorithm. A brief introduction of this scheme, named K2G, is described as following. HF orbitals in Eq. 2 can be rewritten as

where

and $\tilde{\mu}$ is the $\mu$-th atomic orbital in $n$-th ”replica”. The Fock matrix with the supercell atomic orbital basis is expressed as

Diagonalizing the Fock matrix, we obtain the real HF orbitals expanded with the supercell atomic basis functions.

The quantum subspace expansion approach has experimentally and theoretically proven to be one of the most useful techniques on near-term noisy intermediate-scale quantum devices McCKimCar17 ; ColRamDah18 ; SanWanGen18 ; PauAbhChu19 . In the original implementation of QSE, the reference state is prepared with the UCCSD ansatz and excited states are obtained by diagonalizing the Hamiltonian sampled on the linear response space with single excitations McCKimCar17 . It is straightforward to include higher-order excitation operators in QSE but at the expense of a sharp increase of measurements.

Given a set of linear-response excitation operators, $\{\hat{T}_{u}\}$, the configuration state space is defined as

The QSE wave function can be expressed as a linear combination of configuration state functions

The energy and wave function of the ground and excited states can be obtained by solving a generalized eigenvalue problem in the configuration state space

with eigenvectors $\mathbf{C}$, and a diagonal matrix of eigenvalues $\mathbf{E}$. The Hamiltonian matrix elements projected onto the configuration state space are given by

The overlap matrix elements are given by

which is required because the configuration states are not necessarily orthogonal to each other.

The QSE approach is kind of inspired by the linear response method. In order to accurately describe the target state, either a well-defined reference state or the inclusion of high-order excitations is necessary to obtain a converged result. The inclusion of high-order excitations indicates a polynomially increasing computational cost, which is prohibitive in medium- and large-size calculations. The VQE algorithm, especially ADAPT-VQE, is able to efficiently generate a reasonable reference state to approximate the target state even starting from a single-configurational Hartree-Fock state. Here, we combine VQE with the QSE approach by truncating excitation operators up to second order and apply it to study the ground state and excited states.

In Figure 1, we study the ground-state potential energy curve and the absolute energy error with HF, UCCSD-VQE, UCCGSD-VQE, ADAPT(3)-SD, ADAPT(3) and FCI by varying $R$(H-H). HF fails to reproduce the exact energy curve with the mean error (ME) of 32.37 kcal/mol because the correlation effect totally misses in HF. UCCSD-VQE performs much better than HF with ME of 20.0 kcal/mol, but still significantly deviates from the exact FCI result. UCCSD is able to treat most weakly correlated systems while it suffers from the well-known problem of describing the strong electronic correlation effect. For example, as $R$(H-H) increases from 0.5 to 2.0 Å, the absolute energy error of UCCSD-VQE increase from 7.19 to 38.16 kcal/mol. Note that the error of UCCSD-VQE in calculations of the 1D hydrogen chain is much larger than those in molecular calculations. For example, the largest absolute energy error is only 9.19 kcal/mol in calculations of a challenging system N₂ with the STO-3G basis set LeeHugHea19 . ADAPT-SD shares the same operator pool with UCCSD. Here, ADAPT(3)-SD produces almost exactly the same result as UCCSD-VQE with the comparable ME of 19.87 kcal/mol. However, it should be kept in mind that ADAPT-SD and UCCSD are two different wave function ansatzes as discussed later.

UCCGSD-VQE gives a more promising ground state potential curve with the ME of only 1.58 kcal/mol. This agrees with the conclusion revealed in the previous study that the UCCGSD ansatz is far more robust and accurate than the simple UCCSD ansatz LeeHugHea19 . Analogous to ADAPT(3)-SD and UCCSD-VQE, ADAPT(3) and UCCGSD-VQE also share the same operator pool and give quite close results. However, ADAPT-VQE requires much fewer parameters to achieve the accuracy of UCCGSD-VQE. This is mainly because ADAPT-VQE approximates the exact wave function with a compact sequence of unitary transformation acting on the reference state. For example, the number of parameters in UCCGSD-VQE is fixed to be 236 while ADAPT(3) requires at most 144 parameters with $R$(H-H) varying from 0.5 to 2.0 Å. This is achieved at the cost of computing pre-estimated gradients in Eq. (23).

In addition, for ADAPT(3), the absolute energy error of the 1D hydrogen chain is 2-3 magnitude larger than those of small molecules (LiH, BeH₂ and H₆) presented in Ref. GriEcoBar19, . As discussed in Section II.3, these large deviations of the energy results from the complex wave function and Hamiltonian used in the VQE algorithm. In order to analyze errors in detail, Table 1 shows the maximum absolute residual error (MARE) of ACSE-Re and ACSE-Im for HF, UCCSD-VQE, UCCGSD-VQE and ADAPT(3). MAREs of ACSE-Im in UCCGSD-VQE and ADAPT(3) are very close to those in HF and UCCSD-VQE while MEs of the energy in UCCGSD-VQE and ADAPT(3) are much smaller than those in HF and UCCSD-VQE. This reveals that a majority of the absolute energy error in HF and UCCSD-VQE originates from larger residual errors of ACSE-Re. In UCCGSD-VQE and ADAPT(3), the variational algorithm efficiently minimizes the MARE of ACSE-Re ($<$0.1 kcal/mol) while it is not able to simultaneously minimize the MARE of ACSE-Im, which is left to be as large as $\sim$5.2 kcal/mol. Therefore, the error of UCCGSD-VQE and ADAPT(3) are largely attributed to the MARE of ACSE-Im. In order to improve the accuracy of the quantum variational algorithm for periodic systems, it is necessary to remove the residual error of ACSE-Im.

In Figure 2(a), we show the absolute energy error of various VQE-K2G approaches as a function of $R$(H-H), which systematically remove the residual error of ACSE-Im through an unitary transformation of HF orbitals at sampling $k$-points. Analogous to molecular simulations, the potential energy curve computed with $\Gamma$-ADAPT(3) agrees quite well with the exact curve since the variationally optimized wave function strictly satisfies ACSE with the norm of residual errors less than $1\times 10^{-3}$. ADAPT-VQE is an adaptive algorithm and the wave function ansatz is self-consistently grown. Therefore, when the potential energy curve is scanned, the energy discontinuity may appear if the convergence thresh is not small enough. For example, errors in the absolute energy computed with $\Gamma$-ADAPT(3)-SD dramatically fluctuates as shown in Figure 2(a). Here, we also compute the energy with a tighter convergence thresh $\epsilon$=$2\times 10^{-4}$ in $\Gamma$-ADAPT(m)-SD, named $\Gamma$-ADAPT(X)-SD, for comparison. $\Gamma$-ADAPT(X)-SD gives well converged results with ME of 0.07 kcal/mol, which is even comparable to ME of 0.05 kcal/mol in $\Gamma$-ADAPT(3).

In comparison with UCCSD-VQE, $\Gamma$-UCCSD-VQE gives a much better description of the potential energy curve of 1D hydrogen chain. The maximum deviation of the energy for $\Gamma$-UCCSD-VQE at $R=2.0$ Å is only 2.75 kcal/mol. Different from what is shown in Figure 1, $\Gamma$-ADAPT(3)-SD and $\Gamma$-ADAPT(X)-SD give more accurate results than $\Gamma$-UCCSD-VQE. In the ADAPT-VQE method, high-order excitations are ultimately included after a sequence of low-order exponentiated excitation operators acting on the wave function. A tighter convergence thresh implies that a larger $N_{k}$ in Eq. 18 is required to converge the wave function ansatz. As $N_{k}$ increases, the ADAPT-VQE wave function is in principle expanded in a larger configuration state space, which will improve the accuracy of the wave function ansatz. Therefore, the absolute energy error of $\Gamma$-ADAPT(m)-SD is reduced as $m$ increases. This is very similar with the $k$-UpCCGSD approach, in which a larger $k$ was often used to improve the accuracy LeeHugHea19 . UCCSD truncates excitation operators up to the second order and its wave function is optimized in a fixed parameter space. The inclusion of higher-order excitation operators is the most straightforward way to improve the accuracy of $\Gamma$-UCCSD while at a steeply increasing computational cost.

It is worthy to mention that the accuracy of ADAPT(m)-SD can not be improved by increasing $m$. The main reason is the limitation of crystal momenta conservation in excitation operators when HF orbitals at sampling k-points are used. Given excitation operators in Eq. (24), it is impossible to decompose these high-order excitations into the product of low-order excitation operators with crystal momentum conservation. For example, a quadruple excitation operator,

conserves crystal momentum

Eq. 26 is only sufficient but not necessary condition of Eq. (40). Therefore, ADAPT(m)-SD optimizes the wave function in a limited subspace of configuration states, which can not be simply overcome by increasing $m$. Analogous to UCCSD, the explicit inclusion of higher-order excitations is necessary to improve the accuracy of ADAPT-SD.

The quantum expansion subspace method is an alternative approach to improve the estimation of the energy over the reference state. Here, VQE/QSE can be considered as a Post-VQE method to improve the VQE algorithm. In Figure 2(b), we present the absolute energy error of the VQE/QSE approach as a function of $R$(H-H). Here, the reference states are prepared with ADAPT(3)-SD, ADAPT(1), ADAPT(2) and ADAPT(3). The Hamiltonian matrix elements are sampled in the linear-response space of the VQE wave function with single and double excitations. All energy curves computed with ADAPT(m)/QSE exactly reproduce the FCI result with errors less than $1\times 10^{-3}$ kcal/mol. This demonstrates that ADAPT(m) is able to obtain a quite reasonable reference state even with a very large convergence thresh, such as $\epsilon$=$1\times 10^{-1}$. Although the mean error of 1.30 kcal/mol in ADAPT(3)-SD/QSE gets slightly beyond the chemical accuracy, this is almost 3-4 magnitude larger than MEs of ADAPT(m)/QSE. Since the Hamiltonian matrix elements are sampled in the same linear-response space, the reference state prepared with ADAPT(3)-SD is even worst than ADAPT(1). While the number of parameters in both ADAPT(3)-SD and ADAPT(1) is 30. For QSE-SD, the Hamiltonian is sampled in a much smaller linear response space of dimension being 153 (625 for QSE). This leads to a remarkable deviation of the potential energy curve for ADAPT(3)-SD/QSE-SD with the ME of 3.75 kcal/mol. Therefore, the accuracy of the QSE approach strongly depend on the reference state and the truncation of excitation operators. Sometimes higher-order excitations should be included in order to obtain a converged result.

As the substantial success of the VQE algorithm in simulating ground-state properties, there is a broad interest in applying it to study excited states. For the 1D hydrogen chain, the potential energy surface of different electronic states is important to understand the rich and fascinating phase diagram in quantum material physics. In Figure 3, we present the first excited state energy and the absolute energy error of 1D hydrogen chain as a function of $R$(H-H). The ground state wave function is obtained with ADAPT(m) or ADAPT(m)-SD and the first excited state is computed with the QSE approach. Regardless of the convergence thresh $\epsilon$ used in ground-state ADAPT(m) calculations, we obtain exactly the same curves as FCI. For ADAPT(m), mean errors of the first excited state energy are $5.91\times 10^{-4}$, $5.05\times 10^{-4}$, $2.04\times 10^{-4}$ kcal/mol for $\epsilon$=$10^{-1}$, $\epsilon$=$10^{-2}$ and $\epsilon$=$10^{-3}$, respectively. The performance of ADAPT(3)-SD/QSE in excited-state calculations is quite similar with that in the ground state calculation. Although the absolute energy error of ADAPT(3)-SD/QSE is acceptable, we recommend to prepare the ground state with ADAPT(m) since it has been validated to be more robust and accurate.