A Quantum Algorithm to Calculate Band Structure at the EOM Level of Theory

Many properties of a chemical system can be predicted by the electronic structure theory. Over the past decades, density functional theory (DFT)^{1, 2, 3, 4} has become one of the most widely used electronic structure method. However, its accuracy is difficult to be systematically improved since the exchange-correlation part of the exact energy functional in DFT is unknown. For band structure calculation, GW approximation⁵ can be used to improve the accuracy of DFT results with local functionals. However, GW calculation usually has a high computational cost and its error is also hard to assess due to different self-consistency protocols adopted in practical GW calculations⁶.

Wave function based methods, such as configuration interaction (CI), Møller-Plesset perturbation theory, and coupled cluster theory (CC), offer a systematic way to approach the exact solution of the Schrödinger equation. For example, although typical CC calculations are truncated at single and double excitations (CCSD), we can in principle obtain accurate wave function by including all possible excitations. Recently, equation-of-motion coupled-cluster (EOM-CC) theory has been implemented in studying excitations and quasiparticle band structures in solid systems^{7, 8, 9, 10, 11}. Numerical calculations indicate that EOM-CCSD produces quasiparticle energies which are usually more accurate than those from the GW approximation¹². It can be proved that the EOM-CCSDT Green’s function includes all diagrams contained in the GW approximation, along with many high-order vertex corrections. Although EOM-CC can reach a very high accuracy by incorporating high-order terms, it must be addressed that the computational cost also grows fast and quickly goes beyond the capability of current supercomputers.

Recent advances in quantum computing provides possibilities to obtain the exact solution to the many-electron Schrödinger equation within a polynomial time complexity. Quantum algorithms for electronic structure calculations include quantum phase estimation (QPE) and variational quantum eigensolver (VQE)^{13, 14, 15, 16, 17, 18, 19, 20, 21, 22, 23, 24, 25, 26, 27, 28, 29, 30, 31, 32}. The former can efficiently calculate the energy spectrum of a specific Hamiltonian. However, the high gate fidelity and error-correction devices required by QPE is typically far beyond the capability of the near-term noisy intermediate-scale quantum (NISQ) computers. By embedding efficient quantum state preparation and measurement in a classical optimization process, VQE yields a much shallower quantum circuit compared to QPE and it is thus more suitable for NISQ devices. The first experiment of VQE to study chemical systems was carried out on a photonic quantum processor to study the ground state of HeH^{+ 20}. In recent years, VQE has been demonstrated on most of the leading quantum computer architectures such as trapped-ion system^{21, 22, 23} and superconducting-qubit-based platform^{24, 25, 26}.

The key ingredient of VQE is a parametric wave function ansätze. For example, one possibility is assuming that the wave function can be written in a unitary coupled-cluster (UCC) form^{33, 34, 35}. Recently, a more efficient and accurate wave function ansätze, adaptive derivative-assembled pseudo-trotter (ADAPT)³⁶, is proposed, which relies on an iterative construction of the wave function. ADAPT algorithm provides a very accurate solution to the ground-state wave function, making it a favorable reference state for studying excitations under the EOM framework. Although VQE has been used to extract the bands of silicon from a tight-binding Hamiltonian³⁷, quantum algorithms for band structure calculations at a high ab-initio level is not yet available. It is thus desirable to extend the EOM theory to the ADAPT wave function and predict band structures in an EOM-ADAPT way on a quantum computer.

Notice that the development of an EOM-ADAPT quantum algorithm for band structure calculation is not straightforward. First, like most quantum algorithms used in chemical simulations, the ADAPT algorithm is designed for molecular systems. In a previous study³⁸, we have pointed out that a direct application of the ADAPT algorithm to periodic solid systems suffers from a residual error introduced by complex wave functions. A simple way to solve this problem is transforming complex Hartree-Fock orbitals at sampling $k$-points in a unit cell to real orbitals at the $\Gamma$ point in a supercell, which is termed as the K2G transformation³⁸. However, this leads to a significant growth in the number of operator terms since the restriction of crystal momentum conservation is removed. At the same time, it is not convenient for band structure calculations especially when a non-orthogonal $k$-grid is used. Another difficulty in EOM-ADAPT band structure calculation is that the Hamiltonian similarity transformation technique which is widely used in EOM-CCSD ^{39, 40} to satisfy the killer condition is not applicable as detailed below. Therefore, extra efforts are required to avoid deviations between different EOM functionals, which can otherwise lead to additional errors in the predicted band gap.

