A discontinuous Galerkin scheme for full-potential electronic structure calculations

TL;DR

This paper introduces an efficient discontinuous Galerkin scheme for full-potential electronic structure calculations, combining different basis functions in atomic and interstitial regions with spectral convergence.

Contribution

It develops a novel DG-based numerical scheme with spectral convergence for periodic electronic structure calculations, inspired by (L)APW methods.

Findings

Scheme achieves spectral convergence rate

Provides rigorous a priori error analysis

Demonstrates effectiveness through numerical simulations

Abstract

In this paper, we construct an efficient numerical scheme for full-potential electronic structure calculations of periodic systems. In this scheme, the computational domain is decomposed into a set of atomic spheres and an interstitial region, and different basis functions are used in different regions: radial basis functions times spherical harmonics in the atomic spheres and plane waves in the interstitial region. These parts are then patched together by discontinuous Galerkin (DG) method. Our scheme has the same philosophy as the widely used (L)APW methods in materials science, but possesses systematically spectral convergence rate. We provide a rigorous a priori error analysis of the DG approximations for the linear eigenvalue problems, and present some numerical simulations in electronic structure calculations.

Equations

$$H_{\Phi }{\varphi }_i={\lambda }_i{\varphi }_i,{\lambda }_1\le {\lambda }_2\le \cdots \le {\lambda }_{N_{\mathrm{e}}},$$

$$H_{\Phi }=-\frac{1}{2}\Delta +V_{\mathrm{ext}}+V_{\mathrm{H}}[{\rho }_{\Phi }]+V_{\mathrm{xc}}[{\rho }_{\Phi }].$$

$$|\lambda -{\lambda }_i^{\mathrm{DG}}|+\Vert u_i-u_i^{\mathrm{DG}}{\Vert }_{\mathrm{DG}}\le C_s{\varrho }^{\frac{3}{2}+ϵ-s}\forall s\in {\mathbb{R}}^{+},$$

