The Partition of Unity Finite Element Method for the Schrödinger Equation

1 Introduction

According to the density functional theory (DFT) [25, 21, 17], we can derive material properties from ab initio quantum mechanics calculations by solving the system of the Khon–Sham equations [20].

This approach requires repeatedly solving the steady-state Schrödinger and Poisson equations on a parallelepiped unit cell with Bloch-periodic boundary conditions [3]. The Schrödinger equation is an eigenvalue problem for the Hamiltonian differential operator describing the energy of a given multi-electron atomic or molecular system. The eigenvalues are the discrete energy levels, while the eigenfunctions are the wavefunctions of the corresponding quantum states [28, 24].

DFT represents a quantum mechanical computational methodology that facilitates the examination of the electronic structures and properties of molecular and solid-state systems. The foundational principle of DFT resides in the conceptualization of electronic systems via electron density profiles, thereby circumventing the direct resolution of the complex many-body Schrödinger equation. Within the framework of DFT, the more difficult many-body problem is effectively transformed into a tractable single-electron effective potential equation, analogous to a three-dimensional Schrödinger equation. The method delineated herein seeks to harness DFT within this scope.

The many-body Schrödinger equation delineates the interactions among the electrons in a system, presenting a formidable computational challenge, particularly for substantial systems. DFT streamlines this challenge by postulating the electron density as a fundamental descriptor. Electron density conveys the positional probability density of electrons within a spatial domain. The seminal Hohenberg–Kohn theorems assert that the ground-state electron density coherently determines the ground-state energy of a system. Consequently, it is feasible to ascertain the ground-state electron density and subsequently employ it to deduce diverse system properties, obviating the need to explicitly solve the many-body Schrödinger equation.

In the DFT paradigm, the derivation of the effective potential equation incorporates the electron-electron interactions as an effective field influencing each electron in isolation. This effective field amalgamates the external potential attributable to the atomic nuclei, and the exchange-correlation potential, which encapsulates the complexities of electron-electron repulsions. The objective is to pinpoint the electron density that minimizes the total energy of the system, an endeavor accomplished through the resolution of the effective potential equation.

DFT combined with the effective potential equation makes elucidating a system’s electronic structure and properties possible. By resolving the effective potential equation, one can determine the electron density and, thereby, calculate various properties, including the total energy, electronic charge distribution and the configuration of molecular orbitals. These properties yield profound insights into molecular and material behaviors and reactivities.

In essence, the proposed method leverages DFT by simplifying the many-body Schrödinger equation to a manageable single-electron, three-dimensional equation. The solution of this equation grants access to the electron density and ancillary properties of the system, thereby enabling a comprehensive exploration of the electronic structure and properties of molecular entities.

It is by far recognized that solving the Schrödinger equation is the most time-consuming step in such calculations. The plane wave (PW) pseudopotential method is currently the most advanced technology implemented in application software to accomplish this task. The PW method employs a Fourier basis set and calls for the resolution of a discrete standard eigenproblem.

