Ab initio Quantum Simulation of Strongly Correlated Materials with Quantum Embedding

Quantum computing has shown great potential in various quantum chemical applications such as drug discovery, material design, and catalyst optimization. Although significant progress has been made in quantum simulation of simple molecules, ab initio simulation of solid-state materials on quantum computers is still in its early stage, mostly owing to the fact that the system size quickly becomes prohibitively large when approaching the thermodynamic limit. In this work, we introduce an orbital-based multi-fragment approach on top of the periodic density matrix embedding theory, resulting in a significantly smaller problem size for the current near-term quantum computer. We demonstrate the accuracy and efficiency of our method compared with the conventional methodologies and experiments on solid-state systems with complex electronic structures. These include spin polarized states of a hydrogen chain (1D-H), the equation of states of a boron nitride layer (h-BN) as well as the magnetic ordering in nickel oxide (NiO), a prototypical strongly correlated solid. Our results suggest that quantum embedding combined with a chemically intuitive fragmentation can greatly advance quantum simulation of realistic materials, thereby paving the way for solving important yet classically hard industrial problems on near-term quantum devices.


I. INTRODUCTION
Quantum simulation for the ground state problem of chemical systems has been among the most promising applications of quantum computers [1][2][3][4][5][6]. Although the practical quantum advantage of quantum chemistry is still on the way, the development of efficient quantum algorithms for the noisy intermediate-scale quantum (NISQ) devices has been a rapidly growing field of quantum technology [7]. Hybrid quantum-classical algorithms like variational quantum eigensolver (VQE) [3] and its variants have been viewed as the promising candidates that are compatible with NISQ devices for quantum simulation of chemical systems. Considerable progress centered around variational algorithms has been made in quantum simulations of molecular systems via both numerically [8][9][10][11][12] and experimental explorations [3][4][5][6][13][14][15], among which the experimental capability is no more than 20 qubits.
Despite the progress in molecular systems, algorithmic development in quantum simulation of materials is still in the early stage. Such studies could greatly enhance the capability of a quantum computer in solving both fundamental and industrial challenges, such as the understanding of unconventional superconductivity [16] and the optimization of heterogeneous catalysis [17,18]. Nonetheless, quantum simulation of periodic systems remains challenging in the foreseeable future. The primary challenge to ab initio simulation of periodic systems on a quantum computer has its root in the necessity of approaching the thermodynamic limit (TDL) for quantitative predictions. This results in an extra dimension of computational variable compared to their molecular counterparts, thus demanding numerous qubits. So far, most of the current studies have been limited to toy models of simple solids like hydrogen chain with a minimal basis set [19][20][21][22][23][24].
To tackle this problem and enable NISQ devices for quantum simulation of solid-state materials, various multiscale hybrid quantum-classical methods have been proposed [2]. Quantum embedding theory has been the most widely used framework for such a multiscale simulation of chemical systems. In principle, a quantum embedding theory divides a large system into smaller fragments, and enables an accurate treatment using a high-level computational method. Via integrating with a quantum solver like quantum phase estimation or VQE, quantum embedding allows a representation of a fragment on a NISQ device and the quantum subsystem is usually coupled to a classical environment. Fruitful progress has been made on quantum simulation of lattice models and small molecules within the framework of quantum embedding, for instance, dynamical mean-field theory (DMFT) [25][26][27][28][29], density matrix embedding theory (DMET) [30][31][32][33][34], and projection-based embedding [35]. So far the pioneering attempts of quantum simulations on periodic solids have been focused on the active space ansatz where the effective Hamiltonian is constructed using an active subspace of orbitals embedded in a Hartree-Fock (HF) or density functional theory (DFT) mean-field [36][37][38]. Although such an active space approach is a natural choice for systems where an active site can be straightforwardly selected, for example, spin defect in solids [39][40][41], a general framework for ab initio quantum simulation of solid-state materials is still missing. We are particularly interested in DMET as a promising platform to derive an efficient hybrid quantum-classical algorithm for two reasons. DMET intrinsically allows high-level treatments for multiple fragments or correlated sites at the same time unlike the active space or projection-based embedding which have been mostly developed for one single fragment or one correlated site (for instance, solid-state defect [42]). Additionally, DMET is a computationally cheaper alternative to the successful DMFT technique [43] in order to accelerate the convergence to TDL.
Inspired by recent developments in ab initio DMET and especially its multifragment extension for periodic solids [43][44][45][46][47][48], we introduce a hybrid quantum-classical algorithm within the DMET framework where the unit cell is partitioned into multiple subsets of correlating orbitals. The embedding Hamiltonian is then solved in a hybrid manner, where the strongly correlated subsets are solved by a quantum solver and the weakly correlated subsets are treated classically. This orbital-based multifragment approach is motivated by the observation that strong correlation in materials can be attributed to a few correlated orbitals, for example, 3d orbitals in transition metal oxides or cuprates [49]. The strongly correlated subsets are often small and well suited given the current limitation of quantum resources. We investigate the performance of our method using prototypical systems from weakly to strongly correlated electronic structure, including the spin polarization of one-dimensional hydrogen chain (1D-H), equation of states (EOS) of twodimensional hexagonal boron nitride (h-BN), as well as magnetic ordering of three-dimensional nickel oxide (NiO). We show that the quantum resource required for the largest simulation on NiO can be reduced to 20 qubits from requirement of 9984 qubits using a non-embedding algorithm. Overall, our results suggest that it is possible to perform quantitative prediction on realistic solids using quantum computers. This work encourages further research in ab initio simulation of large and complex periodic systems on near-term quantum devices.

