# Towards chemical accuracy using a multi-mesh adaptive finite element method
in all-electron density functional theory

Yang Kuang<EMAIL_ADDRESS>Yedan Shen  Guanghui Hu  School of Mathematics
and Statistics & Center for Mathematics and Interdisciplinary Science (CMIS),
Guangdong University of Technology, China  School of Mathematics and
Information Science, Guangzhou University, China  Hunan Key Laboratory for
Computation and Simulation in Science and Engineering, Xiangtan University,
China.  Department of Mathematics & Guangdong-Hong Kong-Macao Joint Laboratory
for Data-Driven Fluid Mechanics and En- gineering Applications, University of
Macau, China  Zhuhai UM Science and Technology Research Institute, Zhuhai,
China.

###### Abstract

Chemical accuracy serves as an important metric for assessing the
effectiveness of the numerical method in Kohn–Sham density functional theory.
It is found that to achieve chemical accuracy, not only the Kohn–Sham
wavefunctions but also the Hartree potential, should be approximated
accurately. Under the adaptive finite element framework, this can be
implemented by constructing the _a posteriori_ error indicator based on
approximations of the aforementioned two quantities. However, this way results
in a large amount of computational cost. To reduce the computational cost, we
propose a novel multi-mesh adaptive method, in which the Kohn–Sham equation
and the Poisson equation are solved in two different meshes on the same
computational domain, respectively. With the proposed method, chemical
accuracy can be achieved with less computational consumption compared with the
adaptive method on a single mesh, as demonstrated in a number of numerical
experiments.

###### keywords:

Chemical accuracy; Kohn–Sham equation; Adaptive finite element method; Multi-
mesh method

##  1 Introduction

The Kohn–Sham density functional theory plays an important role in the area of
quantum physics, condensed matter, and computational chemistry. For a quantum
system containing $N_{nuc}$ nuclei and $N_{ele}$ electrons, the ground state
can be obtained by solving the lowest $N_{occ}$ eigenpairs
${(\varepsilon_{i},\psi_{i})}_{i=1,\dots,N_{occ}}$ from the Kohn–Sham equation

$\left(-\frac{1}{2}\nabla^{2}+V_{\mathrm{ext}}(\mathbf{r})+V_{\mathrm{Har}}(\mathbf{r})+V_{\mathrm{xc}}(\mathbf{r})\right)\psi_{i}(\mathbf{r})=\varepsilon_{i}\psi_{i}(\mathbf{r}),~{}i=1,\dots,N_{occ},$
(1)

where $N_{occ}=N_{ele}/2$ if $N_{ele}$ is even and $(N_{ele}+1)/2$ if
$N_{ele}$ is odd is the number of occupation number.
$V_{\mathrm{ext}}(\mathbf{r})=-\sum_{I}Z_{I}/|\mathbf{R}_{I}-\mathbf{r}|$
represents the external potential due to the nuclei located at
$\\{\mathbf{R}_{I}\\}_{I=1,\dots,N_{nuc}}$ with charge
$\\{Z_{I}\\}_{I=1,\dots,N_{nuc}}$. $V_{\mathrm{Har}}(\mathbf{r})$ stands for
the Hartree potential which describes the interaction among electrons.
$V_{\mathrm{xc}}(\mathbf{r})$ is the exchange-correlation potential used to
absorb the difference between the Kohn–Sham one-body non-interactive system
and the real many-body interactive system.

Various numerical methods have been presented to solve the Kohn-Sham density
functional theory [27, 1, 29, 2, 4, 18, 21, 5]. To verify the effectiveness of
these methods, the ability to achieve chemical accuracy has been a commonly
accepted criterion. The chemical accuracy is generally considered as 1
kcal/mole accuracy ($\approx 1.59\times 10^{-3}$ Hartree/particle) with
respect to the total energy [20]. In this work, we consider the chemical
accuracy as the total energy to achieve the accuracy within $1\times 10^{-3}$
Hartree/atom compared with the energy obtained from reference software.

Due to the rapid variations of the wavefunctions around the nuclear positions,
a uniform discretization of the Kohn–Sham equation in the computational domain
would result in an unaffordable implementation, especially in all-electron
calculations. The adaptive finite element (AFE) methods which enable different
mesh sizes are favored in solving the Kohn–Sham equation. Furthermore, the
ability of handling the non-periodic boundary condition and complex
computational domain makes AFE methods more and more attractive in electronic
structure calculation in recent decades [23, 2, 18, 3, 15, 5, 17].

There have been lots of efforts for efficiently achieving chemical accuracy by
the adaptive finite element methods, such as using the high-order finite
elements [13, 18, 21, 15, 5, 9], high-order numerical quadrature [21, 5], etc.
In addition, an accurate approximation to the Hartree potential is shown to be
key in achieving chemical accuracy [23, 18, 21, 5]. Under the finite element
discretization framework, the Hartree potential is generally obtained from
solving a Poisson equation on the same computational domain as the Kohn–Sham
equation. It is noted that the decay of Hartree potential in the infinity
behaves as $1/r$ which is much slower than the exponential decay of the
wavefunctions, hence the homogeneous boundary condition in the finite
computational domain for the Hartree potential may introduce large errors in
the ground state energy [5]. To efficiently reduce such errors, the multipole
expansion of the Hartree potential [12] can be used to generate an
inhomogeneous boundary condition for the Poisson equation. It is noted that to
keep using the homogeneous boundary condition, there also have been several
works in the literature. The neutralization techniques have been studied [13,
19] where an additional source term is introduced in the Poisson equation
leading an $r^{-2}$ or even faster decay of the potential such that the
homogeneous boundary condition makes sense. Another approach is adopting a
sufficiently large computational domain for the Poisson equation as in [18],
wherein the diameter for the computational domain of the Poisson equation is
around 100 times larger than that of the Kohn–Sham equation. In such a way,
the zero boundary condition is adequate for the Hartree potential.

As stated in [21, 5], one key issue in accurately approximating the Hartree
potential is to choose an appropriate finite element space for the Poisson
equation. Notably, the electron density, which represents the sum of squared
wavefunctions and constitutes the right hand side of the Poisson equation,
belongs to a finite element space twice the polynomial degree of that used for
the wavefunctions. This suggests a preference for a larger finite element
space in the quest to find the Hartree potential. Various approaches,
including using a finer mesh with the same polynomial degree as the Kohn–Sham
equation [23] or doubling the polynomial degree [21, 5], have shown
improvements in the total energy calculations. However, it should be pointed
out that these approaches are far from perfect, due to the reason that the
approximate space for the Hartree potential in each approach is inherited from
the one for wavefunctions. This space may not ideally suit the Hartree
potential’s distinct behaviors in comparison to wavefunctions. A potential
remedy for creating a suitable space accommodating both the Hartree potential
and the wavefunctions is the use of an a posteriori error indicator [25],
considering contributions from both the Poisson and Kohn-Sham equations. While
this approach can yield accurate results, it also introduces significant
computational complexity, resulting from the substantial number of degrees of
freedom defined on a single mesh, driven by differences not only in decay
behavior but also in the magnitude of the error indicator between the two
quantities. In light of these observations, a practical solution, that is, the
multi-mesh adaptive method, to simultaneously address accuracy and complexity
concerns emerges through a combination of the aforementioned approaches. This
involves the adaptive design of separate approximate spaces for wavefunctions
and the Hartree potential. With such an idea, the numerical accuracy issue can
be handled well since the approximate space for each variable is tailored
according to its own feature, while the complexity is effectively reduced by
replacing a large system of linear equations with two smaller ones.

Distinct from the method that solves all the equations on the same mesh
(hereafter we call such method single-mesh method), the multi-mesh adaptive
method (abbreviated as multi-mesh method) is to solve different variables with
different regularity on different meshes in the same computational domain. The
quality of mesh associated with a certain variable can be improved from the
mesh adaption process similar to the mesh adaption in the single-mesh method.
In this study, the implementation of the multi-mesh method is based on the
framework proposed in [14]. Within this framework, mesh adaptation is achieved
by locally refining or coarsening the element, a process known as
$h$-adaptation. A similar idea can be found in [8, 22] where $p$-adaption of
the finite element space is also considered. The multi-mesh method has been
successfully applied to various fields and models, such as in simulating the
two- and three-dimensional dendritic growth [11, 6], the wetting and spreading
problems [7], the photonic band structure optimization [28], the eigenvalue
problems [10], etc. With the multi-mesh method, chemical accuracy is able to
be achieved since both wavefunctions and the Hartree potential are solved on
different but qualified spaces.

