Discontinuous Galerkin method with Voronoi partitioning for Quantum Simulation of Chemistry

Molecular orbitals based on the linear combination of Gaussian type orbitals are arguably the most employed discretization in quantum chemistry simulations, both on quantum and classical devices. To circumvent a potentially dense two-body interaction tensor and obtain lower asymptotic costs for quantum simulations of chemistry, the discontinuous Galerkin (DG) procedure using a rectangular partitioning strategy was recently piloted [McClean et al, New J. Phys. 22, 093015, 2020]. The DG approach interpolates in a controllable way between a compact description of the two-body interaction tensor through molecular orbitals and a diagonal characterization through primitive basis sets, such as a planewave dual basis set. The DG procedure gives rise to a block-diagonal representation of the two-body interaction with reduced number of two-electron repulsion integrals, which in turn reduces the cost of quantum simulations. In the present work we extend this approach to be applicable to molecular and crystalline systems of arbitrary geometry. We take advantage of the flexibility of the planewave dual basis set, and combine the discontinuous Galerkin procedure with a general partitioning strategy based on the Voronoi decomposition. We numerically investigate the performance, at the mean-field and correlated levels, with quasi-1D, 2D and 3D partitions using hydrogen chains, H$_4$, CH$_4$ as examples, respectively. We also apply the method to graphene as a prototypical example of crystalline systems.


Introduction
In [20], the computational domain is always partitioned into a set of rectangular subdomains, which is ideally suited for describing quasi-1D systems such as hydrogen chains. The rectangular partitioning strategy is also consistent with the use of the DG basis functions for density functional theory (DFT) based calculations [17,11], and will be referred to as the DG-Rectangular method (or DG-R for short). We find that when atoms are placed near the boundary of a subdomain, the error can increase compared to when the atoms appear in the center of a subdomain. This error eventually diminishes as the basis set approaches the complete basis set limit. However, this can increase the overhead of quantum simulation with a relative small basis set. One drawback of the DG-R method is that even for simple molecular configurations, it is not always possible to ensure that all atoms are located at the interior of the subdomain.
The main contribution of this paper is to present a simple solution to the problem above. Our method is based on the observation that the planewave dual basis set can be identified with a discrete set of grid points that can be partitioned into arbitrary subsets, without increasing the complexity of the algorithm or the computational cost. This allows us to partition the computational domain into subdomains of polygons of arbitrary shapes, so that all atoms are located at the interior of a subdomain. This is achieved using a partitioning strategy based on the Voronoi decomposition. We refer to this procedure as the discontinuous Galerkin formalism with Voronoi partitioning (DG-Voronoi or DG-V for short). For quasi one-dimensional systems with a linear geometry, the DG-R and DG-V partitioning strategies coincide with each other. However, the DG-V approach is much more flexible and naturally generalizes to higher-dimensional molecular geometries.
We start by introducing the discretization of the electronic-structure problem in a planewave dual basis, and define precisely what is meant by a diagonal representation of the ERI tensor and strictly localized functions. The discontinuous Galerkin approach for an arbitrary partitioning is introduced as a general framework for maintaining these properties by using the planewave dual basis set as a primitive building block. We then introduce the Voronoi partitioning strategy, which will be used in practical application to obtain a partitioning of the supercell for molecular and crystalline systems in the discontinuous Galerkin procedure. We implement the DG-V method (together with a new implementation of the DG-R method) using the PySCF software package [21]. The performance of the DG-V method is demonstrated using molecular systems including H 4 , hydrogen chain, CH 4 , as well as a periodic graphene system.

Theory
The DG-R formalism described in Ref. [20] interpolates between a primitive basis set and an active space basis set. As primitive basis set we here investigate the use of planewave dual functions. The reason for this is twofold: First, the planewave dual functions are grid based nascent delta functions, which trivializes the characterization of the basis projection matrix. Second, the periodic boundary condition module allows an implementation relying solely on PySCF-a program suit for electronic structure calculations in Python [21].

