A Nonnested Augmented Subspace Method for Kohn-Sham Equation

Density functional theory (DFT) is one of the major breakthroughs in quantum physics and quantum chemistry. In the framework of density functional theory, the many-body problem can be simplified as the motion of electrons without interaction in the effective potential field, which includes the Coulomb potential of the external potential field, the Hartree potential generated by the interaction between electrons, and the exchange-correlation potential for the nonclassical Coulomb interaction. For the exchange-correlation potential, it has always been a difficulty in density functional theory. At present, there is no exact analytical expression for the exchange-correlation potential. Generally, it is approximately described by the local density approximation (LDA), the local spin density approximation (LSDA) and the generalized gradient approximation (GGA), etc. The Kohn-Sham equation is one of the most important models in the calculation of electronic structure, which transforms the wavefunctions (in $\mathbb{R}^{3N}$ space) describing $N$ particle states into density function (in $\mathbb{R}^{3}$ space), so as to reduce the degrees of freedom of equation and reduce the computational work of numerical simulation. The idea of Kohn-Sham equation can be traced back to Thomas-Fermi model in 1927. The model gave a primary description of the electronic structure of atoms. The strict theoretical analysis was described by Hohenberg and Kohn in [26], which proved the correctness and feasibility of using single electron orbit to replace the wavefunctions describing multiple electrons.

So far, lots of numerical methods for solving Kohn-Sham equation have been developed. For instance, plane-wave method [10, 35] is the most popular method in the computational quantum chemistry community. Since the basis function is independent from the ionic position, plane-wave method has advantage in calculating intermolecular force. Combined with the pseudopotential method, plane-wave method plays an important role in the study of the ground and excited states calculations, and geometry optimization of the electronic structures. Although the plane-wave method is the most popular one in the computational quantum chemistry community, it is inefficient in solving non-periodic systems like molecules, nano-clusters, or materials systems with defects, etc. Furthermore, the plane-wave method uses the global basis which significantly affect the scalability on parallel computing platforms. The atomic-orbital-type basis sets [21, 22] are also widely used for simulating materials systems such as molecules and clusters. However, they are well-suited only for isolated systems with special boundary conditions. It is difficult to develop a systematic basis-set for all materials systems. Thus over the past decades, more and more attentions are attracted to develop efficient and scalable grid-based methods such as the finite element method for electronic structure calculations. The advantages of grid-based methods include that it can use unstructured meshes and local basis sets, hence owing high scalability on parallel computing platforms. So far, the applications of the finite element method in solving Kohn-Sham equation have been studied systematically. We refer to [6, 8, 11, 33, 36, 39, 40, 42, 43, 47, 48, 49] and references therein for a comprehensive overview.

In this paper, we propose a nonnested augmented subspace method to solve the Kohn-Sham equation based on the moving mesh adaptive technique and the augmented subspace method. Since the Coulomb potential and Hartree potential have strong singularities, the adaptive mesh refinement is a competitive strategy to improve solving efficiency. Adaptive methods mainly include $h$-adaptive method, $p$-adaptive method and $r$-adaptive method. The $h$-adaptive method uses a posteriori error indicator for local refinement, $p$-adaptive method uses higher-order polynomials on local mesh elements, $r$-adaptive method (moving mesh method) uses control function or metric tensor (also known as metric matrix) for mesh redistribution. In this paper, we adopt the $r$-adaptive method to generate a series of nonnested adaptive finite element spaces. The $r$-adaptive method uses the control function or metric matrix to guide the mesh movement, which can improve the accuracy by moving the grid nodes to the area with large errors. This method can be traced back to Alexander’s work in [3], and then Miller’s work [37, 38] promoted the development of $r$-adaptive method. In 2001, Li and Tang etc [31, 32] proposed the moving mesh method based on harmonic mapping. Later, Di and Li etc [18] applied the moving mesh method based on harmonic mapping to solve the incompressible Navier-Stokes equation. More results about moving mesh method can be found in [19, 20, 53, 54] and the references therein.