An anticipated advantage of the multi-mesh method, when compared to the
single-mesh method, is the potential to reduce computational costs in
achieving chemical accuracy. This can be explained by a sample as displayed in
Figure 1. The first mesh can be viewed as the mesh for solving the Kohn–Sham
equation, with a dense grid distribution near the singularity and a sparse
distribution for regions far from it. The middle mesh is used to solve the
Hartree potential, and due to its smooth behavior and low decay speed, the
mesh grid distribution remains uniform. Finally, the last mesh, formed without
employing the multi-mesh strategy, should be able to capture both the
singularity and the low decay of the Hartree potential. The number of mesh
grids for these three meshes is 37, 49, and 57, respectively. Assume that the
linear finite element method is adopted. In the multi-mesh method, a
generalized eigenvalue problem $Ax=\lambda Bx$ with
$A,B\in\mathbb{R}^{37\times 37}$ and a linear system $Sx=b$ with
$S\in\mathbb{R}^{49\times 49}$. While for the finite element method using one
mesh, the eigenvalue problem and the linear system are of the same size
$A,B,S\in\mathbb{R}^{57\times 57}$. Apparently, the computational cost for
solving these two problems in the multi-mesh method is much less than the
single-mesh method in this sample.

$\oplus$$>$ Figure 1: A sample for multi-mesh method. The left two images show
the meshes used in the multi-mesh method. The right one shows the mesh used in
the single-mesh method. The number of the mesh grids for the three meshes are
37, 49, and 57, respectively.

Nevertheless, the application of the multi-mesh method is accompanied by
several challenges. The first of these challenges is the effective management
of mesh grids, requiring the capacity to perform local grid refinements or
coarsening in response to specific problem demands. This management process
must seamlessly facilitate the update of solutions from the previous mesh to
the new one. The second challenge centers on the critical necessity for
efficient and accurate multi-mesh communication, particularly when computing
integrals. The establishment of smooth communication is indispensable for
maintaining consistency and accuracy across different meshes. To address these
challenges, we introduce the hierarchical geometry tree data structure. This
tree-based approach facilitates the natural and efficient manner in the
management of the mesh grids. Furthermore, in order to prevent any loss of
accuracy during communications across various meshes, we adopt a strategy that
maximizes the utilization of quadrature points in numerical integrals, relying
solely on interpolation without the need for a projection process. This
approach ensures the preservation of numerical accuracy throughout the
communication process.

By employing the multi-mesh approach, we establish a multi-mesh adaptive
finite element framework for the all-electron density functional theory. With
the additional treat on the Hartree potential, systematical convergence of the
total energy to chemical accuracy is observed and the computational cost is
less compared with the single-mesh method given the same accuracy. The mesh
adaption is based on the residual-based a posteriori error estimation of the
Kohn–Sham equation and the Poisson equation whose solution is the Hartree
potential. The ability to achieve chemical accuracy is verified by a series of
numerical examples of atoms and molecules. By a detailed comparison of
computational time between the single-mesh method and the multi-mesh method,
the efficiency of the multi-mesh adaptive method is justified.

The organization of this paper is listed below. In Section 2, we introduce the
Kohn–Sham equation, the Poisson equation, and the finite element
discretizations of these two equations. Then in Section 3, the single-mesh
adaptive finite element framework for Kohn–Sham density function theory is
reviewed and tested on the Helium atom example. Followed by Section 4 in which
the multi-mesh adaptive method is demonstrated and examined also on the Helium
example. In Section 5, the efficiency and accuracy of the presented method are
discussed in detail from numerical examples. Finally, this paper ends with the
conclusion in Section 6.

##  2 The Kohn–Sham equation and finite element discretization

###  2.1 Kohn–Sham density functional theory

The Kohn–Sham equation for a $p$-electron system is read as the following
eigenvalue problem