Planewave dual basis set and the active space basis set
The many-body Hamiltonian for an isolated system with N electrons and M nuclei is given bŷ where we have assumed atomic units, the Born-Oppenheimer approximation such that the positions of nuclei R I are constants, r i represent the positions of electrons, and Z I is the charge of the I-th nuclei. We assume the system is charge neutral, i.e. M I=1 Z I = N . The practical simulation of electronic structure problems often uses the supercell model. For simplicity, the computational domain Ω is defined to be an orthorhombic cell [0, L 1 ] × [0, L 2 ] × [0, L 3 ], with periodic boundary conditions. Its volume is denoted by |Ω| = L 1 L 2 L 3 . The supercell defines a Bravais lattice The supercell setup is necessary for simulating solid state systems. Even for molecular systems, the supercell model is also widely used, particularly in the context of density functional theory calculations [19,15]. In the second quantized formulation, the quantum many-body Hamiltonian in the supercell model in a planewave dual basis set takes the formĤ = Np µ,ν=1 where and v (p) µν = v per (r µ − r ν ).
The nuclei-nuclei interaction takes the form Here {χ µ (r)} Np µ=1 denotes the planewave dual basis set,b † µ ,b µ are the associated creation and annihilation operator, andn µ =b † µbµ is the number operator. Note that the Coulomb kernel should be replaced by a properly shifted Ewald potential denoted by v per (r − r ). In the following discussions, all basis sets will be represented as linear combination of the planewave dual basis set. Hence the superscript (p) indicates that the planewave dual basis set is a "primitive" basis set. We refer readers to Appendix A for more details. Here the number of primitive basis functions N p is equal to the number of grid points N g as in Eq. (A.1).
From Eq. (A.8), the electron repulsion integral (ERI) tensor in the planewave dual basis set can be viewed as a diagonal matrix. Therefore the Hamiltonian of the form Eq. (3) is called a "diagonal Hamiltonian". For diagonal Hamiltonians, the number of nonzero terms scales quadratically with respect to the number of basis functions as O(N 2 p ). For quantum simulation, the planewave dual basis set allows the efficient implementation using "fermionic fast Fourier transform" (FFFT), and "fermionic swap networks" to reduce the asymptotic complexity. We refer readers to [3,13] for more details.
The planewave basis set can be systematically improved by increasing the number of grid points N p . The main disadvantage is that N p can be very large. Even in the pseudopotential formulation, the ratio N p /N can often be 10 3 ∼ 10 4 . Therefore the direct usage of the planewave dual basis set for correlated electronic structure calculations, as well as practical simulation on quantum computers can be prohibitively expensive. As a result, quantum chemistry calculations commonly utilize physically motivated basis sets that result in a more compact description (e.g. the Pople basis sets such as X-YZg split-valence basis sets [6] or the correlation consistent basis sets [7] in the LCAO approach [10]). Let {ϕ p (r)} Na p=1 be a set of orthonormal single-particle functions; we refer to them as "active space" orbitals (for instance, canonical Hartree-Fock orbitals, or natural orbitals), which can be expanded using a primitive basis set as where Φ ∈ C Np×Na is a coefficient matrix with orthogonal columns, with entries given by the nodal representation as Φ µp = ϕ p (r µ ), and |δΩ| := |Ω|/N g the size of the volume element. For flexibility and accuracy, we will allow quite general definitions of the active space basis set in this work. It will pertain to traditional Hartree-Fock canonical orbitals as well as Gaussian basis sets such as the Dunning cc-pVDZ basis set [7]. The main advantage of using an active space basis is that we may have N a N p , and the ration N a /N is usually 10 ∼ 100.
Taking as granted the construction of the matrix Φ, we define a rotated set of creation and annihilation operators in the active space aŝ where Φ µp denotes the complex conjugate and we may project the Hamiltonian aŝ which we refer to as the active space Hamiltonian.
Generally we see that the primitive basis sets have favorable scaling in number of terms in the Hamiltonian (O(N 2 p ) vs O(N 4 a )), which often corresponds to better scaling algorithms. While N p and N a have the same asymptotic scaling, for modest sized calculations it is often observed that N p N a in order to achieve comparable accuracy. Here we will seek a way to split the difference between these two regimes by forming a more compact basis that partially retains the diagonal properties, i.e. the resulting Hamiltonian is block-diagonal.

Discontinuous Galerkin basis set
The core idea of Ref. [20] is to construct the block-diagonal basis by fitting spatially connected blocks of the primitive basis set to the active space basis set, while preserving the properties of the primitive basis set. The DG approach can therefore be seen as an interpolation between the primitive basis set and the active space basis set.
The DG approach systematically compress the active space basis set {ϕ p (r)} Na p=1 into a set of orthonormal basis functions partitioned into spatial elements-the DG elements. We enforce that the basis functions associated with different DG elements have mutually disjoint supports. Assume that the index set Ω = {1, . . . , N p } can be partitioned into The basis-projection matrix Φ can then be partitioned into where U κ is a matrix with orthonormal columns corresponding to the leading n κ singular values up to some truncation tolerance τ , yields a compressed basis Again, the coefficient matrix can be obtained from the nodal values as (U κ ) µj = φ κ,j (r µ ). Note that the matrix Φ representing the nodal values of the active space basis set does not need to be orthogonal. This is in particular the case, when the active space basis is chosen to be a Gaussian type atomic orbital basis. Yet, the compressed DG basis φ κ,j is always orthogonal as U κ is orthonormal for all κ and as DG basis functions of different DG elements have disjoint (discrete) support sets.
The basis set is adaptively compressed with respect to the given set of basis functions, and are locally supported (in a discrete sense) only on a single index set κ. We refer to this basis set {φ κ,j } as the DG basis set. Note that each DG basis function φ κ,j is a linear combination of primitive basis functions which are themselves continuous, so φ κ,j is also technically continuous in real space. In fact, φ κ,j might not be locally supported in real space if each primitive basis function χ p is delocalized. When the primitive basis functions are localized, φ κ,j can be very close to a discontinuous function (see Figure 3). When computing the projected Hamiltonian, we do not need to evaluate the surface terms in the standard DG formalism [1,17]. Furthermore, in the absence of truncation, we clearly have span{ϕ p } ⊆ span{φ κ,j }.