Another technique adopted in this paper is the multilevel correction method (augmented subspace method) [27, 34, 51, 52]. The traditional multilevel correction method can only be performed on the nested multilevel space sequence. In this paper, we develop a nonnested multilevel correction technique for Kohn-Sham equation on a nonnested mesh sequence generated by the moving mesh technique. Different from the classical self-consistent field iterative algorithm which requires to solve the Kohn- Sham equation directly in each adaptive finite element space, our algorithm transforms the Kohn-Sham equation into some linear boundary value problems of the same scale on each level of the adaptive refinement mesh sequence, and then the wavefunctions derived from the linear boundary value problems are corrected by solving a small-scale Kohn-Sham equation defined in a low-dimensional augmented subspace. Since the new algorithm avoids solving large-scale Kohm-Sham equation directly, a significant improvement for the efficiency can be obtained. In addition, the adaptive mesh is produced using the moving mesh technique according to the singularity of the wavefunctions, which can dramatically generate a high accuracy with less computational work. In our algorithm, the approximate wavefunctions are used as the adaptive function, and the modified Hessian matrix of the density function is used as the metric matrix to redistribute the mesh. Through the moving mesh adaptive technique, the redistributed mesh is almost optimal. By combining the moving mesh technique and the nonnested augmented subspace method, the solving efficiency for Kohn-Sham equation can be significantly improved.

The remainder of this paper is as follows. Section 2 introduces Kohn-Sham equation and its finite element discretization method. Section 3 introduces the nonnested augmented subspace method for Kohn-Sham equation by using the moving mesh technique and the multilevel correction method. In Section 4, we propose the convergence analysis and the estimate of computational work for the presented algorithm. In Section 5, a number of numerical examples are demonstrated to verify the efficiency and accuracy of the proposed method.

To describe the Kohn-Sham equation and its finite element discretization, we first introduce some notation. Following [1], we use $W^{s,p}(\Omega)$ to denote Sobolev spaces, and use $\|\cdot\|_{s,p,\Omega}$ and $|\cdot|_{s,p,\Omega}$ to denote the associated norms and seminorms, respectively. In case $p=2$, we denote $H^{s}(\Omega)=W^{s,2}(\Omega)$ and $H_{0}^{1}(\Omega)=\{v\in H^{1}(\Omega):v|_{\partial\Omega}=0\}$, where $v|_{\partial\Omega}=0$ is in the sense of trace, and denote $\|\cdot\|_{s,\Omega}=\|\cdot\|_{s,2,\Omega}$. For convenience, the symbols $x_{1}\lesssim y_{1}$, $x_{2}\gtrsim y_{2}$ and $x_{3}\approx y_{3}$ are used to represent $x_{1}\leq C_{1}y_{1}$, $x_{2}\geq c_{2}y_{2}$, and $c_{3}x_{3}\leq y_{3}\leq C_{3}y_{3}$, respectively, where $C_{1}$, $c_{2}$, $c_{3}$, and $C_{3}$ denote some mesh-independent constants. In this paper, $\Omega$ is dropped from the subscript of the norm for simplicity.

Let $\mathcal{H}:=(H_{0}^{1}(\Omega))^{N}$ be the Hilbert space with the inner product

where $\Omega\subset\mathbb{R}^{3}$. For any $\Psi\in\mathcal{H}$ and a subdomain $\omega\subset\Omega$, we define $\rho_{\Psi}=\sum_{i=1}^{N}|\psi_{i}|^{2}$ and

Let $\mathcal{Q}$ be a subspace of $\mathcal{H}$ with orthonormality constraints:

where $\Phi^{T}\Psi=\big{(}\int_{\Omega}\phi_{i}\psi_{j}dx\big{)}_{i,j=1}^{N}\in\mathbb{R}^{N\times N}$.

We consider a molecular system consisting of $M$ nuclei with charges $\{Z_{1}$, $\cdots$, $Z_{M}\}$ and locations $\{R_{1}$, $\cdots$, $R_{M}\}$, respectively, and $N$ electrons in the non-relativistic and spin-unpolarized setting. The general form of Kohn-Sham energy functional can be demonstrated as follows

for $\Psi=(\psi_{1},\psi_{2},\cdots,\psi_{N})\in\mathcal{H}$. Here, $V_{ext}$ is the Coulomb potential defined by $V_{ext}=-\sum_{k=1}^{M}Z_{k}/|x-R_{k}|$, $D(\rho_{\Phi},\rho_{\Phi})$ is the Hartree energy defined by

and $E_{xc}(t)$ is some real function over $[0,\infty)$ denoting the exchange-correlation energy.

The ground state of the system is obtained by solving the minimization problem:

and we refer to [10, 16] for the existence of a minimizer under some conditions.

The Euler-Lagrange equation corresponding to the minimization problem (2.5) is the well-known Kohn-Sham equation: Find $(\Lambda,\Psi)\in\mathbb{R}^{N}\times\mathcal{H}$ such that

where $H_{\Psi}$ is the Kohn-Sham Hamiltonian operator defined by

with $V_{Har}(\rho_{\Psi})=r^{-1}*\rho_{\Psi}$, $V_{xc}(\rho_{\Psi})=E^{\prime}_{xc}(\rho_{\Psi})$, $\Lambda=(\lambda_{1},\cdots,\lambda_{N})$ and $\lambda_{i}=(H_{\Psi}\psi_{i},\psi_{i})$. The variational form of the Kohn-Sham equation can be described as follows: Find $(\Lambda,\Psi)\in\mathbb{R}^{N}\times\mathcal{H}$ such that

Now, let us define the finite element discretization of (2.8). First we generate a shape regular decomposition $\mathcal{T}_{h}$ of the computing domain $\Omega$. Then, based on the mesh $\mathcal{T}_{h}$, we construct the linear finite element space denoted by $S_{h}\subset H_{0}^{1}(\Omega)$. Define $V_{h}=(S_{h})^{N}\subset\mathcal{H}$. Then the discrete form of (2.8) can be described as follows: Find $(\bar{\Lambda}_{h},\bar{\Phi}_{h})\in\mathbb{R}^{N}\times V_{h}$ such that

with $\bar{\lambda}_{i,h}=(H_{\bar{\Psi}_{h}}\bar{\psi}_{i,h},\bar{\psi}_{i,h})$.

For simplicity and generality, some assumptions are given for the error analysis of the proposed algorithm in this paper. We would like to mention that the following assumptions are satisfied by many practical physical models. For detail of the assumptions, please refer to [9, 10, 14, 15].

• •

  Assumption A: $|V_{xc}(t)|+|tV_{xc}^{\prime}(t)|\in\mathcal{P}(p_{1},(c_{1},c_{2}))$ for some $p_{1}\in[0,2]$, where $\mathcal{P}\big{(}p,(c_{1},c_{2})\big{)}$ denotes the following function set:

  $\displaystyle\mathcal{P}(p,(c_{1},c_{2}))=\Big{\{}f:\exists a_{1},a_{2}\in\mathbb{R}\text{ such that }c_{1}t^{p}+a_{1}\leq f(t)\leq c_{2}t^{p}+a_{2},\ \forall t\geq 0\Big{\}}$

  (2.10)

  with $c_{1}\in\mathbb{R}$ and $c_{2}$, $p\in[0,\infty)$.

• •

  Assumption B: There exists a constant $\alpha\in(0,1]$ such that $|V_{xc}^{\prime}(t)|+|tV_{xc}^{\prime\prime}(t)|\lesssim 1+t^{\alpha-1}$.

• •

  Assumption C: Let $(\Lambda,\Psi)$ be a solution of (2.8). When the mesh size $h$ is small enough, the discrete solution $(\bar{\Lambda}_{h},\bar{\Psi}_{h})\in\mathbb{R}^{N}\times V_{h}$ satisfies the following estimate

  	$\displaystyle\|\Psi-\bar{\Psi}_{h}\|_{1}$	$\displaystyle\lesssim$	$\displaystyle\delta_{h}(\Psi),$		(2.11)
  	$\displaystyle\|\Psi-\bar{\Psi}_{h}\|_{0}+|\Lambda-\bar{\Lambda}_{h}|$	$\displaystyle\lesssim$	$\displaystyle r(V_{h})\|\Psi-\bar{\Psi}_{h}\|_{1},$		(2.12)

  where $r(V_{h})\ll 1$ and

  $$\delta_{h}(\Psi):=\inf\limits_{\Phi\in V_{h}}\|\Psi-\Phi\|_{1}.$$

In this section, we design the nonnested augmented subspace method for Kohn-Sham equation based on the moving mesh adaptive technique and the augmented subspace method. In the first two subsections, we introduce some computing techniques including the solving process for exchange-correlation potential and Hartree potential, and the detailed process of moving mesh technique. Then next, we combine these techniques to generate the nonnested augmented subspace method for Kohn-Sham equation.