In this study, we present an EOM algorithm for quasiparticle band structure calculations on quantum computers. First, we introduce a modified ADAPT-VQE ansätze termed ADAPT-C for preparing the accurate ground state. In contrast to the K2G transformation scheme, the Hartree-Fock orbitals at sampling $k$-points are explicitly treated and the residual errors are minimized by introducing a complementary operator pool. Then, the quasiparticle states are obtained through a quantum subspace expansion (QSE)³¹ style technique derived from the projected-operator EOM-IP and EOM-EA formalism³⁹. The converged ground-state wave function obtained with ADAPT-C is served as the reference state, and the excitation operators are transformed using the ground-state projector to satisfy the killer condition. Finally, the quasiparticle bands are extracted by identifying correct quasiparticle excitations from the solved eigenstates. We name this method as EOM-ADAPT-C, which, as a demonstration, is successfully applied to simulate band structures of silicon and diamond.

ADAPT-C algorithm. Given the Hartree-Fock orbitals as linear combinations of Bloch atomic orbitals at sampling ${k}$ points, the Hamiltonian can be written in the second quantized form as

with

$h_{q\bm{k}_{q}}^{p\bm{k}_{p}}$ and $g_{r\bm{k}_{r}s\bm{k}_{s}}^{p\bm{k}_{p}q\bm{k}_{q}}$ are one- and two-electron integrals, respectively, which can be calculated with classical routines. Primed summations enforce crystal momentum conservation, i.e. the difference between summations of crystal momenta of creation and annihilation operators should be a reciprocal lattice vector. The creation and annihilation operators in Hamiltonian can be mapped to weighted Pauli strings and thus part of the quantum circuit using the Jordan-Wigner or Bravyi-Kitaev transformations^{41, 13, 42, 43}.

The ground-state wave function can be obtained by solving the Schrödinger equation

via a variational algorithm, where the energy expectation value is minimized with a specific wave function ansätze. The UCC ansätze is a widely used physically motivated ansätze for quantum simulations, where the wave function takes the form

$|\Phi_{0}\rangle$ is a reference state, typically the Hartree-Fock ground state. If truncated to the second order (UCCSD), the cluster operators are

It has been suggested that the performance of UCCSD can be improved by replacing occupied and virtual orbitals by general orbitals (UCCGSD)^{44, 45, 46, 47}.

Here, we use $\{i,j,k,\dots\}$, $\{a,b,c,\dots\}$, and $\{p,q,r,\dots\}$ to indicate occupied, virtual, and general orbitals, respectively.

The exponential unitary operator in the UCC ansätze cannot be converted to a quantum circuit directly. In practice, it is decomposed into one- and two-qubit quantum gates using an approximated scheme named Trotter-Suzuki decomposition^{48, 49}:

The corresponding UCC wave function has an expression of

with $\{\hat{\tau}_{i}\}$ being anti-Hermitian operators:

Here, general excitation operators $\{\hat{\tau}_{q\bm{k}_{q}}^{p\bm{k}_{p}},\hat{\tau}_{r\bm{k}_{r}s\bm{k}_{s}}^{p\bm{k}_{p}q\bm{k}_{q}}\}$ should also satisfy the crystal momentum conservation rule. The Trotterization procedure is to be truncated at finite order and the ordering of the Trotterized cluster operators is not unique, which may introduce a considerable Trotterization error^{48, 49}.

The recently proposed ADAPT ansätze³⁶ is capable of bypassing the Trotterization error by generating a compact sequence of unitary operators. Instead of optimizing a large number of UCC parameters simultaneously, a pseudo-Trotter wave function with finite variational parameters is constructed

where the exponential terms are taken from the operator pool $\mathcal{O}$ defined by anti-Hermitian operators $\hat{\tau_{i}}\in\mathcal{O}\equiv\{\hat{\tau}_{q\bm{k}_{q}}^{p\bm{k}_{p}},\hat{\tau}_{r\bm{k}_{r}s\bm{k}_{s}}^{p\bm{k}_{p}q\bm{k}_{q}}\}$. At the $k$th iteration, the energy is minimized through the variational procedure:

At the $(k+1)$th iteration, among all operators in $\mathcal{O}$, the one with the largest residual gradient is selected as $\hat{\tau}^{(k+1)}$. The residual gradient of operator $\hat{\tau}_{i}$ is defined as