$\left\\{\begin{array}[]{lr}\hat{H}\psi_{l}(\mathbf{r})=\varepsilon_{l}\psi_{l}(\mathbf{r}),&l=1,2,\dots,p,\\\
\displaystyle\int_{\mathbb{R}^{3}}\psi_{l}\psi_{l^{\prime}}\,\mathrm{d}\mathbf{r}=\delta_{ll^{\prime}},&l,l^{\prime}=1,2,\dots,p,\end{array}\right.$
(2)

where $(\varepsilon_{l},\psi_{l})$ is the $l$-th eigenpair,
$\delta_{ll^{\prime}}$ is the Kronecker delta, and $\hat{H}$ stands for the
Hamiltonian operator. We denote
$\rho(\mathbf{r})=\sum_{l=1}^{p}|\psi_{l}|^{2}(\mathbf{r})$ as the electron
density. Specifically, $\hat{H}$ consists of the following four terms

$\hat{H}([\rho];\mathbf{r})=-\frac{1}{2}\nabla_{\mathbf{r}}^{2}+V_{\mathrm{ext}}(\mathbf{r})+V_{\mathrm{Har}}([\rho];\mathbf{r})+V_{\mathrm{xc}}([\rho];\mathbf{r}),$
(3)

where the notation $V([\rho];\mathbf{r})$ implies that $V$ is a functional of
the electron density $\rho$. The first term $-\nabla^{2}/2$ in $\hat{H}$ is
the kinetic operator. The second term in $\hat{H}$ describes the Coulomb
external potential due to the nuclei which takes the form

$V_{\mathrm{ext}}(\mathbf{r})=-\sum_{j=1}^{M}\frac{Z_{j}}{\lvert\mathbf{r}-\mathbf{R}_{j}\rvert},$

where $M$ is the number of nuclei. The third term is the Hartree potential
describing the Coulomb repulsion among the electrons

$V_{\mathrm{Har}}([\rho];\mathbf{r})=\int_{\mathbb{R}^{3}}\frac{\rho(\mathbf{r^{\prime}})}{\lvert\mathbf{r}-\mathbf{r^{\prime}}\lvert}\,\mathrm{d}\mathbf{r^{\prime}}.$
(4)

The last term $V_{\mathrm{xc}}$ stands for the exchange-correlation potential,
which is caused by the Pauli exclusion principle and other non-classical
Coulomb interactions. Note that the analytical expression for the exchange-
correlation term is unknown and therefore an approximation is needed.
Specifically, the local density approximation (LDA) from the library Libxc
[16] with the slater exchange potential and the Vosko-Wilk-Nusair (VWN4) [26]
is adopted in this work.

Note that direct evaluation of the Hartree potential (4) requires
computational cost $\mathcal{O}(N^{2})$ with $N$ being the number of grid
points on the computational domain $\Omega$. For simplicity, we denote
$\phi=V_{\mathrm{Har}}(\mathbf{r})$ hereafter. In this paper, the Hartree
potential is obtained by solving the Poisson equation

$\left\\{\begin{aligned}
-\nabla^{2}\phi(\mathbf{r})=4\pi\rho(\mathbf{r}),&~{}\text{in }\Omega,\\\
\phi(\mathbf{r})=\phi_{\partial\Omega}(\mathbf{r}),&~{}\text{on
}\partial\Omega.\end{aligned}\right.$ (5)

The boundary value $\phi_{\partial\Omega}(\mathbf{r})$ is evaluated by the
multipole expansion method. Specifically, the following approximation is used

$\displaystyle\left.\phi(\mathbf{r})\right|_{\mathbf{r}\in\partial\Omega}\approx$
$\displaystyle\frac{1}{\left|\mathbf{r}-\mathbf{r}^{\prime\prime}\right|}\int_{\Omega}\rho\left(\mathbf{r}^{\prime}\right)\,\mathrm{d}\mathbf{r}^{\prime}+\sum_{i=1,2,3}p_{i}\cdot\frac{r^{i}-r^{\prime\prime,i}}{\left|\mathbf{r}-\mathbf{r}^{\prime\prime}\right|^{3}}+\sum_{i,j=1,2,3}q_{ij}\cdot\frac{3\left(r^{i}-r^{\prime\prime,i}\right)\left(r^{j}-r^{\prime\prime,j}\right)-\delta_{ij}\left|\mathbf{r}-\mathbf{r}^{\prime\prime}\right|^{2}}{\left|\mathbf{r}-\mathbf{r}^{\prime\prime}\right|^{5}},$

where

$p_{i}=\int_{\Omega}\rho\left(\mathbf{r}^{\prime}\right)\left(r^{\prime,i}-r^{\prime\prime,i}\right)\,\mathrm{d}\mathbf{r}^{\prime},\quad
q_{ij}=\int_{\Omega}\frac{1}{2}\rho\left(\mathbf{r}^{\prime}\right)\left(r^{\prime,i}-r^{\prime\prime,i}\right)\left(r^{\prime,j}-r^{\prime\prime,j}\right)\,\mathrm{d}\mathbf{r}^{\prime}.$

In the above expressions, $\mathbf{r}^{\prime\prime}$ stands for an arbitrary
point in $\Omega$. In the simulations, we choose it to be

$\mathbf{r}^{\prime\prime}=\frac{\int\mathbf{r}\rho(\mathbf{r})\,d\mathbf{r}}{\int\rho(\mathbf{r})\,d\mathbf{r}}.$

###  2.2 Finite element discretization

In practical simulations, a bounded polyhedral domain
$\Omega\subset\mathbb{R}^{3}$ is served as the computational domain. Assume
that the finite element space $V_{h}$ is constructed on $\Omega$ and the
finite element basis of $V_{h}$ is denoted as
$\\{\varphi_{1},\dots,\varphi_{n}\\}$ with $n$ being the dimension of the
space. Then the wavefunction $\psi_{i}$ can be approximated as $\psi_{i}^{h}$
on $V_{h}$ via

$\psi_{i}^{h}(\mathbf{r})=\sum_{k=1}^{n}X_{k,i}\varphi_{k},\quad
X_{k,i}=\psi_{i}\left(\mathbf{r}_{k}\right),\quad X\in\mathbb{R}^{n\times p},$
(6)

where $\mathbf{r}_{k}$ denotes the node corresponding to the $i$-th basis
function. As a result, to find the approximation of the wavefunctions
$\\{\psi_{i}\\}$ on $V_{h}$ is to find $X$, i.e., the value of $\psi_{i}$ at
each node $\mathbf{r}_{k}$.

In the finite element space $V_{h}$, the variation form for the KS equation
(2) becomes: find $(\psi_{i}^{h},\varepsilon_{i}^{h})_{i=1,\dots,p}\in
V_{h}\times\mathbb{R}$, such that

$\int_{\Omega}\varphi\hat{H}\psi_{i}^{h}\,\mathrm{d}\mathbf{r}=\varepsilon_{i}^{h}\int_{\Omega}\varphi\psi_{i}^{h}\,\mathrm{d}\mathbf{r},\quad\forall\varphi\in
V_{h}.$

By letting $\varphi$ be the finite element basis function and inserting (6) to
the above variational form, we have the following discrete eigenvalue problem

$AX=\varepsilon MY.$ (7)

Here $A$ and $M$ are symmetric matrices with the entries

$\displaystyle A_{i,j}$
$\displaystyle=\int_{\Omega}\frac{1}{2}\nabla\varphi_{j}\cdot\nabla\varphi_{i}+\big{(}V_{\mathrm{ext}}+\phi+V_{\mathrm{xc}}\big{)}\varphi_{j}\varphi_{i}\,\mathrm{d}\mathbf{r},$
(8) $\displaystyle M_{i,j}$
$\displaystyle=\int_{\Omega}\varphi_{j}\varphi_{i}\,\mathrm{d}\mathbf{r}.$ (9)

Similarly, we can obtain the linear system for the Poisson equation (5) on
$V_{h}$:

$S\Phi=\mathbf{b},$ (10)

where $S$ is the stiff matrix with the entry

$S_{i,j}=\int_{\Omega}\frac{1}{2}\nabla\varphi_{j}\cdot\nabla\varphi_{i}\,\mathrm{d}\mathbf{r},$
(11)

and the right hand side $\mathbf{b}$ and the discretized $\phi^{h}$

$\mathbf{b}_{i}=\int_{\Omega}4\pi\rho(\mathbf{r})\varphi_{i}\,\mathrm{d}\mathbf{r},\quad\phi^{h}(\mathbf{r})=\sum_{k=1}^{n}\Phi_{k}\varphi_{k},\quad\Phi\in\mathbb{R}^{n}.$

Owing to the singularity in the Hamiltonian (3), a uniform finite element
discretization would lead to a large number of mesh grids to achieve chemical
accuracy. Hence the adaptive mesh method is necessary to efficiently solve the
Kohn–Sham equation, which will be discussed in the next section.

##  3 The adaptive finite element method

In this section, we will begin by introducing the adaptive finite element
method that relies on the _a posteriori_ error estimates [25]. Subsequently,
we will present and compare two error indicators used to facilitate mesh
adaptation for solving the Kohn–Sham equation (2). The sole distinction
between these two error estimators is their utilization of the error involved
in solving the Hartree potential.

###  3.1 The adaptive method based on the residual _a posteriori_ error
estimation

The adaptive mesh techniques offer enhanced numerical accuracy while requiring
fewer mesh grids compared to uniform mesh methods. In this study, our primary
focus lies on the $h$-adaptive methods, which involve local refinement and/or
coarsening of mesh grids. An important aspect of $h$-adaptive methods is the
determination of an error indicator. Generally speaking, error indicators
identify regions within the domain that necessitate local refinement or
coarsening, and they are typically derived from _a posteriori_ error
estimations [25]. When there is only one orbital in the system, it is natural
to generate the indicator based on information from that specific orbital.
However, when there are multiple orbitals in the system, generating the
indicator solely from an individual orbital is no longer advisable. This is
because every orbital in the system is expected to be well-resolved using the
mesh grids after mesh adaptation. To address this, we adopt the strategy
proposed in [2] for indicator generation. First, indicators are individually
generated for each orbital using a specific method. Then, normalization is
applied to each indicator. The final indicator is obtained by combining these
normalized indicators.

Specifically, based on the _a posteriori_ error estimates [25], the residual-
based _a posterirori_ error indicator for the KS equation (2) in the element
$\mathcal{T}_{K}$ can be defined as

$\eta_{K,\mathrm{KS}}=\left(h_{K}^{2}\sum_{l=1}^{p}\big{\|}\mathbb{R}_{K,\mathrm{KS}}(\psi_{l})\big{\|}_{K}^{2}+\sum_{e\subset\partial\mathcal{T}_{K}}\frac{1}{2}h_{e}\sum_{l=1}^{p}\big{\|}\mathbb{J}_{e}(\psi_{l})\big{\|}_{e}^{2}\right)^{\frac{1}{2}},$
(12)

where $h_{K}$ represents the largest length of the edges of the element
$\mathcal{T}_{K}$, $h_{e}$ stands for the largest length of the common face
$e$ of $\mathcal{T}_{K}$ and $\mathcal{T}_{J}$,
$\mathbb{R}_{K,\mathrm{KS}}(\psi)$ and $\mathbb{J}_{e}(\psi)$ are the residual
and jump term on the element $\mathcal{T}_{K}$, whose formulations are written
as

$\left\\{\begin{aligned}
\mathbb{R}_{K,\mathrm{KS}}(\psi_{l})&=\hat{H}\psi_{l}-\varepsilon_{l}M\psi_{l},\\\
\mathbb{J}_{e}(\psi_{l})&=(\nabla\psi_{l}\bigm{|}_{\mathcal{T}_{K}}-\nabla\psi_{l}\bigm{|}_{\mathcal{T}_{J}})\cdot\mathbf{n}_{e},\end{aligned}\right.$

where $\mathbf{n}_{e}$ is the out normal vector on the face $e$ w.r.t the
element $\mathcal{T}_{K}$. The definition of the norms are

$\left\|f(x)\right\|_{K}=\left(\int_{K}(f(x))^{2}dx\right)^{\frac{1}{2}},~{}~{}\left\|f(x)\right\|_{e}=\left(\int_{e}(f(x))^{2}dx\right)^{\frac{1}{2}}.$

The indicator (12) involves the error arising from the Kohn–Sham equation (2)
and is usually directly adopted to guide the mesh adaption.

An aspect that the indicator (12) may overlook is the error associated with
the Hartree potential. It is important to emphasize the accurate approximation
of the Hartree potential, as it is a crucial component of the Hamiltonian (3)
and contributes to the Hartree potential energy. As the Hartree potential is
obtained by solving the Poisson equation (5), a similar approach can be
employed to generate an error indicator specifically for the Hartree
potential, just as done for the Kohn–Sham equation. Specifically, the error
indicator for the Hartree potential can be defined as

$\eta_{K,\mathrm{Har}}=\left(h_{K}^{2}\big{\|}\mathbb{R}_{K,\mathrm{Har}}(\phi)\big{\|}_{K}^{2}+\sum_{e\subset\partial\mathcal{T}_{K}}\frac{1}{2}h_{e}\big{\|}\mathbb{J}_{e}(\phi)\big{\|}_{e}^{2}\right)^{\frac{1}{2}},$
(13)

where the residual part is defined as
$\mathbb{R}_{K,\mathrm{Har}}(\phi)=\nabla^{2}\phi+4\pi\rho$. Based on the
analysis above, the second indicator for solving the Kohn–Sham equation can be
designed as

$\eta_{K,\mathrm{KS+Har}}=\sqrt{\eta_{K,\mathrm{Har}}^{2}+\eta_{K,\mathrm{KS}}^{2}}.$
(14)

The aim of using this error indicator is to generate a mesh on which both the
wavefunctions and the Hartree potential are approximated well.

The mesh adaptation process, along with the solution of the Kohn–Sham
equation, can be carried out using either the first (12) or the second error
indicator (14). It is important to note that the adaptive algorithms utilizing
these two different indicators are essentially identical, except for the
choice of the indicator. Moreover, all simulations are performed on a single
mesh simultaneously, i.e., the Poisson equation (5) for the Hartree potential
is discretized and solved on the same finite element space for the Kohn–Sham
equation. Consequently, we present the adaptive finite element method for
solving the Kohn–Sham equation, outlined in Algorithm 1, and refer to this
approach as the _single-mesh adaptive method_. Briefly, the adaptive algorithm
consists of an outer iteration and an inner iteration. In the inner iteration,
it is an SCF method for solving the nonlinear Kohn–Sham equation. Meanwhile,
it is the mesh adaption procedure in the outer iteration. In the next
subsection, we would like to show an example using the single mesh adaptive
method with the first indicator (12).

Input: Initial mesh $\mathcal{T}^{(0)}$, initial electron density $\rho$,
initial energy $E=0$, energy tolerance $tol_{1}$, density tolerance $tol_{2}$.

while  _$|E-E_{old}| <tol_{1}$_ do

$E_{old}=E$;

while  _$\|\rho-\rho_{old}\| <tol_{2}$_ do

$\rho_{old}=\rho$;

Calculate the Hartree potential $\phi$;

Generate the Hamiltonian matrix $A$;

Solve the the eigenvalue problem $AX=\varepsilon MX$;

Update the electron density $\rho$;

Calculate the ground state energy $E$;

Mesh adaption based on the error indicator.

Output: Total ground state energy $E$

Algorithm 1 Single-mesh adaptive method.

###  3.2 An issue for numerical solutions towards chemical accuracy

####  3.2.1 Single-mesh adaptive method using indicator (12)

To assess the convergence and behavior of the single mesh adaptive finite
element method with the error indicator (12), we conduct a series of eight
experiments for the helium atom, each involving a different number of mesh
grids, controlled by varying adaption tolerance values. In this example, the
referenced value is obtained from the state-of-art software NWChem using the
_aug-cc-pv6z_ basis set [24]. The hardware and software configurations of the
experiment are introduced in the beginning of Section 5.

The results are displayed in the Figure 2. As the tolerance for error
indicator decreases, the number of degrees of freedom increases. The
convergence of the total energy would serve as evidence of the superiority of
the adaptive method compared to the method employing a uniform mesh partition.
However, it is evident that the total energy does not converge to the
reference result observed from the left of Figure 2. Even in the mesh
generated from the smallest tolerance, which consists of 1,737,569 mesh grids,
a discrepancy of 0.01 Hartree in the energy is still observed, indicating that
the adaptive method does not achieve the desired convergence.

Figure 2: Convergence history of the total enegy (left), absolute errors of
energies and eigenvalue (right) for the adaptive method using error indicator
(12).

As discussed earlier, the error in the total energy primarily stems from the
inaccurate approximation of the Hartree energy. This is confirmed by comparing
the Hartree energy with the reference value $E_{\mathrm{Har,ref}}=1.9961$
Hartree, as depicted on the right side of Figure 2. It is observed that
although the error in the Hartree energy diminishes with respect to the number
of Dofs, on the finest mesh, the Hartree energy is found to be 0.01 Hartree
lower than the reference value, resulting in a 0.01 Hartree error in the total
energy. Furthermore, an imprecise approximation of the Hartree potential can
also introduce a 0.01 Hartree error in the eigenvalue, as depicted on the
right side of Figure 2. As a comparison, it can be observed that the errors in
the one-electron energy and the exchange-correlation energy systematically
decrease as the number of Dofs increases.

The reason behind the inaccurate approximation of the Hartree potential can be
elucidated by referring to Figure 3 (upper row), which displays the results
obtained from the finest mesh consisting of 1,765,990 mesh grids. In the left
two columns, the global and local mesh distributions are depicted. It is
evident that the mesh grid density is considerably high around the singularity
located at the origin, with the smallest mesh size being approximately 0.001
Bohr. Conversely, in regions distant from the singularity, the mesh grid
distribution is sparse. Specifically, near the boundary, the largest mesh size
can exceed 2 Bohr. Such a mesh distribution proves to be suitable for
representing wavefunctions and electron density, as illustrated in the third
column of Figure 3 (upper row), since both quantities exhibit exponential
decay, and the region where the density exceeds 0.1 is confined to a small
interval of $[-1,1]^{3}$ for the helium atom. However, in the fourth column of
Figure 3 (upper row), it is apparent that the Hartree potential exhibits a
much slower decay behavior, with potential values surpassing 0.1 throughout
the entire computational domain of $[-10,10]^{3}$. Furthermore, the contours
near the boundary oscillate and are not smooth. As a result, a mesh size as
large as 2 Bohr is evidently insufficient for capturing the variations in the
Hartree potential accurately. In order to enhance the numerical accuracy of
the Kohn–Sham equation, it is essential to obtain a more precise approximation
of the Hartree potential.

Figure 3: Top row: results from a single mesh adaptive method using the
indicator (12), with total 1,765,990 mesh grids, i.e., Global mesh and zoomed-
in mesh on the sliced $X$-$Y$ plane (left two), density profile (third one),
and Hartree potential profile (fourth one). Bottom row: corresponding results
using the indicator (14), with total 2,847,807 mesh grids.

####  3.2.2 Single mesh adaptive method using indicator (14)

To achieve a more accurate Hartree potential, a direct approach is to
incorporate the corresponding error indicator (13) with the indicator of the
Kohn–Sham equation (12). This combination results in the second indicator
(14). Consequently, we apply the single mesh adaptive method again to the
Helium example, but this time utilizing the second indicator (14).

The convergence of energies and eigenvalues is depicted in Figure 4. Notably,
there is a systematic convergence observed in both the total energy and the
Hartree potential energy. Moreover, the eigenvalue demonstrates convergence
towards the reference value. On the finest mesh, which consists of 2,874,807
grid points, chemical accuracy is achieved. Sliced mesh representations and
profiles of the electron density and Hartree potential are shown in Figure 3
(lower row). A noticeable difference from the mesh in Figure 3 (upper row) is
the presence of smaller mesh sizes in regions far away from the origin. As a
result, the Hartree potential is better captured, as illustrated in the fourth
column of Figure 3 (lower row), from which we can find smoother contour lines
than that in Figure 3 (upper row) using the first indicator (12).

Figure 4: Convergence history of the total energy (left), absolute errors of
energies and eigenvalue (right) for the adaptive method using error indicator
(14).

While the Hartree potential is approximated accurately and chemical accuracy
is obtained, the computational cost increases significantly. The number of
mesh grids on the finest mesh reaches 2,874,807, surpassing the number of
grids used in the previous section. This increment in the number of mesh grids
is primarily caused by the need to capture the variations of the Hartree
potential. Such an increase in the computational grid introduces additional
complexity and computational overhead, leading to longer computation time and
higher memory requirements. This becomes especially problematic when dealing
with larger system sizes or when conducting complex simulations. To address
this issue and reduce the computational cost, an alternative approach, known
as the multi-mesh adaptive method, will be introduced in the next section.

##  4 The multi-mesh adaptive method

In the preceding section, we explored the solution of the discretized
Kohn–Sham equation on a single mesh, employing either the residual-based _a
posteriori_ error indicator (12) or (14). The results indicated that employing
the first indicator (12) alone may fall short of achieving chemical accuracy.
On the other hand, the second indicator (14) demonstrated the ability to
attain chemical accuracy; however, it came at the cost of significantly
increased computational requirements.

To tackle this challenge, we introduce an approach known as the multi-mesh
adaptive method [14]. The primary objective of this method is to strike a
balance between achieving chemical accuracy and managing computational costs
effectively. By utilizing multiple meshes instead of a single mesh, the multi-
mesh adaptive method offers greater flexibility in adapting the mesh
resolution to capture variations for different quantities of interest.

Specifically, in this method, we will utilize two distinct meshes. Both two
meshes will be adapted during the simulation. The first mesh is specifically
designed for solving the Kohn–Sham equation, with the primary objective of
accurately capturing the variations in the wavefunctions. This mesh is
tailored to ensure high resolution in regions where the wavefunctions exhibit
significant changes. Conversely, the second mesh is dedicated to solving the
Poisson equation (5). Its primary purpose is to capture the variations in the
Hartree potential effectively. This second mesh is strategically designed to
provide optimal resolution and precision in regions where the Hartree
potential exhibits substantial changes. By employing two separate meshes with
specific focus areas, we can ensure that each equation is solved with the
appropriate level of accuracy and capture the variations unique to the
wavefunctions and the Hartree potential, respectively.

To implement the multi-mesh method, careful handling of two components is
crucial. The first component involves effectively managing the mesh grids,
allowing for flexible local refinement or coarsening as needed. This
management mechanism should also facilitate efficient solution updates from
the old mesh to the new mesh. The second component focuses on ensuring
efficient and accurate communication between different meshes, particularly in
the calculation of integrals. Efficient communication protocols play a vital
role in maintaining consistency and accuracy across the various meshes. These
requirements can be fulfilled by the hierarchical geometry tree data
structure, which will be introduced in detail in the following.

###  4.1 The hierarchical geometry tree

A well-designed data structure for the mesh grids is needed for an effective
management mechanism. In the presented algorithm, the hierarchical geometry
tree (HGT) [14, 2] is utilized.

Firstly, the mesh structure is described hierarchically, which means that the
mesh information is given from the lowest dimension (0-D, the points) to the
highest dimension (3-D, the tetrahedron) hierarchically. An element such as a
point for 0-D, an edge for 1-D, a triangle for 2-D, a tetrahedron for 3-D is
called a geometry. In the hierarchical description of a tetrahedron, all
geometries have a belonging-to relationship. For example, if an edge is one of
the edges of a triangle, this edge belongs to this triangle. With this
hierarchical structure, the geometry information of a tetrahedron can be
referred to flexibly, and the refinement and coarsening of a mesh can also be
implemented efficiently.

Secondly, the mesh is stored and managed by a tree data structure. The
validity of using the tree data structure is due to the strategy of element
refinement and coarsening. Specifically, for a tetrahedron element (the left
of Figure 5), the refinement of this tetrahedron is dividing it into eight
equally small tetrahedrons via connecting the midpoints on each edge (the
right of Figure 5). As a result, a belonging-to relationship can be
established, for example, any small tetrahedron that is called “child” belongs
to the large “parent” tetrahedron. Meanwhile, the coarsening of any child
tetrahedron in the right of Figure 5 is releasing all the children to obtain
the parent tetrahedron. By this refinement and coarsening strategy, the tree
data structure is able to be established.

Figure 5: The parent tetrahedron (left) and the eight child tetrahedrons
(right).

For a better understanding of the tree data structure, a two-dimensional
illustration is presented. Similar to the 3D case, the refinement of a 2D
element, i.e. the triangle, is to divide the triangle into four equal
triangles, as displayed in the left two columns in Figure 6. By refining a
triangle $\mathcal{T}_{0}$, four sub-triangles
$\\{\mathcal{T}_{0,0},\mathcal{T}_{0,1},\mathcal{T}_{0,2},\mathcal{T}_{0,3}\\}$
are generated. Furthermore, via refining $\mathcal{T}_{0,0}$ the sub-triangles
$\\{\mathcal{T}_{0,0,0},\mathcal{T}_{0,0,1},\mathcal{T}_{0,0,2},\mathcal{T}_{0,0,3}\\}$
are obtained, which is demonstrated in Figure 6. To manage this procedure. the
quadtree data structure in which each internal node has exactly four children
is utilized, as described in Figure 7. When the local refinement and
coarsening techniques are adopted, only some triangles are to refined, as
indicated in the right of Figure 6 and the bottom of Figure 7 where only the
triangle $\mathcal{T}_{0,0}$ is refined. In the quadtree, we call
$\mathcal{T}_{0}$ the root node, and those nodes without further descents like
$\mathcal{T}_{0,1}$ and $\mathcal{T}_{0,0,1}$ the leaf nodes. Suppose that
there is a set of root nodes $\\{\mathcal{T}_{i}\\},i=0,1,2,\dots$, which form
a mesh for a domain $\Omega$, then from the above definition we know that a
set of all leaf nodes also form a mesh. For example, the set
$\\{\mathcal{T}_{0}\\}$ forms the mesh in the left of Figure 6, the set
$\\{\mathcal{T}_{0,0},\mathcal{T}_{0,1},\mathcal{T}_{0,2},\mathcal{T}_{0,3}\\}$
forms the middle mesh in Figure 6, and the set
$\\{\mathcal{T}_{0,0,0},\mathcal{T}_{0,0,1},\mathcal{T}_{0,0,2},\mathcal{T}_{0,0,3},\mathcal{T}_{0,1},\mathcal{T}_{0,2},\mathcal{T}_{0,3}\\}$
forms the right mesh in Figure 6.

ABC$\mathcal{T}_{0}$ABCDEF$\mathcal{T}_{0,0}$$\mathcal{T}_{0,1}$$\mathcal{T}_{0,2}$$\mathcal{T}_{0,3}$ABCDEFGHI$\mathcal{T}_{0,0,0}$$\mathcal{T}_{0,0,1}$$\mathcal{T}_{0,0,2}$$\mathcal{T}_{0,0,3}$$\mathcal{T}_{0,1}$$\mathcal{T}_{0,2}$$\mathcal{T}_{0,3}$
Figure 6: The refinement of a triangle $\mathcal{T}_{0}$ and its sub-triangle
$\mathcal{T}_{0,0}$.

[$\mathcal{T}_{0}$ [$\mathcal{T}_{0,0}$ [$\mathcal{T}_{0,0,0}$]
[$\mathcal{T}_{0,0,1}$] [$\mathcal{T}_{0,0,2}$] [$\mathcal{T}_{0,0,3}$] ]
[$\mathcal{T}_{0,1}$] [$\mathcal{T}_{0,2}$] [$\mathcal{T}_{0,3}$] ]

Figure 7: The quadtree data structure for the mesh in the right of Figure 6.

With the hierarchical description of the geometry and the tree data structure,
the mesh is effectively managed by the HGT. However, building the finite
element space directly on a mesh results in non-conforming finite element
because of the hanging points in the direct neighbors of those refined
triangles. These hanging points can be handled mainly in two ways. If there
are more than one hanging point on a triangle, this triangle will be refined.
This process is called the semi-regularization procedure [14]. If there exists
only one hanging point on a triangle, for example, the triangle $\triangle
DEF$ in the right of Figure 6 has a hanging point $I$ arose in the refinement
of $\mathcal{T}_{0,0}$. To deal with these kind of hanging points, the twin-
triangle geometry is introduced, which is demonstrated in Figure 8. The twin-
triangle actually consists of two triangles $\triangle FDI$ and $\triangle
FIE$. Meanwhile, the degrees of freedom and the interpolation information of
the twin-triangle $\triangle DEF$ inherit from the two triangles $\triangle
FDI$ and $\triangle FIE$. To make the finite element space conforming, the
following strategy is adopted to construct the basis functions in the twin-
triangle geometry: For each basis function, its value is 1 at its
corresponding interpolation point and 0 at other interpolation points. The
support of the basis function whose interpolation point is a common point of
two sub-triangles in a twin-triangle geometry such as $F$ and $I$ in Figure 8,
is the whole twin-triangle geometry. For the non-common points like $D$, the
support of its basis function is only the triangle $\triangle FDI$. With such
a strategy, a conforming finite element space can be built in a mesh where the
local refinement is implemented.

$D$$E$$I$$F$ Figure 8: Twin-triangle.

In virtue of the tree data structure, multiple meshes are allowed to be
described by the same HGT. Naturally, a problem arises: is it possible to make
the information communicates among these meshes without any loss? The answer
is affirmative, and the reason relies on the belonging-to relationship between
any two nodes in the HGT. Briefly speaking, there exist only three kinds of
relationship between any two nodes: equal, belonging-to, and no-overlap. The
communication among equal or no-overlap elements is trivial. Consequently, we
only need to take care of the second kind of relationship. We introduce the
implementation details in the following.

In the presented method, the Hamiltonian matrix and energy require information
from both the KS mesh and the Hartree mesh, which mainly involve numerical
integration. We take the evaluation of the Hartree potential energy for
example, which could be written as

$E_{\mathrm{Har}}=\frac{1}{2}\int_{\Omega}\phi(\mathbf{r})\left(\sum_{l=1}^{p}\psi_{l}(\mathbf{r})^{2}\right)\mathrm{d}\mathbf{r}.$
(15)

The integration should be carefully treated since the Hartree potential
$\phi(\mathbf{r})$ and the wavefunction $\psi_{l}(\mathbf{r})$ belong to
different finite element spaces built on different meshes. An intuitive
illustration is given to show how the presented method calculates the integral
(15). Assume the finite element space $V_{\mathcal{T}_{\mathrm{KS}}}$ for the
wavefunction is built on the mesh $\mathcal{T}_{\mathrm{KS}}$, and the space
$V_{\mathcal{T}_{\mathrm{Har}}}$ is built on the mesh
$V_{\mathcal{T}_{\mathrm{Har}}}$, see Figure 9. For the common elements such
as $\triangle CDF$ which belongs to both $\mathcal{T}_{\mathrm{KS}}$ and
$\mathcal{T}_{\mathrm{Har}}$, the numerical integration can be calculated
directly. While for the remained elements, for example, in
$\mathcal{T}_{\mathrm{KS}}$, the triangle $\triangle DAE$ is refined, while in
$\mathcal{T}_{\mathrm{Har}}$ this triangle is kept, and a similar case for the
triangle $\triangle FEB$, special treatment is needed to avoid the loss of
accuracy.

$A$$B$$C$$D$$E$$F$$I$$H$$G$$A$$B$$C$$D$$E$$F$$J$$K$$L$ Figure 9: Left:
$\mathcal{T}_{\mathrm{KS}}$. Right: $\mathcal{T}_{\mathrm{Har}}$.

In order to prevent the loss of accuracy, we employ a strategy that maximizes
the utilization of quadrature points in numerical integrals. Specifically, the
numerical integration on the element $\triangle ADE$ is divided into the
integrations on its four refined sub-triangles,

$\displaystyle\frac{1}{2}\int_{\triangle
ADE}\phi(\mathbf{r})\left(\sum_{l=1}^{p}\psi_{l}(\mathbf{f})^{2}\right)\mathrm{d}\mathbf{r}\approx\sum_{K=1}^{4}\sum_{j=1}^{q}\mbox{area}(K)J_{j}^{K}w_{j}^{K}\phi(\mathbf{r}_{j}^{K})\left(\sum_{l=1}^{p}\psi_{l}(\mathbf{r}_{j}^{K})^{2}\right),$

where the element $K$ represent the four sub-triangles of $\triangle ADE$,
$\mathbf{r}_{j}^{K}$ is the $j$-th quadrature point of $K$, $J_{j}^{K}$ is the
jacobian at $\mathbf{r}_{j}^{K}$, and $w_{l}^{K}$ is the associated weight in
the numerical quadrature. The values of the Hartree potential on these
quadrature points can be obtained by numerical interpolation. Similarly, the
integral on the triangle $\triangle FEB$ is evaluated by doing the numerical
integral on its sub-triangles. In this way, the accuracy of the integral will
not be affected in the communication of different meshes.

With the HGT, the solution update from the old mesh to the new mesh after the
local refinement and coarsening can also be implemented efficiently. As we
mentioned before, a set of all leaf nodes forms a mesh. Although different
meshes correspond to different sets of leaf nodes, all meshes are from the
same set of root nodes. Then the belonging-to relationship of the geometries
between two meshes can be analyzed easily. As a result, the solution update
can be implemented efficiently according to the relationship.

###  4.2 The multi-mesh adaptive algorithm

In virtue of the HGT data structure, the multi-mesh method can be efficiently
implemented. To solve the Kohn–Sham equation, we propose the multi-mesh
adaptive algorithm illustrated in the flowchart presented in Figure 10. The
Kohn–Sham equation (2) is discretized and solved on the finite element space
built on the mesh $\mathcal{T}_{\mathrm{KS}}$, and the Hartree potential (5)
is solved on the mesh $\mathcal{T}_{\mathrm{Har}}$. The KS mesh
$\mathcal{T}_{\mathrm{KS}}$ is adapted using the error indicator (12) and the
Har mesh $\mathcal{T}_{\mathrm{Har}}$ is adapted using the error indicator
(13). Both of these two mesh adaptions occur after the end of the SCF
iteration. The mesh adaptations continue until the final total energy reaches
convergence, at which point the mesh adaptation ceases, and the results are
outputted.

Initial mesh
$\mathcal{T}_{\mathrm{KS}}=\mathcal{T}_{0},\mathcal{T}_{\mathrm{Har}}=\mathcal{T}_{0}$
Initial guess $X_{0}$ Calculation on $\mathcal{T}_{\mathrm{KS}}$Calculation on
$\mathcal{T}_{\mathrm{Har}}$Calculation on $\mathcal{T}_{\mathrm{KS}}$ and
$\mathcal{T}_{\mathrm{Har}}$ SCFmesh adaption Calculate $\phi$ Construct $A$
Solve $AX=\varepsilon BX$ Update $\rho$ $\|\rho-\rho_{old}\|<tol_{1}$
$|E-E_{old}|<tol_{2}$ Adapt mesh $\mathcal{T}_{\mathrm{KS}}$ with
$\eta_{K,\mathrm{KS}}$ Adapt mesh $\mathcal{T}_{\mathrm{Har}}$ with
$\eta_{K,\mathrm{Har}}$ Output NoYesNoNoYes Figure 10: Flowchart of the multi-
mesh adaptive algorithm for the KS equation.

In order to assess the efficiency and make a comparison between the multi-mesh
algorithm Figure 10 and the single mesh algorithm Algorithm 1, we categorize
the calculations into three distinct types. Namely, the calculations on the
Kohn–Sham mesh $\mathcal{T}_{\mathrm{KS}}$ with a white background in Figure
10, the calculations conducted in the Hartree mesh
$\mathcal{T}_{\mathrm{Har}}$ with a blue background, and the calculations that
require information from both $\mathcal{T}_{\mathrm{KS}}$
and$\mathcal{T}_{\mathrm{Har}}$ with a red background. When comparing the
multi-mesh algorithm to the single mesh algorithm, it becomes evident that
additional calculations fall into the third category. These calculations
necessitate communication between the Kohn–Sham mesh and the Hartree mesh.
Specifically, such communication is required during the generation of the
Hamiltonian, the updating of the electron density (which forms the right-hand
side of the Poisson equation (5) for the Hartree potential), and the
computation of the total energy. We will conduct a detailed comparison of the
time consumed by these individual parts to evaluate the overall performance of
the algorithms in the next section.

By employing the multi-mesh method and incorporating adaptive algorithms
guided by appropriate error indicators, we can achieve accurate and efficient
solutions for the Kohn–Sham equation. The flexibility of mesh adaptation
ensures that the mesh resolution aligns with the specific requirements of each
equation, ultimately leading to improved convergence and reliable results.
This can be verified by the Helium example.

##  5 Numerical Experiments

In this section, we examine the convergence and efficiency of the multi-mesh
adaptive method through a series of numerical examples. All the simulations
are performed on a workstation “Moss” with two AMD EPYC 7713 64-Core
Processors (at 2.0GHz$\times$64, 512M cache) and 900GB of RAM, and the total
number of cores is 128. The software is the C++ library AFEABIC [2] under
Ubuntu 20.04.

###  5.1 Case study: for chemical accuracy and computational efficiency

The multi-mesh adaptive method Figure 10 is first tested for three cases: the
helium atom, the LiH molecule, and the BeH2 molecule. Our examination unfolds
in three phases: firstly, we assess the method’s systematic convergence in the
context of the helium atom example; next, we determine whether it can attain
chemical accuracy across all three cases; and lastly, we undertake a
comparative analysis of the computational costs required to achieve chemical
accuracy between the multi-mesh method and the single-mesh method.

####  5.1.1 Helium atom

Similar to the adaptive method using error indicator (14), the systematical
convergence of the energies and eigenvalue are observed in Figure 11, and the
chemical accuracy is obtained when the number of degrees of freedom in
$\mathcal{T}_{\mathrm{KS}}$ achieves 1,262,599, which is far less than that
requires in the single-mesh algorithm. Furthermore, there are 2,185,920
degrees of freedom in the Hartree mesh $\mathcal{T}_{\mathrm{KS}}$, which is
also less than the number of mesh grids used to achieve the chemical accuracy
using the single mesh adaptive method.

Figure 11: Convergence history of the total energy (left), absolute errors of
energies, and eigenvalue (right) for the multi-mesh adaptive method.

Figure 12: Global and zoomed-in meshes on the sliced $X$-$Y$ plane (top:
$\mathcal{T}_{\mathrm{KS}}$ (1,2262,599 mesh grids); bottom:
$\mathcal{T}_{\mathrm{Har}}$ (2,185,920 mesh grids)), density profile (top
right), and Hartree potential profile (bottom right) for the multi-mesh
adaptive method.

The final meshes are displayed in Figure 12. In the top of Figure 12 is the
sliced mesh of $\mathcal{T}_{\mathrm{KS}}$ which is quite similar with Figure
3 (upper row). The meshes of $\mathcal{T}_{\mathrm{Har}}$ are demonstrated in
the bottom of Figure 12, which has a similar mesh grid distribution as that in
Figure 3 (lower row) as shown in a global vision, while it has a sparser grid
distribution comparing with Figure 3 (lower row) since the Hartree potential
behaves less singular than the Kohn–Sham wavefunctions. Although there is a
larger number of mesh grids in $\mathcal{T}_{\mathrm{Har}}$ than in
$\mathcal{T}_{\mathrm{KS}}$, it should be mentioned that on
$\mathcal{T}_{\mathrm{Har}}$ we only need to solve a linear system, while in
$\mathcal{T}_{\mathrm{KS}}$, a nonlinear eigenvalue problem should be solved.
Consequently, the major computational cost is from the calculations on
$\mathcal{T}_{\mathrm{KS}}$, and the increment of the computational cost due
to the large number of mesh grids on $\mathcal{T}_{\mathrm{Har}}$ is slight.
As a result, the computational cost of the multi-mesh method is of the same
order of magnitude of the adaptive method, while it can achieve chemical
accuracy.

As discussed in previous sections, both the single-mesh and multi-mesh methods
are able to achieve chemical accuracy when the Hartree potential is considered
in constructing the error indicator. To demonstrate the efficiency of the
multi-mesh method, we compare the serial time for these two methods to achieve
chemical accuracy. For the sake of fairness, we start the two algorithms from
the same mesh and initial guess, and the tolerance for guiding the mesh
adaption is also the same. With these settings, it is found that the
convergence in SCF is quite similar and the number of mesh adaptions is the
same, as displayed in Figure 13. From the figure, it suggests that the
introduction of Hartree mesh will not affect the convergence of the SCF
procedure in KS mesh. Therefore, the comparison of the CPU time is fair and
effective.

Figure 13: SCF convergence history (left) and mesh adaption history (right)
for Helium.

A summary of the comparison is illustrated in Table 1. Both methods achieved
chemical accuracy, and the number of mesh grids in the single-mesh method is
the largest among all meshes. The CPU time is compared with respect to two
parts: the SCF part and the mesh adaption part. From Table 1, the multi-mesh
method is proved to be more efficient than the single-mesh method in achieving
the chemical accuracy of the helium atom.

Table 1: Comparison on He with respect to results and serial computational time. $t_{all}$ represents the total CPU time. $t_{SCF}$ stands for the time for SCF and $t_{MA}$ represents the time for mesh adaption.  | $E_{tot}$ | $\Delta E$ | $N_{KS}$ | $N_{Har}$ | $t_{all}$ | $t_{SCF}$ | $t_{MA}$
---|---|---|---|---|---|---|---
single-mesh | -2.8341 | 0.0007 | 2,874,807 | 2,874,807 | 10290.45 | 7585.62 | 2704.83
multi-mesh | -2.8340 | 0.0008 | 1,262,599 | 2,185,920 | 7063.97 | 5135.69 | 1928.28

A more detailed comparison is presented in Figure 14. The method can be
divided into six parts as indicated in Figure 14. The four images represent
CPU time on four meshes during the mesh adaption process. Obviously, the cost
for solving the Hartree potential and the eigenvalue value problem is less in
the multi-mesh method, as expected. Furthermore, the cost of constructing a
matrix in the multi-mesh method is also less than that in the single-mesh
method on the finest mesh. It is noted that this part contributes the largest
portion of the total cost since in the helium example the eigensolver only
needs to solve one eigenpair. In larger systems, solving the eigenvalue
problem will be the most time-consuming part.

Figure 14: CPU time for single-mesh method and multi-mesh method on six parts:
MA (mesh adaption), SH (solve Hartree potential), CM (construct matrix), SE
(solve eigenvalue problem), UD (update electron density) and CE (calculate
energy). Four mesh adaptations are needed to achieve chemical accuracy in both
methods starting from the same initial setup.

####  5.1.2 LiH molecule

A similar comparison for the LiH molecule is employed. The referenced total
energy is $E_{\mathrm{ref}}=-7.9195$ Hartree. The results are displayed in
Table 2, revealing the attainment of chemical accuracy. From this table, the
multi-mesh approach demonstrates its efficiency compared to the single-mesh
method.

Table 2: Comparison on LiH with respect to results and serial computational time. The referenced total energy is $E_{\mathrm{ref}}=-7.9195$ Hartree.  | $E_{tot}$ | $\Delta E/2$ | $N_{KS}$ | $N_{Har}$ | $t_{all}$ | $t_{SCF}$ | $t_{MA}$
---|---|---|---|---|---|---|---
single-mesh | -7.9178 | 0.0008 | 6,029,109 | 6,029,109 | 18613.60 | 13279.79 | 5333.82
multi-mesh | -7.9181 | 0.0007 | 3,786,064 | 3,665,015 | 14025.13 | 10245.13 | 4262.25

The detailed comparison regarding CPU time is illustrated in Figure 15. As
expected, in the LiH example, the solution of the eigenvalue problem
significantly contributes to the CPU time due to the increased number of
required eigenvalues. Notably, the cost of the third mesh is higher than that
of the final mesh. This is because, during the third mesh adaptation step, the
mesh tends to stabilize with minimal changes, leading to fewer iteration steps
in the final mesh, as shown in Figure 13. Consequently, the CPU time on the
third mesh becomes the most significant part contributing to the total CPU
time. Moreover, since the number of Dofs in the multi-mesh method is smaller
than that in the single-mesh method, the SCF time for the multi-mesh approach
is also reduced. Overall, the multi-mesh method proves to be more efficient in
terms of computational time than the single-mesh method in this case.

Figure 15: CPU time for single-mesh method and multi-mesh method for LiH
molecule. Four mesh adaptations are needed to achieve chemical accuracy in
both methods starting from the same initial setup.

####  5.1.3 BeH2 molecule

The comparison between the multi-mesh method and the single-mesh method is
also presented in the BeH2 example, as displayed in Table 3 and Figure 16. A
similar conclusion to the previous examples on the accuracy and efficiency can
be delivered from the results. In this case, given the same mesh adaption
tolerance, the multi-mesh method succeeds in achieving the chemical accuracy,
while the single-mesh method fails. This observation illustrates the superior
accuracy through the multi-mesh method compared to the single-mesh method.

Table 3: Comparison on BeH2 with respect to results and serial computational time. The referenced value is $E_{\mathrm{ref}}=-15.6606$ Hartree.  | $E_{tot}$ | $\Delta E/3$ | $N_{KS}$ | $N_{Har}$ | $t_{all}$ | $t_{SCF}$ | $t_{MA}$
---|---|---|---|---|---|---|---
single-mesh | -15.6573 | 0.0011 | 7,338,638 | 7,338,638 | 42445.26 | 32331.87 | 10113.40
multi-mesh | -15.6580 | 0.0009 | 4,852,528 | 4,276,863 | 32133.92 | 23925.93 | 8207.99
Figure 16: CPU time for single-mesh method and multi-mesh method for BeH2
molecule. Four mesh adaptations are needed to achieve chemical accuracy in
both methods starting from a same initial setup.

###  5.2 Collection of results for more general electronic structures

In this subsection, we list the results of atoms and molecules using multi-
mesh adaptive methods. The referenced values are generated from the state-of-
art software NWChem using _aug-cc-pv5z_ basis sets except for the helium atom
in which the _aug-cc-pv6z_ basis set is adopted. According to Table 4, we
observe that chemical accuracy is achieved in the first four examples.
However, as the number of atoms increases, the accuracy diminishes. This
degradation can be attributed to the limitations of computational resources.
In such scenarios, we try to fully utilize our computational resources for
simulations using both the single-mesh and multi-mesh adaptive methods, which
can be implemented by dynamically adjusting the mesh adaption tolerance. Hence
a comparative assessment between the two approaches can be delivered.

Table 4: List of examples. The referenced value is obtained from NWChem. $N_{s,\mathrm{KS}}$ represents the number of mesh grids for the single mesh method. $E_{s,\mathrm{KS}}$ stands for the energy obtained from the single-mesh method. The last row refers to the comparison with chemical accuracy.  molecule | He | LiH |  BeH2 |  CH4 |  CH3OH |  CH3CH2OH |  C6H6
---|---|---|---|---|---|---|---
$N_{s,\mathrm{KS}}$ | 2,874,807 | 6,029,109 | 7,338,638 | 8,121,252 | 7,552,976 | 7,720,132 | 7,210,851
$N_{\mathrm{KS}}$ | 1,262,599 | 3,786,064 | 4,852,528 | 5,440,384 | 6,226,903 | 6,864,714 | 6,036,918
$N_{\mathrm{Har}}$ | 2,185,920 | 3,665,015 | 4,276,863 | 5,322,827 | 4,063,275 | 3,593,622 | 4,798,747
$E_{\mathrm{ref}}$ | -2.8348 | -7.9196 | -15.6606 | -40.1198 | -114.8503 | -153.8158 | -230.1916
$E_{s,\mathrm{KS}}$ | -2.8341 | -7.9178 | -15.6573 | -40.1105 | -114.8033 | -153.7418 | -229.9824
$E_{\mathrm{KS}}$ | -2.8340 | -7.9181 | -15.6580 | -40.1147 | -114.8205 | -153.7734 | -230.0663
$|\Delta E_{s,\mathrm{KS}}/E_{\mathrm{ref}}|$ | 0.0002 | 0.0002 | 0.0002 | 0.0002 | 0.0004 | 0.0005 | 0.0009
$|\Delta E_{\mathrm{KS}}/E_{\mathrm{ref}}|$ | 0.0003 | 0.0002 | 0.0002 | 0.0001 | 0.0003 | 0.0003 | 0.0005
$|\Delta E_{s,\mathrm{KS}}/n_{atom}|$ | 0.0007 | 0.0009 | 0.0011 | 0.0019 | 0.0078 | 0.0082 | 0.0174
$|\Delta E_{\mathrm{KS}}/n_{atom}|$ | 0.0008 | 0.0008 | 0.0009 | 0.0010 | 0.0050 | 0.0047 | 0.0104

Figure 17: Sliced KS meshes (top) and isosurfaces (bottom) for the following
molecules: LiH, BeH2, CH4, CH3OH, and CH3CH2OH. Contours are depicted for the
two-dimensional structures (LiH and BeH2). The computational domain is set as
$[-10,10]^{3}$, with the exception of the CH4 example, for which it is
extended to $[-20,20]^{3}$.

The results are summarized in Table 4 and we display several meshes and
isosurfaces for these examples in Figure 17. From the table and the figure, we
can see that

  * 1.

For the last four molecules, besides a similar observation for the CPU time in
Figures 14, 15, and 16, another observation can be made, that with comparable
computational resource, our multi-mesh method can generally produce a much
better result. It means that the multi-mesh method would be an answer to the
question of how to generate the most accurate result, by fully using the given
computational resource.

  * 2.

When maximizing our computational resources, the multi-mesh approach yields
more precise results. To be specific, the multi-mesh method attains chemical
accuracy for the first four examples, as indicated in the last row of Table 4.
In contrast, the single-mesh method only achieves chemical accuracy in the
first two examples.

  * 3.

In the multi-mesh method, the number of mesh grids for the Hartree mesh
exceeds the number for the KS mesh only in the case of helium. As the system
size scales up, the number of Hartree mesh grids becomes progressively less
than the number of KS mesh grids. Consequently, the computational time needed
for the Hartree potential becomes less significant as the system size
increases. There is an exception in the cases of CH4 and C6H6, which arises
because we opted for a larger computational domain, specifically,
$[-20,20]^{3}$.

##  6 Conclusion

In this paper, we introduce a multi-mesh adaptive finite element framework for
Kohn-Sham density functional theory, aiming to achieve chemical accuracy. We
investigate the impact of the Hartree potential approximation on total energy,
finding that chemical accuracy cannot be attained without a well-approximated
Hartree potential. While the single-mesh adaptive method, considering the
Hartree potential, achieves chemical accuracy, it comes with significant
computational costs. To address this, we propose the multi-mesh adaptive
method, which offers a more efficient route to achieving chemical accuracy.

We demonstrate the effectiveness and accuracy of the multi-mesh adaptive
method through various numerical examples. However, we observe that even for
the helium atom, achieving chemical accuracy demands a substantial number of
mesh grids, mainly due to the use of linear Lagrange finite elements. One
strategy is to adopt higher-order elements under the multi-mesh framework.
Additionally, we consider future work on conducting numerical analysis for the
multi-mesh adaptive method. Moreover, we recognize the need for accelerating
the multi-mesh adaptive approach, particularly through parallel implementation
and the development of more effective error indicators for mesh adaptation
methods. These enhancements will further refine the efficiency and accuracy of
the multi-mesh adaptive framework.

## Acknowledgement

The work of Y. Kuang was supported in part by the National Natural Science
Foundation of China (Grant No.12201130) and the Guangzhou Municipal Science
and Technology Bureau (Grant No.2023A04J1321). The work of Y. Shen was
supported by Hunan Key Laboratory for Computation and Simulation in Science
and Engineering (Grand No. LCSSE202307). The work of G. Hu was supported by
the National Natural Science Foundation of China (Grant Nos.11922120 and
11871489), the FDCT of MacaoSAR (Grant No.0082/2020/A2), the MYRG of the
University of Macau (Grant No. MYRG2020-00265-FST), and the Guangdong-Hong
Kong-Macao Joint Laboratory for Data Driven Fluid Mechanics and Engineering
Applications (Grant No. 2020B1212030001).

## References

  * [1] Reinhart Ahlrichs, Michael Bär, Marco Häser, Hans Horn, and Christoph Kölmel.  Electronic structure calculations on workstation computers: The program system turbomole.  Chemical Physics Letters, 162(3):165–169, 1989.
  * [2] Gang Bao, Guanghui Hu, and Di Liu.  An $h$-adaptive finite element solver for the calculations of the electronic structures.  Journal of Computational Physics, 231(14):4967–4979, 2012.
  * [3] Huajie Chen, Xiaoying Dai, Xingao Gong, Lianhua He, and Aihui Zhou.  Adaptive finite element approximations for Kohn–Sham models.  Multiscale Modeling & Simulation, 12(4):1828–1869, 2014.
  * [4] Or Cohen, Leeor Kronik, and Achi Brandt.  Locally refined multigrid solution of the all-electron Kohn–Sham equation.  Journal of Chemical Theory and Computation, 9(11):4744–4760, 2013\.
  * [5] Denis Davydov, Toby D Young, and Paul Steinmann.  On the adaptive finite element analysis of the Kohn–Sham equations: methods, algorithms, and implementation.  International Journal for Numerical Methods in Engineering, 106(11):863–888, 2016.
  * [6] Yana Di and Ruo Li.  Computation of dendritic growth with level set model using a multi-mesh adaptive finite element method.  Journal of Scientific Computing, 39(3):441–453, 2009.
  * [7] Yana Di and Xiao-Ping Wang.  Precursor simulations in spreading using a multi-mesh adaptive finite element method.  Journal of Computational Physics, 228(5):1380–1390, 2009.
  * [8] Lenka Dubcova, Pavel Solín, J Cerveny, and P Kus.  Space and time adaptive two-mesh $hp$-finite element method for transient microwave heating problems.  Electromagnetics, 30(1-2):23–40, 2010.
  * [9] Jun Fang, Xingyu Gao, and Aihui Zhou.  A Kohn–Sham equation solver based on hexahedral finite elements.  Journal of Computational Physics, 231(8):3166–3180, 2012.
  * [10] Stefano Giani and Pavel Solín.  Solving elliptic eigenproblems with adaptive multimesh $hp$-FEM.  Journal of Computational and Applied Mathematics, 394:113528, 2021\.
  * [11] Xianliang Hu, Ruo Li, and Tao Tang.  A multi-mesh adaptive finite element approximation to phase field models.  Communications in Computational Physics, 5(5):1012–1029, 2009.
  * [12] John David Jackson.  Classical electrodynamics.  American Association of Physics Teachers, 1999.
  * [13] Lauri Lehtovaara, Ville Havu, and Martti Puska.  All-electron density functional theory and time-dependent density functional theory with high-order finite elements.  The Journal of Chemical Physics, 131(5):054103, 2009.
  * [14] Ruo Li.  On multi-mesh $h$-adaptive methods.  Journal of Scientific Computing, 24(3):321–341, 2005.
  * [15] Yvon Maday.  $h$-$P$ finite element approximation for full-potential electronic structure calculations.  In Partial Differential Equations: Theory, Control and Approximation, pages 349–377. Springer, 2014.
  * [16] Miguel A L Marques, Micael J T Oliveira, and Tobias Burnus.  Libxc: A library of exchange and correlation functionals for density functional theory.  Computer Physics Communications, 183(10):2272–2281, 2012.
  * [17] Phani Motamarri, Sambit Das, Shiva Rudraraju, Krishnendu Ghosh, Denis Davydov, and Vikram Gavini.  DFT-FE–A massively parallel adaptive finite-element code for large-scale density functional theory calculations.  Computer Physics Communications, 246:106853, 2020.
  * [18] Phani Motamarri, Michael R Nowak, Kenneth Leiter, Jaroslaw Knap, and Vikram Gavini.  Higher-order adaptive finite-element methods for Kohn–Sham density functional theory.  Journal of Computational Physics, 253:308–343, 2013.
  * [19] JE Pask, N Sukumar, and SE Mousavi.  Linear scaling solution of the all-electron coulomb problem insolids.  International Journal for Multiscale Computational Engineering, 10(1), 2012.
  * [20] John A Pople.  Nobel lecture: Quantum chemical models.  Reviews of Modern Physics, 71(5):1267, 1999.
  * [21] Volker Schauer and Christian Linder.  All-electron Kohn–Sham density functional theory on hierarchic finite element spaces.  Journal of Computational Physics, 250:644–664, 2013.
  * [22] Pavel Solín, J Cerveny, Lenka Dubcova, and David Andrs.  Monolithic discretization of linear thermoelasticity problems via adaptive multimesh $hp$-FEM.  Journal of Computational and Applied Mathematics, 234(7):2350–2357, 2010.
  * [23] Eiji Tsuchida and Masaru Tsukada.  Electronic-structure calculations based on the finite-element method.  Physical Review B, 52(8):5573, 1995.
  * [24] Marat Valiev, Eric J Bylaska, et al.  NWChem: A comprehensive and scalable open-source solution for large scale molecular simulations.  Computer Physics Communications, 181(9):1477–1489, 2010.
  * [25] Rüdiger Verfürth.  A posteriori error estimation techniques for finite element methods.  OUP Oxford, 2013.
  * [26] Seymour H Vosko, Leslie Wilk, and Marwan Nusair.  Accurate spin-dependent electron liquid correlation energies for local spin density calculations: a critical analysis.  Canadian Journal of Physics, 58(8):1200–1211, 1980.
  * [27] Erich Wimmer, Henry Krakauer, Michael Weinert, and AJ Freeman.  Full-potential self-consistent linearized-augmented-plane-wave method for calculating the electronic structure of molecules and surfaces: O2 molecule.  Physical Review B, 24(2):864, 1981.
  * [28] Shengyang Wu, Xianliang Hu, and Shengfeng Zhu.  A multi-mesh finite element method for phase-field based photonic band structure optimization.  Journal of Computational Physics, 357:324–337, 2018.
  * [29] Chao Yang, Juan C Meza, and Lin-Wang Wang.  A trust region direct constrained minimization algorithm for the Kohn–Sham equation.  SIAM Journal on Scientific Computing, 29(5):1854–1875, 2007.