Resolving all-order method convergence problems for atomic physics applications
aa r X i v : . [ phy s i c s . a t o m - ph ] A p r Resolving all-order method convergence problems for atomic physics applications
H. Gharibnejad, E. Eliav, M. S. Safronova, and A. Derevianko Department of Physics, University of Nevada, Reno, Nevada 89557, USA. Department of Chemistry, Tel Aviv University, Tel Aviv, Israel. Department of Physics, University of Delaware, Newark, Delaware 19716, USA.
The development of the relativistic all-order method where all single, double, and partial tripleexcitations of the Dirac-Hartree-Fock wave function are included to all orders of perturbation theoryled to many important results for study of fundamental symmetries, development of atomic clocks,ultracold atom physics, and others, as well as provided recommended values of many atomic prop-erties critically evaluated for their accuracy for large number of monovalent systems. This approachrequires iterative solutions of the linearized coupled-cluster equations leading to convergence issuesin some cases where correlation corrections are particularly large or lead to an oscillating pattern.Moreover, these issues also lead to similar problems in the CI+all-order method for many-particlesystems. In this work, we have resolved most of the known convergence problems by applyingtwo different convergence stabilizer methods, reduced linear equation (RLE) and direct inversion ofiterative subspace (DIIS). Examples are presented for B, Al, Zn + , and Yb + . Solving these conver-gence problems greatly expands the number of atomic species that can be treated with the all-ordermethods and is anticipated to facilitate many interesting future applications. PACS numbers: 31.15.bw, 31.15.ac, 06.30.Ft, 31.15.ag
I. INTRODUCTION
The coupled-cluster (CC) method has been success-fully applied to solve quantum many-body problems inquantum chemistry [1, 2] as well as computationalatomic [3] and nuclear physics [4]. A relativistic lin-earized variant of the coupled-cluster method (which isnumerically symmetric and is generally referred to as“all-order method”) was developed for atomic physicsapplications in Refs. [5–7]. It is one of the most accu-rate methods currently being used in the atomic struc-ture calculations. This approach was extremely success-ful and led to accurate predictions of energies, transi-tion amplitudes, hyperfine constants, polarizabilities, C and C coefficients, isotope shifts, and other propertiesof monovalent atoms, as well as the calculation of parity-nonconserving (PNC) amplitudes in Cs, Fr, and Ra + (see[8–13] and references therein). Further development ofthe all-order approach, that included triple excitationsand non-linear terms yielded the most precise evaluationof the PNC amplitude in Cs [14, 15] and consequent re-analysis of Cs experiment [16]. This work provided themost accurate low-energy test of the electroweak sector ofthe Standard Model to date, placed constraints on a va-riety of new physics scenarios beyond the SM, and, whencombined with the results of high-energy collider experi-ments, confirmed the energy dependence (or “running”)of the electroweak force over an energy range spanningfour orders of magnitude (from ∼
10 MeV to ∼
100 GeV).All-order method was also used for development of ultra-precise atomic clocks [17–21], ultracold atom and quan-tum information studies [22–26] and many other appli-cations. We refer the reader to review [12] for detailsof the all-order method and its applications. The all-order method is also used as a part of the CI+all-orderapproach for study of more complicated systems [27]. The all-order method requires iterative solutions of thelinearized coupled-cluster equations leading to conver-gence issues in some cases when correlation correctionsare very large or produce an oscillating iterative pattern.The initial guess of the solution is based on the low-orderperturbation theory. Therefore, if high-order correlationcorrections are large, initial guess is very poor leadingto very slow convergence or failure of the straightforwarditerative scheme. In addition, initial non-linear CC equa-tions may have more than one solution, so a convergenceto non-physical solutions may occur. Several such prob-lems have been identified over the years and led to failureto apply all-order approach for many important applica-tions. For example, all or almost all of the low-lying nd and nf states of B, Al, Zn + , Cd + , Hg + , and Yb + do notconverge if standard Jacobi-type iterative procedure isapplied. In the case of Yb + , even core equations do notconverge. Convergence problems also cause complete fail-ure of the all-order approach for super-heavy elements,such as element 113 (eka-Tl). All these convergence is-sues in monovalent systems lead to the same problems inthe application of the CI+all-order approach [27] to thecorresponding divalent systems, such as Al + , Hg, Yb,etc. since this method required prior solution of LCCSDequations for one-particle orbitals. There are several in-teresting present applications of these atoms and ionsthat require high-precision calculations possible with all-order techniques. For example, several of these systemsare used or proposed for optical clocks [17, 28–31] requir-ing precise knowledge of the blackbody radiation (BBR)shift which is hard to accurately measure. BBR shift isa leading source of uncertainties for many of the atomicclock schemes. Yb is used for an ongoing PNC experi-ment [32] as well as studies of degenerate quantum gases[33, 34] owing to a number of available isotopes. Thebest available Yb PNC amplitude value is only accurateto 20% [35].The convergence issues that arise in the solutions ofeigenvalue equations have a long history in general quan-tum chemistry and several methods have been developedto address them [36–42]. Most of these methods are basedon the fundamental idea of effective reduction of originallarge functional space and solution of the projected to thereduced (Krylov) subspace of the simplified equations.This idea was implemented for the first time in a quan-tum chemical application by Lanczos [43], who facilitateda partial diagonalization of a large matrix by transform-ing to a much smaller Krylov subspace, followed by amatrix triangularization procedure. In the present work,we consider two such convergence techniques, namely, re-duced linear equation (RLE) [36, 37] and the direct in-version of iterative subspace (DIIS) [38, 39]. Both ofthe methods use approximate solutions obtained fromfew last iterations as Krylov reduced functional subspaceonto which the linearized equations are projected and inwhich the projected system of equations is solved. Theconvergence of the methods is based on the constructionof error vectors. Different choices of the error vectorslead to different implementation of the methods. Amongthe most popular error vectors are: 1) the difference ofsubsequent iterations and 2) “true” error vector (e.g. dif-ference between exact solution and it’s approximation).In our work, both the convergence methods use the samebest least squares approximation to the true error vectorand thus are rather relative. Moreover, our variant ofDIIS can be regarded as a “symmetric” version of RLE(see below). However, while DIIS method is chosen tominimize the error vector in the least-squares sense, theRLE differs from it by requiring that this vector withinthe basis vanishes. We formulate here implementationsof the RLE and DIIS methods for our variant of thecoupled-cluster equations and test these stabilizer meth-ods on several specific examples, in which we were ableto resolve the convergence problems listed above. Wealso studied the effectiveness of these two techniques insolving specific types of the convergence problems as wellas accelerating convergence in all other cases. Accelera-tion of convergence is particulary important for furtherCI+all-order use since it requires solving all-order equa-tions for a large number of one-particle orbitals.Below, we briefly outline the essence of the convergencestabilization procedures. In the coupled-cluster method,the desired exact wave function | ψ i is obtained by ap-plying (a yet unknown operator) exp( T ) on some refer-ence wave function | φ i , for example, the Dirac-Hartree-Fock (DHF) wavefunction. For a closed-shell systemwith N electrons, the cluster operator T = P T p (where p = 1 , , ..., N ) has the form: T p = 1 /p ! X mn..ab.. ρ mn...ab... a † m a † n ...a a a b ... (1)Here, orbitals m, n... are single-particle excited states; a, b... are core states which are occupied in | φ i ; ρ ’s arecluster amplitudes (also called excitation coefficients) FIG. 1: (Color online) The failure of the LCCSD straight-forward iteration procedure for the 3 s state in boron. Thecalculated correlation energy is plotted as a function of theiteration number. The dashed (red) line indicates the valueof the experimental correlation energy. and a † and a are creation and annihilation operatorswith respect to the quasi-vacuum state | φ i . Finally, p is the number of core electrons excited when applying T p to | φ i . In the linearized coupled-cluster single-double(LCCSD) method, only T and T are retained and non-linear terms in the expansion of exp( T ) are truncated.The LCCSD equations are conventionally solved by an it-erative scheme, symbolically written as ρ ( n +1) = F ( ρ ( n ) ), F being specified later in Section II. In this paper, thistype of straightforward iteration procedure is referred toas the conventional iterations scheme (CIS).Both RLE and DIIS convergence stabilization proce-dures form ρ ( n +1) solution as the linear combination ofcluster amplitudes ( ρ ( n ) , ρ ( n − , ..., ρ ( n − l ) ) accumulatedfrom l previous CIS iterations. Further details of theLCCSD method and RLE and DIIS schemes are dis-cussed in Sections II and III.An example of failed conventional iteration procedureis shown in Fig. 1, where we plot the LCCSD correla-tion energy, δE , as a function of a number of valenceLCCSD iterations for the 3 s state of boron. The experi-mental correlation energy ( − . s correlation energy diverges from the experi-mental values dramatically and begins to oscillate aftera number of iterations. The convergence criteria is setto terminate the iteration procedure when the relativedifference between two consecutive iterations is reducedbelow 0.00001. The convergence is not reached even after500 iterations. As demonstrated below, this problem iscompletely resolved by the use of either RLE or DIIS pro-cedures and convergence to the above criteria is reachedwithin 30 iterations.This paper is organized as follows: in Section II, wedescribe the LCCSD method and the conventional iter-ation procedure (CIS) of solving the LCCSD equations.In Section III, we formulate RLE and DIIS schemes forLCCSD equations. In Section IV, we analyze perfor-mance of the RLE and DIIS schemes for various cases.Finally, in Section V, we draw the conclusions. II. LINEARIZED SD COUPLED-CLUSTERMETHOD
In the present implementation of the CC method, theexact valence wave function | ψ v i is obtained from thelowest-order DHF state, | φ v i = a † v | c i , (2)by applying a wave operator Ω = N [exp( T )] [3]: | ψ v i = Ω | φ v i , (3)where | c i is the core DHF state and N [ ... ] designates thenormal product of operators with respect to a closed-shellcore. Taking into account only the T and T terms inEq. (1), and truncating Ω past the linear terms in theexpansion of the exponential leads to the LCCSD ansatzfor the wave operatorΩ ≃ X ma ρ ma a † m a a + 12 X mnab ρ mnab a † m a † n a b a a + X m = v ρ mv a † m a v + X mna ρ mnva a † m a † n a a a v =1 + S c + D c + S v + D v . (4)Here, S c and D c ( S v , D v ) are the core (valence) singleand double terms, respectively.To find the cluster amplitudes (or excitation coeffi-cients) ρ , we need to specify the Hamiltonian. In ourapproach, we use the Hamiltonian [7] H = H + G : H = X i ε i N [ a † i a i ] + 12 X ijkl g ijkl N [ a † i a † j a l a k ] , (5)where H is the one-electron lowest-order DHF Hamilto-nian and G is the residual Coulomb interaction. Indices i , j , k , and l range over all possible single-particle orbitals,and g ijkl are the two-body Coulomb matrix elements. Aset of coupled equations for the cluster operators ( T ) n :( T c ) = S c , ( T v ) = S v , ( T c ) = D c , and ( T v ) = D v may be found from the Bloch equation [3]. For monova-lent systems [44]:( ε v − H )( T c ) n = { QG Ω } connected, n , (6)( ε v + δE v − H )( T v ) n = { QG Ω } connected, n , (7) where δE v = h φ v | G Ω | φ v i is the valence correlation en-ergy and Q = 1 − | φ v ih φ v | is the projection operator.Note that Eq. (6) contains only the core cluster opera-tors, while Eq. (7) contains both core and valence clusteroperators. The core equations (6) are solved first, andthe resulting CC core amplitudes are subsequently frozenand used in the valence equations (7).The summations over the magnetic quantum numbers m in Eqs. (6) and (7) are performed analytically. Afterthe angular reduction, the equation for the reduced singlecore cluster amplitudes ρ ( ma ) takes form [9, 45]:( ε a − ε m ) ρ ( ma ) = (8) δ κ m κ a { X nb δ κ n κ b s [ j b ][ j a ] Z ( mban ) ρ ( nb ) − X k X ncb ( − j a + j b + j c + j n [ j a ][ k ] Z k ( cbna ) ρ k ( nmcb )+ X k X rnb ( − j a + j b + j r + j n [ j a ][ k ] Z k ( mbrn ) ρ k ( rnab ) } . Here, [ j ] = 2 j + 1, κ is the relativistic angular momen-tum quantum number, ρ ( ma ) and ρ k ( mnab ) are reducedsingle and double cluster amplitudes, X k ( mnab ) are re-duced two-body Coulomb matrix elements, and Z k ( mnab ) = X k ( mnab )+ X k ′ [ k ] (cid:18) j m j a kj n j b k ′ (cid:19) X k ′ ( mnba ) . The equations for the reduced double core cluster am-plitudes ρ k ( mnab ) are given by:( ε ab − ε mn ) ρ k ( mnab ) = X k ( mnab ) (9)+ X cd X l,k ′ A X l ( cdab ) ρ k ′ ( mncd )+ X rs X l,k ′ A X l ( mnrs ) ρ k ′ ( rsab )+ "X r X k ( mnrb ) ρ ( ra ) δ κ r κ a + X c X k ( cnab ) ρ ( mc ) δ κ m κ c − X rc ( − j c + j r + k [ k ] Z k ( cnrb ) e ρ k ( mrac ) + (cid:20) a ↔ bm ↔ n (cid:21) , where A i are angular coefficients given in [45] and ε ij = ε i + ε j . The valence equations have exactly the same formas the core equations with the replacement of index a bythe valence index v everywhere and an addition of thevalence correlation energy δE v into the energy differenceon the left-hand side, i.e., ( ε a − ε m ) ρ ( ma ) −→ ( ε v − ε m + δE v ) ρ ( mv ).Implementation of the RLE and DIIS procedures re-quires rewriting the equations for the cluster amplitudesin a specific vector form. We introduce the vector nota-tion: t = (cid:18) ρ ( ma ) ρ k ( mnab ) (cid:19) , where ρ ( ma ) and ρ k ( mnab ) are to be understood ascolumns composed of all amplitudes for single and doubleexcitations, respectively, i.e., for all possible values of m , n , a , b , and k indexes. Then, the core equations given byEqs. (8) and (9) may be combined as D · t = − a − ∆ · t , (10)where a = − (cid:18) X k ( mnab ) (cid:19) , D = (cid:18) ε a − ε m ε ab − ε mn (cid:19) , and ∆ · t includes all terms on the right-hand sides of theEqs. (8) and (9) except for X k ( mnab ), which is includedin a .Valence equations may be written in the same way with t = (cid:18) ρ ( mv ) ρ k ( mnvb ) (cid:19) and a = − (cid:18) X k ( mnvb ) (cid:19) , D = (cid:18) ε v − ε m + δE v ε vb − ε mn + δE v (cid:19) . The main difference between the core and valence equa-tions for the implementation of the RLE and DIIS is thedependence of the valence array D on the iteration num-ber, since δE v is recalculated after every iteration. In thecore case, D remains constant.Solving Eq. (10) for t gives t = − D − ( a + ∆ · t ) . (11)The above equation can be solved iteratively as t ( m +1) = − D − ( a + ∆ · t ( m ) ) . (12)The iteration usually starts by inserting t (0) = 0 on theright hand side of Eq. (12) and finding t (1) . As we demon-strated in Fig. 1, convergence of this straightforward iter-ative scheme is occasionally very slow or fails altogether.The convergence methods that we develop in the nextsection will alleviate such problems and lead to fasterconvergence rates. III. RLE AND DIIS
In this section, we formulate implementation of RLEand DIIS methods for the LCCSD equations (11) dis-cussed in the previous section. Both methods are two-step procedures. In the first step, a few iterative solutions t ( i ) of Eq. (12) are found (same as the CIS). In the secondstep, a linear combination of these t ( i ) is used to find the next best solution of Eq. (12). The new answer is thenused for another initialization of the CIS and the twosteps are repeated until convergence to specified criteriais reached. In this section, we present the general RLEand DIIS formulas and derive their explicit form for theLCCSD equations.After accumulating m + 1 iteratively found solu-tions, t (1) , t (2) , ..., t ( m +1) , next best approximation canbe found as their linear combination, t [ m +1] = m X i =1 σ i t ( i ) = σ · T . (13)The quantities σ i are the weights that have to be de-termined by solving a system of equations constructedfrom previously found m + 1 CIS solutions. We note that t ( m +1) is not included in the linear combination (13),but is used to find σ i coefficients. Therefore, we use thenotation t [ m +1] instead of t ( m +1) to distinguish betweenthe m + 1 th solution found through the use of RLE/DIISmethods and the initial CIS result, respectively.Both direct inversion of iterative space (DIIS) andreduced linear equation (RLE) methods seek to mini-mize the error between the iteratively found solutionsof Eq. (11) and the exact answer. The error minimiza-tion is the basis for finding the appropriate σ i to formthe approximate solution t [ m +1] . Both methods also usea least square approach to the error minimization. Sincethe exact answer is unknown, approximations are used inthe minimization process. The approximate solution, asmentioned before, is constructed as a linear combinationof a series of iteratively found solutions. The differencebetween the DIIS and the RLE methods is in the assump-tions they make in order to minimize the errors. Furtherdetails of the difference between the two methods andderivations of the DIIS/RLE formulas can be found inthe Appendix A.We rewrite Eq. (10) as a + (∆ + D ) t = a + Bt = 0 . (14)The DIIS formula for determining σ i is given by Eq. (A7): T T B T a + T T B T BT σ = 0 . (15)The RLE formula for determining σ i is given by Eq. (A9): T T ( a + BT σ ) = 0 . (16)Both Eqs. (15) and (16) can be written as a system of m equations: α + R σ = 0 . (17)Solving the above system of equations for σ can be easilydone with standard linear algebra methods. The result-ing coefficients σ i are substituted into Eq. (13) to obtainbest new approximate solution t [ m +1] .Next, we write R and α of Eq. (17) in their explicitforms for DIIS and RLE methods. Substituting ∆ − D for B into DIIS equation (15) yields for the σ i t T ( i ) (∆ + D ) T a + t T ( i ) (∆ + D ) T (∆ + D ) t ( i ) σ i = 0 . Using Eq. (12), we find that ∆ · t ( i ) = − ( D · t ( i +1) + a ). Replacing the dot products involving ∆ with onesinvolving D yields explicit form of DIIS matrix for coreorbitals R ij = X k D kk a k (cid:16) t ( i +1) k + t ( j +1) k − t ( i ) k − t ( j ) k (cid:17) + X k ( a k ) + X k D kk (cid:16) t ( i ) k t ( j ) k + t ( i +1) k t ( j +1) k − t ( i +1) k t ( j ) k − t ( i ) k t ( j +1) k (cid:17) α i = X k a k D kk (cid:16) t ( i ) k − t ( i +1) k (cid:17) − X k ( a k ) . (18)The RLE equations for core orbitals are obtained byrepeating the same steps as for the DIIS approach butstarting from Eq.(16). The resulting RLE equations for R and α are R ij = X kl t ( i ) k (∆ kl + D kl ) t ( j ) l = X k [ t ( i ) k D kk t ( j ) k − t ( i ) k D kk t ( j +1) k − a k t ( j ) k ] ,α i = X k t ( i ) k a k . (19)We noted in the previous section that D depends onthe correlation energy, δE v , in the case of the valenceequations leading to the dependence of D on the iterationnumber. Therefore, the substitution D → D ( i ) mustbe made to rewrite the DIIS and RLE equations abovefor the valence orbitals. To derive the final form of theequations, we have to introduce a somewhat arbitrary dotproduct and normalization definitions. The explicit formof the core RLE equations is obtained by substituting theexpressions for D , a , and t from the previous section intoEq. (19): R ij = X ma ( ε a − ε m ) ρ ( i ) ( ma ) [ ρ ( j ) ( ma ) − ρ ( j +1) ( ma )]+ X L X mnab L ] ( ε ab − ε mn ) ρ ( i ) L ( mnab ) × [ ρ ( j ) L ( mnab ) − ρ ( j +1) L ( mnab )] − α i ,α i = − X L X mnab L ] X L ( mnab ) ρ ( i ) L ( mnab ) . (20) (cid:1)(cid:0)(cid:2)(cid:3)(cid:4)(cid:5) (cid:6)(cid:7) (cid:8)(cid:9)(cid:10)(cid:11)(cid:12)(cid:13)(cid:14)(cid:15)(cid:16)(cid:17)(cid:18) (cid:19) (cid:20) (cid:21) (cid:22)(cid:23) (cid:24)(cid:25) (cid:26)(cid:27) (cid:28)(cid:29) (cid:30)(cid:31) !" FIG. 2: (Color online) Comparison of the RLE5, DIIS5, andDIIS7 schemes for the 3 s state of boron. The correlationenergy is given in a.u. RLE equations for valence case are given by R ij = X ma (cid:16) ε a − ε m + δE ( j ) v (cid:17) ρ ( i ) ( ma ) × [ ρ ( j ) ( ma ) − ρ ( j +1) ( ma )]+ X L X mnab L ] 1[ j v ] (cid:16) ε ab − ε mn + δE ( j ) v (cid:17) ρ ( i ) L ( mnab ) × [ ρ ( j ) L ( mnab ) − ρ ( j +1) L ( mnab )] − α i ,α i = − X L X mnab L ] X L ( mnab ) ρ ( i ) L ( mnab ) , (21)where [ L ] = 2 L + 1.The implementation of the RLE and DIIS methodsproceeds as follows. In step one, our code makes a limitednumber, m + 1, of LCCSD iterations using the CIS. Thisis done to find m +1 series of single and double cluster am-plitudes, ρ ( i ) ( ma ) and ρ ( i ) L ( mnab ) that are then saved. Instep two, a separate subroutine applies the DIIS or RLEequations to these stored cluster amplitudes to find theappropriate R and α matrices. The m -dimensional lin-ear equation (17) is solved for σ . The next best solutionof the LCCSD equations is then found by substituting σ into Eq. (13). These two steps are repeated until con-vergence is reached according to a specified criteria. Inthe next section, we discuss the results of the applicationof the DIIS and RLE procedures to the solution of theLCCSD equations in the cases that do not converge orconverge to non-physical answers with the conventionaliteration scheme. IV. RESULTS AND DISCUSSION
In this section, we study and compare the convergencecharacteristics of the DIIS and the RLE methods. Weinclude a number of test cases in four different systems,
TABLE I: Convergence tests of the LCCSD equations withCIS, RLE, and DIIS methods for B and Al. CIS is the con-ventional iterations scheme (no convergence stabilizer). RLE5designates RLE convergence method with 5 pre-stored itera-tions. Last column gives resulting correlation energy in a.u.*Cases where maximum number of iterations allowed duringrun was reached.Atom State Method δE v ( a.u. )B Core CIS 21 YesRLE5 9 YesDIIS5 15 YesB 2 p / CIS 18 Yes -0.0293907RLE5 13 Yes -0.0293906DIIS5 13 Yes -0.0293908B 3 s CIS 70* No -0.0091643RLE4 70* No -0.0070438DIIS4 70* No -0.0089454RLE5 30 Yes -0.0089491DIIS5 22 Yes -0.0089472DIIS7 31 Yes -0.0089488B 3 p / CIS 44 Yes -0.0056292RLE5 23 Yes -0.0056294DIIS5 19 Yes -0.0056284RLE8 25 Yes -0.0056295DIIS8 24 Yes -0.0056293B 3 d / CIS 85* No -0.0884489DIIS7 71 Yes -0.0007535RLE8 66 Yes -0.0007536DIIS8 51 Yes -0.0007533Al Core CIS 13 YesRLE5 8 YesAl 3 p / CIS 16 Yes -0.0245810RLE5 11 Yes -0.0245798DIIS5 14 Yes -0.0245811Al 4 s CIS 20 Yes -0.0079907RLE5 13 Yes -0.0079906DIIS5 18 Yes -0.0079905Al 3 d / CIS 70* No -0.0209573DIIS6 89 Yes -0.0208637RLE8 86 Yes -0.0208662DIIS8 49 Yes -0.0208662Al 4 d / CIS 70* No -0.0213727RLE8 300* No 0.0022280DIIS8 81 Yes 0.0022478DIIS9 91 Yes 0.0022477Al 4 d / CIS 70* No -1.3775005RLE8 165 Yes 0.0022737DIIS8 97 Yes 0.0022710DIIS9 73 Yes 0.0022712
B, Al, Zn + , and Yb + which have a large number ofstates that do not converge with the conventional iter-ative scheme (CIS). We also test the ability of the RLEand the DIIS to accelerate convergence in the cases wherethe CIS does converge. The main purpose of these testsis to provide general guidelines of how to accelerate or toachieve convergence using the RLE and the DIIS meth- ods. The conclusions and observations of this sectionmay be extrapolated to other systems for both all-orderand CI+all-order approaches.The summary of B and Al convergence tests is given inTable I. We find that convergence patterns for two fine-structure multiplet states, for example 3 p / and 3 p / states, are generally very similar. Therefore, we list only np / and nd / states with the exception of the 4 d statesof Al. Tests were performed for both states of the mul-tiplet as an additional check, since similar results areexpected. The results are given for the 2 p / , 3 s , 3 p / ,and 3 d / states of B and the 3 p / , 4 s , 3 d / , 4 d / and 4 d / states of Al. The resulting LCCSD correla-tion energy is listed in the last column of the table ina.u. The convergence method is specified in the third col-umn. CIS refers to the initial straightforward iterationscheme. RLE5 designates the RLE convergence methodwith 5 pre-stored CIS iterations. Similarly, DIIS8 refersto the DIIS convergence scheme with 8 pre-stored itera-tions. The fourth column indicates the iteration numberat the end of the run. Cases where the maximum num-ber of iterations allowed during the run was reached aremarked with asterisk. In these cases, convergence did notoccur. The convergence criteria was set to terminate theiteration procedure when the relative difference betweentwo consecutive iterations is reduced below 0.00001. Thesame convergence criteria is used in all valence test runs.Only the core and the 2 p states of B converge with theCIS. In the case of Al, all nd states do not converge withthe CIS. All of the cases in Table I converge with theDIIS8. We note that we did not list the RLE5 and theDIIS5 results for many of the nd states because conver-gence was not achieved. In cases where all methods leadto convergence, both the RLE and the DIIS have accel-erated convergence rates relative to the CIS. (cid:144)(cid:145)(cid:146)(cid:147)(cid:148)(cid:149) (cid:150)(cid:151) (cid:152)(cid:153)(cid:154)(cid:155)(cid:156)(cid:157)(cid:158)(cid:159)(cid:160)¡¢ £⁄ ¥ƒ §¤ '“ «‹ ›fifl(cid:176)–†‡·(cid:181)¶•‚„”»…‰(cid:190)¿(cid:192)`´ˆ˜¯ ˘˙¨(cid:201)˚¸(cid:204)˝˛ˇ—(cid:209)(cid:210)(cid:211)(cid:212)(cid:213)(cid:214)(cid:215)(cid:216)(cid:217)(cid:218)(cid:219)(cid:220)(cid:221)(cid:222)(cid:223)(cid:224)Æ(cid:226)ª(cid:228)(cid:229)(cid:230)(cid:231)ŁØŒº(cid:236)(cid:237)(cid:238)(cid:239)(cid:240)æ(cid:242)(cid:243) (cid:244)ı(cid:246)(cid:247)łøœß(cid:252)(cid:253)(cid:254)(cid:255)(cid:18)(cid:20)(cid:0)(cid:1)(cid:2)(cid:3)(cid:4)(cid:5)(cid:6)(cid:7)(cid:10)(cid:8)(cid:9) (cid:11)(cid:12)(cid:13)(cid:14)(cid:15) (cid:16)(cid:17)(cid:19)(cid:21)(cid:22)(cid:23)(cid:24)(cid:25)(cid:26)(cid:27)(cid:28) (cid:29)(cid:30)(cid:31) !" FIG. 3: (Color online) Comparison of the DIIS5, RLE8, andDIIS8 schemes for the 3 d / state of boron. The correlationenergy is given in a.u. We may draw two general conclusions from our tests:1. If a particular LCCSD run converges with the CIS,then the RLE5 appears to be the most efficientmethod in accelerating the convergence.2. If a particular LCCSD run does not converge withthe CIS, the DIIS8 or the DIIS9 appear to be themost efficient in attaining and accelerating the con-vergence.We note that these two rules are not absolute, but theyserve to be good initial guidelines. Our further tests onother (much heavier) systems confirmed these guidelines.We note that the RLE5 is not sufficient to achieve conver-gence for most of the divergent cases. The only exceptionin Table I is the 3 s state of B. However, while the CISnever converges to our standard criteria for the 3 s state,it nearly converges to correct result before exhibiting di-verging and oscillating pattern of Fig. 1. In this case,accumulation of only 5 iterations is sufficient. However,in the case of the nd states, the CIS is never close to con-verging to a correct number and subsequently the RLE5does not work. Occasionally, DIIS9 may achieve conver-gence where DIIS8 would not. Using even larger numberof stored iterations does not improve convergence or ef-ficiency. DIIS10-DIIS12 runs for the 4 d states convergedto non-physical answers in two instances, but to correctresults in all other cases. Number of iterations variedsignificantly from case to case. The results of all con-verged runs listed in Table I are consistent within theconvergence criteria, as expected.We implemented the DIIS/RLE strategies for two sep-arately developed LCCSD codes. The calculations werecarried out using two different finite basis sets, the B-splines of Ref. [46] and the dual-kinetic-basis sets ofRef. [47]. Even though the basis sets and the conver-gence criteria used for each code made slight differencesin the values, the general observations on the convergencepatterns remains the same.We illustrate different convergence patterns of the RLEand the DIIS methods for the 3 s and 3 d / states of boronin Figs. 2 and 3. In Fig. 2, the values of the correla-tion energies for the 3 s states of boron are obtained fromdifferent schemes that are listed on the graph. RLE5,DIIS5, and DIIS7 results after N = 20 interactions areindistinguishable at this plot scale and are not shown.These schemes converge after 30, 22, and 31 iterations,respectively. While the CIS results appear close to con-verged value, convergence was never reached and corre-lation energy began to oscillate as illustrated in Fig. 1.The RLE5 and the DIIS5 results are identical to the CISones for the first four iterations. The 5th value is differ-ent for three of the schemes as this ( m + 1)-th value (seeEq. (13)) is replaced by the RLE or the DIIS predictionsfor RLE5 and DIIS5. We observe that these predictionsare significantly closer to converged result than the 5thCIS iteration. After that, the RLE5 and DIIS5 resultsare sharply adjusted at N = 10 when the second call tothe RLE/DIIS stabilizer codes is made. The DIIS7 be-havior is similar to the one just described, except thatit accumulates 7 iterations before the DIIS procedure is TABLE II: Convergence tests of the LCCSD equations withCIS, RLE, and DIIS methods for Zn + and Yb + . CIS isthe conventional iterative scheme (no convergence stabilizer).DIIS8 designates the DIIS convergence method with 8 pre-stored iterations. Last column gives resulting correlation en-ergy in a.u. *Cases where maximum number of iterationsallowed during run was reached.Atom State Method δE v ( a.u. )Zn + p / CIS 39 Yes -0.0066119DIIS9 12 Yes -0.0066088Zn + d / CIS 70* No -0.0977215RLE5 67 Yes -0.0045508DIIS8 18 Yes -0.0045511Zn + d / CIS 70* No -0.1149058RLE5 93 Yes -0.0045266DIIS8 18 Yes -0.0045267Zn + d / CIS 70* No -0.1936330RLE5 28 Yes -0.0018929DIIS8 18 Yes -0.0018942Zn + f / CIS 200* No -0.0008697DIIS8 200* No -0.0007949DIIS9 153 Yes -0.0007948Zn + f / CIS 70* No -0.0007970DIIS8 173 Yes -0.0007891DIIS9 195 Yes -0.0007891Yb + Core CIS NoRLE5 12 YesDIIS5 12 Yes invoked and now the 7th value gets much closer to thefinal answer.In Fig. 3, the values of the correlation energies obtainedfrom CIS, DIIS5, RLE8, and DIIS8 are plotted for the3 d / states of boron. The RLE8 and the DIIS8 resultsafter N = 35 appear identical on the graph at this scaleand are not shown. The RLE8 and the DIIS8 converge toour criteria after 66 and 51 iterations, respectively. Verysimilar behavior of the RLE8 and the DIIS8 is observed,with the RLE8 energy oscillations being slightly largerafter the RLE subroutine pass. However, other tests showthat the RLE8 in general converges slower, sometimesdramatically so, than the DIIS8. The CIS values divergecompletely and increase rapidly. The DIIS5 seems to beconverging at N = 35, but does not in fact reach selectedcriteria even after 100 iterations.The summary of the selected Zn + and Yb + conver-gence tests is presented in Table II. The results are givenfor the 5 p / , 4 d / , 4 d / , and 5 d / , 4 f / , and 4 f / states of Zn + and the Yb + core. The 4 s and 4 p j states ofZn + and the 6 s , 6 p j , 7 s , and 5 d j states of Yb + convergewith the CIS, so we have omitted these results from thetable. However, it is worth pointing out that RLE5 ac-celerates convergence for all these states compared to theCIS. Table II demonstrates that the DIIS reduces num- TABLE III: Comparison of B, Al, Zn + , and Yb + removal energies (in cm − ) with experiment [48]. Rows labeled “Dif.” giverelative difference with experimental values in %. B p / p / s p / p / d / d / Expt. -66928 -66913 -26888 -18316 -18314 -12160 -12160SD -67049 -67035 -27105 -18497 -18495 -12494 -12494Dif. -0.18% -0.18% -0.81% -0.99% -0.99% -2.7% -2.7% Al p / s d / d / d / d / Expt. -48278 -22931 -15843 -15842 -6045 -6041SD -48271 -23069 -17295 -17289 -6652 -6647Dif. 0.02% -0.60% -8.4% -8.4% -9.1% -9.1% Zn + s p / p / p / d / d / d / f / Expt. -144691 -96027 -95157 -43360 -47950 -47902 -26913 -27606SD -144684 -96221 -95352 -43421 -47929 -47880 -26898 -27633SDpT -145232 -96559 -95679 -43492 -47994 -47946 -26929 -27633Dif. (SD) 0.14% 0.20% 0.19% 0.24% 0.11% 0.11% 0.09% -0.02%Dif. (SDpT) -0.23% -0.15% -0.15% 0.08% -0.02% -0.03% -0.02% -0.01% Yb + s p / p / s d / d / f / f / Expt. -98207 -71145 -67815 -43903 -75246 -73874 -27704 -27627SD -98961 -71016 -67480 -44060 -76141 -74700 -28080 -28062SDpT -99107 -71084 -67592 -44115 -77764 -76317Dif. (SD) -0.77% 0.18% 0.49% -0.36% -1.19% -1.12% -1.36% -1.57%Dif. (SDpT) -0.91% 0.09% 0.33% -0.48% -3.2% -3.2%
DEFGHI JK LMNOPQRSTUV W X Y Z [ \ ] ^ _‘ ab cdefghijklmnopqrstuvwxyz{|}~(cid:127) (cid:128)(cid:129)(cid:130)(cid:131)(cid:132)(cid:133)(cid:134)(cid:135)(cid:136)(cid:137)(cid:138)(cid:139)(cid:140)(cid:141)(cid:142)(cid:143)(cid:144)(cid:145)(cid:146)(cid:147)(cid:148)(cid:149)(cid:150)(cid:151)(cid:152)(cid:153)(cid:154)(cid:155) (cid:156)(cid:157)(cid:158)(cid:159)(cid:160)¡¢£⁄¥ƒ§¤ '“« ‹›fifl
FIG. 4: (Color online) Comparison of the CIS, RLE5, andDIIS5 schemes for the Yb + core. The correlation energy isgiven in a.u. RLE5 and DIIS5 data appear identical at thisscale and are shown as a single curve. ber of iterations for the 5 p / states by a factor of 3 orbetter. Zn + and Yb + tests confirm our conclusions (1)and (2), on the previous page. We were unable to achieveconvergence for higher 7 p j states in Yb + . This problemis not present in Zn + , where LCCSD for the 5 p statesconverges even with CIS as shown in Table II. Perhapsother convergence approaches are needed to resolve thisissue.The case of Yb + core is particularly interesting, since core iterations generally converge well with the CIS. Yb + core is an exception, however, most likely due to verylarge 4 f shell contributions that lead to oscillation ofthe correlation energy. We plot the CIS, the RLE5, andthe DIIS5 results for the Yb + core correlation energy inFig. 4. The RLE5 and the DIIS5 results appear identicalat this scale and are shown as a single curve. Both meth-ods are successful at fixing the CIS’s oscillation problem.The comparison of the B, Al, Zn + , and Yb + removalenergies with experiment [48] is given in Table III. Rowslabeled “Dif.” give relative difference with experimentalvalues in %. The energies here are given in cm − . Most ofthese states did not converge with the CIS, so it is impor-tant to establish the accuracy of this approach for suchcases. Breit interactions and contributions from higherpartial waves are also included. The B and Al ioniza-tion potentials, B 2 p / , and Al 4 s SD energies are inagreement with experiment. We consider only monova-lent states for all of these systems. The SD approxima-tion does not account for mixing with the core-excitedstates such as 3 s p in Al. Therefore, larger disagree-ment with experiment is expected in cases where mixingwith these hole-two-particle states is large. A particularexample is the 3 d and 4 d states of Al. The lower 3 s nd levels heavily mix with the 3 s p D levels. However,the mixing coefficient for this configuration never exceeds30%. As a result, these levels are distributed over sev-eral lower nd levels, resulting in two sets of levels beinglisted as 3 s d D [48, 49] ([ y D ] and [ D ]). In Table III,we compare the 4 d results with the second sets of levels([ D ]).We also included partial valence triples perturbatively(LCCSDpT) to investigate if the LCCSDpT methodwould improve the theory-experiment agreement for Zn + and Yb + . This method is described in detail in [9]. Sincetriple equations are not explicitly iterated in this ap-proach, implementation of the RLE and the DIIS methodis exactly the same as in the SD code. Convergence testsof the LCCSDpT method exhibit essentially the samepattern as the tests of the LCCSD method discussedabove, and a similar number of iterations was generallyrequired for LCCSD and LCCSDpT calculations for thesame states run with the same parameters.As shown in Table III, we find an excellent agreementof all Zn + data with experiment. Inclusion of perturba-tive triples somewhat improves the agreement with ex-periment for most states. The accuracy decreases forYb + , as expected, owing to much softer and heavier coreand strong mixing of monovalent states with one-hole-two-particle states in this system. Nevertheless, for Yb + the average accuracy for removal energies is at the levelof 1% (see Table III.) V. CONCLUSION
We have successfully implemented the RLE and DIISconvergence techniques in the LCCSD and LCCSDpTmethods for high-precision atomic many-body calcula-tions. Most of the convergence problems were resolvedusing these methods. Acceleration of convergence wasdemonstrated for all cases where all-order equations con-verge with straightforward iteration scheme. Numeroustests were performed to establish general recommenda-tions for the RLE/DIIS use for various purposes. Wefind that if particular case converges with CIS, RLE5appears to be the most efficient in achieving and accel-erating convergence. If particular case does not convergewith CIS, DIIS8 or DIIS9 appear to be the most efficientin accelerating convergence. Solving these convergenceproblems greatly expands the number of atomic speciesthat can be treated with the all-order methods and isanticipated to facilitate many interesting future applica-tions for studies of fundamental symmetries as well asatomic clock and ultracold atom research.
Appendix A
Derivations of general formulas in this Appendixmainly follows Appendix of Ref. [36]. Consider solving ageneral linear equation of the form a + Bt = 0 , (A1)which is a system of linear equations of dimension k withvector t being the exact solution that we would like tofind. We make the best approximation to the exact so-lution by using m ( < k ) nonorthogonal and linearly inde-pendent vectors T = ( t (1) , t (2) , ..., t ( m ) ), where each t ( i ) is a k -dimensional vector. We find this best approxima-tion as a linear combination of t ( i ) ’s: t [ m +1] = m X i =1 σ i t ( i ) = σ · T . (A2)Here, σ i are the weights of the optimized solution thatneeds to be determined. We note that t ( m +1) is not in-cluded in the linear combination (A2), but is used toconstruct matrices such as shown in Eqs.(20) and (21).Therefore, t ( m +1) needs to be also found through the CIS.Therefore, we use the notation t [ m +1] instead of t ( m +1) todistinguish between the m + 1 th solution found throughthe use of RLE/DIIS methods and the CIS result, respec-tively.First, we try to derive an ideal equation to find σ i ’sas if we know the exact solution to Eq. (A1). To findthe best approximation, we need to minimize the errorbetween the approximate and the exact answers. To thisend, we use the least square optimization approach. Theerror is e = t − t [ m +1] . The least square optimization of E = e T e with respect to σ then yields ∂E∂σ = − T T ( t − T · σ ) = 0 . (A3)After solving for σ and substituting it in Eq.(A2), we get t [ m +1] = T ( T T T ) − T T t . (A4)However, not knowing what the exact solution t is,the above formula is of little use. The DIIS and RLEare based on replacing t with approximations. Substi-tuting t [ m +1] instead of t in Eq.(A1), will make Eq.(A1)inhomogeneous: a + Bt [ m +1] = a + BT · σ = ǫ . (A5)where ǫ is a vector with constant elements.The difference between the RLE and DIIS methods isin their choice of error to minimize, e . The DIIS takesthe error to be ǫ of Eq. (A5). Then to get the best ap-proximation, we need to minimize E = ǫ T ǫ with respectto σ : ∂E∂σ = 2( − BT ) T ( a + BT σ ) = 0 . (A6)Therefore, the coefficients σ that lead to the best approx-imation satisfy the DIIS equation: T T B T a + T T B T BT σ = 0 . (A7)The RLE requires that the best least square approxi-mation ǫ [ m +1] to ǫ vanishes in the space of T . Followingthe structure of Eq. (A4): ǫ [ m +1] = T ( T T T ) − T T ǫ = T ( T T T ) − T T ( a + BT σ ) = 0 . (A8)0Since T is made of linearly independent vectors, ǫ [ m +1] is only zero if: T T ( a + BT σ ) = 0 . (A9)Eqs. (A7) and (A9) for DIIS and RLE, respectively,correspond to Eqs. (15) and (16) in the paper. Acknowledgments
The work of H.G. and A.D. was supported in part bythe US National Science Foundation Grant No. PHY-9-69580. The work of M.S.S. was supported in part byNational Science Foundation Grant No. PHY-07-58088. [1] J. ˇC`ıˇzek, J. Chem. Phys. , 4256 (1966).[2] R. Bartlett, J. Phys. Chem. , 1697 (1989).[3] I. Lindgren and J. Morrison, Atomic Many–Body Theory (Springer–Verlag, Berlin, 1986), 2nd ed.[4] F. Coester and H. K¨ummel, Nuclear Physics , 477(1960).[5] S. A. Blundell, W. R. Johnson, Z. W. Liu, andJ. Sapirstein, Phys. Rev. A , 3768 (1989).[6] S. A. Blundell, W. R. Johnson, Z. W. Liu, andJ. Sapirstein, Phys. Rev. A , 2233 (1989).[7] S. A. Blundell, W. R. Johnson, and J. Sapirstein, Phys.Rev. A , 3407 (1991).[8] S. A. Blundell, J. Sapirstein, and W. R. Johnson, Phys.Rev. D , 1602 (1992).[9] M. S. Safronova, W. R. Johnson, and A. Derevianko,Phys. Rev. A , 4476 (1999).[10] A. Derevianko, W. R. Johnson, M. S. Safronova, andJ. F. Babb, Phys. Rev. Lett. , 3589 (1999).[11] M. S. Safronova and W. R. Johnson, Phys. Rev. A ,022112 (2000).[12] M. S. Safronova and W. R. Johnson, Adv. At. Mol., Opt.Phys. , 191 (2007).[13] R. Pal, D. Jiang, M. S. Safronova, and U. I. Safronova,Phys. Rev. A , 062505 (2009).[14] S. G. Porsev, K. Beloy, and A. Derevianko, Phys.Rev. Lett. , 181601 (pages 4) (2009), URL http://link.aps.org/abstract/PRL/v102/e181601 .[15] S. G. Porsev, K. Beloy, and A. Derevianko, Phys. Rev.D , 036008 (2010).[16] C. S. Wood, S. C. Bennett, D. Cho, B. P. Masterson,J. L. Roberts, C. E. Tanner, and C. E. Wieman, Science , 1759 (1997).[17] V. V. Flambaum, V. A. Dzuba, and A. Derevianko, Phys.Rev. Lett. , 220801 (2008).[18] K. Beloy, U. I. Safronova, and A. Derevianko, Phys. Rev.Lett. , 040801 (2006).[19] M. S. Safronova, D. Jiang, and U. I. Safronova, Phys.Rev. A , 022510 (2010).[20] D. Jiang, B. Arora, M. S. Safronova, and C. W. Clark,J. Phys. B , 154020 (2009).[21] M. S. Safronova, D. Jiang, B. Arora, C. W. Clark, M. G.Kozlov, U. I. Safronova, and W. R. Johnson, IEEE Trans.on Ultrasonics, Ferroelectrics, and Frequency Control ,94 (2010).[22] A. Derevianko, Phys. Rev. Lett. , 033002 (2010).[23] M. J. Morrison, V. A. Dzuba, and A. Derevianko, Phys.Rev. A , 013604 (2011).[24] B. Ravaine, A. Derevianko, and P. R. Berman, Phys.Rev. A , 022330 (2006).[25] B. Arora, M. S. Safronova, and C. W. Clark, Phys. Rev.A , 052509 (2007). [26] M. S. Safronova, C. J. Williams, and C. W. Clark, Phys.Rev. A , 040303(R) (2003).[27] M. S. Safronova, M. G. Kozlov, W. R. Johnson, andD. Jiang, Phys. Rev. A , 012516 (2009).[28] C. W. Chou, D. B. Hume, J. C. J. Koelemeij, D. J.Wineland, and T. Rosenband, Phys. Rev. Lett. ,070802 (2010).[29] T. Rosenband, D. B. Hume, P. O. Schmidt, C. W. Chou,A. Brusch, L. Lorini, W. H. Oskay, R. E. Drullinger,T. M. Fortier, J. E. Stalnaker, et al., Science , 1808(2008).[30] V. A. Dzuba and A. Derevianko, J. Phys. B , 074011(2010).[31] C. Tamm, S. Weyers, B. Lipphardt, and E. Peik, Phys.Rev. A , 043403 (2009).[32] K. Tsigutkin, D. Dounas-Frazer, A. Family, J. E. Stal-naker, V. V. Yashchuk, and D. Budker, Phys. Rev. Lett. , 071601 (2009).[33] V. V. Ivanov, A. Khramov, A. H. Hansen, W. H. Dowd,F. Muenchow, A. O. Jamison, and S. Gupta, ArXiv e-prints (2011), 1101.5142.[34] S. Tassy, N. Nemitz, F. Baumer, C. H¨ohl, A. Bat¨ar, andA. G¨orlitz, J. Phys. B , 205309 (2010).[35] S. G. Porsev, Yu. G. Rakhlina, and M. G. Kozlov, Pis’maZh. Eksp. Teor. Fiz. , 449 (1995), [JETP Lett. , 1284(1981).[37] G. W. Trucks, J. Noga, and R. J. Bartlett, ChemicalPhysics Letters , 548 (1988).[38] P. Pulay, Chem. Phys. Lett. , 393 (1980).[39] P. Pulay, Journal of Computational Chemistry , 556(1982).[40] N. S. Mosyagin, E. Eliav, and U. Kaldor, Journal ofPhysics B Atomic Molecular Physics , 339 (2001),arXiv:physics/0009060.[41] R. J. Harrison, Journal of Computational Chemistry ,328 (2004).[42] E. Eliav, M. J. Vilkas, Y. Ishikawa, and U. Kaldor, J.Chem. Phys. , 224113 (2005).[43] C. Lanczos, Journal of Research of the National Bureauof Standards , 225 (1950).[44] A. Derevianko and E. D. Emmons, Phys. Rev. A ,012503 (2002).[45] M. S. Safronova, Ph.D. thesis, University of Notre Dame(2000).[46] W. R. Johnson, S. A. Blundell, and J. Sapirstein, Phys.Rev. A , 307 (1988).[47] K. Beloy and A. Derevianko, Comp. Phys. Comm. ,310 (2008).[48] J. E. Sansonetti, W. C. Martin, and S. L. Young, Hand- book of Basic Atomic Spectroscopic Data (NIST, 2006),http://physics.nist.gov/PhysRefData/Handbook/.[49] V. Kaufman and W. C. Martin, Phys. Chem. Ref. Data20