II. RESULTS
A. An orbital-based multi-fragment embedding for solids In this section, we introduce the key steps in the ab initio periodic DMET algorithm used in this work. We discuss how this algorithm can be integrated with an orbital-based partition of the unit cell to enable the simultaneous use of both quantum and classical solvers, resulting in a hybrid quantum-classical algorithm for realistic solids. For periodic systems, each atomic Bloch orbital, or a k-adapted atomic orbitals (AOs), is simply the Fourier transform of a contracted Gaussian-type orbital (cGTO): where R is a real space lattice vector and k is a crystal momentum. The crystalline orbitals are linear combinations of atomic Bloch orbitals, ψ k whereâ kp † p (â kp p ) is the fermionic creation (annihilation) operator to p-th orbital in k p sampled in momentum space (kspace), with p, q, r, s denoting the indices of crystalline orbitals. Ab initio simulation on solids usually requires the computation on a supercell with a very large size. Thus, it is a common practice to employ the transnational symmetry by solving the Hamiltonian in the momentum space (k-space). In Eq. (2), we note that k p − k q = 0 and k p +k r −k s −k q = 0 owing to the translational symmetry (this is preferred as the conservation of crystal momentum). In the following, the localized orbitals are denoted by i, j, k and l while the orbitals in embedding space are denoted byĩ,j,k andl.
Total energy Schematic diagram of the hybrid quantum-classical periodic DMET algorithm for the ab initio simulation of solids. In this cartoon, we assume the unit cell is divided to two fragments a and b without loss of generality. Due to the transnational symmetry in the periodic system, we only perform the partition covering one unit cell and solve them in high-level solver without repeating the calculations on all the other unit cells.
The quantum-classical hybrid flowchart is illustrated in Fig. 1. The algorithm begins with the procedure of projecting the non-orthogonal atomic Bloch orbitals to the orthogonal localized ones, i.e. periodic intrinsic atomic orbitals (IAO) [50] and projected atomic orbitals [51], and partitioning the unit cell into several correlated fragments. In previous DMET algorithms for periodic systems, the whole unit cell is usually selected as one unique fragment embedded in an environment composed of unit cells [46,47]. Due to the translational symmetry, only one embedding problem needs to be solved. It has also been extended to multi-fragments in multi-layer compounds or inhomogeneous systems, where a sub-structure containing several atoms is chosen as the fragment [48]. In this work, we introduce an orbital-based fragmentation of which each fragment may contain several correlated orbitals from different atoms. Within this scheme, it is possible to perform the partition with arbitrary number of fragments. For the sake of simplicity, we assume the unit cell is divided into two fragments, a and b, to be solved with quantum and classical solvers, respectively.
The bath of each fragment is constructed solely by treating the rest of the supercell as the environment in R-space. For fragment a, it is built by performing singular value decomposition on the off-diagonal block of density matrix between the fragment and the environment, which is expressed as, where B is a rectangular matrix used to transform the local orbitals to the bath orbitals. The embedding orbitals are composed of the fragment and the bath orbitals. The number of the bath orbitals is no larger than that of the fragment. Regardless of the size of the supercell, the size of each embedding problem is at most 2N with N is the number of local orbitals within the fragment. Then we obtain the embedding Hamiltonian similarly to molecular systems through a projector operator (which is described in detail in the literature [33,45,47]) as followŝ whereâ † i andâj are operators that create and annihilate electrons in the embedding orbitalĩ andj, respectively. hĩj and gĩjkl are the one-and two-body component of the embedding Hamiltonian which are evaluated classically; and µ is the global chemical potential, which is a parameter to be optimized in DMET procedures to ensure the conservation of electron numbers.
Once constructed, a correlated wave function method is typically employed to solve the embedding Hamiltonian. The key difference is we only utilize the quantum resources for the embedding problem corresponding to the strongly correlated group while a classical solver, for instance, coupled cluster singles and doubles (CCSD) is used for the weakly correlated group. The size of the strongly correlated group orbitals can be managed to accommodate the current as well as any future NISQ devices. In our implementation, the embedding Hamiltonian H emb for fragment a is passed to the quantum device and solved by the VQE with a unrestricted unitary CCSD (UCCSD) ansatz [52][53][54]. The central idea of using VQE as the DMET solver is to construct and measure a parameterized embedding quantum state Ψ( θ) emb on a quantum computer (see Fig. 1). For an UCCSD ansatz, the parameter θ represents the amplitude for an excitation. After obtaining the expectation value Ψ( θ) emb |Ĥ emb |Ψ( θ) emb of embedding Hamiltonian, the amplitudes are updated on a classical computer variationally and iteratively until reach the terminal condition. The quantum state constructed from UCCSD is expressed as |Ψ( θ) emb = exp (T ( θ) −T † ( θ))|Ψ 0 , where |Ψ 0 is the Hartree-Fock ground state of the embedding in the periodic system andT represents the summation of single excitations and double-excitations and defined as, where the occupied orbitals are denoted bym andñ while the virtual orbitals are denoted byã andb; σ and τ are the corresponding spin indices. For an unrestricted UCCSD, occupation and excitation amplitude of the spin-up and spin-down orbitals are allowed to be different. The converged quantum circuit with variationally optimal parameters will be used to measure the one-and two-particle reduced density matrix, i.e. 1 D QC emb and 2 D QC emb , respectively; Similar to the classical algorithm these quantities are required to calculate the energy contribution from each fragment as well as to optimize a global chemical potential.
Once we solve the embedding Hamiltonian, the correlation potential u is employed to minimizes the difference of one-particle density matrix between the correlated density matrix 1 D QC emb and the mean-field (mf) density matrix 1 D mf emb . This is often defined as a self-consistent DMET, denoted in this work as sc-DMET, to distinguish itself with a one-shot DMET (o-DMET) where the mean-field bath is not optimized. In our algorithm, we locally search for the sub-block of the correlation potential or u a for only strongly correlated fragment a while leave the weak correlated The cost function thus can be written as where 1 D classical pq is the mean-field one-particle reduced density matrix corresponding to the total system. Note if there is more than one strongly correlated fragments, the global correlation potential u will be reconstructed by simply concatenating all submatrices u a into the corresponding submatrix. From our numerical experiments below, we have observed that this strategy usually leads to a fast convergence of the correlation potential with a satisfactory accuracy. Note that our proposed algorithm can be directly applied to multiple fragments, thus enabling a flexible framework to combine different quantum and classical solvers for the electronic structure of one crystal.
Once the self-consistent iteration terminates, the physical property of interest, such as total energy of the unit cell can be obtained through a sum over the contributions from different fragments, i.e. E cell = E a + E b , in a democratic manner that has been discussed several times in the literature [55].