However, Galerkin methods with enrichment can be highly competitive with the PW approach to obtain a numerical solution with a required accuracy using the same number of basis functions (or even less), i.e., by significantly reducing the number of degrees of freedom. An incomplete list of research articles exploring this possibility includes the references [2, 32, 22, 36, 26, 33, 34, 11, 27, 8, 18, 12]. The two most important issues with the enriched Galerkin method are the need to solve a generalized rather than a regular eigenvalue problem and the poor conditioning of the associated system matrices. These two crucial points can negatively impact computational costs. The partition of unity finite element method (PUFEM) [14] offers a very effective way to address these points and numerically solve the Schrödinger eigenvalue problem in a parallelepiped unit cell subject to Bloch-periodic boundary conditions. The partition of unity (PUM) [23, 5] is a generalization of the finite element method (FEM) for the spatial discretization of partial differential equations (PDEs); see also [4, 30]. Its application to the Schrödinger eigenvalue problem with Bloch-periodic boundary conditions was recently investigated in [32, 26]. The basic idea is that we cover the three-dimensional domain by overlapping patches so that each patch compactly supports a nonnegative weight function. These weight functions form a partition of unity and have the flat-top property, i.e., it is identically equal to one over some finite subset of its support. We employ 𝑝-th degree orthogonal polynomials on each patch built as the tensor product of univariate rescaled and translated Legendre polynomials. This formulation guarantees the 𝑝-th order completeness. We enrich the PUFEM by incorporating the radial Schrödinger’s eigenfunctions. The discrete variational form results in a generalized eigenvalue problem. Using the variational lumping of [31], we obtain a regular eigenvalue problem for which effective eigensolvers are available. This approach for solving the Schrödinger eigenvalue problem with the harmonic potential is accurate, stable and effective. In fact, for a given required accuracy in the energy eigenvalues, the suggested strategy provides optimal rates of convergence in the energy and the application of orbital enrichment dramatically lower the number of degrees of freedom.

In this work, we provide the proper mathematical and theoretical setting of the Schrödinger eigenvalue problem in a unit cube with Bloch-periodic condition for the partition of unity method. Our approach is based on the spectral approximation theory of compact operators, cf. [6] , and is, in practice, independent of the type of partition of unity functions that we can use in the PUFEM construction. To derive the convergence rates for the eigenvalue and eigenfunction approximation, we need the error bounds for the source problem, which is equivalent to a diffusion-reaction equation (the reaction term being the potential of the Hamilton operator). We obtain these estimates from the basic assumptions on the PUM of the original works of references [23, 5].

The paper is organized as follows. In Section 2, we present the Schrödinger eigenvalue problem on the unit cube with Bloch-periodic boundary conditions. In Section 3, we review the partition of unity finite element method. In Section 4, we prove the convergence of the method and deduce error estimates for the eigenvalue and eigenfunction approximation. In Section 5, we briefly detail the implementation of the method. In Section 6, we assess the numerical behavior of the method and show that it agrees with the theoretical expectations. In Section 7, we offer our final remarks and mention a few possible future research topics along the lines of this work.

2 The Schrödinger Eigenvalue Problem

In the context of atomic and molecular quantum theory, the Schrödinger eigenvalue problem on the domain Ω reads as

where H : = − 1 2 Δ + V eff ( x ) is the Hamilton operator. The eigenvalue 𝜀 and the eigenfunction 𝜓 are the energy level and the corresponding wavefunction of the physical system described by ℋ. From density functional theory [7, 35, 16], the effective potential V eff ⁢ ( x ) in ℋ collects all local and nonlocal interaction terms.

To describe a periodic condensed matter system, we consider the parallelepiped unit-cell Ω ⊂ R d , with primitive lattice vectors a d , d = 1 , 2 , 3 , and the lattice translation vectors R = n 1 ⁢ a 1 + n 2 ⁢ a 2 + n 3 ⁢ a 3 , for all n d ∈ Z , d = 1 , 2 , 3 . Then we assume that the effective potential is periodic in this lattice structure; hence, we have V eff ⁢ ( x + R ) = V eff ⁢ ( x ) . The wavefunction 𝜓 solving Schrödinger’s equation (2.1) must satisfy Bloch’s theorem; see [3, Chapter 8]. Therefore, for any eigenfunction 𝜓, there is a wavevector 𝐤 such that ψ ⁢ ( x + R ) = ψ ⁢ ( x ) ⁢ exp ⁡ ( i ⁢ k ⋅ R ) . For k = 0 (Γ-point), the wavefunction is periodic. For k ≠ 0 , there is a phase shift exp ⁡ ( i ⁢ k ⋅ R ) associated with the translation 𝐑. In this setting, we write problem (2.1) as

We consider the functional space of complex-valued functions

endowed with the usual H 1 norm denoted by ∥ ⋅ ∥ V ; the L 2 -scalar product will be denoted by

We introduce the sesquilinear form A : V × V → C ,