For the block diagonal matrix
the total number of basis functions is thus given by N d := κ∈K n κ . We remark that the number of basis functions n κ can in principle be different across different elements. We can furthermore define a new set of creation and annihilation operators that correspond to the DG basis functions, i.e.
with κ = 1, . . . , N b and j = 1, . . . , n κ . Note that unlike Eq. (8), the basis set rotation in Eq. (14) is restricted to the individual DG-elements. We readily obtain the projected Hamiltonian in the DG basis aŝ The matrix elements are Generally the one-body matrix h (d) can be a full dense matrix, but the two-body interaction tensor v (d) always takes a "block diagonal" form in the following sense (it has a specific sparsity pattern). In principle, the two-body interaction tensor in the DG basis set should take the form 1 2 Compared to Eq. (15), we find that In other words, v can be viewed as a block diagonal matrix with respect to the grouped indices (κκ , λλ ). We remark that the convergence of the DG basis set is independent of the choice of the primitive basis set so long as the primitive basis has sufficient degrees of freedom to form a good approximation to the active space functions of interest. At the end of this adaptive procedure, we expect the number of elements in the Hamiltonian to scale as O(N 2 b n 4 κ ), where we have assumed without loss generality that n κ is a constant and that N d = N b n κ .

Discontinuous Galerkin basis with Voronoi partitioning
In the previous work [20], each element κ always corresponds to a rectangular region, and hence this method will be referred to as the DG-R below. In particular, [20] uses a quasi-1D partitioning which is most suited for systems with an (approximate) linear geometry. To illustrate why a general partitioning strategy is needed, consider a minimal example of an H 4 molecule undergoing the transition from square to linear geometry (see Figure 1 and Section 3.1).
Considering the system to be in linear configuration (α = 90), we obtain four DG-R elements, but as α tends to zero, we transition to the square configuration and obtain merely two DG-R elements containing two atoms each. Therefore it is not possible to apply a quasi-1D rectangular partitioning strategy while maintaining that each atom belongs to a distinct element (note however in this example we may partition the system into two elements, so that each element contains two atoms for every configuration).  In order to find a general partitioning strategy that can be adapted to more complex geometries, we employ the Voronoi decomposition of the considered spatial domain. The Voronoi decomposition [2] partitions the domain into a set of Voronoi cells, i.e. a set of disjoint and domain filling subsets, see e.g. Figure 2(b) (this is in fact the Voronoi partitioning defined via periodic distances). Consider the supercell Ω with M nuclei placed at R 1 , ..., R M . The Voronoi cell κ I -or Thiessen polygon-associated with the I-th nuclear position R I is the set of all points in the supercell Ω whose distance to R I is not greater than the distance to any other nucleus positioned at R J , where J = I, i.e.
While the molecular configuration is always in three-dimensional space, in order to improve efficiency, we may also project all grid points, including the atomic positions to a d-dimensional subspace (e.g. a line with d = 1, or a plane with d = 2) before applying the partitioning. Since we furthermore apply the Voronoi decomposition in a supercells model with periodic boundary conditions, we utilize a periodic Voronoi decomposition. To that end, we tile 3 d − 1 copies of the original supercell around the original supercell Ω (see Figure 2(b) for an example with d = 2). Characterizing the system of tilled supercells by means of displacement vectors {R k } 3 d k=1 ⊂ L, we can adjust the definition of the Voronoi cell in Eq. (19) in the following way which is referred to as the periodic Voronoi decomposition (or simply as the Voronoi decomposition below). We see that the Voronoi decomposition fulfills by definition the premise that each grid point associated with a DG-V element is (periodically) closest to the corresponding atom. Eq. (20) also provides a practical way for implementing the Voronoi decomposition via a sorting operation. In practice we find that computing the Voronoi decomposition only takes a very small amount of time compared to other components, such as forming the two-electron repulsion integral (ERI) tensor.

Controlling the number of basis functions
Since the DG procedure effectively divides one continuous basis functions into several functions that are close to be discontinuous (see Figure 3), the number of DG functions will in general be larger than the total number of active space basis functions. Note that the number of DG basis functions is connected to the low-rank approximation imposed in Eq. (11). For systems with increasing sizes, e.g. hydrogen chains and graphene, we control the number of DG basis functions consistently through a relative-tolerance criterion: Given an exact singular value decomposition of the basis projection matrix Φ κ with singular values σ κ,1 , ..., σ κ,Na ordered non-increasingly, we neglect all singular values that are below a defined threshold ε tol when normalized by the largest singular value σ κ,1 , i.e. we neglect σ κ,i if σ κ,i /σ κ,1 < ε tol . Asymptotically, we find that the number of basis functions per element reaches a constant as the system size increases (see Section 3.2). This is also consistent with the studies of the DG basis functions in the context of DFT calculations [17,11]. However, for small molecules, it is more convenient to directly control the number of DG basis functions per element, i.e. for a given DG-V partitioning, we restrict the number of DG basis functions to an a priori defined percentage of the number of active space basis functions. We also remark that the SVD truncation criterion may also be affected by the representation of the active space basis functions as the input to the DG procedure. For instance, for a fixed given tolerance ε tol , the DG basis functions may depend on whether the active space basis functions are orthonormalized, or how they are orthonormalized. We will discuss these in Section 3.2.