$$L_{\#}^2(\Omega )=\{u\in L_{\mathrm{loc}}^2({\mathbb{R}}^3):u\mathrm{is}\mathcal{R}-\mathrm{periodic}\}.$$

$$u(\mathbf{r})=\sum _{\mathbf{k}\in {\mathcal{R}}^{\ast }}{\hat{u}}_{\mathbf{k}}{e}_{\mathbf{k}}(\mathbf{r})\mathrm{with}{\hat{u}}_{\mathbf{k}}=(u,e_{\mathbf{k}}{)}_{L_{\#}^2(\Omega )}=|\Omega {|}^{-1/2}{\int }_{\Omega }u(\mathbf{r})e^{-i\mathbf{k}\cdot \mathbf{r}}d\mathbf{r}.$$

$$\displaystyle H_{\#}^{s}(\Omega)=\left\{u({\bf r})=\sum_{{\bf k}\in\mathcal{R}^{*}}\hat{u}_{\bf k}e_{\bf k}({\bf r})~{}:~{}\sum_{{\bf k}\in\mathcal{R}^{*}}\big{(}1+|{\bf k}|^{2}\big{)}^{s}|\hat{u}_{\bf k}|^{2}<\infty\right\},$$

$${\mathcal{V}}_K=\left\{v_K(\mathbf{r})=\sum _{\mathbf{k}\in {\mathcal{R}}^{\ast },|\mathbf{k}|\le \frac{2\pi }{D}K}{c}_{\mathbf{k}}{e}_{\mathbf{k}}(\mathbf{r})\right\}.$$

$$\Vert v-{\Pi }_{K}v{\Vert }_{H_{\#}^t(\Omega )}=\underset{v_K\in {\mathcal{V}}_K}{\min }\Vert v-v_K{\Vert }_{H_{\#}^t(\Omega )}\lesssim K^{t-s}\Vert v{\Vert }_{H_{\#}^s(\Omega )}.$$

$$a(u,v)=\lambda (u,v)\forall v\in H_{\#}^1(\Omega ),$$

$$a(u,v)=\frac{1}{2}{\int }_{\Omega }\nabla u\cdot \nabla v+{\int }_{\Omega }Vuv$$

$$\Vert v{\Vert }_{H^1({\Omega }_{\mathrm{in}})}\le C_R\Vert v{\Vert }_{H^1([0,R]\times S^2)}.$$

$$\Delta =\frac{1}{r^2}\frac{\partial }{\partial r}\left(r^2\frac{\partial }{\partial r}\right)+\frac{1}{r^2\sin \theta }\frac{\partial }{\partial \theta }\left(\sin \theta \frac{\partial }{\partial \theta }\right)+\frac{1}{r^2{\sin }^2\theta }\frac{{\partial }^2}{{\partial }^2\varphi },$$

$$\underset{r\to 0}{\lim }{r}^{2}v\frac{\partial v}{\partial r}=0,$$

$$\Vert v{\Vert }_{H^1({\Omega }_{\mathrm{in}})}^2$$

$${\tilde{H}}^s(\Omega )=\left\{v\in H_{\#}^1(\Omega ):v{|}_{{\Omega }_{\mathrm{in}}}\in H^s([0,R]\times S^2),v{|}_{{\Omega }_{\mathrm{out}}}\in H^s({\Omega }_{\mathrm{out}})\right\}$$

$$\Vert v{\Vert }_{{\tilde{H}}^s(\Omega )}:=\Vert v{\Vert }_{H^s({\Omega }_{\mathrm{out}})}+\Vert v{\Vert }_{H^s([0,R]\times S^2)}.$$

$$\displaystyle\mathcal{P}_{K}(\Omega_{\rm out})=\left\{u\in H^{1}(\Omega_{\rm out}),~{}u({\bf r})=\sum_{|{\bf k}|\leq\frac{2\pi}{D}K}c_{\bf k}e_{\bf k}({\bf r})\bigg{|}_{\Omega_{\rm out}}\right\}$$

$$\displaystyle\mathcal{B}_{NL}(\Omega_{\rm in})=\left\{u\in H^{1}(\Omega_{\rm in}),~{}u({\bf r})=\check{u}(r,\theta,\phi)=\sum_{0\leq n\leq N,0\leq l\leq L,|m|\leq l}c_{nlm}\chi_{n}(r)Y_{lm}(\theta,\phi)\bigg{|}_{\Omega_{\rm in}}\right\},$$

$${\mathcal{S}}_{NL}^K(\Omega ):$$

$$=\left\{u\in L_{\#}^2(\Omega ),u{|}_{{\Omega }_{\mathrm{in}}}\in {\mathcal{B}}_{NL}({\Omega }_{\mathrm{in}})\text{mbox}andu{|}_{{\Omega }_{\mathrm{out}}}\in {\mathcal{P}}_K({\Omega }_{\mathrm{out}})\right\}.$$

$$\max \{K,N,L\}\le \varrho \le C\min \{K,N,L\},$$

$$\Vert u-P_{K}u{\Vert }_{H^1({\Omega }_{\mathrm{out}})}=\underset{U_K^{\mathrm{out}}\in {\mathcal{P}}_K({\Omega }_{\mathrm{out}})}{\inf }\Vert u-U_K^{\mathrm{out}}{\Vert }_{H^1({\Omega }_{\mathrm{out}})}.$$

$$\Vert u-P_{K}u{\Vert }_{H^t({\Omega }_{\mathrm{out}})}\le C{K}^{t-s}\Vert u{\Vert }_{H^s({\Omega }_{\mathrm{out}})}.$$

$$u(\mathbf{r})=\sum _{lm}^{\infty }{u}_{lm}(r){\tilde{Y}}_{lm}(\mathbf{r})$$

$$\displaystyle\tilde{u}({\bf r})=\left\{\begin{array}[]{ll}u({\bf r})&\mbox{in }\Omega_{\rm out},\\[4.30554pt] \displaystyle\sum_{lm}^{\infty}\varphi_{lm}(r)\tilde{Y}_{lm}({\bf r})&\mbox{in }\Omega_{\rm in},\end{array}\right.$$

$$\Vert \tilde{u}{\Vert }_{H^s(\Omega )}\le \beta \Vert u{\Vert }_{H^s({\Omega }_{\mathrm{out}})},$$

$${\tilde{u}}_K=\sum _{|\mathbf{k}|\le \frac{2\pi }{D}K}{\tilde{c}}_{\mathbf{k}}{e}_{\mathbf{k}}\mathrm{with}{\tilde{c}}_{\mathbf{k}}={\int }_{\Omega }\tilde{u}(\mathbf{r})e_{\mathbf{k}}^{\ast }(\mathbf{r})d\mathbf{r},$$

$$\underset{U_K^{\mathrm{out}}\in {\mathcal{P}}_K({\Omega }_{\mathrm{out}})}{\inf }\Vert u-U_K^{\mathrm{out}}{\Vert }_{H^t({\Omega }_{\mathrm{out}})}$$

$$\le C{K}^{t-s}\Vert \tilde{u}{\Vert }_{H^s(\Omega )}\le C\beta K^{t-s}\Vert u{\Vert }_{H^s({\Omega }_{\mathrm{out}})},$$

$$\Vert v-P_{N}v{\Vert }_{H^1([0,R])}=\underset{{\psi }_N\in {\Psi }_N}{\inf }\Vert v-{\psi }_N{\Vert }_{H^1([0,R])},$$

$$P_L\phi (\theta ,\varphi )=\sum _{l=0}^L\sum _{m=-l}^l{\hat{\phi }}_{lm}{Y}_{lm}(\theta ,\varphi )\mathrm{with}{\hat{\phi }}_{lm}={\int }_0^{\pi }\sin \theta {\int }_0^{2\pi }\phi (\theta ,\varphi )Y_{lm}^{\ast }(\theta ,\varphi )d\varphi d\theta .$$

Full text

1 Introduction

Electronic structure calculations describe the energies and distributions of electrons, which plays a fundamental role in many different fields: materials science, biochemistry, solid-state physics, and surface physics. Among different electronic structure models, the Kohn-Sham density functional theory (DFT) [34] so far achieves the best compromise between accuracy and computational cost. For an $N_{\rm e}$-electron system with the presence of $M$ nuclei of charge $Z_{k}$ and located at ${\bf R}_{k}\in\mathbb{R}^{3}$ $(k=1,\cdots,M)$, Kohn-Sham DFT gives rise to the following nonlinear eigenvalue problems

[TABLE]

with $\Phi=\{\phi_{1},\cdots,\phi_{N_{\rm e}}\}$ and the Kohn-Sham Hamiltonian

[TABLE]

Here, $\displaystyle V_{\rm ext}({\bf x})=-\sum_{k=1}^{M}\frac{Z_{k}}{|{\bf x}-{\bf R}_{k}|}$ is the external potential generated by nuclear attraction, $\displaystyle V_{\rm H}[\rho_{\Phi}]=\int_{\mathbb{R}^{3}}\frac{\rho_{\Phi}({\bf y})}{|\cdot-{\bf y}|}d{\bf y}$ and $V_{\rm xc}[\rho_{\Phi}]$ are the so-called Hartree potential and exchange-correlation potential, respectively, with the electron density $\displaystyle\rho_{\Phi}({\bf x})=\sum_{i=1}^{N_{\rm e}}|\phi_{i}({\bf x})|^{2}$. A self-consistent field (SCF) iteration algorithm is commonly resorted to for these nonlinear problems. In each iteration, a Hamiltonian $H_{\tilde{\Phi}}$ is constructed from a trial electronic state $\tilde{\Phi}$, and a linear eigenvalue problem is then solved to obtain the low-lying eigenvalues and corresponding eigenfunctions. The loop continues until self-consistency of the electronic states is achieved. The efficiency of the algorithm is mainly determined by the discretization of the Hamiltonian, the self-consistent iteration, and the linear eigensolver. We shall focus ourselves on the discretization method in this paper.

For periodic systems, plane waves with pseudopotentials are natural methods which are simple to implement and give relatively accurate simulations. The pseudopotential approximations [34] replace singular nuclear attraction potential and complicated effects of the motion of core electrons by a smooth potential. They give satisfactory results in most cases, but sometimes fail. The mathematical analysis of the pseudopotential approximations is very rare, and we refer to [8, 14] for two recent works. Moreover, the core electrons have to be considered sometimes and are responsible for many properties. Therefore, the full-potential/all-electron calculations are necessary.

For eigenvalue problems with singular potentials in full-potential calculations, plane waves are inefficient bases for describing the cusps at the nuclei positions [23, 24, 25, 28]. In contrast, it is observed that a significant part of the rapid oscillations can be captured by atomic orbitals such as Gaussians and Slater-type orbitals [27, 34], which have been widely used in quantum chemistry (we refer to [7, 18] for their numerical analysis). Therefore, it would be practically efficient to approximate the wavefunction in a crystal by using combinations of plane waves and appropriate atomic orbitals. Several computational methods using this idea have been developed, for example, augmented plane waves (APW) [38, 40], linearized augmented plane waves (LAPW) [38], and their extensions by including local orbitals (lo), LAPW+lo methods [33, 37, 39]. Exploiting the idea of constructing basis functions for different domains separately, we construct a numerical scheme in this paper. The smoothly varying parts of the wavefunctions away from the atoms are represented by plane waves, the rapidly varying parts near the nuclei are represented by radial basis functions times spherical harmonics, and the approximations inside and outside the spheres are patched together by DG methods.

The DG framework has been widely used in numerical solutions of partial differential equations and investigated theoretically in a lot of works (see, e.g., [3, 6, 11, 26, 41] and references cited therein). For electronic structure calculations, we refer to works by Lin et al. [31, 42], which constructs basis functions adaptively from the local environment and patches them together in global domain by DG methods.

We further present an a priori error analysis of our DG approximations for the linear eigenvalue problems. Thanks to the asymptotic regularity result developed by Flad et al. [22], we can guarantee smoothness of the wavefunctions on the domain $[0,R]\times S^{2}$ in spherical coordinates. Our analysis for DG approximations is also closely related to the technique used in [1, 26, 36, 41]. The main theoretical result in this paper is the following superalgebraic convergence rate under certain assumptions (see Theorem 3.1):

[TABLE]

where $\epsilon>0$ can be arbitrarily small, $\varrho$ denotes the discretization parameters (see (3.1)), and the constant $C_{s}$ depends only on $s$ and the eigenfunctions.

We shall briefly compare our DG method with other existing full-potential/all-electron methods in electronic structure calculations. (a) APW: The augmented plane wave (APW) method [40] introduces basis functions that are plane waves in the interstitial region and radial solutions of Schrödinger equations inside the atomic spheres. A great disadvantage of the APW method is that the basis functions are energy dependent, which results in a nonlinear eigenvalue problem and must be solved separately for each eigenstate by “root tracing” technique [34] or iteration methods. This method is much more complicated to solve than the straightforward linear eigenvalue equations expressed with a fixed basis set, such as plane waves, Gaussians, LAPW (in the following), and our DG schemes. (b) LAPW (+lo): The linearized augmented plane wave (LAPW) method [34, 38] is a linearization of APW, which defines basis functions as linear combinations of a radial solution and its energy derivative evaluated at a chosen fixed energy. This forms a basis set adapted to a particular system that is suitable for calculation of all states in an energy “window”. The accuracy depends heavily on the choice of energy parameter and the width of the energy window under consideration. Although the inclusion of additional variational freedoms (the energy derivatives and sometimes local orbitals (lo) [38]) in the LAPW method facilitates the computation for non-spherical symmetric parts of the potential, there is no proof that it can give solutions of arbitrarily great accuracy for general potentials as our DG scheme. Here we would like to mention a recent work [20] which uses similar ideas as LAPW+lo and may possess systematical convergence. (c) OPW: The orthogonalized plane wave (OPW) method [27] constructs basis functions by orthogonalizing the plane waves to special local functions around each nucleus. The ambiguity of this method arises from inaccuracies of the core wave functions, which are not precise eigenfunctions of the given Hamiltonian. Thus, there is always an uncertainty about the accuracy of OPW results which can not be refined out by more extended calculations. (d) PAW/VPAW: The projector augmented wave (PAW) method [9] replaces the original eigenvalue problem (with singular potential) by a new one with the same eigenvalues but smoother eigenvectors. A slightly different method, called variational projector augmented wave (VPAW), was proposed and analyzed recently [8]. This new method allows for a better convergence with respect to the number of plane waves. But we mention that the PAW method is more of a pseudopotential method.

The rest part of this paper is organized as follows. In Section 2, we set up the model problem and present some regularity results. In Section 3, we introduce a DG discretization scheme, provide a numerical analysis of the convergence and a priori error estimates of the DG approximations. In Section 4, we give some details of the numerical implementations and present some numerical experiments to support our theory. Finally, we give some concluding remarks.

2 Preliminary

Throughout this paper, we shall use $C$ to denote a generic positive constant which may stand for different values at its different occurrences and is independent of finite dimensional subspaces. For convenience, the symbol $\lesssim$ will be used and the notation $A\lesssim B$ means that $A\leq CB$ for some generic positive constant $C$.

Let $\mathcal{R}$ be a discrete periodic lattice of $\mathbb{R}^{3}$, $\Omega$ be the unit cell of the lattice, and $\mathcal{R}^{*}$ be the dual lattice. For simplicity, we take $\Omega=[-\frac{D}{2},\frac{D}{2}]^{3}~{}(D>0)$, $\mathcal{R}=D\mathbb{Z}^{3}$, and $\mathcal{R}^{*}=\frac{2\pi}{D}\mathbb{Z}^{3}$.

For ${\bf k}\in\mathcal{R}^{*}$, we denote by $e_{\bf k}({\bf r})=|\Omega|^{-1/2}e^{i{\bf k\cdot r}}$ the plane wave with wavevector ${\bf k}$. The family $\{e_{\bf k}\}_{{\bf k}\in\mathcal{R}^{*}}$ forms an orthonormal basis set of

[TABLE]

For all $u\in L_{\#}^{2}(\Omega)$, we have

[TABLE]

We introduce the Sobolev spaces of $\mathcal{R}$-periodic functions

[TABLE]

with $s\in\mathbb{R}^{+}$. For $K\in\mathbb{N}^{+}$, we denote the finite dimensional subspace by

[TABLE]

For $v\in H^{s}_{\#}(\Omega)$, the best approximation of $v$ in $\mathcal{V}_{K}$ is $\Pi_{K}v=\sum_{{\bf k}\in\mathcal{R}^{*},|{\bf k}|\leq\frac{2\pi}{D}K}\hat{v}_{\bf k}e_{\bf k}({\bf r})$ for any $H^{t}$-norm ($t\leq s$). The more regularity $v$ has, the faster this truncated series converge to $v$: For real numbers $t$ and $s$ satisfying $t\leq s$, we have that for each $v\in H^{s}_{\#}(\Omega)$,

[TABLE]

As a model problem, we consider the following Schrödinger-type linear eigenvalue problem, which can be viewed as a linearization of (1.1): Find $\lambda\in\mathbb{R}$ and $0\neq u\in H_{\#}^{1}(\Omega)$ such that $\|u\|_{L_{\#}^{2}(\Omega)=1}$ and

[TABLE]

where the bilinear form $a(\cdot,\cdot):H_{\#}^{1}(\Omega)\times H_{\#}^{1}(\Omega)\rightarrow\mathbb{C}$ is given by

[TABLE]

with a $\mathcal{R}$-periodic potential $V\in L^{2}_{\#}(\Omega)$.

To represent the wavefunctions separately in different regions, $\Omega$ is divided into atomic spheres and an interstitial region (see Figure 2.1 (left) for decomposition of a single-atom system, and Figure 2.1 (right) for similar construction of a two-atom system).

For sake of simplicity, we shall restrict our discussions to a single atom located at the origin, the algorithms and analysis of which can be easily generalized to multi-atom systems. Throughout this paper, we shall denote by $\Omega_{\rm out}$ the interstitial region, by $\Omega_{\rm in}$ the sphere centered at the origin with radius $R$, and by $\Gamma$ the spherical surface. We also assume throughout this paper that the potential $V$ equals to $-Z/|{\bf r}|$ in the neighborhood of ${\bf 0}$, and belongs to $C_{\rm loc}^{\infty}(\mathbb{R}^{3}\setminus\mathcal{R})\cap L^{2}_{\#}(\Omega)$.

It was shown in [23, 24, 25] that the exact electron densities are analytic away from the nuclei and satisfy certain cusp conditions at the nuclei. The plane wave approximations can not have as good convergence rate as (2.1) due to the cusps at the nuclear positions. The following lemma concerning the regularity of eigenfunctions of (2.2) is heavily used in our analysis, the proof of which can be referred to [22, Theorem 1,4 and Proposition 1].

Lemma 2.1**.**

If $u$ is an eigenfunction of (2.2), then $u\in H^{s}([0,R]\times S^{2})$ for any $s\in\mathbb{Z}^{+}$.

The following lemma states the relationship between two Sobolev norms.

Lemma 2.2**.**

If $v\in H^{1}(\Omega_{\rm in})\bigcap H^{3}([0,R]\times S^{2})$, then there exits a constant $C_{R}$ depending on $R$ such that

[TABLE]

Proof.

Note that in spherical coordinates

[TABLE]

where the last two terms multiplied by $r^{2}$ is the total angular momentum operator $\Delta_{S^{2}}$ on spherical surface, i.e. the Laplace-Beltrami operator.

Since $v\in H^{1}(\Omega_{\rm in})\bigcap H^{3}([0,R]\times S^{2})$ implies

[TABLE]

we have

[TABLE]

where Green’s formula and (2.4) are used for the second equality. ∎

3 DG discretization

In this section, we construct a DG discretization scheme using radial basis functions times spherical harmonics inside the sphere and plane waves outside. We provide an a priori error analysis of the numerical approximations. Our analysis is composed of three steps: first, we estimate the best approximation errors inside and outside the sphere separately; then we give an error estimate for the DG approximation of the corresponding source problem; finally, we derive an error estimate for the eigenvalue problem. Note that the errors generated by numerical quadratures and linear algebraic solvers are not considered in this paper, which deserve separate investigations.

Note that if $u$ is an eigenfunction of (2.2), then we have from Lemma 2.1 that for any $s>0$, $u|_{\Omega_{\rm in}}\in H^{s}([0,R]\times S^{2})$ in spherical coordinates and $u|_{\Omega_{\rm out}}\in H^{s}(\Omega_{\rm out})$. We can therefore introduce the following space

[TABLE]

with induced norm

[TABLE]

3.1 Approximation space

Denote by $\mathcal{P}_{K}(\Omega_{\rm out})$ the space of functions on $\Omega_{\rm out}$ expanded by plane waves

[TABLE]

and by $\mathcal{B}_{NL}$ the space of functions on $\Omega_{\rm in}$ expanded by radial basis functions times spherical harmonics

[TABLE]

where $\{\chi_{n}\}_{n=0}^{N}$ are basis functions on $[0,R]$. Here, we denote by $\check{u}(r,\theta,\phi)$ the spherical coordinate representations of the function $u({\bf r})$, i.e., $u({\bf r})=\check{u}(r,\theta,\phi)$.

For simplicity, we may assume that the radial basis functions $\{\chi_{n}\}$ are polynomials that span the space of all polynomials of degree no greater than $N$.

Remark 3.1**.**

There are many choices of the radial basis functions $\{\chi_{n}\}$. One can use spectral methods, such as Legendre polynomials, Chebyshev polynomials and Jacobi polynomials, et al. These types of functions form a complete basis set on $[0,R]$, and possess spectral convergence rates for any sufficiently smooth function [15].

Another type of basis set is atomic orbitals [30, 34], such as Gaussians, Slater-type orbitals and numerical solutions of radial Schrödinger equations [30, 34]. These basis functions are closely related to physical problems, and can be very efficient in practice. In some cases, we can rigorously derive their convergence rates [7, 18].

For simplicity of presentations, we shall focus our analysis on the polynomial-type radial basis functions and investigate the atomic orbitals as well in numerical experiments.

Define the finite dimensional space

[TABLE]

Throughout this paper, we may assume that there exists a constant $\varrho$ such that

[TABLE]

which can denote the discretization parameters.

In the following, we shall define some “best” approximations of the function in the interstitial region and atomic spheres respectively.

For the interstitial region, we define the projections $P_{K}:L^{2}_{\#}(\Omega_{\rm out})\rightarrow\mathcal{P}_{K}(\Omega_{\rm out})$ satisfying

[TABLE]

Proposition 3.1**.**

If $u\in H^{s}(\Omega_{\rm out})$, then for $0\leq t<s$, there exists a constant $C$ such that

[TABLE]

Proof.

The proof is similar to that of [19, Proof of Lemma 3.1]. We keep this proof for sake of completeness.

We shall first extend the function $u|_{\Omega_{\rm out}}$ smoothly into the sphere. The wavefunction around the sphere can be represented by

[TABLE]

with $\displaystyle u_{lm}(r)=\int_{0}^{\pi}\sin\theta\int_{0}^{2\pi}u(r,\theta,\phi)Y_{lm}(\theta,\phi)d\phi d\theta$, where we use spherical coordinates ${\bf r}\rightarrow(r,\theta,\phi)$ to express $\tilde{Y}_{lm}({\bf r})=Y_{lm}(\theta,\phi)$ for spherical harmonics on $S^{2}$. Then we can define

[TABLE]

where $\displaystyle\varphi_{lm}(r)=\tau(r)\sum_{n=1}^{s+1}c_{n}u_{lm}(R+\frac{1}{n}(R-r))$ with $\displaystyle\sum_{n=1}^{s+1}(-\frac{1}{n})^{k}c_{n}=1~{}(k=0,1,\cdots,s)$, $\tau\in C^{\infty}([0,R])$ satisfying $\tau=0$ on $[0,\frac{R}{3}]$ and $\tau=1$ on $[\frac{2R}{3},R]$. We observe that $u\in H^{s}(\Omega_{\rm out})$ leads to $\tilde{u}\in H^{s}(\Omega)$ and moreover

[TABLE]

where the constant $\beta$ is only related to $s,~{}R$ and $\|\tau\|_{C^{\infty}([0,R])}$. Let

[TABLE]

we have from (3.6) that

[TABLE]

which completes the proof of (3.2). ∎

For the atomic spheres, we define $P_{N}:H^{1}([0,R])\rightarrow\Psi_{N}\equiv{\rm span}\{\chi_{n}\}_{n=1}^{N}$ satisfying

[TABLE]

and $P_{L}:L^{2}(S^{2})\rightarrow\mathcal{Y}_{L}\equiv{\rm span}\{Y_{lm},0\leq l\leq L,-l\leq m\leq l\}$ satisfying

[TABLE]

For $P_{N}$ and $P_{L}$, we have the following standard estimates (see, e.g., [4, 32])

[TABLE]

for any $0\leq t\leq 1$ and $t<s$. Define the projection $P_{NL}:H^{1}([0,R]\times S^{2})\rightarrow\Psi_{N}\times\mathcal{Y}_{L}$ by $P_{NL}=P_{N}\circ P_{L}$, we have that for $w\in H^{s}([0,R]\times S^{2})$, $0\leq t\leq 1$ and $t<s$,

[TABLE]

Proposition 3.2**.**

If $u\in H^{s}([0,R]\times S^{2})\bigcap H^{1}(\Omega_{\rm in})$, then for $0\leq t\leq 1$ and any $s\geq 3$, there exists a constant $C$ such that

[TABLE]

Proof.

Using (3.7) and Lemma 2.2, we have

[TABLE]

which completes the proof. ∎

3.2 DG approximations of the source problem

We shall discuss the DG discretization for the source problem and our analysis is related to the framework in [1].

For vector-valued ${\bf w}$ and scalar-valued function $u$ which are not continuous on the spherical surface $\Gamma$, we define the jumps by

[TABLE]

and the averages by

[TABLE]

where ${\bf w}^{\pm}$ and $u^{\pm}$ are traces of ${\bf w}$ and $u$ on $\Gamma$ taken from inside and outside the sphere, ${\bf n}^{\pm}$ are the normal unit vectors.

Since $\mathcal{S}^{K}_{NL}(\Omega)$ is a finite dimensional space, there exists a constant $\gamma_{\varrho}$ depending on $\varrho$ such that the following inverse estimate holds

[TABLE]

In our analysis, we assume

[TABLE]

We are not able to justify (3.10) rigorously, however, we provide some numerical experiments in Appendix A to show that it could be true. Then we get from the “interpolation” arguments (see Appendix A) that

[TABLE]

for some $\epsilon\in(0,1)$.

We then define the bilinear form $a^{\rm DG}(\cdot,\cdot):\left(\mathcal{S}^{K}_{NL}(\Omega)\cup H^{1}_{\#}(\Omega)\right)\times\left(\mathcal{S}^{K}_{NL}(\Omega)\cup H^{1}_{\#}(\Omega)\right)\rightarrow\mathbb{C}$ by

[TABLE]

where $\sigma=C_{\sigma}\varrho^{2+2\epsilon}$ is the discontinuity-penalization parameter with a constant $C_{\sigma}$ independent of the discretization.

Note that there are many other types of DG formulations (see, e.g., [1, 3]), and (3.12) is the classical symmetric interior penalty (SIP) method [1, 2, 31].

We further define the broken Sobolev space

[TABLE]

equipped with the following DG-norm

[TABLE]

Lemma 3.1**.**

If $C_{\sigma}$ is sufficiently large, then there exist constants $\alpha,\beta>0$ such that

[TABLE]

Proof.

Using Hölder inequality, Sobolev’s embedding theorem and Young’s inequality, we obtain that

[TABLE]

where $\delta>0$ is arbitrarily small. Hence we have

[TABLE]

with constants $C,~{}b>0$. Moreover, we have

[TABLE]

Using (3.9), (3.15) and (3.16), we can derive (3.14) and complete the proof. ∎

For simplicity, we can take $\beta=0$. Note that $a^{\rm DG}_{\beta}(u,v)=a^{\rm DG}(u,v)+\beta(u,v)$ makes this true for $\beta>0$.

Define the solution operators

[TABLE]

and

[TABLE]

Proposition 3.3**.**

Assume that (3.10) is true and $C_{\sigma}$ is sufficiently large. If $Tf\in H^{s}(\Omega_{\rm out})\oplus H^{s}([0,R]\times S^{2})$ for $f\in L^{2}_{\#}(\Omega)$ and $s\geq 3$, then there exists a constant $C$ such that

[TABLE]

Proof.

Denote $w=Tf$ and $w^{\rm DG}=T^{\rm DG}f$. Define the projection $\mathcal{P}u=P_{K}u|_{\Omega_{\rm out}}+P_{NL}u|_{\Omega_{\rm in}}$, We decompose the error $e=w-w^{\rm DG}$ as $e=\eta+\xi$, where $\eta=w-\mathcal{P}w$ and $\xi=\mathcal{P}w-w^{\rm DG}$. With simple calculations, we can easily obtain that $a^{\rm DG}(w,\xi)=(f,\xi)$, which leads to the property that $a^{\rm DG}(w-w^{\rm DG},\xi)=0$. Using (3.14) and the property, we have

[TABLE]

Thus we deduce that

[TABLE]

where

[TABLE]

Since $V\in L_{\#}^{2}(\Omega)$, we have

[TABLE]

Using the trace inequality, $I_{2}$ can be estimated by

[TABLE]

Similarly, $I_{3}$ can be estimated by

[TABLE]

Collecting (3.18) and the error bounds (3.19) to (3.21), we have

[TABLE]

We obtain from (3.2) and (3.8) that if $u\in H^{s}(\Omega_{\rm out})\oplus H^{s}([0,R]\times S^{2})~{}(s\geq 3)$, then

[TABLE]

which together with

[TABLE]

leads to

[TABLE]

Then we can derive (3.17) by using (3.11). ∎

3.3 DG approximations of the eigenvalue problem

We construct DG methods for eigenvalue problem (2.2): Find $\lambda^{\rm DG}\in\mathbb{R}$ and $u^{\rm DG}\in\mathcal{S}^{K}_{NL}(\Omega)$, such that $\|u^{\rm DG}\|_{L^{2}_{\#}(\Omega)}=1$ and

[TABLE]

Note that (2.2) and (3.22) are equivalent to $\lambda Tu=u$ and $\lambda^{\rm DG}T^{\rm DG}u^{\rm DG}=u^{\rm DG}$, respectively.

Denote by $\sigma(T)$ the spectrum and $\rho(T)$ the resolvent set of the solution operator $T$. For any $z\in\mathbb{C}$ in $\rho(T)$, we define the resolvent operator $R_{z}(T)=(z-T)^{-1}$. Let $\lambda^{-1}$ be an eigenvalue of $T$ and $\gamma$ be a circle in the complex plane that is centered at $\lambda^{-1}$ and does not enclose any other point of $\sigma(T)$.

Define the following operators with contour integrations:

[TABLE]

If $\varrho$ is sufficiently large, then $\mathscr{E}$ and $\mathscr{E}^{\rm DG}$ are the spectral projectors of $T$ and $T^{\rm DG}$ relative to $\lambda^{-1}$, respectively (see [36]).

Define the distances

[TABLE]

Using proposition 3.3 and similar arguments as those in [1], we have the following convergence results (including non-pollution and completeness) for DG eigenvalues and eigenspaces.

Remark 3.2**.**

Let $A\subset\mathbb{R}$ be an open set containing $\sigma(T)$. If $C_{\sigma}$ and $\varrho$ are sufficiently large, then $\sigma(T^{\rm DG})\subset A$. Moreover, for all $z\in\sigma(T)$, we have

[TABLE]

In addition, we have

[TABLE]

where $\mathcal{R}$ denotes the range.

Now we can derive the following a priori error estimate for DG approximations.

Theorem 3.1**.**

Assume that (3.10) is true and $C_{\sigma}$ is sufficiently large. Let $\lambda$ be an eigenvalue of (2.2) with algebraic multiplicity m. Then for $\varrho$ sufficiently large, there exist $m$ eigenpairs $(\lambda_{i}^{\rm DG},u_{i}^{\rm DG})~{}(i=1,2,\cdots m)$ of (3.22) such that

[TABLE]

where the constant $C_{s}$ depends only on $\lambda,u_{i}$ and $s$.

Proof.

Note that for $f\in L_{\#}^{2}(\Omega)$, $Tf\in\widetilde{H}^{2}(\Omega)$ and $\|Tf\|_{\widetilde{H}^{2}(\Omega)}\leq C\|f\|_{L_{\#}^{2}(\Omega)}$ (see [21], p.257, Thm. 9 and (8.137)).

For $f\in H_{\#\delta}(\Omega)$, we have from Proposition 3.3 that

[TABLE]

which implies

[TABLE]

Using (3.24) and [36, Theorem 1], we have the convergence of the eigenvalues and

[TABLE]

Then it is only necessary for us to estimate the right-hand side of (3.25).

Using proposition 3.3, the regularity result Lemma 2.1 and the fact $Tv=\lambda^{-1}v$ for $v\in\mathcal{R}(\mathcal{E})$, we have that for any $s\geq 3$,

[TABLE]

where $C_{s}$ is a constant depending only on $\mathcal{R}(\mathscr{E}),\lambda$ and $s$ .

It is apparent from (3.25) and (3.26) that

[TABLE]

Let $m$ and $m_{\varrho}$ be the dimensions of $\mathcal{R}(\mathcal{E})$ and $\mathcal{R}(\mathcal{E}^{\rm DG})$, respectively. Then, (3.27) indicates that, for $\varrho$ large enough, $m=m_{\varrho}$ (see [29, p.200]) and there exist $m$ eigenfunctions $u_{i}^{\rm DG}\in\mathcal{R}(\mathcal{E}^{\rm DG})$ and $m$ eigenpairs $(\lambda_{i}^{\rm DG},u_{i}^{\rm DG})~{}(i=1,2,\cdots m)$ satisfying (3.22). Moreover, according to the definition of distance $\mathscr{D}(X,Y)$, we can find $u_{i}\in\mathcal{R}(\mathcal{E})$ and $\|u_{i}\|_{L^{2}_{\#}(\Omega)}=1$ such that

[TABLE]

This completes the proof of error estimates for eigenfunctions.

For eigenvalues, we obtain by a simple calculation that

[TABLE]

with the consistency error

[TABLE]

Using (3.26) to (3.30), we obtain

[TABLE]

where the constant $C_{s}$ depends only on $\lambda,u_{i}$ and $s$. ∎

Remark 3.3**.**

We emphasize that our result works not only for the case of single eigenvalue ($m=1$), but also for general cases of multiple eigenvalue ($m>1$).

Remark 3.4**.**

It is shown in many cases that the convergence rate of finite dimensional approximations under a weaker norm is faster than that under a stronger norm (see, e.g., [5, 16]). By making this assumption for our DG approximations, for example,

[TABLE]

it may be true from (3.3) and (3.30) that the eigenvalue approximations have better convergence rate than that of eigenfunctions.

Remark 3.5**.**

Within the framework of Kohn-Sham density functional theory, one has to solve the nonlinear eigenvalue problem (1.1) with a SCF iteration. Using our DG discretizations, the linear eigenvalue problem (3.22) is solved at each iteration step and complex mixing schemes such as Roothaan, level-shifting and DIIS algorithms (see, e.g., [12, 30]) are used to achieve convergence.

If the exchange-correlation potential $V_{\rm xc}$ is sufficiently smooth and the trial state (from previous DG approximations) $\tilde{\Phi}\in(\mathcal{S}^{K}_{NL}(\Omega))^{N_{\rm e}}$, then we have from similar arguments as those in [22] that the eigenfunctions $\{\phi_{i}\}_{i=1,\cdots,N_{\rm e}}$ of $H_{\tilde{\Phi}}$ belong to $C_{\rm loc}^{\infty}(\mathbb{R}^{3}\setminus\mathcal{R})\cup C^{\infty}([0,R]\times S^{2})$. This regularity together with the analysis in Theorem 3.1 gives spectral convergence rates for DG approximations of the (linear) eigenvalue problem in each SCF iteration step.

Note that we have not obtained a priori error estimates for approximations of nonlinear eigenvalue problems but only for linearized equations in SCF iterations. We refer to [13, 17] for numerical analysis of nonlinear eigenvalue problems.

4 Numerical experiments

In this section, we will present some details for implementing our DG scheme, and some numerical experiments in electronic structure calculations.

4.1 Hamiltonian matrix elements

With our DG scheme, we can discretize the continuous eigenvalue problem into a (finite dimensional) matrix generalized eigenvalue problem

[TABLE]

where $\hat{u}_{i}$ are eigenvectors that correspond to the DG approximations $u_{i}^{\rm DG}$. We shall explain in the following how the matrix elements of $H$ and $M$ are generated.

For basis functions ${\bf p}$ and ${\bf q}$, We divide the integrals for overlap matrix $M_{\bf pq}$ and stiff matrix $H_{\bf pq}$ into three parts. The scattering identity (see, e.g., [35])

[TABLE]

with $k=|{\bf k}|$, is heavily used to bridge the gap between plane waves and spherical harmonics.

For ${\bf p,q}\in\mathcal{P}_{K}(\Omega_{\rm out})$, we have

[TABLE]

where $U({\bf k})$ is the Fourier transform of the step function with 0 inside the sphere and 1 outside

[TABLE]

with $k=|{\bf k}|$ and $j_{l}$ the $l$th spherical Bessel function. Similarly, we have from (3.12) that

[TABLE]

where

[TABLE]

with $k=|{\bf k}|$ and the potential inside the sphere expanded by $V({\bf r})=\sum_{lm}v_{lm}(r)\tilde{Y}_{lm}({\bf r})$. The discontinuity and penalization term $\mathfrak{D}^{a}_{\bf pq}$ in (4.4) is given by

[TABLE]

with $k_{\bf p}=|{\bf k}_{\bf p}|$ and $k_{\bf q}=|{\bf k}_{\bf q}|$. Note that the first term of (4.5) is obtained by fast Fourier transform (FFT) and the second term is calculated by numerical integrations.

For ${\bf p,q}\in\mathcal{B}_{NL}(\Omega_{\rm in})$, we have from the orthogonality of $Y_{lm}$ on the surface that

[TABLE]

and

[TABLE]

where the potential inside the sphere is expanded by $V({\bf r})=\sum_{\hat{l}\hat{m}}v_{\hat{l}\hat{m}}(r)\tilde{Y}_{\hat{l}\hat{m}}({\bf r})$ and $G$ is the integral of three spherical harmonics that can be written in terms of Gaunt coefficients (see, e.g., [34]). The discontinuity and penalization term $\mathfrak{D}^{b}_{\bf pq}$ is

[TABLE]

For ${\bf p}\in\mathcal{P}_{K}(\Omega_{\rm out}),~{}{\bf q}\in\mathcal{B}_{NL}(\Omega_{\rm in})$, we have

[TABLE]

and

[TABLE]

Since we use a symmetric DG scheme, the elements for ${\bf p}\in\mathcal{B}_{NL}(\Omega_{\rm in}),~{}{\bf q}\in\mathcal{P}_{K}(\Omega_{\rm out})$ can be obtained immediately.

Combining (4.2) – (4.1), we can obtain the matrices $H$ and $M$, and further solve the matrix eigenvalue problems by linear eigensolvers.

4.2 Numerical results

All the numerical results are presented by atomic units (a.u.). When we test the convergence with respect to one parameter (say, $K$, $N$ or $L$), the other two parameters are fixed and chosen to be sufficiently large.

Example 1. (linear problem for a single-atom system) Consider the linear eigenvalue problem: Find $\lambda\in\mathbb{R}$ and $u\in H^{1}_{\#}(\Omega)$ such that

[TABLE]

with $\Omega=[-5,5]^{3}$ and $\displaystyle V({\bf r})=-\frac{4\pi}{|\Omega|}\sum_{{\bf k}\in\mathcal{R}^{*},{\bf k}\neq 0}\frac{e^{i{\bf k}\cdot{\bf r}}}{|{\bf k}|^{2}}$ . Note that the potential $V$ can be viewed as a periodized version of the potential $\displaystyle-\frac{1}{|\bf r|}$ for a hydrogen atom. It is periodic and sufficiently smooth everywhere except at the origin. Then due to Lemma 2.1, the error estimates in Theorem 3.1 hold.

We first compare the numerical errors of the lowest eigenvalue approximations by plane waves and our DG methods (see Figure 4.4), from which we observe that the DG approximations converge much faster. We compare the eigenfunctions along the $x$-axis obtained by plane waves and DG methods (see Figure 4.4). We observe that the DG approximations can capture the cusp at the nuclear position while that plane waves can not. For a more precise comparison, when the required accuracy is $10^{-1}$ (for the first eigenvalue), the DG method needs around $50$ degrees of freedom (DOFs) while the plane waves method need about $60$ DOFs; when the required accuracy is $10^{-2}$, the DG method needs around 300 DOFs while the plane waves method needs more than 1100 DOFs.

We further show the convergence rates of the eigenvalue errors with respect to plane wave truncations $K$ (see Figure 4.4), and observe exponential decay for different sizes of atomic spheres. We find a slightly faster convergence rate of the numerical errors (in Figure 4.4) with a bigger size of atomic sphere. The reason is that the eigenfunctions are less varying outside a larger atomic sphere. However, we see that the choice of $R$ does not affect the numerical simulations significantly. In practical simulations, we could choose relatively large atomic spheres as long as they do not overlap.

We also present the numerical errors with respect to the orders of radial basis functions (see Figure 4.4), and compare the polynomials with Slater-type atomic orbitals

[TABLE]

with $\eta$ a fixed parameter. We observe that although a high accuracy can be obtained by Slater-type atomic orbitals with very few degrees of freedom, a better systematically convergence rate is achieved by polynomials.

Example 2. (linear problem for a two-atom system) Consider the linear eigenvalue problem for a two-atom system: Find $\lambda\in\mathbb{R}$ and $u\in H^{1}_{\#}(\Omega)$ such that

[TABLE]

where $\Omega=[-5,5]^{3}$ and $\displaystyle V_{j}({\bf r})=-\frac{4\pi}{|\Omega|}\sum_{{\bf k}\in\mathcal{R}^{*},{\bf k}\neq 0}\frac{1}{|{\bf k}|^{2}}e^{i{\bf k}\cdot({\bf r-R_{j}})}~{}(j=1,2)$ with ${\bf R_{1}}$ and ${\bf R_{2}}$ the positions of atoms.

We first compare the numerical errors of the the lowest eigenvalue approximations by plane waves and our DG methods. We observe a much better convergence rate with respect to $K$ in Figure 4.10 and a more accurate capture of the eigenfunction cusp by our DG approximation in Figure 4.10. For a more precise comparison, when the required accuracy is $10^{-1}$, the DG method needs around 100 DOFs while the plane waves method need about 30 DOFs; when the required accuracy is $10^{-2}$, the DG method needs around 400 DOFs while the plane waves method needs more than 1000 DOFs.

We then show the convergence rates of the eigenvalue errors with respect to plane wave truncations $K$ (see Figure 4.4), and angular momentum truncation $L$ (see Figure 4.10) We observe exponential decay with respect to both $K$ and $L$ for different sizes of atomic spheres. In this example, when the radii of atomic spheres are larger than 0.5, increasing the size of spheres does not improve the convergence rate significantly (see Figure 4.10). From the comparisons of different radii, we get the same conclusion as Example 1 that the choice of $R$ is not important in our scheme, and one can simply choose relative large radii such that the atomic spheres do not overlap.

We further show the convergence rates of numerical errors for the lowest 3 eigenvalues, and observe exponential decay of the numerical errors with respect to $K$ (see Figure 4.10). This supports our theory.

Finally, we test the effect of penalty parameter $C_{\sigma}$. In our DG scheme, the penalty constant $C_{\sigma}$ plays an important role to guarantee stability. The errors with respect to different choices $C_{\sigma}$ are shown in Figure 4.10. We observe that the DG method can be stable and accurate in a large range of values beyond a certain threshold value. Similar discussions for the penalty parameter can also be found in [31].

Example 3. (simulation of a helium atom) Consider the following nonlinear eigenvalue problem: Find $\lambda\in\mathbb{R}$ and $u\in H^{1}_{\#}(\Omega)$ such that

[TABLE]

with the external potential $\displaystyle V_{\rm ext}({\bf r})=-\frac{8\pi}{|\Omega|}\sum_{{\bf k}\in\mathcal{R}^{*},{\bf k}\neq 0}\frac{e^{i{\bf k}\cdot{\bf r}}}{|{\bf k}|^{2}}$, the Hartree potential $\displaystyle V_{\rm H}[\rho]=\int_{\mathbb{R}^{3}}\frac{\rho({\bf y})}{|\cdot-{\bf y}|}d{\bf y}$ and $\rho=2u^{2}$. Note that the exchange-correlation potential $V_{\rm xc}[\rho]$ is ignored here from a standard Kohn-Sham DFT model. We present the eigenfunction along the $x$-axis (see Figure 4.12) and the convergence rates of numerical errors for different sizes of atomic spheres (see Figure 4.12). We observe exponential decay of the numerical errors with respect to $K$ even for the nonlinear problem.

5 Concluding remarks

In this paper, we construct a discontinuous Galerkin scheme for full-potential electronic structure calculations. It exploits the idea of augmented plane wave method which approximates the wavefunction in some ways “the best of two worlds”. The smoothly varying parts of the wavefunctions between the atoms are represented by plane waves, the rapidly varying parts near the nuclei are represented by radial atomic functions times spherical harmonics inside a sphere around each nucleus, and these two parts are patched together by discontinuous Galerkin scheme. We demonstrate a priori error estimate of this approximation to illustrate the accuracy and efficiency of this scheme, and provide some numerical experiments to support the theory.

Besides the accuracy and efficiency we have shown in this paper, the discontinuous Galerkin scheme is also flexible and economical for adaptive procedures since the nonconformity results assuredly in limiting the contamination only to the subdomain where refinement is needed. The a posteriori error analysis and the adaptive algorithm will be addressed in our future works.

Acknowledgement

We are grateful to Reinhold Schneider in Technische Universität Berlin for his help in this work and for inspiring discussions.

Appendix Appendix A Inverse estimates on the surface

In this appendix, we shall provide the numerical tests to support the inverse estimate assumption (3.10) and the “interpolation” arguments to obtain (3.11).

Let

[TABLE]

Consider the largest eigenvalue $\lambda_{\varrho,{\rm max}}$ of the following discrete eigenvalue problem on the spherical surface $\Gamma$: Find $\lambda_{\varrho}\in\mathbb{R}$ and $u_{\varrho}\in\widetilde{V}_{\varrho}(\Gamma)$ such that

[TABLE]

We perform numerical simulations for (A.1) and present the scalings of $\lambda_{\varrho,{\rm max}}$ (with respect to the discretizations $\varrho$) in Figure A.1 for different sizes of atomic spheres.

We observe from the numerics that $\lambda_{\varrho,{\rm max}}=C_{R}\varrho^{4}$ for all different radii $R$, which together with the fact

[TABLE]

implies

[TABLE]

which supports the inverse estimate assumption (3.10).

We then use the “interpolation” between two spaces $L^{2}(\Gamma)$ and $H^{1}(\Gamma)$. For any $u\in L^{2}(\Gamma)$ and $t>0$, define

[TABLE]

and the $H^{\frac{1}{2}}$-norm through “interpolation” [10]: $\displaystyle\|u\|_{H^{\frac{1}{2}}(\Gamma)}=\left(\int_{0}^{\infty}\frac{K^{2}(t,u)}{t^{2}}dt\right)^{\frac{1}{2}}$. For any $0<\alpha<1$, we have

[TABLE]

To see (A.4), we have $K(t,u)\leq t\|u\|_{H^{1}(\Gamma)}$ by taking $v=u$ in (A.3) and $K(t,u)\leq\|u\|_{L^{2}(\Gamma)}$ by choosing $v=0$. Using these two inequalities, we can derive

[TABLE]

Combining (A.2) and (A.4), we can obtain the inverse estimate (3.11) .