Orbital-enriched Flat-top Partition of Unity Method for the Schrödinger Eigenproblem
Clelia Albrecht, Constanze Klaar, John E. Pask, Marc Alexander Schweitzer, N. Sukumar, Albert Ziegenhagel
OOrbital-enriched flat-top partition of unity method for theSchr ¨odinger eigenproblem
Clelia Albrecht a , Constanze Klaar a , John Ernest Pask b , Marc Alexander Schweitzer a,c, ∗ ,N. Sukumar d, ∗ , Albert Ziegenhagel a a Fraunhofer Institute for Algorithms and Scientific Computing SCAI, Schloss Birlinghoven, 53757 Sankt Augustin,Germany b Physics Division, Lawrence Livermore National Laboratory, 7000 East Avenue, Livermore, CA 94550, U.S.A. c Institute for Numerical Simulation, Rheinische Friedrich-Wilhelms-Universit¨at, Wegelerstr. 6, 53115 Bonn, Germany d Department of Civil and Environmental Engineering, University of California, Davis, CA 95616, U.S.A.
Abstract
Quantum mechanical calculations require the repeated solution of a Schr¨odinger equation for thewavefunctions of the system, from which materials properties follow. Recent work has shownthe e ff ectiveness of enriched finite element type Galerkin methods at significantly reducing thedegrees of freedom required to obtain accurate solutions. However, time to solution has been ad-versely a ff ected by the need to solve a generalized rather than standard eigenvalue problem and theill-conditioning of associated systems matrices. In this work, we address both issues by proposinga stable and e ffi cient orbital-enriched partition of unity method to solve the Schr¨odinger boundary-value problem in a parallelepiped unit cell subject to Bloch-periodic boundary conditions. Inthe proposed partition of unity method, the three-dimensional domain is covered by overlappingpatches, with a compactly-supported weight function associated with each patch. A key ingredientin our approach is the use of non-negative weight functions that possess the flat-top property, i.e.,each weight function is identically equal to unity over some finite subset of its support. This flat-top property provides a pathway to devise a stable approximation over the whole domain. On eachpatch, we use p -th degree orthogonal (Legendre) polynomials that ensure p -th order completeness,and in addition include eigenfunctions of the radial Schr¨odinger equation. Furthermore, we adopta variational lumping approach to construct a (block-)diagonal overlap matrix that yields a stan-dard eigenvalue problem for which there exist e ffi cient eigensolvers. The accuracy, stability, ande ffi ciency of the proposed method is demonstrated for the Schr ¨odinger equation with a harmonicpotential as well as a localized Gaussian potential. We show that the proposed approach deliversoptimal rates of convergence in the energy, and the use of orbital enrichment significantly reducesthe number of degrees of freedom for a given desired accuracy in the energy eigenvalues while thestability of the enriched approach is fully maintained. Keywords: quantum mechanics, partition of unity method, Bloch boundary conditions,variational mass lumping, enrichment functions, stability ∗ Corresponding authors
Email addresses: [email protected] (Marc Alexander Schweitzer), [email protected] (N. Sukumar)
Preprint submitted to Journal Name March 2, 2018 a r X i v : . [ phy s i c s . c o m p - ph ] M a r . Introduction The Kohn–Sham (KS) equations of density functional theory (DFT) are the dominant theoreti-cal formulation in quantum mechanical simulations of condensed matter (solids and liquids). TheKS equations require the repeated solution of the steady-state Schr¨odinger and Poisson equationson a parallelepiped unit cell with Bloch-boundary conditions [1]. The solution of the Schr¨odingerequation is the most time-consuming part in KS-DFT calculations. The current state-of-the-artapproach to solve the equations of KS-DFT is the planewave (PW) pseudopotential method thatuses a Fourier basis set, and requires the solution of a discrete standard eigenproblem. It has beenappreciated in recent years that enriched Galerkin methods [2–13] can be very competitive withPW methods in attaining the desired accuracy with comparable or far fewer degrees of freedom(basis functions). While early formulations [2, 6, 7, 10] employed direct enrichment, more recentapproaches have employed discontinuous Galerkin [3, 4, 11] or partition of unity finite elementformulations [2, 5, 9, 13] in order to strictly localize orbital enrichments, thus facilitating flexibleapproximation and e ffi cient parallel implementation. In [5] it was shown that an enriched partitionof unity finite element method (PUFEM) [14, 15] requires an order of magnitude fewer basis func-tions than current state-of-the-art PW based methods to attain the desired 1 mHa / atom accuracyin total energy calculations. However, the ill-conditioning of the resulting system matrices andthe need to solve a generalized rather than standard eigenvalue problem were key issues identi-fied as adversely a ff ecting time to solution in practice. In this work, we use a flat-top partitionof unity method (PUM) [16–18] to address these issues in the approximation of the Schr¨odingerequation. Our flat-top PUM produces well-conditioned system (Hamiltonian and overlap) matri-ces and yields a standard eigenvalue problem via variational lumping. The approximation qualityof our flat-top PUM is comparable to that reported in [2], but it overcomes the two main short-comings of the PUFEM that arise in the solution of the Schr¨odinger eigenproblem. In addition tothe electronic Schr¨odinger equation, the flat-top PUM with Bloch boundary conditions also holdspromise in areas such as acoustic scattering [19], elastodynamics [20] and electromagnetics [21],where large-scale eigenproblems are solved and useful a priori information is available.In condensed matter calculations, the Schr¨odinger equation is solved in a unit cell (paral-lelepiped domain Ω ) subject to Bloch-periodic boundary conditions (see Figure 1). In the flat-toppartition of unity method, the domain Ω is covered by overlapping patches (see Figure 2) andeach patch i is associated with a weight function ϕ i ( x ) with support ω i such that (cid:80) i ϕ i ( x ) = ϕ i ( x ) ≡ ω FT i ⊂ ω i (see Figure 3). The local basis set V i on each patch consists of polynomialsand / or non-polynomial (orbital enrichment) functions, and the global approximation is formed bylinear combinations of the products of ϕ i ( x ) and functions from V i . We perform local orthogonal-ization to ensure that all functions on a patch are linearly independent and thereby attain globalstability [17], and adopt the variational lumping scheme [18] to realize a standard eigenproblem.The remainder of the paper is organized as follows. In the next section, we state the strongand weak forms of the Schr¨odinger eigenproblem. In Section 3, we introduce the partition ofunity method, where we present the proposed flat-top PUM in Section 3.1. The key steps inthe local orthogonalization procedure to construct a stable global approximation are discussedin Section 3.2, and we describe the variational lumping scheme in Section 3.3. Numerical examplesfor the Schr ¨odinger equation are presented in Section 4, where we show that the system matricesare well-conditioned and that the use of orbital-enrichment provides a very e ffi cient solution vis-`a-vis solely using polynomials over each patch. In addition, we also provide comparisons in the2igenspectrum when using the consistent overlap matrix versus the lumped overlap matrix, andthe results reveal that the variational lumping scheme does not adversely a ff ect the accuracy of theenergy eigenvalues. We close with a few concluding remarks in Section 5.
2. The Schr¨odinger eigenproblem
The stationary Schr¨odinger equation reads as H ψ ( x ) (cid:66) − ∇ ψ ( x ) + V e ff ( x ) ψ ( x ) = εψ ( x ) in Ω . (1)We consider a parallelepiped unit-cell Ω ⊂ R with primitive lattice vectors a d ( d = , ,
3) todescribe a periodic condensed matter system (see Figure 1). Thus, the e ff ective potential and thecharge density are periodic, i.e., they satisfy V e ff ( x + R ) = V e ff ( x ) , ρ ( x + R ) = ρ ( x ) , (2)whereas the solution of Schr¨odinger’s equation ψ , the so-called wavefunction, satisfies Bloch’stheorem ψ ( x + R ) = ψ ( x ) exp( i k · R ) (3)for any lattice translation vector R = n a + n a + n a with n d ∈ Z ( d = , ,
3) and wavevector k .Note that for k = Γ -point) the wavefunction is also periodic whereas for all other wavevectorsthere is a phase-shift exp( i k · R ) associated with a translation by R . In this situation the generalproblem (1) becomes − ∇ ψ ( x ) + V e ff ( x ) ψ ( x ) = εψ ( x ) in Ω ,ψ ( x + a d ) = exp( i k · a d ) ψ ( x ) on Γ d , ∇ ψ ( x + a d ) · ˆ n ( x ) = exp( i k · a d ) ∇ ψ ( x ) · ˆ n ( x ) on Γ d , (4)where ( ψ, ε ) denotes an eigenpair consisting of the respective wavefunction ψ and its associatedenergy ε , ˆ n ( x ) is the outward unit normal at x and Γ d are the bounding faces of the domain Ω ,see Figure 1. Even though the boundary conditions and thus the wavefunction are complex-valued,the eigenvalues ε ∈ R due to the self-adjointness of the Hamiltonian H .In this paper, we are concerned with the numerical approximation of (4) by Galerkin methodsand thus need to rewrite (4) in its respective weak form. To this end, we consider the functionspace V (cid:66) { v ∈ H ( Ω , C ) : v ( x + a d ) = v ( x ) exp( i k · a d ) on Γ d , d = , , } (5)and test (4) with v ∈ V to attain a ( v , ψ ) = ε (cid:104) v , ψ (cid:105) L ( Ω , C ) for all v ∈ V (6)with a ( v , ψ ) (cid:66) (cid:90) Ω (cid:16) ∇ v ( x ) ∇ ψ ( x ) + v ( x ) V e ff ( x ) ψ ( x ) (cid:17) d x (7)and (cid:104) v , ψ (cid:105) L ( Ω , C ) (cid:66) (cid:90) Ω v ( x ) ψ ( x ) d x (8)3 a a Γ a Γ Figure 1: Sketch of a parallelepiped unit cell Ω spanned by primitive lattice vectors a d with respective boundarysegments Γ d ( d = , , after integration by parts [2]. Thus, choosing a finite-dimensional space V M ⊂ V with basisfunctions φ i to discretize (6) we obtain the generalized eigenproblem H ˜ ψ = ε S ˜ ψ, (9)where H = ( H i j ) ∈ C M × M denotes the discrete Hamiltonian and S = ( S i j ) ∈ C M × M is the so-calledoverlap (or consistent mass) matrix H i j (cid:66) a ( φ j , φ i ) , and S i j (cid:66) (cid:104) φ j , φ i (cid:105) L ( Ω , C ) . (10)The approximate eigenfunction ψ ∈ V M is given by ψ ( x ) (cid:66) M (cid:88) i = ψ i φ i ( x ) , (11)where ˜ ψ = ( ψ i ) ∈ C M denotes the associated coe ffi cient vector. Throughout this paper we employa particular partition of unity method [16, 22, 23] to construct the respective finite-dimensionalspace V M and its basis functions φ i .
3. Partition of unity method
The partition of unity method (PUM) is a generalization of the finite element method (FEM)that is typically employed for the spatial discretization of a partial di ff erential equation (PDE), seefor example [24, 25]. The notion of a PUM was coined in [14, 15] and is based on the specialFEM developed in [26]. The abstract ingredients of a PUM are a partition of unity (PU) { ϕ i : i = , . . . , N } and a collection of local approximation spaces V i ( ω i ) (cid:66) span (cid:104) ϑ mi (cid:105) dim( V i ) m = defined on thepatches ω i (cid:66) supp( ϕ i ) for i = , . . . , N which overlap and whose union covers the computationaldomain Ω ⊂ R d . With these two ingredients we define the PUM space V PU (cid:66) N (cid:88) i = ϕ i V i = span (cid:104) ϕ i ϑ mi (cid:105) ; (12)that is, the basis functions of a PUM space are simply defined as the products of the PU functions ϕ i and the local approximation functions ϑ mi . The PU functions provide the locality and global4 igure 2: Schematic of a sequence of uniformly refined covers which come from the scaling of uniform grid cells intwo dimensions. Depicted is a single cover patch ω i = (cid:81) dl = ( o li − α h , o li + α h ) (gray) with 2 h = / , / , / (left to right)and its center o i . regularity of the product functions whereas the functions ϑ mi equip V PU with its approximationpower. Note that the local approximation spaces V i can be chosen in a problem-dependent fashionand are independent of each other, i.e., any local basis ϑ mi may be employed on any patch ω i . Ingeneral the local approximation space V i (cid:66) span (cid:104) ϑ ni (cid:105) associated with a particular patch ω i of aPUM space V PU consists of two parts: A smooth approximation space, for example polynomials V P i ( ω i ) (cid:66) span (cid:104) π si (cid:105) , and a problem-dependent enrichment part V E i ( ω i ) (cid:66) span (cid:104) ψ ti (cid:105) , i.e.,span (cid:104) ϑ ni (cid:105) = V i ( ω i ) = V P i ( ω i ) + V E i ( ω i ) = span (cid:104) π si , ψ ti (cid:105) . (13)For the ease of notation, we make the following conventions. For an arbitrary function u ∈ V PU with the representation u ( x ) = N (cid:88) i = V i ) (cid:88) m = u mi ϕ i ( x ) ϑ mi ( x ) = N (cid:88) i = ϕ i ( x ) dim( V i ) (cid:88) m = u mi ϑ mi ( x ) = : N (cid:88) i = ϕ i ( x ) u i ( x )we denote its associated overall coe ffi cient vector by˜ u = ( u ( i , m ) ) ∈ R dim( V PU ) with dim( V PU ) (cid:66) N (cid:88) i = dim( V i ) . The main di ff erence between our PUM approach [22, 27] and most other generalized or ex-tended FEM techniques, compare [24, 25], is that we employ a so-called flat-top PU instead ofclassical linear Lagrange FE functions, compare Figure 3. For the construction of such a flat-topPU let us first define a cover C Ω (cid:66) { ω i } of the domain Ω with the help of a uniform regular meshof mesh-width 2 h by an isotropic scaling of the mesh-cells C i = d (cid:89) l = ( o li − h , o li + h ) , .
25 0 . .
75 101 0 0 .
25 0 . .
75 101
Figure 3: Partition of unity comprised of linear FE functions (left) and a flat-top PU (right) in one dimension. .
33 0 .
67 101 ω ω ω .
33 0 .
67 101 ω ω , T ω ω ω , T Figure 4: Comparison of one-dimensional flat-top PU weight functions for Dirichlet or Neumann boundary conditions(left) and periodic boundary conditions (right). To realize periodic weight functions we copy patches ω i with ω i ∩ ∂ Ω (cid:44) ∅ periodically (see also Figure 5 and [23]). T ( ω i ) 1 2 3 45 6 7 8910 1011 1213 14 15 161 2 3 45 6 7 8910 11 1213 14 15 16 T , ( ω i ) T , ( ω i ) T , ( ω i ) Figure 5: Schematic of the implementation of the periodicity of a PU function ϕ i (16) associated with a patch ω i (red)by generalization of the notion of neighboring or overlapping patches ω i ∩ ω j (cid:44) ∅ (green) at the boundary ∂ Ω viaperiodic copies of patches ω i ∩ ∂ Ω (cid:44) ∅ (red), see [23] for details. ω i as ω i (cid:66) d (cid:89) l = ( o li − α h , o li + α h ) , with α > , (14)see Figure 2. To obtain a PU on a cover C Ω with N (cid:66) card( C Ω ) we define a weight function W i : Ω → R with supp( W i ) = ω i for each cover patch ω i by W i ( x ) = (cid:40) W ◦ T i ( x ) , x ∈ ω i , otherwise (15)with the a ffi ne transforms T i : ω i → [ − , d and W : [ − , d → R any non-negative compactlysupported function, such as a B-spline. By normalizing these weight functions we obtain thefunctions ϕ i ( x ) (cid:66) W i ( x ) (cid:80) l ∈ C i W l ( x ) , (16)where we define the local neighborhood C i (cid:66) { l : ω l ∩ ω i (cid:44) ∅} of a patch ω i in a slightly moregeneral way to account for the required periodicity of the basis functions, see Figures 4 and 5.To impose Bloch-periodic boundary conditions, we multiply the entries of the resulting systemmatrices in matrix blocks that correspond to Bloch-boundary patches with the appropriate Blochphase factor. This results in a periodic PU and Bloch-periodic basis, analogous to the procedurein [2]. Note that the PU (16) is non-negative since the employed weight functions are non-negativeand that the ϕ i satisfy the flat-top property for any α ∈ (1 , α ∈ (1 , α = α = The approximation power of an enriched PUM is mostly obtained by the choice of high qualityenrichment space V E (enrichment functions ψ ti ) on the patches ω i for the problem at hand. Forinstance, the use of generalized harmonic polynomials or local spaces based on planewaves havebeen employed successfully for smooth problems in [14, 30, 31], whereas for crack propagationproblems the use of discontinuous and singular enrichment functions is appropriate, see e.g. [25,32, 33]. In the context of the Schr¨odinger and Poisson equations of Kohn–Sham density functionaltheory, enrichment functions ψ ti constructed from isolated-atom solutions were employed in [5,9, 10], as in previous work [2] in the context of model densities and potentials. Our generalflat-top PUM allows for the use of arbitrary problem-dependent enrichment functions and thuswe anticipate that our PUM has essentially the same convergence properties as for example thePUFEM of [2, 5, 10] when we employ the same enrichment functions.The fundamental di ff erence of our approach from the PUFEM of [2, 5, 10] is that we employa flat-top PU [17, 27, 34] to overcome the two major challenges encountered in the PUFEM: ill-conditioning of the overlap matrix and the need for the solution of a generalized eigenproblem.Let us first consider the ill-conditioning of the overlap matrix, i.e., the L -stability of (12),which can be encountered in all enriched approximations [24, 25, 35–37]. For the smooth space V P i
7n (13) we employ a local basis π si on ω i , i.e., π si = p s ◦ T i and { p s } denotes a stable basis on [ − , d ,for instance Legendre-polynomials. The enrichment functions ψ ti however are often given as globalfunctions η t on the computational domain Ω since they are designed to capture special behavior ofthe global solution at a particular location. Therefore, the restrictions ψ ti (cid:66) η t | ω i of the enrichmentfunctions η t to a particular patch ω i may be ill-conditioned or even linearly dependent on ω i , even ifthe enrichment functions η t are well-conditioned on the global domain. Furthermore, the couplingbetween the spaces V P i and V E i on the patch ω i must be considered. The set of functions { π si , ψ ti } will degenerate from being a basis of V i to yielding a generating system if the restricted enrichmentfunctions ψ ti = η t | ω i can be well-approximated by polynomials π si on the patch ω i . Moreover, inthe general case when (12) employs an arbitrary non-flat-top PU, e.g., a PU built from FE, wealso need to consider the interactions of the overlapping local approximation spaces V i and V j (13)on neighboring patches ω i ∩ ω j (cid:44) ∅ which essentially introduces global constraints and thereforerenders the use of an arbitrary non-flat-top PU computationally infeasible. If we, however, restrictourselves to the use of a flat-top PU, all issues that may result in ill-conditioning can be e ffi cientlyresolved solely on the local patches ω i , see [17] for further details.To attain a stable basis of V PU for arbitrary local approximation spaces V i with a flat-top PU { ϕ i } , we essentially need to be able to identify a local enrichment space V D i such that we canrewrite (13) as V i = V P i + V E i = V P i ⊕ V D i . To this end, we choose an inner product on the patch ω i and construct V D i such that V P i ⊥ V D i holds. The respective stability transformation P : V PU = N (cid:88) i = ϕ i ( V P i + V E i ) → N (cid:88) i = ϕ i ( V P i ⊕ V D i ) (17)can in fact be computed e ffi ciently on-the-fly by partial orthogonalization with respect to the cho-sen inner product and only involves local operations on the patches [17]. The particular stabilitytransformation on a patch ω i depends on the choice of the inner product and works in four stepsstarting from a local matrix M i = (cid:32) M i ; P , P M i ; P , E M i ; E , P M i ; E , E (cid:33) arising from Gram matrix of all basis functions π si of V P i and ψ ti of V E i in the chosen inner prod-uct. First, all basis functions of V P i and V E i are scaled appropriately with respect to the employedinner product, so that the induced norm of all functions is identical. Then, we compute eigen-value decompositions of M i ; P , P and M i ; E , E individually and use thresholding of small eigenvaluesto eliminate any instabilities within V P i and V E i , so that we can transform M i to attain the form M (1) i = (cid:32) I P M (1) i ; P , E M (1) i ; E , P I E (cid:33) . Then, the overlap of the two spaces V P i and V E i is removed by partial orthogonalization (Schurcomplement) which yields M (2) i = (cid:32) I P M (2) i ; D , D (cid:33) . M (2) i ; D , D we obtain an orthonormal basis of V i with respect to the employed inner product. And so the problem of ill-conditioning which isencountered in PUFEM can be resolved in a flat-top PUM while fully maintaining the improvedapproximation quality due to enrichment. The second drawback of the PUFEM identified in [5] is related to the fact that the Galerkindiscretization of (4) yields a Hermitian generalized eigenproblem (9) where H and S are in gen-eral sparse but not diagonal matrices, which renders the solution of (9) more challenging since thecomputation of S − is implicitly required, i.e., one e ff ectively needs to solve the standard eigen-problem S − H ˜ u = ε ˜ u . In the context of classical (linear) FEM one therefore often replaces the overlap (or consistentmass) matrix S by an approximation ¯ S whose inverse ¯ S − can be computed much more e ffi ciently,e.g., if ¯ S is diagonal, by so-called mass-lumping [38–40]. Various constructions for ¯ S exist in theFEM context [41], yet these approaches are in general not directly applicable to a PUM or anyenriched approximation scheme due to the non-polynomial and non-interpolatory character of theemployed basis functions.Fortunately, there is a natural and rather simple approach to the construction of an appropriateapproximation of the overlap matrix for the PUM in general [18]. It is applicable to any non-negative PU { ϕ i } and arbitrary local approximation spaces (13) (e.g., higher order polynomials,discontinuous, and singular enrichments) and utilizes only the special product structure of thePUM basis functions ϕ i ϑ ni .To introduce this variational lumping scheme for our PUM, recall that the consistent overlapmatrix is given by S = ( S ( i , n ) , ( j , m ) ) , S ( i , n ) , ( j , m ) = (cid:104) ϕ j ϑ mj , ϕ i ϑ ni (cid:105) L ( Ω , C ) , (18)where our flat-top PU functions ϕ i and ϕ j are non-negative. Thus, the entries of the overlap matrixare given by the inner products in L ( Ω , C ) of our PUM basis functions. In [18] it was shownthat the global L ( Ω , C ) inner product in (18) can be replaced by local weighted L ( ω i , C ; ϕ i ) innerproducts without diminishing the global convergence behavior, since there holds the equivalence (cid:104) f , ϕ i ϑ ni (cid:105) L ( Ω , C ) = (cid:90) Ω f ϕ i ϑ ni dx = (cid:90) Ω ∩ ω i f ϕ i ϑ ni dx = : (cid:104) f , ϑ ni (cid:105) L ( Ω ∩ ω i , C ; ϕ i ) for an arbitrary function f ∈ L ( Ω , C ). With the help of these local inner products, we thereforeobtain an approximate overlap matrix¯ S = ( ¯ S ( i , n ) , ( j , m ) ) , ¯ S ( i , n ) , ( j , m ) = (cid:40) i (cid:44) j (cid:104) ϑ mi , ϑ ni (cid:105) L ( Ω ∩ ω i , C ; ϕ i ) i = j (19)which is block-diagonal and symmetric positive definite for any choice of the local approximationspaces V i , see [18] for further details. Therefore, an approximate solution of (9) can be obtainedvery e ffi ciently via H ˜ u = ε ¯ S ˜ u . (20)In fact, choosing the local weighted L ( ω i , C ; ϕ i ) inner products also in the stability transforma-tion (17) yields ¯ S = I , so that (20) becomes H ˜ u = ε ˜ u and the expensive solution of the generalizedeigenproblem encountered in the PUFEM can be completely circumvented in our flat-top PUM.9 . Numerical results In this section we present some numerical results obtained with our flat-top PUM. We con-sider two benchmark problems defined on three-dimensional unit cells [2] to validate our imple-mentation and assess the accuracy and e ffi ciency in terms of the degrees of freedom (DOFs) ofthe flat-top PUM with respect to the Schr¨odinger eigenproblem with periodic and Bloch-periodicboundary conditions. The two main objects of interest of our study, however, are the conditioningof the consistent and lumped overlap matrices (before and after the stabilization (17)) and the e ff ectour lumping scheme (19) has on the quality of the approximation. Since absolute runtime is notthe focus of this study, we employed a simple tensor-product 6 × × As in [2], the first benchmark problem we consider is the Schr ¨odinger equation (4) with har-monic potential V ( | x − τ | ) = | x − τ | , (21)under periodic boundary conditions ( k = Ω , we take a cuboid cell with primitive lattice vectors a (cid:66) a (1 , , a (cid:66) a (0 , . , a (cid:66) a (0 , , . a = τ (cid:66) a + a + a .On each patch, the polynomial approximation space V P i contains Legendre polynomials up tototal degree p = , ,
3, while the enrichment space V E i consists of periodic lattice sums of theinfinite box eigenfunctions ψ (cid:48) nlm ( x ) = R nl ( r ) Y lm ( θ, φ ), ψ nlm = (cid:88) R ψ (cid:48) nlm ( | x − τ − R | ) , (22)over lattice translation vectors R = i a + i a + i a , to enrich each of the ten lowest states. Here, Y lm are the analytically given spherical harmonics and R nl is the radial part of the eigenfunction,which is obtained by solving the radial Schr¨odinger equation [2]. Thus, the radial componentsare given numerically, i.e., at discrete points. To ensure a compact support of R nl we multiply itsdiscrete values by a C cut-o ff function [2] h ( r , r ) (cid:66) + r r − r r + r r − r r , r ≤ r , r > r , (23)where we choose r = × ×
64 mesh whichis accurate to 7 digits (and yields λ ref1 = . (cid:80) i = λ ref i = . α = .
1. Here, we anticipate that the flat-topPUM converges with rates that are comparable to classical finite element methods of the respective10 τ τ
Figure 6: Sketch of refinement by increasing the enrichment radius r e in two dimensions. We only enrich a patch ifits center point lies within the ball of radius r e centered at the potential center τ . Here, this ball and enriched patchesare colored in red. Depicted are three di ff erent values for r e : 0 . . . × order. Then, we use a fixed uniform cover and use enriched local approximation spaces, where weincrease the enrichment support radius from 0 to 4 au (see Figure 6). In the first case we do notuse any enrichments at all, in the latter we enrich each patch whose center lies within the supportradius of the enrichment center τ . The results are shown in Figure 7. From the depicted plots wecan see that the convergence behavior of the non-enriched method is just as expected, while it isalso evident that using enrichment functions is very beneficial in terms of accuracy: significantlyless than 500 DOFs (e.g., 320 DOFs with linear polynomials on a 4 × × × × − Ha, which is a typicalrequirement in quantum mechanical calculations. Apart from reducing DOFs, enrichment alsoreduces the overall wall clock time of the method: in our experiments, the computation (includingsetup, assembling of matrices and solving the EVP) on a 3 × × .
03 seconds on 8 cores, while the non-enriched computation with quadratic polynomials ona 16 × ×
16 cover took 192 .
67 seconds on 8 cores. Both methods yield about the same absoluteerror for the lowest eigenvalue λ (see Figure 7), in both cases less than 10 − Ha. Note that there isa distinct di ff erence in the PUM functions ϕ i ϑ mi on even (e.g., 4 × ×
4) and odd (e.g., 3 × × τ is contained in the overlap of eight neighboring patches. Thus, theresults summarized in Figure 7 show that our approach is also fully robust with respect to therelative position of τ .To properly evaluate the benefits of the lumped method, we first conduct the same experimentas above, only this time using a lumped overlap matrix and summarize the results in Figure 9.We observe that the convergence behavior is practically identical to that of the consistent method,however, with a slightly larger absolute error. Apart from being almost as accurate as the con-sistent method, using the lumped overlap matrix has its own benefits: First of all we can reducethe generalized eigenproblem (9) to a standard eigenproblem as described in Section 3.3, whichsaves a lot of computational time, especially for larger problems (see Figure 8). In the context ofself-consistent Kohn–Sham calculations this is especially important, as it enables the use of highlye ffi cient algorithms to refine eigenvectors in each self-consistent-field iteration [11, 42]. Compu-tational time for the overall method is also again saved by using enrichment functions: while thecomputation on a 16 × ×
16 cover using quadratic polynomials without enrichments took 151 . λ was achieved by110 − − − − − − − degrees of freedom a b s o l u t ee rr o r( H a ) p = p = p = p = + e ( ↑ r e )( p = + e ( ↑ r e )( p = + e ( ↑ r e ) − − − − − − − degrees of freedom a b s o l u t ee rr o r( H a ) p = p = p = p = + e ( ↑ r e )( p = + e ( ↑ r e )( p = + e ( ↑ r e ) Figure 7: Convergence history of the lowest eigenvalue λ for the harmonic oscillator potential (21) attained fordi ff erent refinement schemes. We consider a purely polynomial approximation ( p = , ,
3) on a sequence of uniformlyrefined covers, and a refinement by increasing the enrichment radius (compare Figure 6) with a single enrichmentfunction on a fixed uniform cover (top: 3 × ×
3, 7 × ×
7; bottom: 4 × ×
4, 8 × ×
8) that is labeled by p = , , + e ↑ r e . All results were obtained with the consistent overlap matrix. − − degrees of freedom ti m e standard eigenproblem solvergeneralized eigenproblem solver degrees of freedom ti m e standard eigenproblem solvergeneralized eigenproblem solver Figure 8: Wall clock time measured for the solution of the generalized and standard eigenproblems for the harmonicoscillator potential (21), using the default SLEPc eigenvalue solver, polynomials of degree p = , , n × n × n cover. All these data points are in one plot, sorted by number of DOFs (left: n = , , , , , n = , , ,
32 on 64 cores).Table 1: Measured condition numbers with and without stabilization (17) with respect to the local weighted L innerproduct L ( ω i , C ; ϕ i ) for the consistent and the lumped overlap matrix obtained with p = × × r e and ten enrichment functions per enriched patch. consistent overlap matrix lumped overlap matrix r e D ofs without (17) with (17) without (17) with (17)0 . ,
860 4 · · · · . ,
870 2 · · · · . ,
930 4 · · · · . ,
150 5 · · · · . ,
590 6 · · · · . ,
270 1 · · · · . ,
090 2 · · · · . ,
890 6 · · · · . ,
210 3 · · · · using a 3 × × .
87 seconds on 8 cores.Furthermore, we can overcome the problem of large condition numbers that arises in thePUFEM [5] using the stable transformation introduced in Section 3.2. The measured conditionnumbers for both consistent and lumped overlap matrices are displayed before and after applyingthe stable transformation (17) in Table 1. Here we see that for both types of overlap matrices,the condition number deteriorates quickly with increasing enrichment support radius without ourstabilization (17). The same observation was made in [5] for the PUFEM. However, after applyingthe stable transformation, we can dramatically reduce the condition number to the range of classi-cal FE overlap matrices for the consistent case. Combining the stabilization (17) with our lumpingapproach (19) yields the optimal behavior ¯ S = I so that (20) becomes a well-conditioned stan-dard eigenproblem. Thus, our flat-top PUM overcomes the two main drawbacks of the PUFEMand enables the stable use of enrichment functions on every patch and arbitrary potentials so thathigh-fidelity approximations can be attained with extremely small numbers of DOFs in a robustfashion. For example, the PUFEM employed in [2] yields an accuracy of 10 − Ha for this examplewith 948 DOFs using enriched cubic FE where the condition number of the overlap matrix is of the130 − − − − − − − degrees of freedom a b s o l u t ee rr o r( H a ) p = p = p = p = + e ( ↑ r e )( p = + e ( ↑ r e )( p = + e ( ↑ r e ) − − − − − − − degrees of freedom a b s o l u t ee rr o r( H a ) p = p = p = p = + e ( ↑ r e )( p = + e ( ↑ r e )( p = + e ( ↑ r e ) Figure 9: Convergence history of the lowest eigenvalue λ for the harmonic oscillator potential (21) attained fordi ff erent refinement schemes. We consider a purely polynomial approximation ( p = , ,
3) on a sequence of uniformlyrefined covers, and a refinement by increasing the enrichment radius (compare Figure 6) with a single enrichmentfunction on a fixed uniform cover (top: 3 × ×
3, 7 × ×
7; bottom: 4 × ×
4, 8 × ×
8) that is labeled by p = , , + e ↑ r e . All results were obtained with the lumped overlap matrix. − − − − − − − degrees of freedom a b s o l u t ee rr o r( H a ) p = p = p = (a) α = . − − − − degrees of freedom a b s o l u t ee rr o r( H a ) p = p = p = (b) α = . − − − − degrees of freedom a b s o l u t ee rr o r( H a ) p = p = p = (c) α = . Figure 10: Convergence history of the sum of the ten lowest eigenvalues for the harmonic oscillator potential (21)attained for di ff erent values of α = . , . , . ω i with a lumped overlap matrix. The dashed line indicates the accuracy of the employedreference solution. order O (10 ) whereas our flat-top PUM with stability transform and lumping using a comparablesetup provides an error of 5 · − Ha with 810 DOFs and an optimal condition number of 1 for thelumped overlap matrix.Lastly, we examine the influence of the flat-top region on the eigenvalue convergence withthe lumped overlap matrix. To this end, we consider various values of the stretch factor (14) α = . , . , .
3. Figure 10 shows the convergence of the sum of the ten lowest eigenvalues fordi ff erent values of α . For these tests, we enrich each of the ten lowest states on every patch of ourcover and use the stability transform (17). From the depicted plots it becomes clear that lumpingprovides better results for smaller stretch factors α . The absolute value of the error grows withincreasing α as well as the pre-asymptotic range of the convergence behavior. For fine enoughcovers, i.e., in the asymptotic range, the convergence behavior is very much in agreement with theconsistent scheme. These observations can be explained by the structure of the error introduced bythe lumping scheme [18]. Recall that our PUM degenerates to a DG approach with α = α →
1. Thus, it is expected thatthe approximation quality of the computed eigenvalues will improve with smaller α . The a priori determination of an optimal stretch factor α for a particular choice of the cover and employed localapproximation spaces V i , however, is the subject of current research. As our second benchmark problem, we consider a periodic Gaussian potential [2]. This externalpotential is defined via V ( x ) = (cid:88) R V g ( | x − τ − R | ) , (24)with V g ( r ) = −
10 exp (cid:32) − r . (cid:33) , (25)where we sum over the lattice translation vectors R = i a + i a + i a , i d = − , . . . , d = , , Ω and the potential center τ are defined as in Section 4.1. For this experiment,15e solve (4) subject to Bloch-periodic boundary conditions and choose k = (0 . , . , . α = . R nl of the infinite box eigenfunctionsis multiplied by h ( r , r =
10 au in (23). Also, as before, we are measuring the absoluteerror with respect to a cubic finite element reference solution computed on a 64 × ×
64 mesh (with λ ref1 = − . (cid:80) i = λ ref i = − . × × − Ha for the lowest state.
5. Concluding remarks
In this paper, we addressed two key issues in quantum mechanical calculations using enrichedpartition of unity finite element methods: the need to solve a generalized rather than standardeigenvalue problem, and the ill-conditioning of the associated system matrices. To address these,we developed a stable and e ffi cient orbital-enriched flat-top partition of unity method to solve therequired Schr¨odinger equation subject to Bloch-periodic boundary conditions in a general paral-lelepiped domain. To this end, we employed a stable transformation and variational mass lumpingin a flat-top partition of unity formulation. Compared to PUFEM approaches used previously,the present method yields well- (or even optimally-) conditioned system matrices and a standardeigenvalue problem rather than a generalized one, while maintaining accuracy and systematic con-vergence to benchmark results.With the above key issues resolved, future work will focus on the incorporation of nonlocalpseudopotentials, as occur in practical calculations involving all but the lightest elements, system-atic determination of maximal α consistent with chemical accuracy with variational mass lumping,and e ffi cient numerical cubature. A corresponding implementation for the Poisson equation willthen enable full Kohn–Sham calculations. Acknowledgments
This work was performed, in part, under the auspices of the U.S. Department of Energy byLawrence Livermore National Laboratory under Contract DE-AC52-07NA27344.160 − − − − degrees of freedom a b s o l u t ee rr o r( H a ) p = p = p = p = + e ( ↑ r e )( p = + e ( ↑ r e )( p = + e ( ↑ r e ) − − − − degrees of freedom a b s o l u t ee rr o r( H a ) p = p = p = p = + e ( ↑ r e )( p = + e ( ↑ r e )( p = + e ( ↑ r e ) Figure 11: Convergence history of the lowest eigenvalue λ for the periodic Gaussian potential (24) using di ff erentrefinement schemes. We consider a purely polynomial approximation ( p = , ,
3) on a sequence of uniformly refinedcovers, and a refinement by increasing the enrichment radius (compare Figure 6) with a single enrichment function ona fixed uniform cover (top: 3 × ×
3, 7 × ×
7; bottom: 4 × ×
4, 8 × ×
8) that is labeled by p = , , + e ↑ r e . Allresults were obtained with the lumped overlap matrix. eferences [1] N. W. Ashcroft, N. D. Mermin, Solid State Physics, Holt, Rinehart and Winston, New York,1976.[2] N. Sukumar, J. E. Pask, Classical and enriched finite element formulations for Bloch-periodicboundary conditions, International Journal for Numerical Methods in Engineering 77 (8)(2009) 1121–1138.[3] L. Lin, J. Lu, L. Ying, W. E, Adaptive local basis set for Kohn-Sham density functional theoryin a discontinuous Galerkin framework I: Total energy calculation, Journal of ComputationalPhysics 231 (4) (2012) 2140–2154.[4] G. Zhang, L. Lin, W. Hu, C. Yang, J. E. Pask, Adaptive local basis set for Kohn-Sham densityfunctional theory in a discontinuous Galerkin framework II: Force, vibration, and moleculardynamics calculations, Journal of Computational Physics 335 (2017) 426–443.[5] J. E. Pask, N. Sukumar, Partition of unity finite element method for quantum mechanicalmaterials calculations, Extreme Mechanics Letters 11 (2017) 8–17.[6] S. Yamakawa, S. Hyodo, Electronic state calculation of hydrogen in metal clusters basedon Gaussian-fem mixed basis function, Journal of Alloys and Compounds 356-357 (2003)231–235.[7] S. Yamakawa, S. Hyodo, Gaussian finite-element mixed-basis method for electronic structurecalculations, Phys. Rev. B 71 (3) (2005) 035113.[8] J. S. Chen, W. Hu, M. Puso, Orbital hp -cloud for solving Schr¨odinger equation in quantummechanics, Comput. Methods Appl. Mech. Engrg. 196 (2007) 3693–3705.[9] J. E. Pask, N. Sukumar, M. Guney, W. Hu, Partition-of-unity finite-element method for largescale quantum molecular dynamics on massively parallel computational platforms, Tech.Rep. LLNL-TR-470692, Department of Energy LDRD Grant 08-ERD-052 (March 2011).[10] J. E. Pask, N. Sukumar, S. E. Mousavi, Linear scaling solution of the all-electron Coulombproblem in solids, International Journal for Multiscale Computational Engineering 10 (1)(2012) 83–99.[11] A. S. Banerjee, L. Lin, W. Hu., C. Yang, J. E. Pask, Chebyshev polynomial filtered subspaceiteration in the discontinuous Galerkin method for large-scale electronic structure calcula-tions, Journal of Chemical Physics 145 (15) (2016) 154101.[12] B. Kanungo, V. Gavini, Large-scale all-electron density functional theory calculations usingan enriched finite-element basis, Phys. Rev. B 95 (2017) 035112.[13] D. Davydov, T. Gerasimov, J.-P. Pelteret, P. Steinmann, Convergence study of the h -adaptivepum and the hp -adaptive fem applied to eigenvalue problems in quantum mechanics, Ad-vanced Modeling and Simulation in Engineering Sciences 4 (1) (2017) 7.1814] I. Babuˇska, J. M. Melenk, The partition of unity finite element method: Basic theory andapplications, Comput. Meth. Appl. Mech. Engrg. 139 (1996) 289–314.[15] I. Babuˇska, J. M. Melenk, The partition of unity method, Int. J. Numer. Meth. Engrg. 40(1997) 727–758.[16] M. Griebel, M. A. Schweitzer, A particle-partition of unity method—Part II: E ffi cient coverconstruction and reliable integration, SIAM J. Sci. Comput. 23 (5) (2002) 1655–1682.[17] M. A. Schweitzer, Stable enrichment and local preconditioning in the particle–partition ofunity method, Numer. Math. 118 (1) (2011) 137–170.[18] M. A. Schweitzer, Variational mass lumping in the partition of unity method, SIAM Journalon Scientific Computing 35 (2) (2013) A1073–A1097.[19] F. Ihlenburg, Finite Element Analysis of Acoustic Scattering, Vol. 132 of Applied Mathemat-ical Sciences, Springer-Verlag, New York, NY, 1998.[20] P. A. Deymier (Ed.), Acoustic Metamaterials and Phononic Crystals, Spring Series in Solid-State Sciences, Springer, New York, NY, 2013.[21] J. D. Joannopoulos, S. G. Johnson, J. N. Winn, R. D. Meade, Photonic Crystals: Molding theFlow of Light, 2nd Edition, Princeton University Press, Princeton, NJ, 2008.[22] M. A. Schweitzer, A Parallel Multilevel Partition of Unity Method for Elliptic Partial Dif-ferential Equations, Vol. 29 of Lecture Notes in Computational Science and Engineering,Springer, 2003.[23] C. Klaar, A partition of unity method for quantum mechanical material calculations, Master’sthesis, Rheinische Friedrich-Wilhelms-Universitt Bonn, Institut fr Numerische Simulation(2016).[24] T.-P. Fries, T. Belytschko, The extended / generalized finite element method: An overview ofthe method and its applications, Int. J. Numer. Meth. Engrg. 84 (3) (2010) 253–304.[25] M. A. Schweitzer, Generalizations of the finite element method, Cent. Eur. J. Math. 10 (1)(2012) 3–24.[26] I. Babuˇska, G. Caloz, J. E. Osborn, Special finite element methods for a class of second orderelliptic problems with rough coe ffi cients, SIAM J. Numer. Anal. 31 (1994) 945–981.[27] M. A. Schweitzer, Generalized finite element and meshfree methods, Habilitation, Univerist¨atBonn (2008).[28] M. Griebel, M. A. Schweitzer, A particle-partition of unity method—Part VII: Adaptivity, in:M. Griebel, M. A. Schweitzer (Eds.), Meshfree Methods for Partial Di ff erential EquationsIII, Vol. 57 of Lecture Notes in Computational Science and Engineering, Springer, 2006, pp.121–148. 1929] M. A. Schweitzer, An adaptive hp-version of the multilevel particle–partition of unity method,Comput. Meth. Appl. Mech. Engrg. 198 (2009) 1260–1272.[30] T. Strouboulis, I. Babuˇska, R. Hidajat, The generalized finite element method for Helmholtzequation: Theory, computation, and open problems, Comput. Meth. Appl. Mech. Engrg.195 (37-40) (2006) 4711–4731.[31] T. Strouboulis, R. Hidajat, I. Babuˇska, The generalized finite element method for Helmholtzequation. Part II: E ff ect of choice of handbook functions, error due to absorbing boundaryconditions and its assessment, Comput. Meth. Appl. Mech. Engrg. 197 (5) (2008) 364–380.[32] C. A. M. Duarte, J. T. Oden, hp clouds – A meshless method to solve boundary value prob-lems, Numer. Meth. for PDE 12 (1996) 673–705.[33] T. Belytschko, Y. Y. Lu, L. Gu, Crack propagation by element-free Galerkin methods, Engrg.Frac. Mech. 51 (1995) 295–315.[34] M. A. Schweitzer, A parallel multilevel partition of unity method for elliptic partial di ffff