Numerical Results
In this section we present in silico applications of the DG-V formalism. To that end, we have implemented the DG-V procedure in the PySCF package (version 1.6.4) [21] and apply it to the following test systems: In Section 3.1, we investigate the flexibility of the Voronoi partitioning with respect to the molecular geometry for the H 4 molecule. We study the convergence with respect to the grid size of the underlying planewave dual basis set, as well as the small "eggbox" effect due to the use of the grid. We present energy computations at the mean-field, and correlated level of theory comparing the DG-V formalism with active space basis set calculations. Furthermore, we demonstrate that even for such a small molecule, the GTO basis sets can become nearly linearly dependent, especially when diffuse basis functions are added. In comparison, the DG basis set is an orthonormal set by construction, and therefore does not suffer from the ill-conditioning issue of the overlapping matrix.
In Section 3.2, we investigate the potential benefit of using the DG-V procedure for reducing the cost of quantum simulation of chemistry. Compared to GTOs, the use of the DG-V basis set can already lead to a smaller number of nonzero two-electron repulsion integrals for a small hydrogen chain consisting of tens of atoms. Furthermore, the number of two-electron integral elements is five orders of magnitude less compared to that of the underlying planewave dual basis set.
In Section 3.3, we extend the numerical investigations of ab initio quantum-chemical simulations to methane with a three dimensional Voronoi partitioning. We present energy computations at the mean-field, and correlated level of theory comparing the DG-V formalism with active space basis set calculations.
In Section 3.4, we apply the DG-V formalism to graphene to investigate the performance for a crystalline system. We study the accuracy of Hartree-Fock calculations and MP2 calculations with DG-V basis for the supercells of 4 atoms, 8 atoms, and 16 atoms.
Additional numerical results are presented in the supplementary material.

H 4 transition from square to linear geometry
In order to investigate the flexibility of the DG-V procedure when applying it to molecular systems, we seek a model system that allows for a continuous variation of the degree of quasidegeneracy, while being small enough to enable computations using various and large basis sets at an affordable cost. Moreover, we seek a molecular geometry that cannot be covered by means of a rectangular partitioning procedure. A minimal system that fits these requirements is H 4 transitioning from square to linear geometry, see Figure 1. The degree of quasi degeneracy is here steerable thought the parameter α; as α tends to zero the systems becomes metallic, i.e. the HOMO-LUMO gap vanishes. Following Ref. [12], we set a = 2.0 (a.u.). Our selection of the active space basis set includes Pople basis sets w.o./w. diffuse basis functions (3-21G, 3-21++G, 6-31G, 6-31++G, 6-311G, and 6-311++G), as well as a correlation consistent basis set w.o./w. diffuse basis functions (cc-pVDZ, and augcc-pVDZ). These basis sets are not orthogonal, while the DG-V basis functions will be orthonormal as outlined in Section 2.2. In order to investigate the sensitivity of the Hartree-Fock energy with respect to the truncation of the DG-V basis set per element, we apply a truncation threshold relative to the maximal number of DG-V basis functions per element, which is equal to the total number of active space basis functions. This corresponds to a priori restricting the number of DG-V basis functions to be a fixed percentage of the number of active space basis functions, as outlines in Section 2.4.
In the following calculations we center the molecule in a supercell of size 12 × 12 × 12 (a.u.). As the DG-V formalism is based on a real-space grid discretization, we first study the convergence with respect to the grid spacing, i.e. in order to prohibit errors through the use of a too coarse grid, we determine the utilized grid spacing based on the convergence of the Hartree-Fock energy with respect to the underlying grid spacing, see e.g. Figure B1 for an example using 6-311G as active space basis showing the convergence of the mean-field energy, as well as convergence of the correlation energies at the MP2 level of theory. These calculations show that the Hartree-Fock energy is sufficiently converged for a grid spacing of 0.2 (a.u.), which we shall use subsequently. Another way to study the convergence of the grid spacing is to investigate the "egg-box" effect, i.e. the unphysical phenomenon that the energies are affected by the positions of the grid points relative to the nucleus. We compare the egg-box effect in total, mean-field and correlation energy (at the MP2 level) of the DG-V formalism and the periodic boundary condition module of PySCF, see Figure 4 for 6-31G as active space basis and Figure B2 for cc-pVDZ as active space basis. We find that the egg-box effect caused by the DG-V procedure is on the order of 3 × 10 −5 Hartree and is negligible. The size of the error is also comparable to that in the periodic boundary condition module of PySCF. In order to numerically investigate the applicability of the DG-V procedure to a metal-insulator type transition, we gradually change the system parameter α, i.e. we perform the following energy calculations at α = 15, 30, 45, 60, 75, 90 degrees. We first note that, compared to the rectangular partitioning of the system, the DG-V procedure provides sufficient freedom to decompose the molecular system in an optimal way, i.e. grouping the grid points that are periodically closest to a particular atom into one Voronoi cell. By construction, the DG-V elements change accordingly as we change the parameter α, so that every atom is always located in the interior of each element.
The appearance of (near-computational) linear dependence in quantum-chemical simulations is a common problem when using Gaussian type atomic orbitals. This is indeed the case for molecular systems, and even more so for crystalline systems. The case is particularly severe when diffuse functions are added to the basis set, which are important for describing, e.g., anions or dipole moments [10]. The DG approach remedies this linear dependence while extending the descriptive power of the basis set. In order to monitor the linear dependence of a given basis set, we compute the condition number of the overlap matrix. When we add diffuse basis functions the condition number of the overlap matrix is amplified, see Figure 5 for an example using (aug-)-cc-pVDZ and Figure B3 for a comparison across different basis sets. This is an expected behaviour as the diffuse functions are extended Gaussian basis functions with small exponents, however, we highlight again that by construction, the condition number of the DG-V overlap matrix is equal to one, and is therefore always well-conditioned. As outlined above, the transition of H 4 from square to linear geometry furthermore exemplifies the flexibility of the DG-V approach applied to quasi two-dimensional molecular systems, describing configurations with small HOMO-LUMO gap as well as configurations with well-separated HOMO and LUMO (as α tends to 90 degrees). In the following we investigate the robustness of the DG-V approach with respect to the imposed truncation along the change of the molecular configuration. See Figure 6 for an example using (aug-)cc-pVDZ and Figure B4, B5, and B6 for a comparison across different basis sets. Clearly, the DG-V basis set is-in its untruncated form-larger than the active-space basis, and yields therefore lower energies (given that the primitive basis set is of sufficient size). As we are interested in reducing the number of basis functions per Voronoi element, we used a truncation of the DG-V basis functions per element relative to the maximal number of DG-V basis functions, i.e. only keeping a certain percentage of the number of DG-V basis functions per Voronoi element. We observe that the Hartree-Fock energy of the considered system is fairly robust with respect to said truncation, see Figure 6. For clarity, we list the minimal and maximal difference in millihartree (i.e. max and min of E HF − E (DG) HF ) along the PES between the Hartree-Fock energy in active-basis representation and DG-V discretization for different level of truncation in Table B1 and B2.
We furthermore investigate the performance of DG-V including electron-correlation effects at the level of MP2 and CCSD. Based on Table B1, we can set an appropriate level of truncation. The corresponding potential energy surfaces are in presented in Figure 7 for the Pople bases 6-31G and 6-31++G, see also Figure B7, B8 and B9 for a comparison across different basis sets. The ratio of DG-V basis functions per Voronoi cell to the number of basis functions in regular discretization together with the respective level of truncation is listed in Table 3.1. It is interesting to observe that while the mean-field energy obtained from the DG-V basis set is only slightly lower than that from the GTOs, the correlation energy from the DG-V basis set is considerably lower. This indicates that the DG-V basis set can be effective for capturing dynamical correlation effects. on a semilogarithmic scale. Note that as the number of the DG-V basis functions increases, the energy from the (aug)-cc-pVDZ basis set is higher than that from the DG-V basis set.

