An adaptive planewave method for electronic structure calculations

We consider a Schrödinger-type linear eigenvalue problem: Find $(\lambda_{i},\varphi_{i})\in\mathbb{R}\times H^{1}_{\#}(\Omega)$, such that $\|\varphi_{i}\|_{L_{\#}^{2}(\Omega)}=1$ and

where the potential $V\in C^{\infty}(\mathbb{R}^{d})$ and is $\mathbb{L}$-periodic. The smoothness assumption on the potential $V$ is reasonable, since the planewave discretization is in practice often combined with smooth pseudo-potentials. It implies in particular that $V$ is bounded below.

The weak form of Eq. 1 reads

with the bilinear form $a(\cdot,\cdot):H^{1}_{\#}(\Omega)\times H^{1}_{\#}(\Omega)\rightarrow\mathbb{R}$ given by

Since the operator $A:=-\Delta+V$ is self-adjoint, bounded below with compact resolvent, there exists a non-decreasing sequence of eigenvalues $\lambda_{1}<\lambda_{2}\leq\cdots\leq\lambda_{j}\rightarrow+\infty$, where the $\lambda_{j}$’s are repeated according to their multiplicity.

Up to shifting the operator by a positive constant, e.g. $\displaystyle\max\big{\{}1-\inf_{\boldsymbol{r}\in\Omega}V(\boldsymbol{r}),0\big{\}}$, we can further assume that $V(\boldsymbol{r})\geq 1~(\forall~\boldsymbol{r}\in\Omega)$, which ensures that all eigenvalues of $A$ are positive. Note that such a shifting has no impact on the eigenfunctions of the problem. Then it holds

and we can define the corresponding energy norm by $\|\cdot\|_{a}:=\sqrt{a(\cdot,\cdot)}$.

Given an energy cut-off $E_{\rm c}>0$, we define the following finite dimensional space

For $v\in H^{s}_{\#}(\Omega)$, the best approximation of $v$ in $X_{E_{\rm c}}(\Omega)$ is

for any $H^{t}_{\#}$-norm ($t\leq s$). The more regularity $v$ has, the faster this truncated series converge to $v$: For $s,t\in\mathbb{R}_{+}$ satisfying $t\leq s$, we have that for each $v\in H^{s}_{\#}(\Omega)$ (see e.g. [19]),

For simplicity of notation, we suppress for now the index for the considered eigenpair in the eigenvalue problem. We rewrite the linear problem Eq. 2 as

and present available a priori and a posteriori error estimates. We note that the following discussions are not restricted to a specific eigenpair, but require a gap between the considered eigenvalue and the surrounding ones.

Given an energy cut-off $E_{\rm c}$, the planewave discretization is a Galerkin approximation of the eigenpairs of Eq. 6 within the finite dimensional subspace $X_{E_{\rm c}}(\Omega)\subset H_{\#}^{1}(\Omega)$: Find $(\lambda_{E_{\rm c}},\varphi_{E_{\rm c}})\in\mathbb{R}\times X_{E_{\rm c}}(\Omega)$, such that $\|\varphi_{E_{\rm c}}\|_{L^{2}_{\#}(\Omega)}=1$ and

An a priori estimate for the solutions $(\lambda_{E_{\rm c}},\varphi_{E_{\rm c}})$ of Eq. 7 can be found in [43, Theorems 2 and 3]. There exists a constant $C>0$ independent of $E_{\rm c}$, such that

Hence, if the eigenfunction $\varphi$ is analytic, which we will assume in the following, then the approximation error decays exponentially as (see [19, Section 5.4, Eq. (5.4.5)])

with the constants $C,c,\tilde{c}>0$ independent of $E_{\rm c}$. The above a priori error estimates not only guarantee that the error goes to zero as the energy cut-off increases, but also offer the speed of convergence of the approximate solutions to the exact one.

We observe from the numerical experiments (see Fig. 2 in Section 5.1) that for the smooth tested potentials, the planewave approximation errors decay exponentially fast with respect to $\sqrt{E_{\rm c}}$, which matches the a priori error estimates in Eqs. 10 and 11.

In Section 3, we will combine the a priori bounds Eqs. 10 and 11 with the a posteriori bounds presented below in a linear regression strategy for determining good energy cut-offs.

To develop an adaptive method, it is very useful to have a posteriori error estimators available. Unlike a priori error estimates, the goal of the a posteriori error estimates is to provide a computable bound on the error between the numerical approximation of the solution and the unknown exact solution. This bound should therefore be computable with the knowledge of the approximate solution and the parameters used for the computation only.

In this section, we review an a posteriori error estimator for the planewave approximation of the linear eigenvalue problem Eq. 1, which gives computable and asymptotically accurate approximate errors. This a posteriori error indicator was first introduced for simple eigenvalues for conforming finite element discretizations in [16], nonconforming methods in [17], and then extended to eigenvalue clusters in [18].

Let us denote $H^{1}_{\#}(\Omega)$ by $H$ and its dual by $H^{\prime}$ for simplicity of notations. Let $(\lambda,\varphi)\in\mathbb{R}\times H$ be an eigenpair of Eq. 6 and $(\lambda_{E_{\rm c}},\varphi_{E_{\rm c}})\in\mathbb{R}\times X_{E_{\rm c}}(\Omega)$ be the corresponding planewave approximation with a given energy cut-off $E_{\rm c}$. We define the residual ${\rm Res}(\lambda_{E_{\rm c}},\varphi_{E_{\rm c}})\in H^{\prime}$ by

The idea of residual-based a posteriori error estimation is to estimate the approximation error by the dual norm of the residual in an appropriate Hilbert space. Note that the dual norm must be well chosen, otherwise the norm of the residual might not accurately represent the error or not even go to zero when the error goes to zero. To construct the a posteriori error estimator for the approximation $(\lambda_{E_{\rm c}},\varphi_{E_{\rm c}})$, we consider the dual norm of the residual

where the dual norm (with respect to the energy norm $\|\cdot\|_{a}$) is defined for $f\in H^{\prime}$ by

For simplicity of notation, we will denote the residual and its dual norm by

in the following. We have from [18] that if $\lambda$ is a simple eigenvalue, there exists $\alpha(E_{\rm c})$ with $\alpha(E_{\rm c})\rightarrow 0$ as $E_{\rm c}\rightarrow\infty$ such that

To evaluate the a posteriori error indicator Eq. 13, one needs to solve the infinite dimensional linear system Eq. 12 to compute the residual. A practical evaluation method is to restrict the trial function space and test function space in Eq. 12 to a finite dimensional subspace $X_{E_{\rm g}}(\Omega)$ (see definition Eq. 4) with some cut-off $E_{\rm g}\gg E_{\rm c}$. We will then denote the computable error indicator by $\eta_{E_{\rm c}}^{E_{\rm g}}$, whose detailed calculations are presented in Appendix A. If $E_{\rm g}$ is proportional to $E_{\rm c}$, then a standard calculation by solving the corresponding finite dimensional linear system will require a cost of $\mathcal{O}(E_{\rm c}^{d/2}\ln E_{\rm c})$, and we refer to Section A.1 for the explanations.

In order to reduce the cost of the error indicator, we can alternatively evaluate the residual via a perturbation-based method. The resulting error indicator will be denoted by $\eta_{E_{\rm c}}^{E_{\rm g},[k]}$ (with $k=1,2$ the order of perturbation), whose cost is reduced by avoiding the resolution of a linear system. We refer to Section A.2 for the derivation of this method and discussion of its computational cost.

In the ab-initio quantum mechanical modeling of many electron systems, the Kohn–Sham density functional theory (DFT) [34, 36] has established itself as one of the most powerful electronic structure methods due to its good balance between accuracy and computational cost.

We consider a many-particle system (with $N_{\rm atom}$ nuclei and $N_{\rm e}$ electrons) in the spin-less and non-relativistic setting. The Kohn–Sham ground state solution can be obtained by solving the following Kohn–Sham equations [40]: Find $(\lambda_{i},\varphi_{i})\in\mathbb{R}\times H_{\#}^{1}(\Omega),~i=1,\cdots,N_{\rm b}$, such that $\int_{\Omega}\varphi_{i}\varphi_{j}=\delta_{ij}$ and

where $N_{\rm b}$ is the number of eigenpairs/bands under consideration, $\{\lambda_{i}\}_{i=1}^{N_{\rm b}}$ are the $N_{\rm b}$ lowest eigenvalues, $\rho=\sum_{i=1}^{N_{\rm b}}f_{i}|\varphi_{i}|^{2}$ is the electron density with $f_{i}$ the occupancy number, and

with $V_{\rm{eff}}[\rho]$ the effective potential, $v_{\rm{ps}}$ the pseudo-potential generated by the nuclei and core electrons, $v_{\rm{H}}[\rho]$ and $v_{\rm{xc}}[\rho]$ the so-called Hartree potential and exchange-correlation potential, respectively. The occupancy numbers $f_{i}$ can be chosen by smearing methods like the Fermi–Dirac or the Methfessel–Paxton scheme.

We shall now give more details on the potentials used to define $V_{\rm eff}$. In the norm-conserving pseudo-potential setting [40], the pseudo-potential $v_{\rm{ps}}$ can be written as

where $v_{\rm loc}:\mathbb{R}^{d}\rightarrow\mathbb{R}$ is the local part, $v_{\rm nl}$ is a non-local operator given by

with $M_{\rm ps}\in\mathbb{Z}_{+}$ and $\xi_{I,j}~(I=1,\cdots,N_{\rm atom},~j=1,\cdots,M_{\rm ps})\in L_{\#}^{2}(\Omega)$ being sufficiently smooth functions in $\mathbb{R}^{d}$. In the context of ultra-soft pseudo-potentials (USPP) or PAW methods [39, 53], the non-local parts are formulated by

where $D_{nm}^{I}$ depends on the eigenfunctions and need to be updated during the SCF loop. For more details on the PAW method, which will be the model used in practice for the numerical experiments in Section 5, see Appendix B.

In periodic systems, $v_{\rm{H}}[\rho]$ represents the $\mathbb{L}$-periodic Coulomb potential generated by the $\mathbb{L}$-periodic electron density $\rho$

Finally, for the exchange-correlation potential, we use a generalized gradient approximation (GGA) [46] in our numerical experiments.

Problem Eq. 19 is a nonlinear eigenvalue problem since the operator $H[\rho]$ depends on the electron density $\rho$, and hence on the eigenfunctions $\{\varphi_{i}\}_{i=1}^{N_{\rm b}}$. A self-consistent field (SCF) iteration algorithm [40] is commonly resorted to for this kind of nonlinear eigenvalue problem. At each iteration, a new Hamiltonian is constructed from a trial electronic state and a linear eigenvalue problem, in the form of Schrödinger-type equation Eq. 1, is then solved to obtain the low-lying eigenvalues and corresponding eigenvectors.

Once the ground state solution of Kohn–Sham equations have be obtained, we can evaluate the band structure energy $E_{\textrm{bands}}\equiv\sum_{i=1}^{N_{b}}f_{i}\lambda_{i}$ as a weighted sum of the eigenvalues by occupation numbers. When we consider the $k$-point sampling, the integration over the Brillouin zone also needs to be included.

In the first iterations of the self-consistent field algorithm, the iterates are far form the exact solution, so that the error coming from the iterative solver is dominant. But close to self-consistency, the choice of the basis set becomes important, as the basis set size will ultimately determine the quality of the final approximation. Therefore, an efficient planewave SCF algorithm should be able to adjust the energy cut-off on the fly, so that no computational resource is wasted while the required accuracy can be reached at the end of the iterations. We refer to [30] for a similar discussion of error balancing for the Gross-Pitaevskii equation.

We now discuss the construction of error indicators and the adaptive algorithm for the nonlinear Kohn–Sham equations Eq. 19. As mentioned in Section 4.1, an SCF iteration algorithm is used to solve the nonlinear eigenvalue problem, in which a linear Schrödinger-type equation is solved at each step. Therefore, we can incorporate Algorithm 1 (for linear problems) into each SCF iteration, and thus adapt the energy cut-off on the fly with some well-chosen tolerance, corresponding to the self-consistency error. The use of the self-consistency error as the tolerance for adapting the cut-off is based on the following: when the solution is far from self-consistency, one can choose a relative small energy cut-off for the linearized eigenvalue problem; when the iteration is close to self-consistency, one has to use a large energy cut-off to obtain accurate eigenpairs.

Let $\rho_{\rm in}^{(m)}$ and $\rho_{\rm out}^{(m)}$ represent the input and output electron densities at the $m$-th step of the SCF algorithm. At the $m$-th step, a linearized eigenvalue problem, with the trial input electron density $\rho_{\rm in}^{(m)}$ is solved

with some energy cut-off $E_{\rm c}^{(m)}$. Then an output electron density $\rho_{\rm out}^{(m)}$ is obtained from the eigenfunctions $\big{\{}\varphi_{i}^{(m)}\big{\}}_{i=1}^{N_{\rm b}}$. The self-consistency error is defined by

The parameter $\alpha$ in Eq. 24 is set to control the ratio between discretization error and self-consistency error in the total error. In our numerical simulations, $\alpha$ is set empirically: we take $\alpha$ to be the ratio between the tolerances of the discretization error and self-consistency error. Note that a good choice of $\alpha$ could be obtained by a careful error analysis of the nonlinear eigenvalue problems, estimating explicit optimal weights (in the total error) of the discretization error and self-consistent error. This will be investigated in future work.

The planewave discretization error indicator is defined with respect to the linear eigenvalue problem Eq. 23, which involves all $N_{\rm b}$ eigenpairs under consideration, as discussed in Remark 2.2. In analogy with definition Eq. 13, we define the discretization error indicator as

where $f_{i}$ is the occupancy number, $\|\cdot\|_{a^{\prime}}$ is the dual norm corresponding to the linear operator

and ${\rm Res}^{(m)}_{i}$ is the residual for the $i$-th eigenpair at the $m$-th step

With a given $E_{\rm g}\gg E_{\rm c}^{(m)}$, we define the computable discretization error indicator by

based on the standard calculation Eq. 33, and for $k=1,2$.

based on a perturbation method Eq. 38.

Combining the SCF algorithm and Algorithm 1 for linear eigenvalue problems, we propose the following adaptive algorithm for nonlinear eigenvalue problems.

In Step 4 of Algorithm 2, by taking $\left(\eta_{\rm SCF}^{(m)}\right)^{2}$ as the tolerance for planewave discretizations, we are actually estimating the total energy error of the full nonlinear eigenvalue problem by the sum $\left(\eta_{E_{\rm c}}^{(m)}\right)^{2}+\left(\eta_{\rm SCF}^{(m)}\right)^{2}$. We note that it could be possible to construct some a posteriori error estimate for the nonlinear eigenvalue problems directly (see [14, 23, 25, 30]). However, we will not use or discuss this type of estimates in this work.

In Step 5 of Algorithm 2, one can exploit a charge mixing scheme to ensure or accelerate the convergence of the SCF algorithm. The charge mixing in our implementation is performed in the Fourier space. Assume that we need to mix two charge densities $\rho_{1}$ and $\rho_{2}$ with cut-offs $E_{\rm c}^{1}$ and $E_{\rm c}^{2}$, respectively, and $E_{\rm c}^{2}>E_{\rm c}^{1}$. We simply fill the rest of planewave components of $\rho_{1}$ by zeros to expand it to the same grid as $\rho_{2}$ and then carry out the mixing. The additional computational cost of this implementation is negligible.

The most important feature of Algorithm 2 is that it is adaptive. It automatically determines a good discretization parameter, in our case the energy cut-off, on the fly to achieve the required accuracy. Moreover, compared to the conventional SCF iterations, Algorithm 2 avoids solving large linear eigenvalue problems at all iteration steps. Since the energy cut-offs are determined on the fly and increased during the SCF iterations, computational cost can be saved from that one only needs to solve small eigenvalue problems while the solution is still away from convergence. The extra cost of this adaptive algorithm comes from calculating the error indicator, which scales as $\mathcal{O}(N_{\rm b}E_{\rm g}^{d/2}\ln E_{\rm g})$. Since the cost for one linear eigensolver is $\mathcal{O}(N_{\rm b}^{2}E_{\rm c}^{d/2})$ and $E_{\rm g}$ is proportional to $E_{\rm c}$, we have that the additional cost for calculating the error indicator is way smaller than the cost saved from the linear eigensolvers used with small energy cut-offs.

A linear eigenvalue problem. Consider the 2D linear eigenvalue problems of the form Eq. 1, with $\Omega=[-5,5]^{2}$ and $V(x,y)=1.5(x^{2}+y^{2})+e^{-(x-1)^{2}-y^{2}}$. We use a large energy cut-off $E_{\rm c}=100.0$ to obtain a sufficiently accurate solution as the reference in this example.

We first compare in Fig. 2 the numerical errors in the lowest two eigenvalues with different a posteriori error indicators discussed in Sections 2.4 and A. We observe in the plots that they match almost perfectly. Moreover, the errors decay exponentially with respect to the energy cut-off which suggests that a linear regression strategy could work well. We also observe that in this case, the first order perturbation method (Eq. 38 with $k=1$) to compute the error indicator already leads to a very accurate approximation of the true error.

We further test the adaptive algorithm (Algorithm 1) on this linear problem. Tables 1 and 2 show the energy cut-offs obtained by using the error indicators Eq. 34 and Eq. 38 respectively, with Strategy A and Strategy B, and tolerances $\varepsilon=10^{-3}$ and $10^{-6}$. We obtain by performing a traditional convergence test that the optimal cut-offs for these two tolerances are 7.1 and 16.0 respectively. Algorithm 1 captures energy cut-offs respectively of 7.5 and 17.12 in a few steps, which suffice for the target accuracy without overshooting too much the optimal cut-offs.

A nonlinear eigenvalue problem. Consider the following Gross–Pitaevskii equation (GPE) originated from modelling Bose–Einstein condensates [47]: Find $(\lambda,\varphi)\in\mathbb{R}\times H^{1}_{\#}(\Omega)$, such that $\|\varphi\|_{L^{2}(\Omega)}=1$ and

where $\mu=1.0$, $\Omega=[-5,5]^{2}$, $V(x,y)=10(x^{2}+y^{2})+e^{-\big{(}(x-1)^{2}+y^{2})\big{)}}$, and $\rho(\boldsymbol{r})=|\varphi(\boldsymbol{r})|^{2}$. Note that this equation can be viewed as a simplified version of Eq. 19, where only the lowest eigenvalue is considered and the nonlinear term is much simpler.

We apply Algorithm 2 to solve Eq. 29 with tolerance $\delta=10^{-6}$. Note that the “best” choice of the parameter $\alpha$ (defined in Eq. 24) is not known in this example, and we perform experiments with $\alpha=1.0$ and 0.1, respectively, to illustrate how the choice of $\alpha$ affects the number of convergence steps of our algorithm. We show in Fig. 3 the discretization and SCF errors with respect to the SCF steps. The energy cut-off is iteratively adjusted with Algorithm 2 and we observe that a good balance is achieved between the SCF and discretization errors. The oscillation in the SCF errors arises from the simple mixing scheme used in the SCF iterations.

We implemented the a posteriori error indicator and the adaptive algorithm (Algorithm 2) in an in-house planewave PAW code [31]. For the validity tests we chose two bulk systems. The first is a diamond structure with two carbon atoms in the unit cell. This is an insulator and we used only the $\Gamma$ point for the $k$-point sampling. The number of occupied bands is $N_{\textrm{b}}=4$. The second system is an FCC structure with four aluminum atoms in the unit cell. It is a metallic system with a vanishing band gap, and we employed a $6\times 6\times 6$ Monkhorst–Pack $k$-mesh. The Brillouin zone grid is fixed over all calculations. The first order Methfessel–Paxton smearing method [41] with a smearing width of $0.3\,\textrm{eV}$ was used and the number of occupied bands is $N_{\textrm{b}}=7$. For the efficiency tests on relatively large systems, we constructed a $4\times 4\times 4$ diamond supercell containing 128 carbon atoms, and the number of occupied bands is $N_{\textrm{b}}=256$.

In the following, we first check the accuracy of the a posteriori error indicators on the linearized Kohn–Sham equations and then test the behavior of the adaptive algorithm for solving the nonlinear equations. In order to take into account all the $N_{\textrm{b}}$ occupied eigenstates, we present error results on the band structure energy $E_{\textrm{bands}}$. Note that the discretization error in the band structure energy and that in the total energy decay in the same rate with respect to $E_{c}$ (see [13, Theorem 4.2 and (4.90)]).

Accuracy of the error indicators for the linearized Kohn–Sham equations. We are primarily concerned with the accuracy of the error indicator for the linearized Kohn–Sham equation since it determines the effectiveness of the adaptive algorithm. In our tests we fixed the density as the superposition of atomic charge densities and solved the linearized Kohn–Sham equations for different energy cut-off values $E_{\rm c}$. The a posteriori error indicators were calculated by using the standard and the perturbation-based methods detailed in Appendix A.

The error indicators $|\eta_{E_{\rm c}}^{E_{\rm g}}|^{2}$, $|\eta_{E_{\rm c}}^{E_{\rm g},[1]}|^{2}$, and $|\eta_{E_{\rm c}}^{E_{\rm g},[2]}|^{2}$ are compared with real errors in Tables 3 and 4 for the diamond structure and the aluminum systems, respectively. It is observed from the tables that the first-order indicator $|\eta_{E_{\rm c}}^{E_{\rm g},[1]}|^{2}$ already provides a very reliable approximation of the real error. The second-order indicator $|\eta_{E_{\rm c}}^{E_{\rm g},[2]}|^{2}$ exhibits a clear improvement over the first-order indicator and derives results very close to those from the standard indicator. And the errors estimated by the standard indicator $|\eta_{E_{\rm c}}^{E_{\rm g}}|^{2}$ closely match the real error.

In terms of cost, however, the standard calculation of the error indicator is remarkably more expensive than the perturbation-based methods, as shown by the timings for the diamond example in Table 5. As mentioned, the standard calculation requires to solve a linear system of size $N_{E_{\rm g}}$ for each of the $N_{\textrm{b}}$ eigenvalues. Since the planewave discretization leads to large-scale dense matrices, the solution of the corresponding linear systems can be very demanding.

Therefore, the perturbation-based calculation of the a posteriori error indicator is the proper choice in practice to balance accuracy versus computational cost. We employed the perturbation-based methods in the following tests. It should be emphasized that the proposed error estimate is accurate both for the tested insulator and metallic systems. This high-quality a posteriori error indicator is critical to a reliable tuning of the energy cut-off.

Adaptive solution of the Kohn–Sham equations with an unknown converged cut-off. We now apply Algorithm 2 to solve the nonlinear Kohn–Sham equations. The density mixing is performed using Pulay’s DIIS method [48]. The tolerance for the discretization error in $E_{\textrm{bands}}$ and that for the SCF error were chosen to be $2\times 10^{-3}$ eV for both the diamond structure and the aluminum systems.

The evolution of the discretization and SCF errors with respect to the SCF step are shown in Fig. 4 for the diamond and the aluminum systems. For a detailed illustration, we take the diamond structure and further present detailed data on the energy cut-off and the two components of errors in Table 6. It is seen that a converged solution is obtained by adaptively tuning the energy cut-off depending on the discretization and SCF errors.

Comparison with convergence-test computations and discussion. We have stated in the introduction that the purpose of our algorithm is to determine the converged energy cut-off adaptively and avoid the conventional convergence-test computations. In this section we compare and discuss the two approaches in detail.

In a typical convergence test on the energy cut-off $E_{\rm c}$, first, one needs to choose three parameters: the smallest and largest cut-offs to be tested and a fixed step size of values in between. The largest value should be greater than the converged cut-off, but the latter is unknown in advance. So the largest cut-off needs to be set appropriately with some redundancy. Then the Kohn–Sham equations are solved under all the chosen cut-offs, and the converged cut-off is determined as the one beyond which the variation in results of the concerned quantity (the energies or the forces) is consistently smaller than the given tolerance.

In the adaptive algorithm proposed here, users need to choose the starting cut-off and the ratio between the tolerances for planewave discretizations and SCF iterations. At each iteration step, the discretization error associated with $E_{\rm c}$ is estimated with the a posteriori error indicator. Then the estimated error can be used directly to determine the computation is converged or not. If not, the known error data are combined with the a priori error analysis to predict the converged $E_{\rm c}$ in step sizes not limited by the fixed cut-off step like the case of convergence tests.

It can be seen that compared with the convergence-test approach, the proposed algorithm is based on theoretically justified bounds: the error estimates. As for the computational cost, two factors should be taken into account here: (i) the extra computational cost of the a posteriori error estimation and (ii) the number of iterations used to reach the converged solution. The time complexity for calculating the a posteriori error indicator is $\mathcal{O}(N_{\textrm{b}}\,N_{E_{\rm g}}\ln N_{E_{\rm g}})$ and that for the iterative diagonalization of eigenvalue problems is $\mathcal{O}(N_{\textrm{b}}^{2}\,N_{E_{\rm c}})$. When the number of bands $N_{\textrm{b}}$ is relatively large, the additional cost for the error estimation is small compared to the cost for solving eigenvalue problems. We expect that the efficiency advantage of the adaptive algorithm over the convergence-test approach grows with the number of bands. The number of iterations to convergence is well controlled in the adaptive algorithm if the tuning of $E_{\rm c}$ proceeds efficiently by combining the a posteriori error with the a priori analysis and using strategies stated in Section 3.

In many practical simulations, the number of occupied bands $N_{\textrm{b}}$ is large, and the convergence tests are computationally demanding. This is the case for systems with defects or vacancies, surface systems, etc., where large supercells must be used. This happens as well for high-temperature simulations, where the cell may be small but $N_{\textrm{b}}$ grows rapidly with the temperature. In the latter case, besides using a large parameter $N_{\textrm{b}}$, the converged cutoff energy $E_{\rm c}$ is usually significantly larger than the one needed in simulations at ambient temperature [10].

In the following comparison of the convergence-test approach and the adaptive algorithm, we use a $4\times 4\times 4$ diamond supercell containing 128 carbon atoms. The tolerance for the discretization error in $E_{\textrm{bands}}$ and for the SCF iteration error was set respectively to 0.128 eV (i.e. 1 meV/atom) and 0.1 meV. Note that although for simple bulk systems it is not necessary to perform convergence tests using large supercells, we still chose the example here because the scale of the system represents well the computational cost of systems with relatively large $N_{\textrm{b}}$.