The $L2$-norm of residual gradients $\vec{R}^{(k)}=(R_{1}^{(k)},R_{2}^{(k)},\dots)$ is used as the convergence criterion

Such a convergence criterion guarantees that the real part of the anti-Hermitian contracted Schrödinger equation (ACSE)⁵⁰ is satisfied. The imaginary part of ACSE is automatically satisfied only for real-valued wave functions. Therefore, there is a residual error originated from the imaginary part of ACSE when the ADAPT algorithm is straightforwardly applied to solids with complex wave function³⁸. In the UCC case, it is suggested to use complex cluster amplitudes to avoid this problem.⁵¹

In ADAPT, it is more natural to introduce a complementary operator pool $\mathcal{O}_{c}$

where $\hat{\tilde{\tau}}_{q\bm{k}_{q}}^{p\bm{k}_{p}}$ and $\hat{\tilde{\tau}}_{r\bm{k}_{r}s\bm{k}_{s}}^{p\bm{k}_{p}q\bm{k}_{q}}$ are anti-Hermitian operators defined as

The wave function is then updated using the extended operator pool

The size of the operator pool is doubled regardless of the Brillouin zone sampling scheme. At the same time, the crystal momentum conservation law is still preserved. We name such a modified ADAPT algorithm with a complementary operator pool as ADAPT-C. In ADAPT-C, variational parameters for operators in both $\mathcal{O}$ and $\mathcal{O}_{c}$ are kept to be real numbers.

EOM theory. The EOM formalism based on the coupled-cluster theory^{52, 53, 54, 55} is a well-established approach to calculate excitation energies in quantum chemistry. In EOM theory, we focus on the excitation operator $\hat{\mathcal{R}}$ connecting the ground state $|\Psi_{0}\rangle$ and a target excited state $|\Psi_{x}\rangle$. From the definition $\hat{\mathcal{R}}=|\Psi_{x}\rangle\langle\Psi_{0}|$, we can obtain an EOM-like equation for $\hat{\mathcal{R}}$

where the excitation energy is defined as $\Delta E_{x}=E_{x}-E_{0}$. Such an equation provides us a way to access excitation energies without the knowledge about the ground and excited states explicitly.

Different types of target states can be studied in the EOM theory by using different excitation operators. Ionization potential (IP) and electron affinity (EA) can be obtained if the operator $\hat{\mathcal{R}}$ includes a different number of creation and annihilation operators and thus not electron conserving. The fundamental band gap of a $N$-electron neutral system is characterized as the energy difference between the IP and EA. Therefore, we can use the EOM-IP energy $\Delta E^{IP}_{x}$ and the EOM-EA energy $\Delta E^{EA}_{x}$ at sampling $k$-points to construct the valence and conduction bands.⁷

Here, we truncate the EOM-IP and EOM-EA excitation operator up to the second order, which means the ionization energies are calculated by diagonalizing the Hamiltonian in the space of 1-hole (1$h$) and 2-hole, 1-particle (2$h$1$p$) states, while electron affinities are obtained in the space of 1-particle (1$p$) and 2-particle, 1-hole (2$p$1$h$) states. The corresponding excitation operators $\hat{\mathcal{R}}_{IP}$ and $\hat{\mathcal{R}}_{EA}$ are expressed as a linear combination of a set of Fermion excitation operators as

where the $r$-coefficients of the operator basis are parameters to be determined. Again, crystal momentum conservation rule $\bm{k}_{p}\mp\bm{k}_{q}-\bm{k}_{s}\pm\bm{k}=\bm{G}_{m}$ should be satisfied^{7, 8, 9, 10, 11}.

In practical EOM calculations, a reference state $|\Psi\rangle$ is chosen. Then, eigenvalues of Eq. (19) can be obtained from stationary values of the following functional

for any reference state $|\Psi\rangle$ that has a nonzero overlap with the exact ground state $|\Psi_{0}\rangle$. Notice that, if the operator basis to represent $\hat{\mathcal{R}}$ is incomplete, the accuracy of the calculated excitation energy depends on the choice of reference state $|\Psi\rangle$. In the Liouvillian superoperator formalism, Eq. (22) corresponds to a simple metric defined in the operator space $(\hat{\rho_{i}}|\hat{\rho_{j}})=\langle\Psi|\hat{\rho_{i}}^{\dagger}\hat{\rho_{j}}|\Psi\rangle$ ³⁹.