B. Spin polarized states of the hydrogen chain
We first perform numerical simulations on the one-dimensional equidistant hydrogen chain (1D-H). As a simple but useful model, 1D-H is regarded as a good starting point to understand more complex systems [57][58][59]. In this part, we use a spin-unrestricted UCCSD ansatz in VQE as quantum computing (QC) solver integrated with DMET scheme to study the spin polarized 1D-H. The initial guess for the spin density is an antiferromagnetic state where the adjacent hydrogen atoms have opposite spins. In our calculations, each unit cell containing two hydrogen atoms is treated as the fragment (defined as a unit cell embedding throughout this article). Fig. 5(a) shows that the difference between the total energies calculated using QC solver (QC-DMET) and CCSD solver (CCSD-DMET) for the embedding system is negligible for this system. The deviation from k-CCSD gets larger with the increasing number of k-points as shown in Fig. S1(b). Noticeably, the results are considerately improved by the self-consistent fitting of the correlation potential, especially around the stretched region where D(H-H) ≥ 1.0Å. In Fig. 2, the EOS of 1D-H calculated by sc-DMET with QC solver is presented. As the distance between two hydrogen atoms increases, the energy curve becomes flat and the spin density on the hydrogen atom gets larger. This indicates that the system acts more like isolated atoms with a weak coupling between the hydrogen atoms. As the system is compressed, the spin polarization vanishes. For the 1x1x11 k-mesh, the phase transition occurs at ca. 0.75 A, which is consistent with the previous studies, 0.85Å [60] and 0.5-0.75Å [61]. It is also noted that the stretched regions with

