Staggered mesh method for correlation energy calculations of solids: Second order Møller-Plesset perturbation theory
aa r X i v : . [ phy s i c s . c o m p - ph ] F e b Staggered mesh method for correlation energycalculations of solids: Second orderMøller-Plesset perturbation theory
Xin Xing, † Xiaoxu Li, † , ‡ and Lin Lin ∗ , † , P , S Department of Mathematics, University of California, Berkeley, CA 94720, USA, School ofMathematical Sciences, Beijing Normal University, No.19, Xinjiekouwai St, HaidianDistrict, Beijing, 100875, P.R.China, Computational Research Division, Lawrence BerkeleyNational Laboratory, Berkeley, Berkeley, CA 94720, USA, and Challenge Institute forQuantum Computation, University of California, Berkeley, CA 94720, USA
E-mail: [email protected]
Abstract
The calculation of the MP2 correlation energy for extended systems can be viewedas a multi-dimensional integral in the thermodynamic limit, and the standard methodfor evaluating the MP2 energy can be viewed as a trapezoidal quadrature scheme.We demonstrate that existing analysis neglects certain contributions due to the non-smoothness of the integrand, and may significantly underestimate finite-size errors.We propose a new staggered mesh method, which uses two staggered Monkhorst-Pack ∗ To whom correspondence should be addressed † Department of Mathematics, University of California, Berkeley, CA 94720, USA ‡ School of Mathematical Sciences, Beijing Normal University, No.19, Xinjiekouwai St, Haidian District,Beijing, 100875, P.R.China P Computational Research Division, Lawrence Berkeley National Laboratory, Berkeley, Berkeley, CA94720, USA S Challenge Institute for Quantum Computation, University of California, Berkeley, CA 94720, USA eshes for occupied and virtual orbitals, respectively, to compute the MP2 energy. Wedemonstrate that the staggered mesh method circumvents a significant error source inthe standard method, in which certain quadrature nodes are always placed on pointswhere the integrand is discontinuous. Numerical results indicate that the staggeredmesh method can be particularly advantageous for quasi-1D systems, as well as quasi-2D and 3D systems with certain symmetries. Correlated wavefunction based methods have long been the standard in quantum chemistryfor accurate solution of the many-electron Schr¨odinger equation in molecular systems. In re-cent years, they are also increasingly used for evaluating energies beyond the mean-field levelin extended systems . In contrast to the zero dimensional molecular systems, propertiesin bulk solids, surfaces and other low-dimensional extended systems need to be calculatedproperly in the thermodynamic limit (TDL). Due to the steep increase of the computationalcost with respect to the system size, reaching convergence in a brute force fashion is oftenbeyond reach, and finite-size corrections must be applied. Common correction methods usedto reduce the finite-size errors in correlation energy calculations include power-law extrapo-lation , structure factor extrapolation , and twist averaging .Unless otherwise stated, throughout the paper, we assume the system extends along allthree dimensions, and a standard Monkhorst-Pack (MP) mesh with N k points sampled inthe first Brillouin zone (BZ) is used. The power law extrapolation typically assumes thatthe finite-size error is proportional to N − / k , N − k , or their linear combinations. The N − / k scaling is due to the fact that the correlation energy may inherit the O ( N − / k ) finite-sizeerror in HF orbital energies . The finite-size errors in the orbital energies can be reducedto O ( N − k ) via the Madelung-constant correction . With this error removed, it has beenargued based on structure factor analysis that the finite-size error in the correlation energyscales as O ( N − k ) due to the omission of certain terms in the structure factor . The struc-2ure factor extrapolation method, as its name suggests, computes the finite-size correctionby extrapolating the omitted structure factor around the Coulomb singularity. The twistaveraging technique calculates and averages the structure factors, and consequently the cor-relation energies using a set of shifted k -point meshes, and is often used as a pre-processingfor power-law extrapolation and structure factor interpolation. The effectiveness of thesecorrection methods can often be strongly system-dependent in practice.In this paper, we focus on the finite-size error of correlation energy calculations andits correction in the simplest scenario, namely the correlation energy from the second or-der Møller-Plesset perturbation theory (MP2) for insulating systems (the MP2 energies formetallic systems may diverge ). In the TDL, the MP2 energy can be expressed as anintegral in the BZ. The numerical evaluation of the MP2 energy then uses a trapezoidalquadrature to replace the integral by a finite sum over the MP mesh. Correspondingly, thefinite-size error in MP2 energy arises from two sources: the error of the integrand, and theerror of the numerical quadrature. The first error comes from the basis set incompletenessand finite-size errors in orbitals and orbital energies, and can be reduced by various existingtechniques .The integrand of the MP2 energy calculation generally has many discontinuous points.In this paper, we demonstrate that existing structure-factor based error analysis neglectscertain contributions due to the discontinuous behavior of the integrand, and underestimatesthe finite-size errors from the numerical quadrature. We show that the error of the numer-ical quadrature comes from the overall non-smoothness of the integrand and the improperplacement of the quadrature nodes. Particularly, the standard MP2 calculation uses thesame MP mesh for both occupied and virtual orbitals. This leads to the sampling of certain q points (the difference between the k points of an occupied-virtual orbital pair) on whichthe integrand is discontinuous. The error due to such improper placement of the quadraturenodes is O ( N − k ).We propose a simple modification to address this problem. Our staggered mesh method3ses one MP mesh for occupied orbitals, and another MP mesh shifted by half mesh size forvirtual orbitals. We show that the integrand is well defined on all q points in the numericalcalculation, thus circumventing the need of structure factor extrapolation. We show thatthe finite-size error of the staggered mesh method is mainly affected by the intrinsic non-smoothness of the integrand in the MP2 calculation.We compare the performance of the staggered mesh method and the standard methodfor a model problem, where the mean-field orbital energies and wavefunctions are obtainedaccurately from a given effective potential. We then demonstrate numerical tests on periodicH dimer, and silicon in the quasi-1D, 2D and 3D bulk settings using the PySCF package.Our results indicate that the use of the staggered mesh can significantly accelerate the con-vergence towards the TDL in two scenarios: 1) quasi-1D systems, where the non-smoothnessof the integrand is removable, 2) quasi-2D or 3D bulk systems with certain symmetries. Let Ω be the unit cell, | Ω | be its volume, and Ω ∗ be the associated BZ. The Bravais latticeis denoted by L and its associated reciprocal lattice is denoted by L ∗ . The MP mesh is usedfor k -point sampling in Ω ∗ and N k denotes the total number of k points. When the MPmesh contains the Γ-point, the system can be identified with a periodic supercell Ω S withvolume (cid:12)(cid:12) Ω S (cid:12)(cid:12) = N k | Ω | . Each molecular orbital can be written as ψ n k ( r ) = 1 √ N k e i k · r u n k ( r ) = 1 | Ω | √ N k X G ∈ L ∗ ˆ u n k ( G ) e i( k + G ) · r , where n is a generic band index, and u n k is periodic with respect to the unit cell. We alsodefine the pair product (of the periodic components) as ̺ n ′ k ′ ,n k ( r ) = u ∗ n ′ k ′ ( r ) u n k ( r ) := 1 | Ω | X G ∈ L ∗ ˆ ̺ n ′ k ′ ,n k ( G ) e i G · r . n ∈ { i, j } refers to the occupied orbital and n ∈ { a, b } refers tothe unoccupied orbital. The two-electron repulsion integral (ERI) tensor in the molecularorbital basis can be written as h i k i , j k j | a k a , b k b i = 1 | Ω S | X ′ G ∈ L ∗ π | q + G | ˆ ̺ i k i ,a k a ( G ) ˆ ̺ j k j ,b k b ( G k a , k b k i , k j − G ) . (1)Here k a − k i =: q , and we have k i + k j − k a − k b =: G k a , k b k i , k j ∈ L ∗ by crystal momentumconservation. The notation P ′ G ∈ L ∗ means that the possible term with q + G = is excluded.The correlation energy per unit cell, up to the levels of coupled cluster singles and doubles(CCSD), is given by E c = 1 N k X ijab X k i k j k a k b (2 h i k i , j k j | a k a , b k b i − h i k i , j k j | b k b , a k a i ) T a k a ,b k b i k i ,j k j , (2)where k i , k j , k a , k b ∈ Ω ∗ . Here T a k a ,b k b i k i ,j k j = t a k a ,b k b i k i ,j k j + t a k a i k i t b k b j k j , and t a k a i k i and t a k a ,b k b i k i ,j k j are singlesand doubles amplitudes obtained from solution of related amplitude equations. In the cou-pled cluster doubles (CCD) theory, we have t a k a i k i = 0, and the MP2 energy is further givenby setting the doubles amplitude to t a k a ,b k b i k i ,j k j = h a k a , b k b | i k i , j k j i ε i k i + ε j k j − ε a k a − ε b k b . (3)Note that Eq. (2) can be rewritten as E c = 1 N k | Ω S | X ijab X k i k j k a k b h i k i , j k j | a k a , b k b i e T a k a ,b k b i k i ,j k j , (4)where we have absorbed the exchange term into the redefined amplitude e T a k a ,b k b i k i ,j k j = (cid:12)(cid:12) Ω S (cid:12)(cid:12) (cid:16) T a k a ,b k b i k i ,j k j − T b k b ,a k a i k i ,j k j (cid:17) , (cid:12)(cid:12) Ω S (cid:12)(cid:12) ensures that each entry e T a k a ,b k b i k i ,j k j does not vanish in the TDL.In order to write down the correlation energy in the TDL, we use the fact that boththe ERI tensor and the T amplitude do not change if we replace any k by k + G for some G ∈ L ∗ . Then fixing k i ∈ Ω ∗ , we may shift k a by some G vector so that the difference q = k a − k i ∈ Ω ∗ . Similarly further fixing k j ∈ Ω ∗ , we may shift k b so that G k a , k b k i , k j = , i.e. k b = k j − q . Note that this requires redefining ˆ ̺ n ′ k ′ ,n k to accommodate the case where k isoutside Ω ∗ . More importantly, such manipulation is only formal and is introduced to simplifythe theoretical analysis. In practical calculations, we may still keep k i , k j , k a , k b ∈ Ω ∗ as instandard implementations. After such modifications, E c in the TDL as N k → ∞ can beconcisely written as a triple integral over BZ (which is a 9-dimensional integral for 3D bulksystems): E TDL c = Z Ω ∗ d q Z Ω ∗ d k i Z Ω ∗ d k j | Ω | (2 π ) X ijab X ′ G ∈ L ∗ π | q + G | ˆ ̺ i k i ,a ( k i + q ) ( G ) ˆ ̺ j k j ,b ( k j − q ) ( − G ) e T a ( k i + q ) ,b ( k j − q ) i k i ,j k j . (5)In this continuous formulation, we may also write P ′ G ∈ L ∗ simply as the regular summation P G ∈ L ∗ since the singularity set { q + G = , q ∈ Ω ∗ , G ∈ L ∗ } = { q = , G = } is only anisolated point. All numerical schemes for evaluating the correlation energy in the TDL amounts to approx-imating the triple integral (5). The quality of the numerical approximation can be affectedby the following error sources: 1) The error introduced by replacing the integral (5) by anumerical quadrature (4), 2) The mean-field orbital energies { ε n k } and orbitals { u n k ( r ) } are not evaluated in the TDL, 3) Basis set incompleteness error, and the truncation of thevirtual orbitals due to the finite size of the basis set, 4) Error in evaluating the T -amplitudes.The last three sources contribute to the errors of the integrand values used in the numerical6uadrature (4).This paper only concerns the first error, i.e. the quadrature error. We assume that mean-field calculations are less expensive than correlation energy calculations, and the finite-sizeerror of the orbitals and orbital energies could be reduced by using other correction methodsand/or a large enough MP mesh if needed. Even when the same MP mesh is used toevaluate mean-field energies and orbitals, after the the Madelung-constant correction to theoccupied orbital energies, the contribution of the finite-size from the orbital energies becomes O ( N − k ) . The error due to the incompleteness of the basis set error and the truncation of thevirtual orbitals are more difficult to assess. Though such errors can be reduced via power-lawextrapolation or explicit correlation methods , we will not consider such improvementsin this paper. We will also only consider the evaluation of the MP2 energy, where the T -amplitudes are given explicitly by orbital energies and ERIs. We will demonstrate belowthat even under such assumptions, the finite-size effect due to the quadrature error remainssignificant.To connect to the commonly used argument in the literature to analyze the quadra-ture error using structure factors, we note that the structure factor S q ( G ) corresponds to apart of the integrand in (5) as S q ( G ) = Z Ω ∗ d k i Z Ω ∗ d k j | Ω | (2 π ) X ijab ˆ ̺ i k i ,a ( k i + q ) ( G ) ˆ ̺ j k j ,b ( k j − q ) ( − G ) e T a ( k i + q ) ,b ( k j − q ) i k i ,j k j . (6)The correlation energy is then E TDL c = Z Ω ∗ d q X ′ G ∈ L ∗ π | q + G | S q ( G ) . (7)We may also combine the information from the structure factors and define the integrand ofEq. (7) as h ( q ) = X ′ G ∈ L ∗ π | q + G | S q ( G ) . (8)7he standard MP2 calculation (4) can be interpreted as two quadrature steps in estimat-ing each S q ( G ) at a finite set of q points and E TDL c as, S q ( G ) ≈ | Ω ∗ | N k X k i , k j ∈K | Ω | (2 π ) X ijab ˆ ̺ i k i ,a ( k i + q ) ( G ) ˆ ̺ j k j ,b ( k j − q ) ( − G ) e T a ( k i + q ) ,b ( k j − q ) i k i ,j k j ! =: e S q ( G ) , q ∈ K q , G ∈ L ∗ , (9) E TDL c ≈ | Ω ∗ | N k X q ∈K q X ′ G ∈ L ∗ π | q + G | e S q ( G ) ! , (10)where K denotes the MP mesh and K q is a same-sized MP mesh containing all q ∈ Ω ∗ defined as the minimum image of k a − k i with k i , k a ∈ K . Furthermore, K q always includesthe Γ-point. These two steps apply the trapezoidal rules with uniform meshes K × K and K q for (6) and (7), respectively.Note that the integrand in (7) is discontinuous at q = and its value at this point isindeterminate due to the term (4 π/ | q | ) S q ( ) with G = . It has been argued that for q + G = , S q ( G ) converges quickly, and hence the error is mainly due to the neglect ofthis discontinuous term from the primed summation in (10), which scales as N − k ∼ (cid:12)(cid:12) Ω S (cid:12)(cid:12) − .However, such an analysis misses two factors:1) Each S q ( G ) in (9) also neglects certain discontinuous terms. Specifically, fixing q and G , the amplitude e T a ( k i + q ) ,b ( k j − q ) i k i ,j k j in the integrand for S q ( G ) in (6) is discontinuous as afunction of ( k i , k j ) when k j − k i − q ∈ L ∗ due to its exchange part, i.e., | Ω S | T b ( k j − q ) ,a ( k i + q ) i k i ,j k j = P ′ G ′ ∈ L ∗ π | k j − k i − q + G ′ | ˆ ̺ ∗ i k i ,b ( k j − q ) ( G ′ ) ˆ ̺ ∗ j k j ,a ( k i + q ) ( − G ′ ) ε i k i + ε j k j − ε b ( k j − q ) − ε a ( k i + q ) . For each pair ( k i , k j ) satisfying the relation k j − k i − q ∈ L ∗ , the exchange term above neglectsthe summation term associated with k j − k i − q + G ′ = , leading to N − k ∼ | Ω S | − error inthe associated volume element corresponding to the multi-index ( k i , k j ). For each q ∈ K q ,there are exactly N k such pairs ( k i , k j ) ∈ K × K . Overall, neglecting the discontinuousterms when evaluating e T a ( k i + q ) ,b ( k j − q ) i k i ,j k j at these quadrature nodes leads to O ( N − k ) error in8omputing each S q ( G ), and hence the additional O ( N − k ) error in computing E TDL c .2) The integrand of Eq. (7), i.e., h ( q ), is not a smooth function, which leads to significantquadrature errors. The error of the trapezoidal rule can be generally analyzed using theEuler-Maclaurin expansion. Let δk denote the mesh size along each direction (i.e., N k ∼ ( δk ) − d for systems that extend along d dimensions). Then for a periodic function withcontinuous derivatives up to m -th order, the quadrature error can be as small as O ( δk m ).However, the integrand for E TDL c already has unbounded second order derivatives. Thereforestandard error analysis predicts that the quadrature error can be O ( δk ) = O ( N − / k ), oreven worse, for three-dimensional systems. If so, the finite-size errors would in general bedominated by such quadrature errors. Fortunately, the points of discontinuity are isolated,and we find that the quadrature error should be O ( δk ) = O ( N − k ) in the worst case for three-dimensional systems. However, the analysis is much more involved than the direct applicationof the Euler-Maclaurin expansion. Instead it generalizes the result of Lyness for a classof punctured trapezoidal rules, and we will report the full analysis in a future publication.Furthermore, for systems with certain symmetries (for instance, three-dimensional systemswith cubic symmetries), the smoothness condition of the integrand can be improved, whichleads to quadrature error that decays faster than O ( N − k ), and such faster decay agrees withthe observations in the literature .Our analysis above is also applicable to quasi-1D and quasi-2D systems, which samples k points on the corresponding 1D axis and 2D plane in Ω ∗ , respectively. Without loss ofgenerality we may assume the MP mesh includes k points of the form k = (0 , , k z ) forquasi-1D systems, and k = (0 , k y , k z ) for quasi-2D systems. The correlation energies of thismodel in the TDL can be written in an integral form similar to Eq. (5), while only changingthe integration domains for k i , k j , and q from Ω ∗ to the corresponding axis/plane in Ω ∗ .The discontinuity of the integrands in (6) and (7) described for 3D systems earlier stillformally applies to this low-dimensional model. Neglecting discontinuous terms in quasi-1Dand quasi-2D systems also leads to O ( N − k ) quadrature error in the structure factor and the9P2 energy. Furthermore, the quadrature error due to the non-smoothness of the integrandis at most O ( N − k ) for quasi-2D systems. The situation for quasi-1D system is qualitativelydifferent. This is because the discontinuous points in quasi-1D systems are removable, i.e., byproperly redefining the integrand values at these isolated points, h ( q ) can become a smoothfunction (see the numerical examples in Figures 1 and 4). Therefore with a properly definedintegrand, the quadrature error for quasi-1D systems can decay super-algebraically (i.e. thequadrature error decays asymptotically faster than O ( δk m ) for any m > e T a ( k i + q ) ,b ( k j − q ) i k i ,j k j and h ( q ) witha quasi-1D model problem (details of the model are given in Section 3). According to thediscussion above, such discontinuous behavior is generic in MP2 calculations. The standardMP2 calculation with any k -point mesh K always places some of its quadrature nodes atsuch points of discontinuity.The discontinuity of h ( q ) at q = is generally not removable in quasi-2D and 3D systems(similarly for the discontinuity of the integrand in (6) for computing S q ( G ) and h ( q )). Forsystems with certain symmetries, lim q → h ( q ) may exist and we can redefine the value of h ( )as this limit. The resulting modified integrand has improved smoothness and the resultingquadrature error can be o ( N − k ). In this case, the overall quadrature error will be dominatedby improper placement of the quadrature nodes. As an example, Figure 2 illustrates thediscontinuity of h ( q ) obtained from two quasi-2D model problems which have an isotropicand an anisotropic Gaussian effective potential fields, respectively. The additional symmetryfrom the isotropic potential leads to the removable discontinuity at q = for h ( q ), while inthe anisotropic case, the values of h ( q ) along the x, y axes are very different near q = , andhence lim q → h ( q ) is not well defined. 10 a) | e T a ( k i + q ) ,b ( k j − q ) i k i ,j k j | with q = (0 , , π ) -3 -2 -1 0 1 2 3-9-8.5-8-7.5-7-6.5-6-5.5-5-4.5-4 10 -5 (b) h ( q ) Figure 1: Illustration of discontinuities in e T a ( k i + q ) ,b ( k j − q ) i k i ,j k j and h ( q ) based on a quasi-1D modelproblem with unit cell [0 , and an anisotropic Gaussian effective potential field describedin (12). All sampled k points are of the form (0 , , k ) with k ∈ [ − π, π ]. There are two linesof discontinuities in e T a ( k i + q ) ,b ( k j − q ) i k i ,j k j : k j − k i − π = 0 and k j − k i − π = − π . -3 -2 -1 0 1 2-3-2-1012 -8.65e-05-8.63e-05-8.58e-05-8.48e-05-8.35e-05-8.19e-05-7.98e-05-7.74e-05-7.46e-05 (a) isotropoic -3 -2 -1 0 1 2-3-2-1012 -9.4e-05 -8.82e-05-8.24e-05-7.65e-05-7.07e-05-6.49e-05-5.91e-05-5.33e-05-4.75e-05 (b) anisotropic Figure 2: Illustration of discontinuities in h ( q ) from two quasi-2D model problems with unitcell [0 , , which have an isotropic and an anisotropic Gaussian effective potential fields,respectively, as described in (12). All sampled q points are of the form (0 , q y , q z ) with q y , q z ∈ [ − π, π ]. 11 .2 Staggered mesh method Based on the analysis above, the potential non-smoothness of the integrand is intrinsic to thedefinition of E TDL c and does not depend on the choice of the quadrature scheme. However,the standard method for MP2 calculations places certain quadrature nodes on points ofdiscontinuity of the integrand, which leads to finite-size errors of size O ( N − k ). We propose asimple modification of the procedure to evaluate the MP2 energy, called the staggered meshmethod . The main idea is to use an MP mesh K occ for occupied momentum vectors k i , k j ,but a different, same-sized MP mesh K vir for virtual momentum vectors k a , k b , where K vir isobtained by shifting K occ with half mesh size in all extended directions to create a staggeredmesh (see Figure 3). The MP2 energy is then computed as E staggered c = 1 N k | Ω S | X ijab X k i , k j ∈K occ X k a , k b ∈K vir h i k i , j k j | a k a , b k b i e T a k a ,b k b i k i ,j k j , (11)with N k = |K occ | = |K vir | and | Ω S | = N k | Ω | .Figure 3: Illustration of the staggered meshes K occ and K vir for a quasi-2D system.This calculation (11) can still be interpreted as a two-step numerical quadrature schemein (9) and (10), but with a different set of quadrature nodes. The induced mesh K q in (10)shifts the Γ-centered MP mesh by half mesh size (recall that K q is the set of all possibleminimum images of k a − k i with k a ∈ K vir , k i ∈ K occ ) and does not contain q = . Recall12hat in (9) for computing S q ( G ), the integrand becomes discontinuous when k j − k i − q ∈ L ∗ .In the staggered mesh method, for each q ∈ K q , all possible values of k j − k i − q (for any k i , k j ∈ K occ ) belong to K q and are always outside L ∗ . As a result, all the defined quadraturenodes in the staggered mesh method do not overlap with any points of discontinuity of theintegrand for computing S q ( G ), h ( q ), or E TDL c . This completely eliminates the first sourceof the error in Section 2.1, i.e., the improper placement of quadrature nodes. Figure 4illustrates the q -point mesh K q and the computed h ( q ) in the standard and the staggeredmesh methods for a quasi-1D model problem. We note that the staggered mesh methodavoids sampling q = for integrating h ( q ). The computed values of h ( q ) are also moreaccurate than those computed by the standard method, as it avoids sampling discontinuouspoints of the integrand in (6).The quadrature error now only comes from the second error source, i.e. lack of smoothnessin corresponding integrands. For systems with certain symmetries, this error is asymptoti-cally o ( N − k ) and hence the staggered mesh method exhibits significant advantages over thestandard method. Another special case is quasi-1D systems, where the point of discontinuityin the integrand becomes removable, and the finite-size error of the staggered mesh methodcan be much smaller. According to the discussion in Section 2.1, there are multiple factors contributing to thefinite-size errors of the MP2 correlation energy. In order to focus on the contribution from thequadrature error, we will first compare the performance of the standard and the staggeredmesh methods for MP2 calculations for a series of model problems with given effectivepotentials in Section 3.1. We then compare the performance of the two methods for periodichydrogen and silicon systems in Section 3.2, using the PySCF software package .In all the following tests, the MP mesh for virtual orbitals includes the Γ point. The13 -5 Figure 4: Illustration of h ( q ) computed by the standard and the staggered mesh methodswith mesh size 1 × ×
10 for a quasi-1D model problem with unit cell [0 , and an anisotropicGaussian effective potential field described in (12). All sampled q points are of the form(0 , , q z ) with q z ∈ [ − π, π ]. The reference curve for h ( q ) is computed based on the standardmethod with mesh size 1 × × h ( q = 0) isremovable.standard method uses the same MP mesh for occupied orbitals. The staggered mesh methodshifts the MP mesh by half mesh size for occupied orbitals. For quasi-1D, quasi-2D, and 3Dsystems, the MP meshes are of size 1 × × N k , 1 × N / k × N / k , and N / k × N / k × N / k ,respectively. To avoid finite-size errors in orbitals and orbital energies computed by mean-field methods,we consider a model system with a (possibly anisotropic) Gaussian effective potential field.More specifically, let the unit cell be [0 , , and use 14 × ×
14 planewave basis functionsto discretize functions in the unit cell. The Gaussian effective potential takes the form V ( r ) = X R ∈ L C exp (cid:18) −
12 ( r + R − r ) ⊤ Σ − ( r + R − r ) (cid:19) , (12)14ith r = (0 . , . , . k in Ω ∗ , we solve the correspondingeffective Kohn-Sham equation to obtain n occ occupied orbitals and n vir virtual orbitals. Thecovariance matrix Σ controls the isotropicity of system. For the isotropic case, we chooseΣ = diag(0 . , . , . ), C = − n occ = 1, n vir = 3. For the anisotropic case, we chooseΣ = diag(0 . , . , . ), C = − n occ = 1, n vir = 1. For such model problems, theselected n vir virtual bands are separated from the remaining virtual bands, which ensuresthat the MP2 correlation energy with a fixed number of virtual bands is a well definedproblem. There is also a direct gap between the occupied and virtual bands in all cases.We first consider the error for estimating the integrand h ( q ) in (8). For quasi-1D systems,we consider the evaluation of h ( q ) at q = (0 , , π ). This particular point is selected because h ( q ) can be directly evaluated by the standard method when N k is even, and by thestaggered mesh method when N k is odd. Similarly, for quasi-2D and 3D systems, we considerthe evaluation of h ( q ) at q = (0 , π, π ) and q = ( π, π, π ), respectively.Figure 5 demonstrates the convergence of h ( q ) with respect to N k using the standardand the staggered mesh methods. For all the systems, we find that the finite-size error ofthe staggered mesh method in estimating h ( q ) is much smaller than that of the standardmethod, regardless of the dimension or the anisotropicity of the system. This suggeststhat the structure factor S q ( G ) is evaluated more accurately as well by the staggered meshmethod. Based on the analysis in Section 2.1, it may be surprising that the staggered meshmethod significantly outperforms the standard method even in the anisotropic case. This isin fact due to the implicit symmetrization effect of the cubic cell, which reduces the errorfor computing h ( q ).Figure 6 demonstrates the convergence of the MP2 correlation energy per unit cell com-puted by the standard and the staggered mesh methods for quasi-1D, quasi-2D, and 3Dmodel systems. For quasi-1D systems, we find that the finite-size errors in the staggeredmesh method decay very rapidly with respect to N k , and the curve is nearly flat. For quasi-2D and 3D model systems, the finite-size errors of the staggered mesh method are also much15
10 15 20-8.9-8.85-8.8-8.75-8.7-8.65 10 -5 (a) Quasi-1D, isotropic -5 (b) Quasi-2D, isotropic -5 (c) 3D, isotropic -5 (d) Quasi-1D, anisotropic -5 (e) Quasi-2D, anisotropic -5 (f) 3D, anisotropic Figure 5: Estimate of h ( q ) at q / q / q using the standard and the staggered mesh methodsfor quasi-1D/quasi-2D/3D model systems with isotropic and anisotropic Gaussian effectivepotential fields. Each of these curve fittings omits the first two data points.16maller for the isotropic systems. However, for the anisotropic systems, the convergencerates of the two methods are comparable, though the error of the staggered mesh methodstill exhibits a smaller preconstant. (a) Quasi-1D, isotropic (b) Quasi-2D, isotropic (c) 3D, isotropic (d) Quasi-1D, anisotropic (e) Quasi-2D, anisotropic (f) 3D, anisotropic Figure 6: MP2 energy per unit cell computed by the standard and the staggered mesh meth-ods for quasi-1D, quasi-2D, and 3D model systems with isotropic and anisotropic Gaussianeffective potential fields. Each of these curve fittings omits the first two data points.17 .2 Hydrogen and silicon systems
We have implemented the staggered mesh method in the PySCF software package. In orderto focus on the quadrature error, we perform our comparisons as follows. For 3D systems,we first compute the orbitals and orbital energies by Hartree-Fock method with a Γ-centeredMP mesh of size 2 N / k × N / k × N / k . Then we choose two uniform sub-meshes of size N / k × N / k × N / k and their associated orbitals and orbital energies to compute the MP2energies for the staggered mesh method. We remark that this procedure is only employedto simplify our preliminary implementation in PySCF, and is not necessary in a productionlevel computation. For the standard method, only one such sub-mesh is used. The setup ofthe tests for quasi-1D and quasi-2D systems is similar. Basis set gth-szv is used in all thetests below.We consider two sets of periodic systems: hydrogen dimer, and silicon. The hydrogendimer is placed at the center of a cubic unit cell of edge length 6 Bohr pointing in the x -direction and has separating distance 1 . O ( N − k ) finite-size error fromthose in HF orbital energies. Thus, the fitting for the staggered mesh method still uses thepower law scaling of N − k . For quasi-2D and 3D cases, the staggered mesh method performssignificant better than the standard one for silicon. This is likely due to the symmetryof the system that leads to improved smoothness in the integrand and thus reduces thecorresponding quadrature error. In comparison, the performance of the two methods becomesimilar for the quasi-2D and 3D hydrogen dimer systems due to the anisotropicity.18 (a) Quasi-1D (b) Quasi-2D (c) 3D Figure 7: MP2 energy per unit cell computed by the standard and the staggered meshmethods for the periodic hydrogen dimer systems. (a) Quasi-1D (b) Quasi-2D (c) 3D
Figure 8: MP2 energy per unit cell computed by the standard and the staggered meshmethods for silicon crystal systems. 19
Conclusion
The convergence of the MP2 correlation energy towards the TDL is a fundamental question inmaterials science. Existing analysis in the literature focuses on the missing contribution of thestructure factor S q ( G ) at q + G = , but neglects contributions from 1) certain quadraturenodes coincide with points of discontinuity of the integrand 2) the quadrature error due tothe intrinsic non-smoothness of the integrand. We demonstrate that such contributions canbe at least equally important and scale as O ( N − k ). We propose the staggered mesh methodthat uses a different set of quadrature nodes for the trapezoidal quadrature, which allowsus to completely avoid the first source of the error. Numerical evidence shows that thestaggered mesh method is particularly advantageous over the standard method for quasi-1D systems and systems with symmetries, which reduces the contribution from the seconderror source. We expect that the new approach can also be useful for correlation energycalculations beyond the MP2 level, such as higher levels of perturbation theories and coupleclustered theories. Acknowledgement
This work was partially supported by the Air Force Office of Scientific Research under awardnumber FA9550-18-1-0095 (X.X.,L.L.), by the Department of Energy under Grant No. DE-SC0017867, and by the National Science Foundation under Grant No. DMS-1652330 (L.L.),and by the China Scholarship Council under File No. 201906040071. (X.L).
References (1) Marsman, M.; Gr¨uneis, A.; Paier, J.; Kresse, G. Second-order Møller–Plesset pertur-bation theory applied to extended systems. I. Within the projector-augmented-waveformalism using a plane wave basis set.
The Journal of Chemical Physics , ,184103. 20igure 9: TOC(2) Gr¨uneis, A.; Marsman, M.; Kresse, G. Second-order Møller–Plesset perturbation theoryapplied to extended systems. II. Structural and energetic properties. The Journal ofChemical Physics , , 074107.(3) M¨uller, C.; Paulus, B. Wavefunction-based electron correlation methods for solids. Phys.Chem. Chem. Phys. , , 7605–7614.(4) Sch¨afer, T.; Ramberger, B.; Kresse, G. Quartic scaling MP2 for solids: A highly paral-lelized algorithm in the plane wave basis. J. Chem. Phys. , , 104101.(5) McClain, J.; Sun, Q.; Chan, G. K. L.; Berkelbach, T. C. Gaussian-based coupled-clustertheory for the ground-state and band structure of solids. J. Chem. Theory Comput. , , 1209–1218.(6) Gruber, T.; Liao, K.; Tsatsoulis, T.; Hummel, F.; Gr¨uneis, A. Applying the coupled-cluster ansatz to solids and surfaces in the thermodynamic limit. Phys. Rev. X , , 021043.(7) Booth, G. H.; Gr¨uneis, A.; Kresse, G.; Alavi, A. Towards an exact description ofelectronic wavefunctions in real solids. Nature , , 365–370.(8) Liao, K.; Gr¨uneis, A. Communication: Finite size correction in periodic coupled clustertheory calculations of solids. J. Chem. Phys. , , 141102.219) Mihm, T. N.; McIsaac, A. R.; Shepherd, J. J. An optimized twist angle to find thetwist-averaged correlation energy applied to the uniform electron gas. The Journal ofChemical Physics , , 191101.(10) Mihm, T. N.; Yang, B.; Shepherd, J. J. Power laws used to extrapolate the coupledcluster correlation energy to the thermodynamic limit. arXiv, 2007.11696 ,(11) Chiesa, S.; Ceperley, D. M.; Martin, R. M.; Holzmann, M. Finite-size error in many-body simulations with long-range interactions. Phys. Rev. Lett. , , 6–9.(12) Lin, C.; Zong, F.; Ceperley, D. M. Twist-averaged boundary conditions in continuumquantum Monte Carlo algorithms. Physical Review E , , 016702.(13) Broqvist, P.; Alkauskas, A.; Pasquarello, A. Hybrid-functional calculations with plane-wave basis sets: Effect of singularity correction on total energies, energy eigenvalues,and defect energy levels. Physical Review B , , 085114.(14) Shepherd, J. J.; Henderson, T. M.; Scuseria, G. E. Coupled cluster channels in thehomogeneous electron gas. The Journal of Chemical Physics , , 124102.(15) Gell-Mann, M.; Brueckner, K. A. Correlation energy of an electron gas at high density. Physical Review , , 364.(16) H¨attig, C.; Klopper, W.; K¨ohn, A.; Tew, D. P. Explicitly correlated electrons inmolecules. Chemical Reviews , , 4–74.(17) Gr¨uneis, A.; Shepherd, J. J.; Alavi, A.; Tew, D. P.; Booth, G. H. Explicitly corre-lated plane waves: Accelerating convergence in periodic wavefunction expansions. TheJournal of Chemical Physics , , 084112.(18) Sun, Q.; Berkelbach, T. C.; Blunt, N. S.; Booth, G. H.; Guo, S.; Li, Z.; Liu, J.;McClain, J. D.; Sayfutyarova, E. R.; Sharma, S., et al. PySCF: the Python-based22imulations of chemistry framework. Wiley Interdisciplinary Reviews: ComputationalMolecular Science , , e1340.(19) Lyness, J. An error functional expansion for N -dimensional quadrature with an inte-grand function singular at a point. Mathematics of Computation , , 1–23.(20) Drummond, N. D.; Needs, R. J.; Sorouri, A.; Foulkes, W. M. C. Finite-size errors incontinuum quantum Monte Carlo calculations. Phys. Rev. B ,78