Sometimes, it is more convenient to choose the commutator metric $(\hat{\rho_{i}}|\hat{\rho_{j}})=\langle\Psi|[\hat{\rho_{i}}^{\dagger},\hat{\rho_{j}}]_{+}|\Psi\rangle$. Then, the corresponding EOM functional becomes

Such an expression can lower the rank of density matrices involved in classical calculations^{55, 39} and reduce the number of Pauli terms in quantum simulations. However, it is not equivalent to the simple metric functional Eq. (22) for a practical ansätze of $\hat{\mathcal{R}}$ such as those specified in Eqs. (20) and (21). If the so-called killer condition $\hat{\mathcal{R}}^{\dagger}|\Psi\rangle=0$ (it means that the ground state cannot be de-excited) is satisfied, the simple and commutator metric expressions becomes equivalent.

In band structure calculations with complex wave function involved, another important issue is if the excitation energy will be a real number. With such a consideration, a symmetric double commutator form of the EOM functional has been introduced⁵⁵

where $[\hat{\mathcal{R}}^{\dagger},\hat{H},\hat{\mathcal{R}}]_{+}=\frac{1}{2}([[\hat{\mathcal{R}}^{\dagger},\hat{H}],\hat{\mathcal{R}}]_{+}+[\hat{\mathcal{R}}^{\dagger},[\hat{H},\hat{\mathcal{R}}]]_{+})$. This EOM functional guarantees real-valued excitation energies. A similar expression has been implemented in the qEOM method to calculate the electron-excitation (EE) energies.⁵⁶ If $\hat{\mathcal{R}}=|\Psi_{x}\rangle\langle\Psi_{0}|$, Eq. (24) can be obtained directly from Eq. (19). For a practical ansätze of $\hat{\mathcal{R}}$, Eq. (24) can be obtained from Eq. (23) if the exact ground state $|\Psi_{0}\rangle$ is used as the reference state.

Since in practical calculations $\hat{\mathcal{R}}$ is typically truncated and $|\Psi\rangle$ is not the exact ground state, a direct use of Eq. (24) may lead to significant errors as will be demonstrated in Si and diamond band structure calculations. To solve this problem, we adopt the projected excitation operator formalism ³⁹ for both EOM-IP and EOM-EA

Since $\hat{\tilde{\mathcal{R}}}|\Psi\rangle=\hat{\mathcal{R}}|\Psi\rangle$, we can substituting $\hat{\mathcal{R}}$ in EOM functionals with $\hat{\tilde{\mathcal{R}}}$. Then, the killer condition becomes satisfied for any reference wavefunction. More importantly, the IP or EA energies obtained from Eq. (24) become equivalent to those from Eqs. (22) and (23).

The problem of such a projected operator formalism in quantum simulation is that the projector $|\Psi\rangle\langle\Psi|$ is difficult to be directly implemented in a quantum circuit. Therefore, a simplified working equation is obtained from either Eqs. (22), (23), or (24), which is given as

In contrast to Eqs. (22) and (23), this working equation always leads to real excitation energies.

The second term of the working equation (A1) is a constant and the first term can be converted into a generalized Hermitian eigenvalue problem in a QSE³¹ style

The matrix elements of $\bm{H}$ and $\bm{S}$ can be evaluated as

where $\hat{\rho}_{i}$ stands for Fermion excitation terms in Eqs. (20) and (21), $|\Psi_{C}\rangle$ is the ADAPT-C ground-state wave function. Calculation of these matrix elements can take advantage of a quantum computer, while eigenvalues and eigenvectors of Eq.(27) can be calculated on a classical computer.

Among all the eigenstates, only the quasiparticle states should be considered for band structure. Therefore, we denote $\vec{R}_{1}=\{r_{p\bm{k}},r^{p\bm{k}}\}$ and define the quasiparticle weight ($QPWT$) as

For example, the $QPWT$ for EOM-IP takes the form

where $\{r_{p\bm{k}}\}$ are coefficients of single ionization operators $\{\hat{a}_{p\bm{k}}\}$. To construct the band structure, we find excitation energies with a large $QPWT$, which indicates the corresponding excitation is dominated by single electron addition (for EOM-EA) or removal (for EOM-IP). The protocol described above for band structure calculation is named EOM-ADAPT-C.