C. Equation of states of 2D h-BN
We investigate the performance of DMET with the orbital-based partition in calculating the EOS of a prototypical 2D system, the hexagonal boron nitride (h-BN) [62] with weak electron correlation. We note that this system has been studied before with a unit cell or multiple unit cells as the fragments [47]. Inspired by the sp 2 hybridization of B and N [63], we divide the unit cell of h-BN into three fragments (see Fig. 3a). The first fragment consists of the 2s, 2p x , and 2p y orbitals of B while the second fragment consists of similar orbitals of N. The third fragment consists of all remaining orbitals. Since each of the first two fragments has only six spin orbitals, this allows us to utilize a QC solver while employing classical solvers for the third fragment. DMET using one unit cell as well as k-CCSD calculations are also performed for the reference. Fig. 3b presents the EOS of h-BN using different meshes for k-point. One can see that DMET using the orbitalpartition scheme combined with the QC solver underestimates the absolute correlation energy compared to k-CCSD or DMET using the unit cell as the fragment. However, the shifted EOS computed by three approaches are virtually on top of each other, suggesting that the orbital-based partition scheme accurately captures the topology of the EOS. The same trend is observed at different k-points. We anticipate that the orbital-based embedding approach combined with the chemical intuition can be helpful in predicting reaction barrier or phase transition in solid-state materials. Importantly, the size of the embedding Hamiltonian within this approach is relatively small yet can be systematically expanded by including basis functions with high angular momentum numbers. This allows a simple way to apply the  3. (a) The unit cell of BN is fragmented into three fragments and their corresponding IAOs. The first fragment is the 2s, 2px, and 2py orbitals of B while the second fragment is the same orbitals of N. The third fragment is the pz orbitals from B or N and all other orbitals; only several of them are shown here. (b) EOS of BN by DMET and k-CCSD at different meshes of k-points: the absolute EOS (left) and the shifted EOS (right). The shifted EOS for each method was obtained by shifting both the minimal energy and the lattice constant corresponding to this minimum to zero. The 5×5×1, 6×6×1, and 7×7×1 are denoted by green, red, blue, respectively. current noisy quantum solvers to realistic chemical systems at medium-or large-size basis sets.