and assume that V eff ⁢ ( x ) is uniformly bounded from below and above, i.e., there exist two strictly positive constants V ∗ and V ∗ such that V ∗ ≤ V eff ⁢ ( x ) ≤ V ∗ for almost every x ∈ Ω . Consequently, the sesquilinear form A ⁢ ( ⋅ , ⋅ ) is coercive and continuous, that is, there exist positive constants 𝛼 and 𝐶 such that

The variational formulation of (2.1) reads as: find ( ε , ψ ) ∈ R × V , with ∥ ψ ∥ 0 , Ω = 1 , such that

From the classic eigenvalue theory [9, 6], equation (2.5) admits an infinite discrete set of strictly positive eigenvalues, which diverge to infinity. Moreover, such eigenvalues can have a multiplicity greater than one with a finite-dimensional eigenspace. The eigenfunctions form a basis of 𝑉 that is orthogonal with respect to the L 2 -inner product (2.3) and the inner product defined by the sesquilinear form (2.4).

The discretization of (2.1) is carried out by introducing a sequence of conforming, finite-dimensional subspaces of 𝑉, denoted by V ℓ PU and numbered with the integer label ℓ → ∞ . This label denotes the refinement level in the case of the ℎ-refinement and the degree of the polynomials in the case of the 𝑝-refinement. The construction of the approximation spaces V ℓ PU will be detailed in the next section. The discrete eigenvalue problem reads as: find ( ε ℓ , ψ ℓ ) ∈ R × V ℓ PU , with ∥ ψ ℓ ∥ 0 , Ω = 1 , such that

We also consider the continuous source problem associated with the eigenvalue problem (2.5): find ψ s ∈ V such that

where we assume that g ∈ L 2 ⁢ ( Ω ) . It is standard to prove the well-posedness of problem (2.7), i.e., existence and uniqueness of its solution, by using the Lax–Milgram lemma [10]. In fact, due to the boundedness assumption on the potential field V eff , the bilinear form 𝒜 in (2.4) is coercive and the right-hand side of (2.7) is a continuous functional. Moreover, due to regularity result [1, 15], there exists r > 0 , depending only on Ω, such that ψ s ∈ H 1 + r ⁢ ( Ω ; C ) . The following stability estimate holds:

Eventually, it is known that r = 1 if the domain Ω is convex, as in our case, so we can deduce that ψ s ∈ H 2 ⁢ ( Ω ; C ) .

Similarly, the discrete source problem reads as: find ψ ℓ s ∈ V ℓ PU such that

The well-posedness of problem (2.9) follows from the fact that all the approximation spaces V ℓ PU are conforming subspaces of 𝑉, and the continuous problem on 𝑉 is well-posed.

3 Partition of Unity Finite Element Method

In this section, we discuss the construction of the approximation space V PU used in the partition of unity finite element method (PUFEM). Since the indication of the refinement level is not important, we remove the sub-index ℓ to simplify the notation. The main ingredients of the PUFEM are

1. a set of overlapping patches { Ω i } i , i = 1 , 2 , … , n , forming a finite, open covering of the computational domain Ω, i.e., Ω ⊆ ⋃ i = 1 n Ω i ;

2. a set of Lipschitz-continuous functions { ϕ i } i = 1 , 2 , … , n such that supp ⁢ ( ϕ i ) ⊆ Ω ̄ i and forming a partition of unity, i.e., ∑ i = 1 n ϕ i ⁢ ( x ) = 1 for all x ∈ Ω ;

3. a set of local approximation spaces V i ⊂ H 1 ⁢ ( Ω i ) .

The global regularity of the approximation space V PU stems from the global regularity of the partition of unity functions ϕ i and the overlapping of the patches Ω i . Since, by construction, the global space V PU is a subspace of 𝑉, its functions automatically satisfy the periodicity condition (2.2).