Notice that, in classical EOM-CC calculations, a similarity transformation technique is widely used, where a transformed Hamiltonian $\hat{\tilde{H}}=e^{-\hat{T}}\hat{H}e^{\hat{T}}$ is considered with the Hartree-Fock wave function $|\Psi_{HF}\rangle$ as the reference state^{54, 40}. Such a scheme guarantees that the killer condition is satisfied $\hat{\mathcal{R}}^{\dagger}|\Psi_{HF}\rangle\equiv 0$, although it is not due to the orthogonality between excited and ground states as it suppose to be. The similarity-transformed Hamiltonian has the same excitation energy spectrum as the bare Hamiltonian. More importantly, due to the commutation relationship $[\hat{\mathcal{R}},e^{\pm\hat{T}}]\equiv 0$, we also have $[\hat{\tilde{H}},\hat{\mathcal{R}}]|\Psi_{HF}\rangle=\Delta E_{x}\hat{\mathcal{R}}|\Psi_{HF}\rangle$, which means that $\hat{\mathcal{R}}$ is the same for the bare and transformed Hamiltonians. In the UCC or ADAPT case with anti-Hermitian operators $\{\hat{\tau}_{i}\}$ introduced, since $[\hat{\mathcal{R}},e^{\hat{\tau}}]\neq 0$, $\hat{\mathcal{R}}$ for the bare and transformed Hamiltonians are different. Such a difference is consistent with the intermediate state representation⁵⁷, where the so-called self-consistent excitation operators $e^{T-T^{\dagger}}\hat{\mathcal{R}}e^{T^{\dagger}-T}$ is applied on the UCC or ADAPT ground state and the similarity-transformed Hamiltonian appears in the working equation. As a result of this difference, the identification of quasiparticle state will become more difficult. At the same time, estimation of the matrix element $\langle\Psi|\hat{\rho}_{u}^{\dagger}e^{T^{\dagger}-T}\hat{H}e^{T-T^{\dagger}}\hat{\rho}_{v}|\Psi\rangle$ for the similarity-transformed Hamiltonian leads to a more complicated quantum circuit. Therefore, the similarity transformation technique is not used in this study.

Numerical results. We first use a one dimensional hydrogen chain system to test the accuracy of ADAPT-C. Then, band structures of diamond and silicon are calculated using the EOM-ADAPT-C algorithm. All calculations are performed using an in-house developed code interfaced with several open-source software packages. One- and two-electron coefficients in the Hamiltonian in Eq. (1) are calculated with PySCF⁵⁸. The mapping from Fermion operators to Pauli operators are performed using the Jordan-Wigner transformation implemented in the OpenFermion code⁵⁹. In ADAPT-C optimization, the limited-memory Broyden-Fletcher-Goldfarb-Shanno (L-BFGS-B) algorithm in SciPy⁶⁰ is used. The GTH-SVZ basis set is used together with the GTH pseudopotential. Full configuration interaction (FCI) results are obtained by direct diagonalization of the system Hamiltonian in the qubit space. The convergence threshold of residual gradient $\varepsilon$ in Eq.(14) is taken to be $1\times 10^{-3}$ Hartree in both potential energy surface and band structure calculations.

A unit cell with two hydrogen atoms are used to build the one-dimensional hydrogen chain. A $1\times 1\times 4$ $k$-point grid is used for Brillouin zone sampling. We use the notation ADAPT{$X$}, ADAPT-C{$X$} and ADAPT-K2G{$X$} to represent different types of ADAPT-VQE calculations where $X\in$ {SD,GSD}. For example, ADAPT-C{GSD} stands for an ADAPT-C calculation with general single and double excitation cluster operators.

As shown in Figure 1, for both ADAPT and ADAPT-C, the generalized {GSD} gives better results than {SD}, which is in consistent with previous UCC results^{38, 47}. More importantly, the ADAPT-C ansätze significantly improves the accuracy over ADAPT by 1 to 3 orders of magnitude, since the former eliminating the residual error introduced by complex orbitals. In the whole range of H-H distance from 0.40 to 2.70 Ångstrom, the chemical accuracy is reached with ADAPT-C{GSD} and the mean absolute error (MAE) throughout the potential energy curve is only 0.06 kcal/mol, which is much smaller than that of ADAPT{GSD} (Table 1).

While the accuracy of the K2G transformation technique which we proposed previously³⁸ is at the same order of magnitude as ADAPT-C, the required computational resource is increased roughly by a factor of $N_{k}$, where $N_{k}$ is the total number of $k$-points. The numbers of terms in the Hamiltonian and the size of the operator pool scale as $O(N_{k}^{3}N_{o}^{4})$ for ADAPT and ADAPT-C, while for the K2G transformation technique this will be increased to $O(N_{k}^{4}N_{o}^{4})$, where $N_{o}$ is the number of orbitals in a unit cell. The difference becomes even more significant if higher excitation operators such as triples and quadruples are included. Therefore, the ADAPT-C ansätze provides an accurate and efficient solution to periodic system ground-state energy calculation.