D. Magnetic ordering in NiO
In this section, we further study the accuracy of our orbital-based partition approach using a high spin 3d solid, NiO, a strongly correlated insulator whose electronic structure and magnetic properties have been extensively studied for decades [64][65][66][67][68][69][70][71][72][73]. It is well-documented that NiO crystallizes in a rock-salt structure with a type II antiferromagnetic ground state, i.e., an AFII state. First, all the magnetic moments localized on Ni 2+ in each (111) plane are ferromagnetically aligned. The ferromagnetic (111) planes are stacked antiferromagnetically along the [111] direction, resulting in the AFII state (see Fig.4a). If they are stacked ferromagnetically (all magnetic moments are aligned in the same direction), then we have an FM state of NiO. The energetic difference between these two states of NiO could be straightforwardly calculated using the nearest-neighbor and the next-nearest-neighbor exchange interaction, J 1 and J 2 , respectively. The exchange interactions have been measured or calculated using various techniques, thereby providing reliable references to benchmark our calculations. We note that DMET and DMFT calculations using the unit cell(s) as the fragments have been reported on the AFII state of NiO in previous work [47].
In order to model the AFII and FM state of NiO, we use a rhombohedral structure with two formula units per unit cell. Within the octahedral crystal-field, the five-fold degenerate 3d orbitals of Ni 2+ are split into e g and t 2g orbitals. Inspired by the fact that the t 2g orbitals is fully filled while the e g orbitals is half filled [71] we partition the unit cell of NiO into three fragments as shown in Fig. 4b. The first fragment contains the e g orbitals of the first Ni 2+ and the p orbitals of the first O 2− . We note that the inclusion of the p orbitals from O 2− in this fragment is to account for the superexchange between Ni 2+ 3d and O 2− p orbitals [64,65]. The second fragment is composed of similar orbitals from the second Ni 2+ and O 2− . The third fragment is composed of all remaining orbitals which include the t 2g orbitals as well as other orbitals with high angular momentum quantum numbers. Similar to h-BN, our partition scheme permits a mixed solver scheme where the QC solver is used for the two small fragments which presumably account for the static correlation while the classical CCSD solver is used for the large fragment to recover the remaining dynamic correlation. The resulting solver is denoted as QC/QC/CCSD. We also perform DMET calculations using the entire unit cell as the fragment to compare with our orbital-based partition scheme.  I presents our computed energy gaps between the FM and AFII states using DMET from two different partition schemes and solvers. First, we find that o-DMET tends to underestimate the gaps between FM and AFII and does not significantly differ from the FM-AFII gap computed by HF. At the 2×2×2 k-point, o-DMET using the unit cell as the fragment even incorrectly predicts FM to be the ground state of NiO. The sc-DMET predictions are in better agreement with the existing computational and experimental FM-AFII gaps from the literature. This highlight the importance of the correlation potential fitting in DMET calculations for the electronic structure of magnetic NiO. Second, we find that the orbital-based partition scheme is overall in excellent agreement with the unit cell embedding despite employing much smaller fragments. At the 2×2×2 k-point, the energy gaps computed by the orbital-based sc-DMET are only around 20.0 meV (or 0.5 kcal/mol) different from that by k-CCSD, validating the accuracy of the orbital-based partition for calculating the FM-AFII gaps of NiO. Since there are precisely twenty spin orbitals within the embedding problem for each small fragment, we can employ VQE with twenty qubits as a DMET solver. The FM-AFII gap by the quantum DMET agrees with the classical results within 10 meV. As far as we are concerned, this is one of the first examples, if not the first, of the ab initio quantum computation on strongly correlated solids containing 3d transition metal using a medium-sized basis set. We emphasize that we do not attempt to prove a quantum advantage within this approach, rather we demonstrate how the orbital-based partition for DMET could be utilized in conjunction with quantum solvers to study realistic and challenging solids. Our proposed approach will benefit from both hardware and algorithmic advances from the quantum computing community.

E. Quantum Resource Estimation
We highlight the computational efficiency of our framework by analyzing the required qubits for each system studied in this work (see Table II). One can easily see that the required number of qubits by a naive quantum algorithm would quickly explode as the dimension of the solids increases from 1D to 2D to 3D. For example, the largest calculation performed in this work on NiO with a 4 × 4 × 4 k-mesh would require 9984 qubits, which is far beyond the current status of any near-term quantum devices.