Let diam ⁡ ( Ω i ) = sup x ′ , x ′′ ∈ Ω i | x ′ − x ′′ | denote the diameter of the 𝑖-th patch. According to [23, 5], we consider the following assumption on the partition of unity.

1. A positive integer 𝑀, and two positive real constants C ∞ and C G exist such that

   1. card ⁢ { Ω i : x ∈ Ω i } ≤ M for every x ∈ Ω ;

   2. ∥ ϕ i ∥ ∞ , Ω i ≤ C ∞ ;

   3. ∥ ∇ ϕ i ∥ ∞ , Ω i ≤ C G diam ⁡ ( Ω i ) .

The property A of the partition of unity and the local approximation properties of the spaces V i determine the approximation result of References [23, 5], which we report below for the sake of completeness and future reference in the paper.

In view of the previous theorem, the space V PU inherits the approximation properties of the local spaces V i . In particular, if the patches are star-shaped with respect to a ball of points and the local spaces contain the polynomials of degree at most 𝑝, we can bound the approximation errors ϵ 0 ⁢ ( i ) and ϵ 1 ⁢ ( i ) in the ℎ-refinement setting using the standard estimates of the polynomial interpolation in Sobolev spaces, cf. [13, 10]. Therefore, for every patch Ω i , it holds that

Hence, using (3.2)–(3.3) from Theorem 1, we find the global approximation estimates for u PU ,

The constant 𝐶 is independent of ℎ, but depends on C ∞ , C G and 𝑀.

According to the standard convergence theory of the conforming finite element methods, cf. Cea’s lemma and the duality argument, we obtain the following bounds on the error for the solution of the source problem (2.9).

If we consider the 𝑝-refinement setting, we can bound the approximation errors ϵ 0 ⁢ ( i ) and ϵ 1 ⁢ ( i ) as follows:

Hence, using (3.2)–(3.3) from Theorem 1, we find the global approximation estimates for u PU ,

Again, according to the standard convergence theory of the conforming finite element methods, we obtain the following bounds with respect to 𝑝 on the error for the solution of the source problem (2.9).

4 Convergence Analysis and Error Estimates

4.3 Filtering PUFEM and Enriched PUFEM

While in the classical setting of PUFEM functions in V i are polynomials in Ω i , we are interested in the enriched version of PUFEM, where the elements of V i may be polynomials in Ω i or other functions, defined in an appropriate way in order to better approximate the exact solution. In order to formalize this concept, we might write the 𝑖-th local approximation space V i as

V i = V i P + V i E ,

where V i P = span ⁡ { π k ( i ) } and V i E = span ⁡ { ψ k ( i ) } . The notation refers to the fact that usually V i P will be the polynomial part of the space and V i E the enrichment.

Figure 1

1D structure of the overlapping patches (the top line), and integration cells (bottom line) corresponding to the various patches. The first and last integration cells are due to the overlap with two patches that are not shown in the figure.

We implement PUFEM and enriched PUFEM by assembling the global stiffness and mass matrices from local contributions. We compute these local contributions on a set of nonoverlapping integration cells that partition the computational domain by taking into account possible patch overlappings. For example, in the 1D case, see Figure 1, we distinguish two different types of integration cells:

1. cells containing a single patch;

2. cells shared by two overlapping patches.

In the former case, the local stiffness and mass matrices only contain the contributions from the basis functions { π k ( i ) } of the associated patch. In the latter case, such matrices depend on union of the basis functions of the overlapping patches. This procedure leads to a “redundant” construction and, eventually, to singular matrices.

The same consideration holds for the enrichment part of the local approximation spaces. Moreover, the enrichment procedure may suffer of numerical instability because the refinement may lead to an ill-conditioned basis for V i since some components in V i E might be better and better approximated by V i P . Similar issues exist also in the 2D and 3D cases; technical details are reported in the final appendix.

We now describe a filtering procedure that we are going to use in order to address these issues. At the end of the process, we are going to have a reduced space V ̃ i that will be presented as

V ̃ i = V ̃ i P ⊕ V ̃ i D ,