Hydrogen chains
Next, we investigate hydrogen chains of various lengths as presented in Ref. [20], aiming at empirical data for cost factors in quantum simulation of chemistry. This quasi onedimensional system is well suited for examining different scaling behavior of the DG-V formalism as a function of the system size. As we systematically increase the number of atoms, we choose a consistent truncation procedure of the DG-V basis functions by means of the absolute size of the singular values. In the particular case of hydrogen chains, we observe that around 90% of singular values are enclosed in the interval [0, 1], and the largest singular value is approximately 10. Hence, in this special case, the truncation of the singular values with respect to the absolute size, can be straightforwardly related to a truncation criterion with respect to relative size. Aligning with the previous simulations we here use the cc-pVDZ basis. Furthermore, as we aim at a cost factor comparison for quantum simulation of chemistry, we need to construct a second quantized Hamiltonian with an orthogonormal basis set. To this end we compare the performance of the DG-V basis with the Löwdin orthonormalized [18] cc-pVDZ atomic orbital basis.
We also present numerical evidence that the number of DG-V basis functions per Voronoi cell as function of the system size converges to a constant that depends on the truncation tolerance. The hydrogen chains H n with n − 1 bond length of 1.7 (a.u.) are centered in a supercell of size (n · 3.6 + 10) × 10 × 10 (a.u.). The underlying grid spacing was chosen to be 0.2 (a.u.). In the case of H 30 this yields a super cell of size 118 × 10 × 10 with 590 × 50 × 50 grid points which corresponds to 1475000 planewave dual functions,   i.e. approximately 49000 basis functions per atom. The large number of planewave dual functions needed is due to that cc-pVDZ is a basis set designed for all-electron calculations. Using this configuration, we find that, even for an aggressive truncation threshold, the result from the DG-V basis can reliably matches the accuracy of the results from Gaussian orbitals, while each DG-V basis function is far more compact; see Figure 8 for a comparison of potential energy surfaces (PES) in DG-V discretizations of different tolerances with the corresponding PES in the active space basis set. In fact, the energies obtained from the DG-V basis are slightly lower than those from the Gaussian basis set. This is because as n κ increases, the span of the Gaussian basis set becomes approximately a subspace of the span of the DG-V basis set, and hence the DG-V basis can possibly yield lower energies due to the variational principle. The results suggest that the overall accuracy of the energy is relatively insensitive even to rather aggressive singular value truncations in the DG-V decomposition procedure, at least for quasi onedimensional systems. We also confirm that the average number of DG-V basis functions per atom for fixed truncation thresholds in the DG-V procedure converges rather rapidly to a constant, see Figure 9.
With the parameter settings above, we investigate the scaling of the number of nonzero two-electron integrals in Figure 10. Analogous to the procedure in Ref. [20], we neglect two-electron repulsion integrals with an absolute value smaller than 10 −6 . In particular, when the SVD truncation tolerance is set to 10 −1 , the DG-V basis set becomes advantageous when the system size is larger than 14 atoms. Comparing with results from Ref. [20], we find that the scaling of the number of nonzero two-electron integrals with respect to the number of atoms is lower in the present work using the Gaussian basis set, and the value of α is reduced from 3.9 to 3.0. A more careful investigation reveals that the difference is affected by orthogonalization strategies of the GTOs. In this work, we apply the Löwdin orthonormalized Gaussian orbitals, while  For a fixed accuracy, it is observed that the average number of DG-V functions per Voronoi call converges as a function of system size. Figure 10: The number of nonzero two-electron integrals in different representations for SVD-truncation tolerances of 10 −1 , 10 −2 and 10 −3 , plotted on a log-log scale. We fit a trendline plotted with black dots to extract the scaling as a function of system size as N α + c for some constant c, and list the exponent α beside each representation in the legend. In the shown graph, we neglect integral contributions smaller that 10 −6 .
Ref. [20] uses the SVD-orthonormalized Gaussian orbitals (the number of nonzero twoelectron integrals is reproduced in Figure C1 using the PySCF based code in this work).