III. DISCUSSION
In this work, we presented a flexible hybrid quantum-classical framework for quantum simulation of solid-state materials. With the multi-fragment partition of the unit cell at the orbital level, we show the possibility of ab initio study of electronic structures with reasonable resources on quantum computers. Our approach is somewhat similar to the active space partition schemes from molecular quantum chemistry, for instance, generalized active space [78,79] method, while we allow high-level treatments for multiple fragments. To investigate the accuracy of our method, we studied the antiferromagnetic electronic structure of the 1D hydrogen chain within PBC where the whole unit cell is chosen to be a correlating group of orbitals. This can be seen as a special case of our approach. The agreement between the EOS produced by our method and other classical approaches like k-CCSD or CCSD-DMET demonstrated the accuracy of the unrestricted UCCSD ansatz used in DMET. We then calculated the EOS of h-BN and compare them against those calculated by k-CCSD. A good quantitative agreement was observed, indicating the accuracy of our method for this 2D material. The strongly correlated electronic structure of NiO was also studied. Not only do our calculations agree with the k-CCSD at the available k-point, but are also well consistent with previous studies from both experimental and theoretical literature.
The spin unrestricted UCCSD ansatz was implemented and employed in the VQE procedures for 1D hydrogen chain and 3D NiO system, in particular, both occupations and orbital coefficients are not restricted to be the same for spin-up and spin-down electrons. Previous studies have focused extensively on the spin restricted wave function and the spin unrestricted ansatz has been more or less overlooked. The spin unrestricted UCCSD solver could be useful for cases when it is more convenient to break the spin symmetry to converge to the ground state. We anticipate that the ansatz improvement within this direction, for instance, an unrestricted k-UpCCGSD ansatz [80] or an unrestricted adaptive variational algorithm [81,82] may further reduce the circuit depth with a satisfactory accuracy.
We anticipate several directions in which the recent technical advances that can be augmented with our framework to improve both accuracy and efficiency when dealing with large systems. One direction relates to the systematically improvable embedding scheme. In this scheme, the conventional bath orbitals in DMET can be systematically augmented with some extra orbitals computed from local second-order Møller-Plesset perturbation (MP2) calculations, resulting in a larger correlation space and higher accuracy [83]. Another promising direction is to explore the use of quantum-classical hybrid quantum Monte Carlo (QMC) solvers, for examples, a quantum-classical hybrid full configuration interaction QMC [84] or auxiliary field QMC [15]. These QMC enhanced methods are based on the fact that quantum computers could provide a better initial trial state with good overlaps with the ground state of the system.
Finally, we propose two applications that our method could potentially apply to. First, the correlating orbitals embedding studied herein could be used to study high-temperature superconducting cuprate on a NISQ device owing to the fact that Cu 2+ 3d and O 2− 2p orbitals govern the strong correlation physics in cuprate similar to Ni 2+ 3d and O 2− 2p orbitals in the case of NiO. We note that a comprehensive DMET study on the parent state of cuprate has been recently reported by Cui el al. [48]. Thus, it would be an interesting complementary study on these materials using the orbital-based multi-fragment embedding proposed in this work. Solid-state heterogeneous catalysis [85,86] like single atom/cluster catalyst [85,86] is another important chemistry application where the 3d orbitals of transition metals play a significant role in determining the reaction barriers, thereby requiring an accurate treatment for electron correlations. In both cases, the application of our methodology requires a priori partition of different groups of active orbitals using chemical intuition. Our recommended practice is to systematically expand each correlating fragment by adding orbitals from the weakly correlated fragment until the property of interest, e.g. relative energy, is converged.
The classical calculations of FCI, CCSD and HF were performed by using PySCF in this work [87,88]. GTH pseudopotentials were employed to replace core electrons [89,90]. The GTH-SZV basis set was used to represent the valence electron of H (1s) while GTH-DZVP-MOLOPT-SR basis sets were used for B, N, O, and Ni. [91] These correspond to 2s2p3s3p3d AOs for B, N, O and 3s3p3d4s4p4d4f5s AOs for Ni. Gaussian density fitting (GDF) was applied to evaluate the two-electron integrals [92]. An even-tempered Gaussian basis set was used as the auxiliary basis set with the exponential factor β = 2.0 for 1D-H and β = 2.3 for 2D-hBN and 3D-NiO [93]. The DMET workflow was implemented based on the package libDMET [47,94]. Yao [95] was used as the backend of QC solver when simulating 1D-H and 2D-hBN. For the larger simulations of NiO, FQE [96] was chosen as the backend of QC solver to speed up the numerical calculations.
The momentum space was sampled using a Γ-centered Monkhorst-Pack grid. For 1D-H, 10Å vacuum region was set to avoid interactions between the neighboring layers due to periodic boundary conditions (PBC). A 1×1×N k k-mesh with N k = 3, 5, 7, 9, 11 was studied. Similarly for h-BN, a 20.0Å vacuum region was used to model the single layer structure. We employed three expanding k-meshes, i.e., 5×5×1, 6×6×1, and 7×7×1. For NiO, both FM and AFII states are calculated using the lattice parameter a = 4.17Å from the experiment [97]. The k-point of 2×2×2, 3×3×3 and 4×4×4 were employed, however, only the latter two were used to extrapolate the energy gaps to the thermodynamic limit (TDL). This two-point TDL extrapolation was reported to give good agreements with experimental measurements on the silicon and carbon crystal [98]. For a non-embedding reference, we were only able to perform the k-sampled periodic CCSD (k-CCSD) calculation at the 2×2×2 k-point due to the current limitation of the computational resources.
The spin density in Fig. 2 is populated in Mulliken way [56] that the off-diagonal terms in the one-particle reduced density matrix is equally attributed to the corresponding orbitals. The spin density of atom x is defined as difference of spin-up and spin-down electron density of atom x, which is expressed as where D is the one-particle reduced density matrix. The TDL extrapolation in Table. I is achieved by performing linear regression between the energy gap and N k − 1 Appendix C: The energy differences between the FM and AFII state for NiO The exchange coupling constants J 1 (nearest-neighbor) and J 2 (next-nearest-neighbor) for NiO are defined as [65] where E AFI , E AFII , and E FM are the energy of the AFI, AFII, and FM state, respectively. Thus, the energy differences between the FM and AFII state can be calculated as  [65] 116.0 SIC-LSDA (ES) [65] 106.4