A force-based gradient descent method for $\mathit{\text{ab initio}}$ atomic structure relaxation

Force-based algorithms for $\mathit{\text{ab initio}}$ atomic structure relaxation, such as conjugate gradient methods, usually get stuck in the line minimization processes along search directions, where expensive $\mathit{\text{ab initio}}$ calculations are triggered frequently to test trial positions before locating the next iterate. We present a force-based gradient descent method, WANBB, that circumvents the deficiency. At each iteration, WANBB enters the line minimization process with a trial stepsize capturing the local curvature of the energy surface. The exit is controlled by an unrestrictive criterion that tends to accept early trials. These two ingredients streamline the line minimization process in WANBB. The numerical simulations on nearly 80 systems with good universality demonstrate the considerable compression of WANBB on the cost for the unaccepted trials compared with conjugate gradient methods. We also observe across the board significant and universal speedups as well as the superior robustness of WANBB over several widely used methods. The latter point is theoretically established. The implementation of WANBB is pretty simple, in that no a priori physical knowledge is required and only two parameters are present without tuning.

Atomic structure relaxation determines the ground-state atomic configuration by searching for the local minimum in the energy landscape. It makes the foundation and potential performance bottleneck for the search of global structure minimum [1][2][3] and high-throughput calculations in material design [4][5][6][7][8]. Conjugate gradient methods (CG) [9], the direct inversion in the iterative subspace (DIIS) [10], and limited-memory Broyden-Fletcher-Goldfarb-Shanno quasi-Newton methods (LBFGS) [11] are three widely used methods to relax the atomic structure following the density functional theory (DFT) calculation. While DIIS and LBFGS could be slightly faster in some cases, CG is more stable in the situations where the initial atomic configuration is far from equilibria. For large systems, hundreds of steps may be required to converge to the typical atomic force tolerance of 0.01 eV/Å. So it will be tremendously helpful even if the ab initio atomic relaxation could be sped up to some degree.
One reason for the slow convergence of the local minimization goes to the narrowly curved energy valley, which prevents the efficient execution of the conventional CG. Thus, various preconditioners have been designed for improving search directions via approximating the Hessian matrices. Zhao et al. [12] notice that the Hessian matrix of large quantum dots can be approximated by that of motifs comprising merely a few atoms. Some other works construct general force fields based on a priori physical knowledge. Surrogate potentials are used to either directly approximate the Hessian matrix [13][14][15], possibly with on-the-flight fittings [16,17], or generate the "preconditioned force" with the implicit inversion of the Hessian matrix [18].
Another reason is the inefficient execution of the line minimization (LM) to which not much attention has been paid in the previous study. Actually, gradient descent method and CG are two typical LM-based methods [9,19,20]. A flowchart of the general LM-based method for the ab initio atomic relaxation is depicted in Fig. 1. The ab initio LM process, marked out by the shaded box in Fig. 1, is carried out along search directions such as steepest descent or CG. During this process, one begins with an initial trial stepsize and then repeatedly invokes an LM algorithm to generate new trials until certain LM criterion is met. Involving moving atoms, each trial step is followed by ab initio calculation, e.g., solving the Kohn-Sham (KS) equation. But only the last accepted trial step triggers the update of search direction. Taking the conventional CG for instance, we test nearly 80 systems and find that each update of search direction is followed by 2.4 calls of LM algorithms on average. In other words, nearly 60% of the total computational overhead is paid for the unaccepted trial steps in the LM process. Moreover, preconditioning with the (fitted) classical force fields does not necessarily waive the cost in the LM process [17]. As a result, the computational cost on the LM process is far more than that on simply determining a search direction if the whole procedure gets frequently trapped in the inner loop.
In this work, we present an LM-based gradient descent method, called WANBB, which avoids trapping in the LM process during atomic structure relaxation. We achieve this by (i) calculating an initial trial stepsize that captures the local curvature of potential energy surface (PES) and (ii) devising an unrestrictive LM criterion that tends to accept early trials. Nearly 80 systems have been tested, including organic molecules, metallic systems, semiconductors, surface-molecule adsorption systems, ABX 3 perovskites, etc. The set of benchmark tests favors the universality of everyday ab initio atomic relaxation. The unaccepted trial steps in the LM process account for only about 1% of the total solved KS equations on average in WANBB, compared with nearly 60% in CG. The considerable compression of the LM process leads to a prominent saving on running time. The average speedup factors of WANBB over CG, DIIS, and LBFGS are about 1.5, 1.2, and 1.2, respectively. The robustness of the atomic relaxation method, i.e., Start Input initial atomic positions R0, force tolerance ε > 0; let k := 0 convergence to equilibria regardless of initial configurations, is also investigated. While CG, DIIS, and LBFGS fail on some systems, WANBB manages to converge across the benchmark. This robustness is theoretically established and helpful for the cases with local lattice distortion or reconstruction. Last but not least, WANBB is pretty simple. Free of a priori physical knowledge, it can work with only two parameters present but no tuning required.
In what follows, we are ready to deliver the algorithmic development. To facilitate narration, we collect some notations beforehand. We denote by E(R) ∈ R and F (R) ∈ R 3×N , respectively, the potential energy and atomic forces evaluated at R. When describing algorithms, we use subscripts for abbreviation; e.g., both E k and F k are evaluated at R k . We denote by I the identity matrix in R 3×N . The operator ·, · represents the Frobenius inner product of two 3-by-N matrices whereas · F yields the Frobenius norm as M F := M, M for any M ∈ R 3×N . We further denote by M 2,∞ the ℓ 2,∞ -norm of matrix M , namely, the maximum among the ℓ 2 -norms of columns in M . In particular, F 2,∞ gives the maximum ℓ 2 -norm among the atomic forces. Throughout iterations, we take the atomic forces F k as the search direction The whole procedure shown in Fig. 1 can fall into the trap of the LM process owing to two factors. One is the poorly chosen initial trial stepsize, which degrades the starting states of the LM process. The other is the restrictive LM criterion, which acts as a stringent stopping rule for the LM algorithm. The ideal situation is that the procedure always passes the test of the LM criterion with just the original initial stepsize. Intuitively, we require both good initial stepsizes and a loose but safe criterion to achieve this target. Here, "safe" means that the convergence property is maintained with the LM criterion.
At first sight, a good initial stepsize should grasp, as much as possible, the curvature of PES around the current configuration. As early as 1988, Barzilai and Borwein [21] propose the following renowned Barzilai-Borwein (BB) stepsizes: for any k ≥ 1, where S k−1 := R k − R k−1 and Y k−1 := F k−1 − F k represent the difference of successive atomic positions and forces, respectively. Remarkably, these two stepsizes satisfy In other words, α BB1 k −1 I and α BB2 k −1 I approximate the Hessian matrix around R k along the last search direction. They capture the local curvature of PES with merely two 3-by-N matrices. The alternating BB (ABB) stepsize [22] is presented as follows which is shown to overwhelm BB1 and BB2 in many contexts (e.g., [22,23]). It is not difficult to prove that the implementation of (A)BB stepsizes requires the same computational and storage complexity as computing the conjugate parameter in CG.
Despite the simplicity, the (A)BB stepsizes turn out to be significant benefits in solving various optimization problems. The most related application goes to solving KS equations as orthogonality constrained optimization problems, where the (A)BB stepsizes play as important local acceleration components; see, e.g., [23][24][25][26][27]. However, the application of the (A)BB stepsizes on the atomic structure relaxation problem remains to be investigated. In this work, we take the (A)BB stepsizes as the initial trial stepsizes to exploit their merits.
Without LM, the energy values produced by gradient descent methods with the (A)BB stepsizes do not necessarily decrease monotonically. This phenomenon reminds us of the possible divergence when starting the algorithms far from equilibria. Nevertheless, the monotone LM criterion, adopted by the conventional CG, can ruin the merits of the (A)BB stepsizes, since frequent calls of LM algorithms are entailed to assure the monotonic decrease of the energy and the local curvature of PES may be lost. Therefore, it is more sensible to design a nonmonotone LM (NLM) criterion, accepting the original (A)BB stepsizes in most cases.
The existing NLM criteria primarily fall into two categories: the "max" type (MNLM) [28] and the "average" type (ANLM) [29]. In the (k + 1)-th iteration, the MNLM criterion asks the LM algorithm to return a factor r k > 0 such that , both c and M are given positive constants, i.e., there is a sufficient reduction with respect to (w.r.t.) the largest one of the past m(k) energies. The ANLM criterion requires r k to fulfill where σ ∈ (0, 1), the monitoring sequence {C k } is updated through 1]. Roughly in the ANLM criterion, a sufficient reduction w.r.t. C k is obligatory.
It is straightforward to deduce that C k is a weighted average of the past energies {E j } k−1 j=0 and the current E k . The construction of the monitoring sequence is completely free of a priori physical knowledge. Clearly, entailing reduction w.r.t.Ē k and C k respectively instead of E k , both MNLM and ANLM criteria are less restrictive than the monotone criterion adopted by the conventional CG becauseĒ k , C k ≥ E k . Moreover, as noted in [29], the MNLM criterion is sensitive to the choice of M in some contexts, and is often outperformed by the ANLM criterion.
However, preliminary simulations show that the ANLM criterion is not loose enough for the (A)BB stepsizes during the atomic relaxation. We take the gradient descent direction armed with the ABB stepsizes and ANLM criterion to relax the silicene surface with one adsorbed H 2 O molecule. From the short relaxation history in Fig. 2, one can observe that the monitoring sequence {C k } (the blue dashed line with square markers) decreases fast, stays close to the energies {E k }, and does not tolerate somewhat large increments. As a result, the combination "ABB + ANLM" solves 3 KS equations for the unaccepted trials during the short history, which are marked out by the black arrows in Fig. 2. For example, the initial trial step at the 8-th iteration gives an energy value of -99.92 eV, which is higher than C 7 = -100.03 eV and thus is rejected by the ANLM criterion.
To this end, we devise a reweighted ANLM (WANLM) criterion, which differs from the ANLM criterion in the weights assigned to the past and current energies. Specifically, we denote the monitoring sequence in WANLM by {B k } and update it via with B 0 := E 0 , P 0 := 1, and µ k ∈ [0, 1]. Then WANLM needs r k to satisfy both Eqs. (2) and (3) with C k , Q k , and η k replaced by B k , P k , and µ k , respectively. a This value is recommended by [29].
The looseness of the WANLM criterion can be understood via the comparison with the ANLM criterion. Due to Eq. (4) and Since the monitoring sequence {C k } decreases monotonically (see the appendix for a rigorous proof), it would then stay close to the energy values {E k }, rendering the ANLM criterion stringent for the (A)BB stepsizes as reflected in Fig. 2. In contrast, from Eq. (5) and µ 0 ∈ [0, 1], the coefficient of E 1 is µ 0 /(1 + µ 0 ) ≤ 0.5. Therefore, the monitoring sequence {B k } in the WANLM criterion is relatively far away from the energies {E k } at the beginning and gradually approaches the latter one, leaving sufficient room to accept the (A)BB stepsizes. As shown in Fig. 2, the first 7 energy values obtained by "ABB + ANLM" and "ABB + WANLM" are the same. But one can find a much larger distance from {E k } to {B k } than {C k }, which enables the WANLM criterion to accept the 8-th initial trial step without invoking extra ab initio calculations and retain the merits of the ABB stepsizes. All the subsequent iterations of "ABB + WANLM" benefit from this characteristic and progress to lower energy values faster. Regarding safety, it can be verified that the sufficient reduction in terms of the atomic force norm holds for the monitoring sequence {B k }. One can then establish the convergence to equilibria,regardless of initial configurations, after noticing the lower boundedness of {B k }. For more details, please refer to the appendix.
The gradient descent method combining the (A)BB stepsizes with the WANLM criterion is referred to as WANBB. We set the initial trial stepsize to be 4.8 × 10 −2 Å 2 /eV at the first iteration, and take the alternating strategy described in Eq. (1) afterwards because it leads to better performance than BB1 and BB2 in our simulations. Before invoking an LM algorithm, we take absolute value in case that S k , Y k < 0 and make truncation for safeguard: As is shown in the appendix, the convergence of WANBB to equilibria can be attributed to both Eqs. (2) and (3). For the sake of implementation, it is more efficient to fullfill only Eq. (2). In WANLM, we set µ k ≡ 0.05, c = 10 −4 . The LM algorithm used by WANBB forms quadratic or cubic interpolations based on recent information of PES; for details see [19,Chapter 3].
We have implemented WANBB in the in-house planewave code CESSP [30][31][32][33]. The exchange-correlation energy is described by the generalized gradient approximation [34]. Electron-ion interactions are treated with the projector augmented-wave (PAW) potentials based on the opensource ABINIT Jollet-Torrent-Holzwarth dataset library in the PAW-XML format [35]. The total energies are calculated using the Monkhorst-Pack mesh [36] with the k-mesh spacing of 0.2 Å −1 . The typical plane-wave energy cutoffs are around 500~600 eV. The KS equations corresponding to the updated atomic positions are solved by the preconditioned selfconsistent field (SCF) iteration [32].
In addition to WANBB, we test the performances of three relaxation methods, including the conventional CG [9,19], LBFGS [11], and DIIS [10]. By the way, DIIS is the socalled quasi-Newton method in some software such as VASP [37,38], which can be derived from Broyden's method [39]. The implementations of CG and DIIS follow those in [11,19], and LBFGS is realized by the transition state tool VTST [40]. CG employs Brent's method [41] as the LM algorithm.
The benchmark for performance test contains nearly 80 systems from various categories, including organic molecules, metallic systems, semiconductors, surface-molecule adsorption systems, and ABX 3 perovskites. The number of atoms ranges from 2 to 429. Some of them are available online, e.g., in Materials Project [42] and Organic Materials Database [43]. In the beginning, these systems are at their ideal crystal positions for defects, heterostructures, and substitutional alloys, or simply place a molecule on top of a surface. These are the typical initial configurations used by researchers in their simulations. More information on the benchmark can be found in Supplementary Material [44]. The atomic structure relaxation is terminated if either the maximum atomic force falls below 0.01 eV/Å or the number of the solved KS equations (#KS) arrives at 1000. The convergence criterion for the SCF iteration is 10 −5 eV.
We evaluate the performances of the relaxation methods in two ways. One is the number of the solved KS equations (#KS) needed to reach the 0.01 eV/Å force tolerance. The other is the running time (CPU) of the atomic structure relaxation. We adopt Dolan and Moré's performance profile [45] for an overall comparison among CG, WANBB, DIIS, and LBFGS. To this end, for any method indicated by m ∈ {CG, WANBB, DIIS, LBFGS}, we denote its #KS and CPU for relaxing the n-th system by #KS n,m , CPU n,m , respectively [46], define ratios The quantities π KS m (ω) and π CPU m (ω) estimate respectively the probabilities of r KS n,m ≤ ω and r CPU n,m ≤ ω over the benchmark [47]. For example, π CPU WANBB (1) yields the portion of the systems where WANBB is the fastest. And π CPU WANBB (2) means the estimated probability that, for one system, the running time of WANBB is not more than twice the shortest running time needed by CG, WANBB, DIIS, and LBFGS. In terms of the #KS, π KS m (ω) can be interpreted in a similar way. We plot π KS m (ω) and π CPU m (ω) as a function of ω for each method in Fig. 3. Basically, the larger the area under the curve, the better the overall performance of the method. It is easy to see that the overall performance of WANBB ranks the best. More specifically, one can tell from the supplementary Tables I and II [44] that the average CPU speedup factors of WANBB over CG, DIIS, and LBFGS are 1.51, 1.21 and 1.16, respectively. As shown in Fig. 3, though sometimes requiring solving more KS equations than LBFGS, WANBB turns out to be faster than LBFGS in most tested systems. LBFGS may make larger atomic displacements due to constant stepsize, which results in underestimating the quality of initial wave functions and slowing down the SCF iteration. Later, we will explain it in more detail using the relaxation of the PuO 2 surface with one adsorbed H 2 O molecule (see Fig. 9).
In addition to the efficiency, the robustness is also worth the whistle. As shown in the supplementary Tables I and II [44], WANBB converges in all the benchmark systems, while CG fails in 1 case due to the breakdown of Brent's method, DIIS and LBFGS fail in 8 cases and 1 case, respectively, due to divergence. This kind of robustness is theoretically guaranteed by the convergence property of WANBB proved in the appendix.
Taking the performance of CG as the baseline, we illustrate the acceleration obtained by WANBB for each system, as depicted in Fig. 4. We have an average CPU speedup factor around 1.51, and a factor above 2.0 on about 20% of systems. According to the supplementary Tables I and II [44], CG takes on average 117.68 calls of solving the KS equations to reach the force tolerance. Thus a speedup factor of 1.5~2 can mean a significant saving of computational overhead. This saving is mainly attributed to the enhanced ab initio LM process of WANBB since the distribution of #KS speedup is consistent with that of CPU speedup in Fig. 4. Moreover, we also count the percentage of the LM steps invoked by WANBB and CG for unaccepted trials and draw the frequency distribution over the benchmark in Fig. 5. The average LM step percentages for unaccepted trials are 59.09% and 1.47% for CG and WANBB, respectively. It is obvious from these statistics that WANBB escapes from the LM process quickly, usually without invoking LM algorithms at all. The following discussion is concentrated on the performance comparison between CG and WANBB for each cate-gory. In Fig. 6, we plot the mean speedup factors by category and the number of atoms, respectively. All description below is supported by Tables I and II in Supplementary Material [44]. Whenever speaking of "speedup" below, we are talking about the CPU speedup if not specified. Organic molecules. This class of systems includes various alkanes and amino acids. The amino acids contain common organic atomic types, e.g., oxygen, nitrogen, sulfur, etc. The number of atoms in this class ranges from 26 to 106. Most of the initial configurations are available in Organic Materials Database [43]. Compared with CG, WANBB achieves a mean speedup factor of 1.81 for this class and speedup factors greater than 2 on over 20% of systems. In the case of Glutamate, the acceleration factor even reaches 3.18, which means a saving of nearly 70% computational overhead.
Metallic systems. We consider Ag 16 cluster, Cu(111) surface with adsorbed organic molecules, and some alloys including Ni 3 Al based intermetallic compound as well as highentropy alloys (HEAs) such as FeCrNiCoAl. The speedup factors on these systems w.r.t. CG are mostly around 1.5~3. In particular, for Cu(111) surface with adsorbed organic molecules, CG needs to solve over 400 KS equations for atomic relaxations, while WANBB cuts the required #KS roughly in half. It is also worth mentioning that WANBB can efficiently treat HEAs with potential strong local lattice distortion [48]. Actually, local lattice distortion is an essential factor in studying HEAs [49], and can be described by the supercell method with ab initio atomic structure relaxation. In Fig. 7, the local lattice distortion obtained by relaxing the body-centered cubic supercell of FeCrNiCoAl is illustrated by the smearing radial distribution function (RDF), in sharp contrast to the peaks at coordination shells in the case of the ideal lattice. On average in our benchmark, WANBB achieves a speedup of 1.54 over CG for HEAs.
Semiconductors & heterojunction systems. We test semiconductors (CdSe, GaAs, Si, etc.) and heterojunction systems (GaAs-InAs). Moreover, we consider semiconductors with defects by constructing the bulk GaAs and Si supercells with some vacancies. In most of these systems, the acceleration factors are above 1.5, and sometimes exceed 3. We also investigate the performances of WANBB compared with CG on Si supercells with increasing system sizes; see Fig. 8. As the system size increases, WANBB is more stable in terms of #KS, and its CPU advantage becomes more prominent. ABX 3 perovskite systems. This type of systems is fairly popular in recent solar cell studies. We consider the cases where "A" = CH 3 NH 3 , CH(NH 2 ) 2 , or Cs, "B" = Pb or Sn, "X" = Br, Cl, or I. Their initial configurations come from Materials Project [42]. Due to the existence of the organic molecules and their relative weak interaction with the inorganic framework, the atomic relaxations of such systems are quite challenging. Taking CH 3 NH 3 SnI 3 for instance, CG solves 581 KS equations, while WANBB entails only about one fifth of the expenditure.
Surface-molecule adsorption systems. Several surfacemolecule adsorption systems have been tested, including a 2D silicene surface with one adsorbed H 2 O molecule (total 21 atoms), a 2D MoS 2 surface with one adsorbed NH 3 molecule (total 16 atoms), 2D PbS(001) and (111) surfaces with one adsorbed oleic acid molecule ((001) total 254 atoms, (111) total 429 atoms), and a 2D PuO 2 surface with one adsorbed H 2 O molecule (total 195 atoms). The results show that CG often undergoes quite a long way until convergence in this category,and even fails in one case due to the breakdown of Brent's method. On the contrary, WANBB converges in all cases and usually costs less time; it is even twice as fast as CG on PbS(001) surface with one adsorbed oleic acid molecule. One interesting observation is that WANBB occasionally takes equal or a few more #KS than CG, even though it is less time-costly. This is also observed when comparing WANBB and LBFGS.
To explain, we plot the cumulative SCF iteration numbers of CG, WANBB, and LBFGS when relaxing the PuO 2 surface with one adsorbed H 2 O molecule; see observe that CG and LBFGS often consume more SCF iterations per KS equation than WANBB. In fact, CG and LBFGS make larger atomic displacements, resulting in degrading initial wave functions and slowing down the SCF iteration. On the other hand, Fig. 9 also reveals that CG and LBFGS produce more efficient search directions, yielding fewer iteration numbers. This weakness, however, does not prevent WANBB from being a more favorable alternative in the context of ab initio simulations, because the cut down on the LM process obviously brings much more benefits.
2D and other systems. The 2D materials in our benchmark include C 96 , Si 32 , and Germanene on hexagonal boron-nitride. The average speedup factors of WANBB over CG on these systems are 1.26. We also include systems of other type, e.g., cerium-iron oxide and potassium chloride. The former one is of practical usage for NO x reduction. WANBB reduces the numbers of solved KS equations by half on both of them.
We recognize that the inefficiency of the LM process is an important obstacle to the high-performance execution of ab initio atomic structure relaxation. With the most widely used CG, nearly 60% of the computational overhead on average is attributed to the unaccepted trials in the LM process. To this end, we propose a force-based gradient descent method called WANBB, where the initial trial (A)BB stepsizes grasp well the local curvature of PES and the unrestrictive WANLM criterion often accepts the initial trials without calling extra ab initio calculations. The robustness of WANBB, i.e., the convergence to equilibria regardless of initial configurations, is theoretically guaranteed. The numerical simulations on the benchmark containing nearly 80 various systems reveal the average 1.47% computational cost for unaccepted trials in WANBB and demonstrate its universal speedups as well as superior robustness over CG, DIIS, and LBFGS.
Our work concentrates on the LM process. Nevertheless, efficient search directions together with subtly designed LM can be another great boon. This topic merits futural study. In this part, we analyze the convergence property of WANBB. Before moving on, we define the level set For the ease of reference, we restate the nonmonotone conditions in WANLM as follows: in the (k + 1)-th iteration, R k+1 satisfies By the arguments in [ Proof. We prove by contradiction. Suppose there exists l 1 > 0 such that E l1 > E 0 . By the condition (A.1), Since B l1−1 is a convex combination of {E 0 , . . . , E l1−1 }, there must exist some l 2 ≤ l 1 − 1 for which E l1 ≤ E l2 holds. We can repeat the above arguments to obtain a nonnegative sequence {l i } which is finite for its strict monotonicity and ends up with 0. Therefore, we get E 0 < E l1 ≤ E l2 ≤ · · · ≤ E 0 , a contradiction.
From the update formula (5), one could obtain that the surrogate sequence {B k } is always lower bounded by the potential energy sequence {E k }.
Proof. The claim is true for k = 0 by definition. For k ≥ 1, by the update formula (5) of B k+1 , one has where the inequality follows from the condition (A.1). The proof is complete.
With Lemmas 1 and 2 in place, we are ready to show the convergence of WANBB. The proof relies on the conditions (A.1) and (A.2) as well as the update formula (5). Specifically, the condition (A.2) ensures a uniform positive lower bound of stepsizes, while Eqs. (A.1) and (5) guide us to a sufficient reduction over {B k }. Theorem 1. Suppose the potential energy is bounded from below and F is Lipschitz continuous, with modulus L > 0, on L, namely, Let {R k } be the atomic position sequence generated by WANBB with µ k ∈ [µ min , µ max ] ⊆ (0, 1]. Then we have That is to say, the atoms will approach equilibria, regardless of initial configurations. Proof. Substracting F k 2 F from both side of Eq. (A.2), we have Combining Lemma 1, the Lipschitz continuity of F in L, and the fact σ ∈ (0, 1), we can further derive A byproduct of Eq. (A.5) is the monotonicity of {B k }. Since E is bounded from below on L and B k ≥ E k for any k ≥ 0 (by Lemma 2), {B k } is lower bounded as well. Consequently, by Eq. (A.5), Following from the update formula of {Q k } and the fact µ max ≥ µ k ≥ µ min > 0, Q k ≥ 1, we can derive which, together with Eq. (A.6), yields the summability of F k . Hence, F k → 0 as k → ∞ as desired.