Methane
We use the CH 4 molecule to demonstrate the DG-V procedure with a 3D partitioning strategy. Similar to the calculation presented in Section 3.1, we here consider the methane molecule centered in a supercell of size 15 × 15 × 15 (a.u.) using a grid spacing of 0.125 (a.u.). We again use Gaussian type atomic orbitals as active space basis set and similar to Section 3.1, the Φ matrix is non-orthogonal. As CH 4 is a relatively small system we can control the DG-V basis-set size manually by a priori defining a truncation threshold for the number of DG-V basis functions per element relative to the maximal number of DG-V basis functions per element, i.e. the total number of active space basis functions.
The application to methane shows again the increased flexibility of the Voronoi decomposition opposed to the previously implemented rectangular decomposition of the considered supercell. For a graphical illustration of the decomposition obtained for methane, we show in Figure 11 three different DG-V elements. The green-colored grid points correspond to the DG-V element around the centered carbon atom (depicted in black). The orange-and blue-colored grid points correspond to two different hydrogen atoms (depicted in blue), respectively. Similar to the calculations performed for H 4 , we apply the DG-V procedure to different atomic orbital sets: 6-31G, 6-31++G, 6-311G, 6-311++G, and the correlation consistent basis set with diffuse basis functions cc-pVDZ, and aug-cc-pVDZ. We consider methane to be in equilibrium geometry and investigate different (diffuse) basis sets. The results are listed in Table 2.
We observe that for all basis sets tested the DG-V procedure yielded improved energy results. Moreover, we observe that in most cases, except for cc-pVDZ and 6-31G, an energy improvement is obtained with a DG-V basis set truncation as aggressive as 30%. Based on the accuracy reported in Table 2, we choose the truncation threshold for the individual active space basis sets and perform correlated calculations at the MP2 level. The energy results are listed in Table 3, together with imposed relative truncation threshold and the number DG-V basis functions per element. We see that in all cases, the DG-V procedure improves the correlation energies by more than 0.01 (a.u.), which indicates the effectiveness of the DG-V basis functions for capturing dynamical correlation effects.

Graphene
In contrast to the single-molecule simulation with a large vacuum package presented in the previous section, the planewave dual basis is most natural for describing periodic systems (e.g., crystalline solids). This allows to conveniently extend quantumsimulation methods to condensed phase systems of interacting electrons. Subsequently, we exemplify this by applying the DG-V approach to graphene (see Figure 12). Aligning with the previous investigations, we start with energy computations at the mean-field level of theory. We consider graphene with a gth-pade psuedopotential provided in PySCF. Using the minimal basis set, gth-szv, we compute graphene for differently sized orthohombic cells containing up 16 atoms.
We consider three differently sized orthorhombic cells characterizing graphene: a small-size cell (containing four atoms), a middle-size cell (containing eight atoms), and a large-size cell (containing 16 atoms). For each system, we use a vacuum padding of 5 angstrom along the z-direction, i.e. the length of the different spuercells along the zdirection is always 10 angstrom to simulate the vacuum. Each of the boxes is discretized with 60 × 60 × 150 grid points. To reduce the number of DG-V basis functions, we use the molecular orbitals, instead of atomic orbitals as active space basis to generate DG-V basis. The computational results are provided in Figure 13. We decrease the relative tolerance of the singular value truncation from 0.4 to 0.05; note that this corresponds to controlling the number of DG-V basis functions through the relative-tolerance criterion outlined in Section 2.4. All calculations are performed using the gth-szv basis set as active space basis in the DG-V procedure. We find that the number of DG-V basis functions per element will increase as the relative truncation tolerance decreases, see Figure 13(c). Meanwhile, the Hartree-Fock energy of the considered systems obtained by means of the DG-V formalism decreases and eventually crosses the Hartree-Fock energy obtain by the periodic boundary condition module in PySCF, see Figure 13(a). The MP2 correlation energy of DG-V basis set converges much faster than the mean-field energy. It is surprising that the MP2 correlation energy of the DG-V procedure is already close to the result obtained by the active space basis set discretization when the relative truncation tolerance is set to be 0.4 and the corresponding number of basis functions per element is around 12. We eventually observe the crossing of the MP2 correlation energies as the relative truncation tolerance becomes smaller than 0.1, see Figure 13(b).