where V ̃ i D is that part of V ̃ i E that is orthogonal to V ̃ i P .

The reduction procedure is based on a threshold tolerance τ > 0 that will be responsible to filtering the elements of V i that cause the ill-conditioning. The eigenvalues responsible for the ill-conditioning are very close to the machine precision, typically of the order of 10 − 14 , and are well-separated from the “good” eigenvalues, which are generally clustered around 1 after a matrix rescaling. Therefore, we can set up the threshold tolerance to a small value just above the machine precision, e.g., in the range [ τ = 10 − 13 , 10 − 12 ] . Hence, as long as 𝜏 is sufficiently small not to affect the “good” eigenvalues and sufficiently big to remove all the “bad” ones, the accuracy of the final approximation is independent of its choice.

The filtering process consists of three main steps: an initial rescaling of the functions forming a basis of V i P and V i E ; a reduction step based on a first internal orthogonalization leading to two reduced approximation spaces V ̃ i P and V ̃ i E ; a final orthogonalization step between the basis of V ̃ i P and V ̃ i E .

Consider the so called overlap (mass) matrix

S i = [ S i , P ⁢ P S i , P ⁢ E S i , E ⁢ P S i , E ⁢ E ] ,

where the entries are the block sub-matrices ( S i , P ⁢ P ) | a ⁢ b = ( π a i , π b i ) , ( S i , P ⁢ E ) | a ⁢ b = ( π a i , ψ b i ) , ( S i , E ⁢ P ) | a ⁢ b = ( ψ a i , π b i ) , and ( S i , E ⁢ E ) | a ⁢ b = ( ψ a i , ψ b i ) for every admissible value of the index pair ( a , b ) .

First Step

We rescale the basis functions π a ( i ) and ψ a ( i ) , and introduce two new sets of basis functions { π ̂ a ( i ) } and { ψ ̂ a ( i ) } that are defined as follows:

π ̂ a ( i ) = | Ω i | 1 / 2 ∥ π a ( i ) ∥ 0 , Ω i ⁢ π a ( i ) and ψ ̂ a ( i ) = | Ω i | 1 / 2 ∥ π a ( i ) ∥ 0 , Ω i ⁢ ψ a ( i ) .

According to [29], we adopt the weighted inner product on the functional space L 2 ⁢ ( Ω i , C , { ϕ i } ) and, in practice, we build the enriched functional spaces through a double orthogonalization and reduction process that is aimed at improving the stability of the discrete problem without significantly affecting its convergence properties.

Second Step

By using these scaled basis functions, we recompute a new overlap matrix S ̂ i and note that its diagonal entries are equal to | Ω i | . Since the overlap matrix S ̂ i is real and symmetric, and so are the diagonal matrix blocks S ̂ i , P ⁢ P and S ̂ i , E ⁢ E , we perform the eigendecomposition

D i , P = O i , P T ⁢ S ̂ i , P ⁢ P ⁢ O i , P and D i , E = O i , E T ⁢ S ̂ i , E ⁢ E ⁢ O i , E ,

where D i , ∗ and O i , ∗ with ∗ ∈ { P , E } , are the eigenvalue and eigenvector matrices of, respectively, S ̂ i , P ⁢ P and S ̂ i , E ⁢ E . Note that the eigenvector matrices O i , ∗ are orthogonal, i.e., O i , ∗ T ⁢ O i , ∗ = I i , ∗ . Using these matrices, we define a basis transformation providing two new sets of basis functions, namely, { π ̃ a ( i ) } and { ψ ̃ a ( i ) } given by the linear combinations of { π ̂ a ( i ) } and { ψ ̂ a ( i ) } , and using the components of the matrices O i , P and O i , E as the coefficients. We now make use of the tolerance 𝜏 and write

