Numerical stability and efficiency of response property calculations in density functional theory

Abstract

Response calculations in density functional theory aim at computing the change in ground-state density induced by an external perturbation. At finite temperature, these are usually performed by computing variations of orbitals, which involve the iterative solution of potentially badly conditioned linear systems, the Sternheimer equations. Since many sets of variations of orbitals yield the same variation of density matrix, this involves a choice of gauge. Taking a numerical analysis point of view, we present the various gauge choices proposed in the literature in a common framework and study their stability. Beyond existing methods, we propose a new approach, based on a Schur complement using extra orbitals from the self-consistent field calculations, to improve the stability and efficiency of the iterative solution of Sternheimer equations. We show the success of this strategy on nontrivial examples of practical interest, such as Heusler transition metal alloy compounds, where savings of around 40% in the number of required cost-determining Hamiltonian applications have been achieved.

Appendix A Choosing the number of extra bands

In this paper, we saw through various numerical examples that using a Schur complement to compute the unoccupied–occupied contributions to the orbitals’ response improves the convergence of the Sternheimer equation. In this appendix, we quantify this acceleration and discuss how this idea can be used to select the number of bands to be computed. Considering the straight convergence curves from Figs. 5, 6 suggest that the convergence of the CG is indeed led by the square root of the condition number of the system matrix (see [43, Section 9]) when using the Schur complement. Key idea will thus be to estimate the condition number of the linear system (50).

1.1 A.1 Numerical analysis

To analyse the condition number of the Schur complement, we consider the following specific setting