Conclusion
In this article we have introduced the discontinuous Galerkin formalism with Voronoi partitioning (DG-V), an extension to the discontinuous Galerkin procedure with rectangular partitioning (DG-R) presented in Ref. [20]. The aim of both DG-V and DG-R is to provide a smooth interpolation between the regime of a primitive basis sets (such as a planewave dual basis set) and an activate space basis set (such as a Gaussian basis set). The Hamiltonian, and in particular the ERI tensor in the DG-V basis set takes a block-diagonal form, which leads to asymptotically efficient quantum algorithms, but the preconstant can be significantly reduced compared to that of the primitive basis set. Compared to DG-R, the DG-V method provides a more flexible partitioning procedure allowing to apply it to molecular and solid state systems of arbitrary geometry. DG-V always produces an orthonormal basis set, and the block structure of the basis set naturally leads to easy and efficient implementations of density fitting techniques, which we will explore in the future.
Although the DG-V procedure dramatically reduces the number of basis functions compared the primitive basis set, when used for classical quantum chemistry calculations at the correlated level (e.g. MP2 or CCSD calculations), we find that the increase of the total number of basis functions compared to the standard active orbital implementation is still significant, at least for systems of relatively small sizes. This is particularly the case since our current implementation in PySCF cannot efficiently take into account the sparsity pattern of the ERI tensor in the DG-V basis for calculations at the correlated level. We expect that in the near future, the DG-V scheme could be more useful when effectively taking advantage of such sparsity pattern, such in the settings of quantum embedding methods [14,22]. In this section we briefly describe the use of the planewave dual basis set for quantum many-body calculations. Associated with the Bravais lattice is a reciprocal lattice It is convenient to introduce two sets of grids, one in the Fourier space and one in the real space. The one in the Fourier space is indexed by a subset of L * , denoted by where n 1 , n 2 , and n 3 are assumed to be even for simplicity. The second is a Cartesian grid in the real space is denoted by The total number of grid points of both grids is equal to For simplicity we shall give each point in X a lexicographic order so that X = {r µ } Ng µ=1 . Consider a smooth single-particle function ψ(r), its value at any point r ∈ X is given by The Fourier coefficients {ψ(G), G ∈ G} can be computed using the samples of ψ(r) at X:ψ (G) = Here |δΩ| := |Ω|/N g the size of the volume element. We now define a set of functions in the real space as In particular, they can be viewed as the numerical δ-functions on the discrete set X, satisfying χ r (r) = |δΩ| − 1 2 δ r,r , r, r ∈ X.
For a given function, the nodal representation (A.6) allows one to identify its function values evaluated at Cartesian grid points X with the expansion coefficients under the planewave dual basis set. This is particularly convenient when computing integrals.
The supercell model also requires us to replace the Coulomb kernel in the free space by the Ewald potential [8], denoted by v per (r−r ), defined in terms of its Fourier coefficients as and ξ is a constant. For a cubic cell of length L, we have ξ = −ν M /L where ν M is the the Madelung constant of the supercell [8]. For simplicity we have absorbed the shift due to the Madelung constant into the definition of the Ewald potential, which can significantly reduce the finite size error of mean-field energies and orbital energies. Note that both the planewave basis and the planewave dual basis are delocalized functions in the supercell. Hence a Galerkin formulation generally leads to dense matrices. For instance, the matrix due to electron-nuclei interaction V per en,µν = − is a dense matrix of size N g × N g . However, the sparsity can be improved using the pseudospectral formulation by replacing an integral by a quadrature using the nodal values of functions for r ∈ X, i.e.
V per en,µν ≈ − |δΩ| which is a diagonal matrix. Here we have used Eq. (A.5). Similarly, the two-electron repulsion integral (ERI) matrix elements can be written as If we reshape the ERI tensor as a matrix by considering the multi-indices (µν) and (λγ), respectively, the resulting matrix is a diagonal matrix. Hence Eq. (A.8) will be referred to as the diagonal representation of ERIs. This gives the Hamiltonian in the planewave dual basis in Eq. (3) in the second quantized form. Note that to reduce the number of planewave dual basis sets, the electron-nuclei interaction can also be replaced by a pseudopotential formulation. In this work, we employ the GTH pseudopotential [9] as implemented in PySCF.
Appendix B. H 4 transition from square to linear geometry In the following section we will present further numerical results that align with the results presented in Section 3.1 but extend to further basis sets.
Appendix B.1. Condition number, grid spacing and egg-box effect In this subsection we present additional data on the (near computational) linear dependence, the convergence of the Hartree-Fock energy as a function of the grid spacing and the "egg-box" effect. As remarked in Section 3.1, we have investigated the convergence behavior of the Hartree-Fock and MP2 energies as functions of the grid spacing. We found that at 0.2 (a.u.), the Hartree-Fock and MP2 energies are sufficiently converged, see Figure B1.
(a) (b) Figure B1: Convergence of mean field (a) and MP2 correlation energy (b) as functions of the grid spacing for the Pople basis sets 6-311G and 6-311++G.
Next we show an additional computation of the "egg box" effect for the DG-V procedure using the cc-pVDZ basis set. We confirm that the egg-box effect caused by the DG-V procedure is on the order of 3 × 10 −5 Hartree and is negligible. The size of the error is also comparable to that in the pbc (periodic boundary condition) module of PySCF. We have furthermore the computed condition number of the overlap matrix as a function of the system parameter α for different basis sets in Figure B3. We see that the adding diffuse basis functions increase in general the condition number, and that the condition number in all cases is significantly higher than one, which is the condition number of DG-V overlap matrix.