Band structures of diamond and silicon are calculated at the experimental lattice constants ($a=5.431$ Å for silicon and $a=3.567$ Å for diamond). At each specific $k$-point $\bm{k}$, a two-step process is executed. (1) Perform an ADAPT-C{GSD} calculation to obtain the ground-state wave function $|\Psi_{C}\rangle$ with a $1\times 1\times 1$ $k$-point sampling grid centered at $\bm{k}$. (2) Perform EOM-ADAPT-C calculations with $\hat{\tilde{\mathcal{R}}}_{IP}|\Psi_{C}\rangle$ and $\hat{\tilde{\mathcal{R}}}_{EA}|\Psi_{C}\rangle$ as the target states. Solve the generalized eigenvalue problem defined in Eq.(28) and compare $QPWT$ to obtain the IP and EA energy spectrum at $\bm{k}$.

As shown in Figure 2, the band structures calculated with EOM-ADAPT-C agree well with classical EOM-CCSD results. The average deviation is within 0.05 eV. If DFT at the generalized gradient approximation PBE level⁶¹ is used to calculate the band structure with the same $k$-point scheme, a significant deduction of the band gap is observed as expected. On the other hand, if Eq. (24) is used directly without adopting the projected operator Eq. (25), the so-called EOM-NP result is significantly deviated from the EOM-ADAPT-C result. For example, the band gap of Si is changed from 1.04 to 0.11 eV.

Notice that noise in quantum device should be considered for real applications. Therefore, a noise-included EOM-ADAPT-C simulation to calculate the IP and EA energies of a one-dimensional hydrogen chain is performed. As shown in Figure S1, noise can indeed introduce non-negligible errors in the EOM-ADAPT-C algorithm. These errors can also be significantly reduced by applying error mitigation techniques as for other quantum algorithms.

Relatively small $k$-grids and the minimal basis set are used in this study, which are only suitable for algorithm demonstration and they are not expected to predict accurate electronic structures. Even if clever Brillouin zone sampling technique and accurate Wannier⁶² basis functions can be used and, at the same time, the requirement to use a high-order truncation in the wave function and excitation operator ansätze for strongly correlated systems is not considered, classical simulation of EOM-ADAPT-C algorithm is still strongly restricted by the computational resource cost which grows exponentially. Actually, running an EOM-ADAPT-C simulation on a classical computer is even more expensive than a classical EOM-CC calculation. For example, if the DZVP basis set with a $3\times 3\times 3$ $k$-point grid is used for C and Si band structure calculations as in the previous classical EOM-CCSD study,⁷ a total number of 1404 qubits will be required, leading to a memory consumption of approximately $2^{1404}\times 16$ bytes in a brute-force simulation, which is far beyond the capacity of any existing supercomputer.

However, if the EOM-ADAPT-C algorithm is run on a quantum computer, the expensive construction of Eq.(28) can be realized in parallel via a polynomial number of measurements on a quantum processor. Therefore, the EOM-ADAPT-C algorithm is expected to be a very useful method to obtain highly accurate band structures in the quantum computer era.

In summary, we have proposed the ADAPT-C approach for ground state quantum simulation of solid states and developed a quantum EOM algorithm for quasiparticle band structure calculations. In the ADAPT-C variational ansätze, we introduce a set of complementary operators to the original anti-Hermitian operator pool. The complementary operators are targeted at describing the imaginary part of the anti-Hermitian contracted Schrödinger equation. Optimizing the energy functional of ADAPT-C wave functions, we reach highly accurate ground-state energy as in the K2G transformation case with a significantly reduced computational cost and improved flexibility. EOM calculations can be performed using the accurate ADAPT-C ground-state wave function as the reference state. A QSE-style EOM equation is derived by adopting projected EOM-IP and EOM-EA operators to make sure that the killer condition is satisfied. As a demonstration, band structures of silicon and diamond are calculated, which is in good agreement with the classical EOM-CCSD result. The EOM-ADAPT-C method developed in this study represents an attractive way for accurate quantum simulation of band structures.

This work was partially supported by the NSFC (21825302), by the Fundamental Research Funds for the Central Universities (WK2060000018), and by the USTC Supercomputing Center.