$$\begin{aligned} {\left\{ \begin{array}{ll} H_\rho \phi _n = \varepsilon _n\phi _n,\quad \varepsilon _1 \leqslant \varepsilon _2 \leqslant \cdots \\ \left\langle \phi _n,\phi _m \right\rangle = \delta _{nm}, \end{array}\right. } \end{aligned}$$

(53)

where $H_\rho \in {{\mathbb {C}}}^{{N_{\textrm{b}}}\times {N_{\textrm{b}}}}$ is typically the discretized self-consistent Hamiltonian of the system, at some k-point. We assume that we have N occupied orbitals that have an occupation number higher than the threshold we fixed and that we have ${N_{\textrm{ex}}}$ extra bands, as explained in Sect. 3.3. In summary, we have at our disposal $N + {N_{\textrm{ex}}}$ bands in total: $\Phi = (\phi _1,\dots ,\phi _N)$ are occupied, fully converged bands and the extra bands $\Phi _{\textrm{ex}}^{\ell } = (\phi ^{\ell }_{N+1},\dots ,\phi ^{\ell }_{N+{N_{\textrm{ex}}}})$ are not necessarily all converged. We added here the exponent $\ell$ as we make the following assumptions:

• for any $\ell \in {{\mathbb {N}}}$, $(\Phi ,\Phi _{\textrm{ex}}^{\ell })$ is an orthonormal family;

• for any $\ell \in {{\mathbb {N}}}$, $(\Phi _{\textrm{ex}}^{\ell })^* H_\rho \Phi _{\textrm{ex}}^{\ell } \in {{\mathbb {C}}}^{{N_{\textrm{ex}}}\times {N_{\textrm{ex}}}}$ is a diagonal matrix whose elements are labelled $\varepsilon _{m}^{\ell } {:}{=}\left\langle \phi _m^{\ell },H_\rho \phi _m^{\ell } \right\rangle$ for $N+1\leqslant m\leqslant N+{N_{\textrm{ex}}}$;

• as $\ell \rightarrow +\infty$, $(\phi ^{\ell }_{m}, \varepsilon ^{\ell }_{m}) \rightarrow (\phi _{m}, \varepsilon _{m})$.

All these assumptions are satisfied, for instance, if the sequence $(\Phi ,\Phi _{\textrm{ex}}^{\ell })_{l\in {{\mathbb {N}}}}$ is generated by any Rayleigh–Ritz-based eigensolver (for instance, the LOBPCG eigensolver [31]), which is the case by default in DFTK. For every $\ell$, we can thus decompose the plane-wave approximation space ${{\mathcal {H}}}= X_{{N_{\textrm{b}}}}$ (with ${N_{\textrm{b}}}\gg N+{N_{\textrm{ex}}}$) in two different ways:

$$\begin{aligned} {{\mathcal {H}}}= \textrm{Ran}(P) \oplus \textrm{Ran}(T) \oplus \textrm{Ran}(R) \quad \text {and}\quad {{\mathcal {H}}}= \textrm{Ran}(P) \oplus \textrm{Ran}(T^{\ell }) \oplus \textrm{Ran}(R^{\ell }), \nonumber \\ \end{aligned}$$

(54)

where

$$\begin{aligned} P {:}{=}\sum _{n=1}^N \phi _n\phi _n^*\quad \text {and}\quad {\left\{ \begin{array}{ll} T {:}{=}\sum _{n=N+1}^{{N_{\textrm{ex}}}} \phi _n\phi _n^*, \quad R {:}{=}1 - P - T \\ T^{\ell } {:}{=}\sum _{n=N+1}^{{N_{\textrm{ex}}}} \phi _n^{\ell }(\phi _n^{\ell })^*, \quad R^{\ell } {:}{=}1 - P - T^{\ell }, \end{array}\right. } \end{aligned}$$

(55)

are all orthogonal projectors. In these two decompositions, $H_\rho$ has the associated block representations:

$$\begin{aligned} H_\rho = \begin{pmatrix} E &{} 0 &{} 0 \\ 0 &{} E_{\textrm{ex}} &{} 0 \\ 0 &{} 0 &{} \ddots \end{pmatrix} \quad \text {and}\quad H_\rho = \begin{pmatrix} E &{} 0 &{} 0 \\ 0 &{} E^{\ell }_{\textrm{ex}} &{} R^{\ell }H_\rho T^{\ell } \\ 0 &{} T^{\ell }H_\rho R^{\ell } &{} R^{\ell }H_\rho R^{\ell } \end{pmatrix} \end{aligned}$$

(56)

where $E {:}{=}\textrm{Diag}(\varepsilon _1,\dots ,\varepsilon _n)$, $E_{\textrm{ex}} {:}{=}\textrm{Diag}(\varepsilon _{N+1},\dots ,\varepsilon _{N+{N_{\textrm{ex}}}})$ and $E_{\textrm{ex}}^{\ell } {:}{=}\textrm{Diag}(\varepsilon ^{\ell }_{N+1}, \dots ,\varepsilon ^{\ell }_{N+{N_{\textrm{ex}}}})$ are diagonal matrices. Moreover, note that as $\Phi ^{\ell }_{\textrm{ex}} \rightarrow \Phi _{\textrm{ex}}$, the residuals $R^{\ell }H_\rho T^{\ell }$ converge to 0.

Now, we fix $n\leqslant N$ and we compute the condition number of the linear system (50). Enforcing the CG to stay at each iteration in $\textrm{Ran}(R^{\ell })$, this condition number is given by the ratio of the largest and smallest nonzero eigenvalues of

$$\begin{aligned} H_n^{\ell } + X_n^{\ell }, \end{aligned}$$

(57)

where

$$\begin{aligned} \begin{aligned} H_n^{\ell } {:}{=}R^{\ell }(H_\rho - \varepsilon _n)R^{\ell } \quad \text {and}\quad X_n^{\ell }&= -R^{\ell }(H_\rho -\varepsilon _n)\Phi _{\textrm{ex}}^{\ell } (E_{\textrm{ex}}^{\ell } - \varepsilon _n)^{-1}\\&(\Phi _{\textrm{ex}}^{\ell })^*(H_\rho - \varepsilon _n)R^{\ell }\\&= -R^{\ell }H_\rho \Phi _{\textrm{ex}}^{\ell }(E_{\textrm{ex}}^{\ell } - \varepsilon _n)^{-1} (\Phi _{\textrm{ex}}^{\ell })^*H_\rho R^{\ell }. \end{aligned} \end{aligned}$$

(58)

Here, $E^{\ell }_{\textrm{ex}}-\varepsilon _n$ is diagonal and thus explicitly invertible if $\ell$ is large enough as $\varepsilon _{N+1}^{\ell }\rightarrow \varepsilon _{N+1}>\varepsilon _N\geqslant \varepsilon _n$. We focus for the moment on the smallest nonzero eigenvalue, that is $\varepsilon ^{\ell }_{N+{N_{\textrm{ex}}}+1} - \varepsilon _n$. The condition number being proportional to the inverse of the smallest eigenvalue, we now derive a lower bound of $\varepsilon ^{\ell }_{N+{N_{\textrm{ex}}}+1} - \varepsilon _n$ in order to get an upper bound on the condition number of (57). When $\ell \rightarrow +\infty$, we have $X_n^{\ell }\rightarrow 0$ (as $RH_\rho \Phi _{\textrm{ex}}=0$) and $H_n^{\ell }\rightarrow H_n{:}{=}R(H_\rho - \varepsilon _n)R$ whose smallest nonzero eigenvalue is $\varepsilon _{N+{N_{\textrm{ex}}}+1} - \varepsilon _n$. We use next a perturbative approach to effectively approximate the condition number of (57).

We use a standard eigenvalue perturbation result, whose proof is recalled for the sake of completeness. It is directly adapted from the general case of self-adjoint bounded below operators with symmetric perturbations studied, for instance, in [11].

Proposition 1

Let $N\in {{\mathbb {N}}}$, $H_0,W\in {{\mathbb {C}}}^{N\times N}$ be Hermitian matrices and $\alpha \geqslant 0$ such that $H_0+\alpha >0$. Then, the eigenvalues of $H{:}{=}H_0+W$ and $H_0$ satisfy

$$\begin{aligned} \left| \nu _i(H) - \nu _i(H_0) \right| \leqslant (\nu _i(H_0)+\alpha )\left\| W \right\| _{H_0,\alpha }, \end{aligned}$$

(59)

where $\left\| W \right\| _{H_0,\alpha }$ is the operator norm of $(H_0+\alpha )^{-1/2}W(H_0+\alpha )^{-1/2}$ and $\nu _i(A)$ is the i-th lowest eigenvalue of the matrix A.

Proof

Let $u\in {{\mathbb {C}}}^N$ and define $v{:}{=}(H_0+\alpha )^{1/2}u$. Then,

$$\begin{aligned} \begin{aligned} \left| \left\langle u,Hu \right\rangle - \left\langle u,H_0u \right\rangle \right| = \left| \left\langle u,Wu \right\rangle \right|&= \left| \left\langle v, (H_0+\alpha )^{-1/2}W(H_0+\alpha )^{-1/2}v \right\rangle \right| \\&\leqslant \left\| W \right\| _{H_0,\alpha }\left\langle v,v \right\rangle = \left\| W \right\| _{H_0,\alpha }\left\langle u,(H_0+\alpha )u \right\rangle . \end{aligned}\nonumber \\ \end{aligned}$$

(60)

Therefore,

$$\begin{aligned}{} & {} (1 - \left\| W \right\| _{H_0,\alpha })\left\langle u,H_0u \right\rangle - \alpha \left\| W \right\| _{H_0,\alpha }\left\langle u,u \right\rangle \leqslant \left\langle u,Hu \right\rangle \leqslant (1 + \left\| W \right\| _{H_0,\alpha })\left\langle u,H_0u \right\rangle \nonumber \\{} & {} \quad + \alpha \left\| W \right\| _{H_0,\alpha }\left\langle u,u \right\rangle . \end{aligned}$$

(61)

The min–max theorem then yields for $i=1,\dots ,N$,

$$\begin{aligned} (1 - \left\| W \right\| _{H_0,\alpha })\nu _i(H_0) - \alpha \left\| W \right\| _{H_0,\alpha }\nonumber \\ \leqslant \nu _i(H) \leqslant (1 + \left\| W \right\| _{H_0,\alpha })\nu _i(H_0) + \alpha \left\| W \right\| _{H_0,\alpha }, \end{aligned}$$

(62)

which gives the desired inequality. $\square$

In our case, we can apply this result to

$$\begin{aligned} H_n^{\ell } + X_n^{\ell } = H_n + (H_n^{\ell } - H_n) + X_n^{\ell }, \end{aligned}$$

(63)

with $H_0 = H_n$, $W = W_n^{\ell } {:}{=}(H_n^{\ell } - H_n) + X_n^{\ell }$ and $\alpha =\varepsilon _{N+{N_{\textrm{ex}}}+1} -\varepsilon _n> 0$. Proposition 1 applied to the $(N+{N_{\textrm{ex}}}+1)$-th eigenvalues then yields

$$\begin{aligned} \varepsilon _{N+{N_{\textrm{ex}}}+1}^{\ell } -\varepsilon _n\geqslant \left( \varepsilon _{N+{N_{\textrm{ex}}}+1}-\varepsilon _n \right) \left( 1 - 2\left\| W_n^{\ell } \right\| _{H_n,\varepsilon _{N+{N_{\textrm{ex}}}+1}-\varepsilon _n} \right) \nonumber \\ \approx \left( \varepsilon _{N+{N_{\textrm{ex}}}+1}-\varepsilon _n \right) , \end{aligned}$$

(64)

where we assume that $2\left\| W_n^{\ell } \right\| _{H_n,\varepsilon _{N+{N_{\textrm{ex}}}+1}-\varepsilon _n}$ is small enough to be negligible with respect to 1, which is the case if the extra states are sufficiently converged. Now, if this bound is valid in theory, in practice we do not have access to $\varepsilon _{N+{N_{\textrm{ex}}}+1}$ as we work with $N+{N_{\textrm{ex}}}$ bands only. However, up to loosing sharpness, we can use that $\varepsilon _{N+{N_{\textrm{ex}}}+1}\geqslant \varepsilon _{N+{N_{\textrm{ex}}}}$ where $\varepsilon _{N+{N_{\textrm{ex}}}}$ can be estimated using the last extra band. Indeed, using, for instance, the Bauer–Fike bound ( [24, Theorem 1] or [42]), we obtain

$$\begin{aligned} \varepsilon _{N+{N_{\textrm{ex}}}} \geqslant \varepsilon _{N+{N_{\textrm{ex}}}}^{\ell } - \left\| r_{N+{N_{\textrm{ex}}}}^{\ell } \right\| , \end{aligned}$$

(65)

where $r_{N+{N_{\textrm{ex}}}}^{\ell }$ is the residual associated with the last extra band. Of course, this estimate is not sharp as we expect the error on the eigenvalue to behave as the square of the residual, but this requires to estimate the gap to the rest of the spectrum, see, for instance, the Kato–Temple bound [24, Theorem 2]. In the end, we have the following lower bound for $\varepsilon _{N+{N_{\textrm{ex}}}+1}^{\ell }-\varepsilon _n$:

$$\begin{aligned} \varepsilon _{N+{N_{\textrm{ex}}}+1}^{\ell } -\varepsilon _n\geqslant \left( \varepsilon _{N+{N_{\textrm{ex}}}}^{\ell }-\varepsilon _n - \left\| r_{N+{N_{\textrm{ex}}}}^{\ell } \right\| \right) \approx \varepsilon _{N+{N_{\textrm{ex}}}}^{\ell }-\varepsilon _n, \end{aligned}$$

(66)

where we assume again that $\left\| r_{N+{N_{\textrm{ex}}}}^{\ell } \right\|$ is small enough with respect to $\varepsilon _{N+{N_{\textrm{ex}}}}^{\ell }-\varepsilon _n$.

We can now derive an upper bound on $\kappa _{n}^{\ell }$, the condition number of (57). It is given by the ratio of its highest eigenvalue and $\varepsilon _{N+{N_{\textrm{ex}}}+1}^{\ell } -\varepsilon _n$. Since the Laplace operator is the higher-order term in the Kohn–Sham Hamiltonian, the highest eigenvalue is, as usually in plane-wave simulations, of order $E_{\textrm{cut}}$. With proper kinetic preconditioning, we can assume that its contribution to the condition number of the linear system is constant with respect to $E_{\textrm{cut}}$ and n so that, finally,

$$\begin{aligned} \kappa _n^{\ell } \lesssim \frac{C}{\varepsilon _{N+{N_{\textrm{ex}}}+1}^{\ell } -\varepsilon _n} \lesssim \frac{C}{\varepsilon _{N+{N_{\textrm{ex}}}}^{\ell }-\varepsilon _n}. \end{aligned}$$

(67)

Therefore, the condition number is bounded from above by $C/(\varepsilon _{N+{N_{\textrm{ex}}}}^{\ell }-\varepsilon _n)$ to first order. The number of CG iterations to solve the linear system (50) with a given accuracy is then proportional to the square root of the condition number of the matrix (57) (see [43]):

$$\begin{aligned} \sqrt{\kappa _n^{\ell }} \lesssim \sqrt{\frac{C}{{\varepsilon _{N+{N_{\textrm{ex}}}}^{\ell }-\varepsilon _n}}}. \end{aligned}$$

(68)

Note that this upper bound is valid provided that the extra bands are converged enough, not necessarily fully, and proper kinetic preconditioning is employed.

Estimate (68) leads, as expected, to the qualitative conclusion that the more extra bands we use, the higher the difference $\varepsilon _{N+{N_{\textrm{ex}}}}^{\ell }-\varepsilon _n$ and the faster the convergence. However, note that it is not possible to evaluate directly the convergence speed as the constant C is a priori unknown, in particular if we use preconditioners.

1.2 A.2 An adaptive strategy to choose the number of extra bands

The main bottleneck of (68) is the estimation of the constant C. However, one can reasonably assume that this constant does not depend too much on n, so that the ratio between the number of iterations to reach convergence between the last occupied band ($n=N$) and the first band $(n=1)$ can be estimated by

$$\begin{aligned} \xi _{{N_{\textrm{ex}}}}^{\ell } {:}{=}\sqrt{\frac{\varepsilon _{N+{N_{\textrm{ex}}}}^{\ell }-\varepsilon _1}{\varepsilon _{N+{N_{\textrm{ex}}}}^{\ell }-\varepsilon _N}}, \end{aligned}$$

(69)

This ratio can be of interest as (68) suggests that the Sternheimer solver converges the fastest for $n=1$ and the slowest for $n=N$.

Fig. 9

Comparison between the ratio $\xi _{{N_{\textrm{ex}}}}^{\ell }$ ($\times$) and the ratios of the experimental number of iterations between the first and last occupied bands, with ($\square$) and without ($\circ$) the Schur complement (50). On the x-axis is the index of the k-point. [Left] Al$_{40}$ [Middle] Fe$_{2}$MnAl spin up channels [Right] Fe$_{2}$MnAl spin down channels

We plot in Fig. 9 the upper bound $\xi _{{N_{\textrm{ex}}}}^{\ell }$ as well as the computed ratios between the number of iterations of the first and last bands for the systems we considered in Sect. 4. These plots show that $\xi _{{N_{\textrm{ex}}}}^{\ell }$ is indeed an upper bound of the actual ratio. This bound does not seem to be sharp however. This is due to the successive approximations we made to obtain this estimate. Plots in Fig. 9 also confirm that if, for every k-point, the ratio of the number of iterations between the first and last occupied bands is assumed to be an accurate indicator of the efficiency of the Sternheimer solver, then using the Schur complement (50) always makes this ratio smaller.

If we want the ratio of the number of iterations between the first and the last occupied bands to be lower than some target ratio $\xi _T$ (for instance, 3), Fig. 9 suggests that the computable ratio $\xi _{{N_{\textrm{ex}}}}^{\ell }$ can help in choosing the number of extra bands to reach this target ratio. We propose in Algorithm 1 an adaptive algorithm to select the number of extra bands as a post-processing step after termination of the SCF. The basic idea is that, given the initial output $(\Phi ,\Phi _{\textrm{ex}}^{\ell })$ with $\ell =0$ of an SCF calculation, one iterates $\Phi _{\textrm{ex}}^\ell \rightarrow \Phi _{\textrm{ex}}^{\ell +1}$ where $\Phi _{\textrm{ex}}^{\ell }$ gathers the extra bands. At each iteration $\ell$, we compute $\xi _{{N_{\textrm{ex}}}}^{\ell }$ and check if it is below the target ratio. If not, we compute more approximated eigenvectors, that we converge until the residual $\left\| r_{N+{N_{\textrm{ex}}}}^{\ell } \right\|$ is negligible with respect to $\varepsilon _{N+{N_{\textrm{ex}}}}^{\ell }-\varepsilon _N$, and so on. To generate such a residual, after adding a random extra band properly orthonormalized, we update the extra bands using a LOBPCG with tolerance

$$\begin{aligned} \textrm{tol} = (\varepsilon _{N+{N_{\textrm{ex}}}-1}^{\ell }-\varepsilon _N)/50. \end{aligned}$$

(70)

Note that we use $\varepsilon _{N+{N_{\textrm{ex}}}-1}^{\ell }$ instead of $\varepsilon _{N+{N_{\textrm{ex}}}}^{\ell }$: this is done for the sake of simplicity, instead of updating the tolerance on the fly with $\varepsilon _{N+{N_{\textrm{ex}}}}^{\ell }$ changing at each iteration of the LOBPCG.

1.3 A.3 Numerical tests

We test this strategy on the systems investigated in Sect. 4, with different values for the target ratio $\xi _T$ in order to see a noticeable improvement for each system. In practice, we suggest this ratio to be between 2 and 3.

We first start with the Al$_{40}$ system. Figure 9 [Left] suggests that the default choice of extra bands already gives satisfying results by reaching a ratio of approximately 2.5 for all k-points but the $\Gamma$-point (for which there is no real issue with the Sternheimer equation, according to Table 1). We thus run Algorithm 1 with a smaller target ratio $\xi _T=2.2$. We use as initial value for ${N_{\textrm{ex}}}$ the default value for each k-point. Results are plotted in Table 3 [Left] and suggests adding 15 extra bands. Running again the simulations from Sect. 4 with 72 fully converged bands and 18 additional, not fully converged, bands yields indeed an improvement in the convergence of the CG when solving the Sternheimer equation with the Schur complement method. Moreover, in Fig. 10, the ratio $\xi _{{N_{\textrm{ex}}}}^{\ell }$ indeed lies below the target ratio $\xi _T=2.2$ and matches this ratio for the k-points that caused difficulties for the Sternheimer equation solver to converge. In terms of computational time, the number of Hamiltonian applications to compute the response has been reduced from $\sim 14,800$ with the default number of extra bands to $\sim 12,800$. However, running the algorithm required $\sim 3,400$ additional Hamiltonian applications, making the total amount of Hamiltonian applications higher than that of the Schur approach with the standard heuristic.

Table 3 Suggested number of extra bands for Al$_{40}$ and Fe$_{2}$MnAl to reach the target ratio $\xi _T$, obtained with Algorithm 1 with default ${N_{\textrm{ex}}}$ as starting point, as well as the number of iterations to reach convergence with the newly suggested ${N_{\textrm{ex}}}$. Note that the ratio between iterations indeed lies below the target ratio $\xi _T$

Fig. 10

Comparison between the ratio $\xi _{{N_{\textrm{ex}}}}^{\ell }$ ($\times$) and the ratios of the measured number of iterations between the first and last occupied bands, with ($\square$) and without ($\circ$) the Schur complement (50). On the x-axis is the index of the k-point. [Left] Al$_{40}$ with 15 additional extra bands. [Left] Fe$_{2}$MnAl spin $\uparrow$ with 8 additional bands. [Right] Fe$_{2}$MnAl spin $\downarrow$ with 8 additional bands

Similarly, for Fe$_{2}$MnAl, we run Algorithm 1 with target ratio $\xi _T=2.5$ as well as initial value the default ${N_{\textrm{ex}}}$ for all the 140 k-points and spin polarizations. We present in Table 3 [Right] the output for both spin polarizations of two particular k-points. The results are similar for the rest of the k-points and the maximum additional extra bands suggested by the algorithm is 8. We thus run the same simulations as in Sect. 4 but this time with 35 fully converged bands and 11 extra, nonnecessarily converged, bands. We indeed see for these two k-points that the target ratio has been reached, and that the number of iterations to converge is smaller than for the default choice we made in Table 2. In Fig. 10, we plot the ratios $\xi _{{N_{\textrm{ex}}}}^{\ell }$ as well as the actual ratios and they almost all lie below the target ratio. Contrarily to Al$_{40}$, we note, however, that the actual measured ratios are not always below the indicator $\xi _{{N_{\textrm{ex}}}}^{\ell }$. In terms of computational time, the number of Hamiltonian applications has been reduced from $\sim 208,000$ with the default choice of ${N_{\textrm{ex}}}$ to $\sim 179,000$. Again, running the algorithm required $\sim 49,000$ additional Hamiltonian applications, making it more expensive than using the default number of extra bands.

It appears that Algorithm 1 can be used to choose the number of extra bands in order to reach a given ratio $\xi _T$. However, using the algorithm as such is not useful in practice as it requires a too high number of Hamiltonian applications, making this strategy less interesting than the Schur approach we proposed with the default choice of extra bands. Strategies to reduce the number of Hamiltonian applications in order to choose an appropriate number of extra bands will be subject of future work.