Appendix B.2. Mean field computations
The Figures B4, B5 and B6 show the truncation sensitivity of the Hartree-Fock energy in the DG-V basis discretization for a number of Pople basis set with added diffuse basis functions. Picking the truncation threshold for each basis set based on the calculations presented in Figures B6, B6 and B6, we can extend the results presented in Section 3.1 on the correlated calculations for the respective basis sets. The results are presented in Figures B7, B8, and B9.
We here also present numerical values showing the maximal and minimal improvement by means of the DG-V procedure. We here highlight that the main limitation of this improvement is restricted by the choice of the active basis.   ) along the PES between the Hartree-Fock energy in regular and DG-V discretization for different level of truncation. The truncation threshold for the number DG-V basis functions relative to the maximal number of DG-V basis functions per element is reported in percentage. on a semilogarithmic scale. Note that as the number of the DG-V basis increases, the energy from the 3-21G and 3-21++G basis set is higher than that from the DG-V basis set. on a semilogarithmic scale. Note that as the number of the DG-V basis increases, the energy from the 6-31G and 6-31++G basis set is higher than that from the DG-V basis set. on a semilogarithmic scale. Note that as the number of the DG-V basis increases, the energy from the 6-311G and 6-311++G basis set is higher than that from the DG-V basis set.  ) along the PES between the Hartree-Fock energy in regular and DG-V discretization for different level of truncation. The truncation threshold for the number DG-V basis functions relative to the maximal number of DG-V basis functions per element is reported in percentage.

Appendix C. Hydrogen chains
In a previous work [20], the number of nonzero two-electron repulsion integrals for hydrogen chains of increasing length was compared for discretizations through Gaussian type orbitals, DG-R and planewave dual basis, aiming at an empirical cost-factor scaling analysis for quantum simulations of chemistry. To that end, the Gaussian type orbitals were orthonormalized by means of a singular values decomposition, introducing a somewhat arbitrary gauge. We here extend this study by employing the Löwdin orthonormalization to the Gaussian type orbitals, and reporting the corresponding scaling for an active-space basis set of molecular orbitals. The considered active space of molecular orbitals (MO) consists of the N occ occupied orbitals and the lowest 3N occ virtual orbitals. We confirm the scaling behavior of the number of nonzero twoelectron repulsion integrals using an SVD orthonormalization reported in Ref. [20] for the considered system, see Figure C1(b). We furthermore observe that the scaling of the number of nonzero two-electron repulsion integrals for the molecular orbitals (MO) is close to O(N 4 ). Hence, although an active space introduces a beneficial offset, extrapolation of the curves in Figure C1 show that the use of molecular orbitals becomes more expensive than DG-V beyond 74 atoms. Comparing with the results obtained using the Löwdin orthonormalization Gaussians, see Figure C1 (a), we observe that the scaling of the number of nonzero two-electron repulsion integrals in a Gaussian atomic orbital discretization is significantly reduced. Nonetheless, we observe a crossing with a DG-V curve between 12 and 14 atoms. For certain types of quantum algorithms such as those based on the linear combination of unitaries approach [4,5], the major part of the cost does not come directly from the bare number of nonzero two-electron repulsion integrals but from the λ-value, i.e., the 1norm of the two-body interaction tensor λ = p,q,r,s |v p,q,r,s |. We find that the Löwdin orthonormalization of the Gaussian type atomic orbital basis has a more significant effect to the scaling of the λ-value than on the number of nonzero two-electron repulsion integrals, see Figure C2(a). In particular, the λ values of DG-V and active basis are almost parallel. Note however, that around 30 atoms, the λ-value for the molecular orbital basis (using an active space) crosses the Löwdin orthonormalized Gaussian basis set, i.e., in terms of the λ-value influenced cost factors in quantum simulations, the use of Löwdin orthonormalized Gaussians becomes more efficient compared to the use of molecular orbitals despite the use of a basis reduction through a active space consideration.
Although the improvement in the λ-value scaling through the DG formalism is not as pronounced as it is for the number of nonzero two-electron repulsion integrals, the DG procedure gives rise to a two-body interaction tensor with a block-diagonal sparsity structure. We expect that this can still significantly reduce the cost of implementing quantum algorithms, such as those based on fermionic swap networks [13,20].  Figure C1: The number of nonzero two-electron integrals in different representations for SVD-truncation tolerances of 10 −1 , 10 −2 and 10 −3 , plotted on a log-log scale. We fit a trendline plotted with black dots to extract the scaling as a function of system size as N α + c for some constant c, and list the exponent α beside each representation in the legend. As predicted, for these system sizes the number of two-electron integrals lies between the primitive and active-space representations, tending closer to the O(N 2 ) scaling of the primitive representation, while requiring fewer function. In both cases we orthonormalized the Gaussian basis. In (a) we used the Löwdin orthonormalization and in (b) SVD orthonormalization introducing an arbitrary gauge.  Figure C2: λ value for Gaussian and DG-V basis, which is related to the cost of quantum algorithms such as linear combination of unitaries (LCU). The value of λ is plotted as a function of system size for different representations. We fit a trendline plotted with black dots from the second point onward to extract the scaling as a function of system size as N α + c for some constant c, and list the exponent α beside each representation in the figure legend. In (a) we used the Löwdin orthonormalization and in (b) SVD orthonormalization introducing an arbitrary gauge.