π ̃ a ( i ) = ∑ b ( O i , P ) a , b ⁢ π ̂ b ( i ) = ∑ b : ( D i , P ) b ⁢ b > τ ( O i , P ) a , b ⁢ π ̂ b ( i ) + ∑ b : ( D i , P ) b ⁢ b ≤ τ ( O i , P ) a , b ⁢ π ̂ b ( i ) = π ̃ a ( i , 1 ) + π ̃ a ( i , 2 ) , ψ ̃ a ( i ) = ∑ b ( O i , P ) a , b ⁢ ψ ̂ b ( i ) = ∑ b : ( D i , P ) b ⁢ b > τ ( O i , E ) a , b ⁢ ψ ̂ b ( i ) + ∑ b : ( D i , P ) b ⁢ b ≤ τ ( O i , E ) a , b ⁢ ψ ̂ b ( i ) = ψ ̃ a ( i , 1 ) + ψ ̃ a ( i , 2 ) .

In these relations, we split the summation over the index 𝑏 to separate the contributions corresponding to the eigenvalues bigger or smaller than the given tolerance factor 𝜏. This eigendecomposition has two crucial properties:

1. { π ̃ a ( i ) } is an orthogonal basis set in V i P and { ψ ̃ a ( i ) } is an orthogonal basis set in V i E ;

2. ∥ π ̃ a ( i , 2 ) ∥ 0 , Ω i 2 ≤ τ and ∥ ψ ̃ a ( i , 2 ) ∥ 0 , Ω i 2 ≤ τ .

Property (i) follows immediately from the definitions introduced so far. For example, for the basis set { π ̃ a ( i ) } , a straightforward calculation proves that

( π ̃ a ( i ) , π ̃ a ′ ( i ) ) = ∑ b ⁢ b ′ ( O i , P ) a , b ⁢ ( O i , P ) a ′ , b ′ ⁢ ( π ̂ b ( i ) , π ̂ b ′ ( i ) ) = [ ( O i , P ) T ⁢ S ̂ ⁢ O i , P ] a ⁢ a ′ = ( D i , P ) a ⁢ a ′ = ( D i , P ) a ⁢ a ⁢ δ a ⁢ a ′ .

Property (ii) readily follows by taking a = a ′ in the calculation above and considering only the summation terms providing π ̃ a ( i , 2 ) . We obtain ∥ π ̃ a ( i , 2 ) ∥ 0 , Ω i 2 = ( D i , P ) a ⁢ a ′ ≤ τ . The same arguments prove the corresponding relations for ψ ̃ a ( i , 2 ) . At this point, we perform the spectral reduction by taking π ̃ a ( i ) ≈ π ̃ a ( i , 1 ) and ψ ̃ a ( i ) ≈ ψ ̃ a ( i , 1 ) . We denote the spaces generated by the linear combinations of these reduced set of basis functions by V ̃ i P and V ̃ i E , respectively. Note that these functions still form two sets of orthogonal functions.

Third Step

In the third and final step, we mutually orthogonalize the functions π ̃ a ( i ) and ψ ̃ a ( i ) so that the decomposition

V ̃ i = V ̃ i P ⊕ V ̃ i D ,

holds, where now V ̃ i D is the space of functions of V ̃ i E that are orthogonal to V ̃ i P .

If we reformulate problem (2.6) by considering the space V \PU defined in (3.1) with V ̃ i instead of V i , we can repeat the convergence analysis of the previous sections by using the same arguments as long as the right approximation property is satisfied. In the next section, we will discuss how the spectral reduction strategy can be effectively implemented.

5 Implementation for the Schrödinger Equation

In our implementation, we use the flat-top partition of unit that was initially proposed in [2]. The design of the partition of unity method that we consider in this paper first requires a set of overlapping patches { Ω i } , i = 1 , 2 , … , n . Let h = 1 / n be the mesh size parameter and consider the univariate partition of the unit interval [ 0 , 1 ] in every space direction given by the points x i s = 2 ⁢ h ⁢ i , s = 1 , … , d , where d = 1 , 2 , 3 . Then we define the 𝑖-th patch in the 𝑠-th direction as