G L G Sleijpen

H A Van Der Vorst

H J J Van Dam

J H Van Lenthe

CiteSeerX

An Improvement of Davidson’s Iteration 
Method: Applications to MRCI and 
MRCEPA Calculations 
H. J. J. VAN DAM* and J. H. VAN LENTHE 
Theoretical Chemistry Group, Debye Institute, Utrecht University, Padualaan 14, 3584 CH Utrecht, 
The Netherlands 
G. L. G. SLEIJPEN and H. A. VAN DER VORST 
Department of Mathematics, Utrecht University, Budapestlaan 6, 3584 CD Utrecht, The Netherlands 
Received 23 November 1994; accepted 12 May 1995 
ABSTRACT 
Davidson’s method is widely used for finding the lowest eigenvalues of large 
matrices. Recently, mathematicians have shown that Davidson’s derivation 
could be improved. They have corrected the derivation yielding a new iteration 
method. In this article this new method is adapted for realistic MRCI and 
MRCEPA calculations. Results show that the new method converges sigruficantly 
faster in H,O and 0, with moderately elongated bonds than Davidson’s 
original method. The new method offers new insights into the rate of convergence 
of Davidson’s original method. 0 1996 by John Wiley & Sons, Inc. 
Introduction 
onfiguration interaction (CI) methods are C well established and are widely applied. In 
CI calculations one aims at the few lowest eigen- 
values and the corresponding eigenvectors of CI 
Hamilton matrices. Because these matrices are 
large and sparse, iterative methods are the meth- 
ods of choice. These methods require a 
matrix-vector multiplication with the CI matrix in 
*Author to whom all correspondence should be addressed. 
Journal of Computational Chemistry, Vol. 17, No. 3, 267-272 (1996) 
0 1996 by John Wiley & Sons, Inc. 
every iteration. Because the CI matrices are so 
large, these matrix-vector multiplications deter- 
mine the cost of the algorithm. Therefore, almost 
every modification that reduces the total number 
of these matrix-vector multiplications can be ap- 
plied to reduce the total cost. 
The eigenvalue problem solver which is most 
often used in CI calculations was proposed by 
Da~idson.’-~ In this method the eigenvectors are 
expanded in a subspace and the subspace is ex- 
tended with an update vector in every iteration. 
Recently, mathematicians have studied Davidson’s 
method. Saad4 pointed out that Davidson’s deduc- 
CCC 01 92-8651 1 96 1030267-06 
VAN DAM ET AL. 
tion of the update vector did not lead to an opti- 
mal choice. Sleijpen et a1.5 have corrected the de- 
duction. Their results lead to extended insight in 
the convergence properties of Davidson’s method 
and suggest improvements. 
In this article we summarize the discussion that 
resulted in the corrected derivation. In addition, 
we describe the application of the improved 
method in multireference CI (MRCI) calculations. 
The results are compared with the results obtained 
with Davidson’s original method. A related 
method has been proposed in J. Olsen, P. Jmgen- 
sen, and J. Simons, Chern. Phys. Lett., 169, 463 
(1990). 
Davidson’s Update Vector 
In his original article, Davidson’ started from 
the observation that the lowest eigenvalue of a 
matrix is at the minimum of the corresponding 
Rayleigh quotient: 
Therefore, if the desired eigenvector is approxi- 
mated by a vector lu), one is interested in an 
update vector IS) that minimizes R(lu) + IS)). To 
this end, Davidson expanded the Rayleigh quo- 
tient to second order in the update vector and 
minimized the expression yielding the linear sys- 
tem 
(A - AI)IS) = -1r) (2) 
where Ir) is the residual vector 
Ir) = (A - A1)lu) 
and A = R(lu)). If the matrix (A - AI) is diago- 
nally dominant, that is6 
(3) 
k # j  
then eq. (2) may be approximated by 
(D - AI)IS) = -Ir) (4) 
Equation (4) is the famous Davidson precondi- 
tioning. The fact that tlus method results from a 
second-order energy expression should explain the 
good convergence properties compared to gradient 
methods. 
However, Saad4 pointed out that Davidson’s 
derivation is not complete. Following Davidson’s 
line of reasoning, the best results would be ob- 
tained if one solves eq. (2) for IS). This would 
result in IS) = -lu). However, after orthogonal- 
ization of IS) on the subspace spanned by update 
vectors from earlier iterations, the null vector re- 
mains. Clearly, this vector cannot serve as an effec- 
tive extension of the subspace. Consequently, the 
update vector IS) should be linearly independent 
to the subspace whereas lu) is a linear combina- 
tion in that subspace. This suggests that Davidson’s 
argument is not complete. In a recent article, 
Sleijpen et aL5 pointed out that while deriving eq. 
(2), Davidson assumed IS) orthogonal to lu). 
However, to assure that IS) is an effective update 
vector, one should determine IS) under the condi- 
tion that IS) is orthogonal to Iu). This may be 
effected by computing the update vector from the 
projection of A onto the subspace orthogonal to 
lu). Sleijpen et al. showed that this is equivalent to 
substituting the projected matrix 
(5) B = (I  - Iu)(uI)A(I - Iu>(uI) 
for A in eq. (21, yielding 
(B - AI)IS) = -1r) (6) 
If the matrix (B - AI) is diagonally dominant, 
(7) 
where D, is the diagonal of B. This equation is 
analogous to eq. (4). Moreover, if the matrix (A - 
AI) is strongly diagonally dominant, then a unit 
vector is a good approximation to the eigenvector 
and to lu). In that case, eq. (7) will only differ from 
eq. (4) in the first element. The difference in the 
resulting update vectors nearly vanishes upon or- 
thogonalizing the update vector on the subspace. 
Therefore, both eqs. (4) and (7) may then be con- 
sidered equivalent. However, eqs. (2) and (6) are 
certainly not equivalent. Due to Saad, we know 
that a calculation based on eq. (2) will not con- 
verge. To investigate the convergence of calcula- 
tions with eq. (61, a calculation on H,O was per- 
formed. The results are shown in Table I (compu- 
tational details are provided later in this article). 
Note that the convergence of the energy is at least 
second order. That is, the number of converged 
digits is doubled in every it eration. This is in 
accordance with the fact that eq. (6) resulted from 
a second-order energy expression. 
For practical applications, we will concentrate 
on eq. (6). Approaches based on this equation will 
be referred to by the name generalized Jacobi- 
Davidsod (GJD). When Davidson’s method con- 
eq. (6) may be approximated to give 
(DB - AI)/S) = -Ir) 
268 VOL. 17, NO. 3 
DAVIDSON'S ITERATION 
TABLE 1. 
The Convergence of the Generalized Jacobi- 
Davidson Method in Extreme for H,O at 2.0. Re. 
Iteration 
Number Energy (Hartreel 
- 75.839984920556 
- 75.950820825635 
- 75.951 032652229 
- 75.951 032652292 
verges slowly due to large off-diagonal elements, 
eq. (6) may allow approaches with better conver- 
gence properties. However, in its current formula- 
tion, the equation is rather impractical because of 
the size of the matrix (B - AI). Therefore, some 
suitable approximations to this matrix must be 
found. 
Application of the Generalized 
Jacobi-Davidson Method to 
MRCI Problems 
In MRCI we distinguish two categories of 
molecular orbitals (MOs): internal orbitals and ex- 
ternal  orbital^.^ The set of internal orbitals contains 
the MOs which are of primary importance in the 
construction of a qualitatively correct wave func- 
tion. From the MOs, two types of n-electron states 
may be constructed: vacuum states and external 
states. The vacuum states have no electrons occu- 
pying external orbitals. All other states are external 
states. Based on these definitions, CI vectors can be 
separated into a vacuum part and an external part. 
We will denote a vacuum part with a subscript v 
and an external part with a subscript e. 
In MRCI calculations with a suitably chosen 
vacuum space, the vacuum coefficients are much 
more important than the external coefficients. This 
suggests that if the projector P is defined as 
this may be approximated by 
P -  [p;. 3. (9) 
Introducing this approximation into the matrix 
(B - AI) from eq. (6) yields 
(B - AI) 
1 (I - P)vvA,,(I - P)vv (I - P)vvAve Aev(I - P)w A ee 
-A[';. c] 
Furthermore, we assume that (B - AI) is diago- 
nally dominant in the external space. This means 
we assume that we may approximate the matrix 
by 
(B - AI) 
1 - [(I  - P)vvAv,(I - P>w - AI,, 0 0 Dee - AIee 
(10) 
where Dee is the diagonal of the external-external 
block. Substituting eq. (10) into eq. (6) results in 
(B - AI)wlS), = -1rh (11) 
(D - A1)eeIS)e = -Ir)e (12) 
Because the linear system in the vacuum space 
is relatively small, we expect that these equations 
may be cost effectively applied to realistic MRCI 
problems. Note that in Table C18 calculations, sim- 
ilar approximations can be used based on a selec- 
tion of the most important configurations. Further- 
more, eqs. (11) and (12) may also be applied to 
MRCEPA9 calculations because MRCEPA may be 
expressed in a form analogous to MRCI. That is, 
the MRCEPA matrix is equivalent to a MRCI ma- 
trix with shifted diagonal elements. 
The linear system in eq. (11) can be treated with 
Gaussian elimination if the matrix is small. How- 
ever, if the vacuum space is large (but small com- 
pared to the total CI space), an iterative method 
such as conjugate may be more effi- 
cient. In the conjugate gradient method, a 
matrix-vector multiplication with the vacuum part 
of the matrix (B - AI) is performed in every itera- 
tion. These multiplications may be computed effi- 
ciently as follows. Assume the vacuum matrix 
(B - AI) is to be multiplied with a vector 14). 
Then, using eq. (51, we have 
(B - AI)lq) = Alq) - ld(ulAIq) - Alu>(ulq> 
+ Ib>(uIAIu)(uIq) - Iq)A (13) 
Writing 
A = (ulAlu) 
JOURNAL OF COMPUTATIONAL CHEMISTRY 269 
VAN DAM ET AL. 
defining 
12) = Alu) 
and using that A is hermitian, eq. (13) becomes 
(B - AOIq)  = AIq) + Iu)((uIq)A - (ZIq)) 
- Iz)(uIq) - Iq)A (14) 
Equation (14) involves a matrix-vector product 
with the original vacuum matrix A and some sim- 
ple vector-vector operations. The vector lz) can be 
obtained without additional cost because it is al- 
ready calculated in the Davidson algorithm in the 
computation of the residual vector. 
Symmetric Dissociation of H,O 
The first test case we consider involves the 
cleavage of both OH bonds at a constant angle in 
the water molecule. The geometries and the basis 
sets of Bauschlicher and Taylor” are used. This 
means that the energy of the molecule was calcu- 
lated at the equilibrium OH bond length R,, at 
1.5.Re and at 2.0.Re, where R, was 1.889 726 
Bohr. Additionally, we also calculate the energy at 
4.0, 6.0, and 8.0 times RE.  The HOH angle is fixed 
at 104.5”. The 1 s  orbital is frozen on the self-con- 
sistent field (SCF) level. The active space in the 
Multi Configuration SCF (MCSCF) calculation also 
defines the reference space in the MRSDCI and 
MRCEPA calculations. The results in Tables I1 and 
I11 are obtained with the 55 Configuration State 
Function (CSF) Complete Active Space (CAS) from 
Ruttink et al? The initial vector is computed by 
solving the reference part of the CI matrix for the 
lowest eigenvalue. The calculations are performed 
with the ATMOL program pa~kage.’~ 
As a first result, we found that the GJD itera- 
tions are almost as efficient as the Davidson itera- 
tions in terms of central processing unit (CPU) 
time. This is in accordance with the assumption 
that solving a linear system in the vacuum space is 
efficient compared with a matrix-vector multipli- 
cation with the complete CI matrix. Furthermore, 
the results in Table I1 show that Davidson’s method 
and the GJD method have comparable conver- 
gence rates at the extreme geometries. To explain 
this, we have made a comparison of the diagonal 
dominance of the matrices at the different geome- 
tries. This comparison is a tentative one because a 
robust measure for diagonal dominance has not 
been defined in the literature. We found that 
Davidson’s method converges fast at the extreme 
geometries because the matrix (B - AI) is rela- 
tively diagonally dominant. However, this diago- 
nal dominance has different causes in the equilib- 
rium geometry and near the dissociation limit. On 
the one hand, at the equilibrium geometry the 
matrix (B - AI) is diagonally dominant due to 
the structure of the CI matrix. This is obvious 
from the resulting normalized CI vector, which has 
one element that is larger than 0.96. Therefore, the 
projections contained in B have no significant ef- 
fect on the diagonal dominance of (B - AI). On the 
other hand, near the dissociation limit the CI ma- 
trix has large off-diagonal elements. This is re- 
flected in the resulting CI vector, which has five 
components exceeding 0.42. However, the matrix 
(B - AI) is diagonally dominant because the off- 
diagonal elements are reduced due to the projec- 
tions contained in B. The GJD method converges 
significantly faster with moderately elongated 
bonds than Davidson’s method. This improved 
convergence is due to a more accurately calculated 
reference part of the update vector. 
TABLE II. 
Comparison of the Convergence of Davidson and the Generalized Jacobi-Davidson for H,O at Various 
Geometries in MRSDCI. 
No. of Iterations 
OH Bond 
Length (Re)  Energy (Hartree) Davidson 
Generalized 
Jaco bi-Davidson 
1 .o 
1.5 
2.0 
4.0 
6.0 
8.0 
- 76.25591 04 
- 76.0706700 
- 75.951 0327 
- 75.91 291 60 
- 75.91 28205 
- 75.91 281 77 
9 
13 
14 
11 
12 
12 
8 
11 
10 
7 
9 
9 
The reference space is a CAS containing the 3a1, 4a1, lb,, 2 4 ,  lb,, and 2b, orbitals. The convergence threshold is 
i - E < 0.1E - 6, where is the current approximation to the eigenvalue and E is the fully converged eigenvalue. 
270 VOL. 17, NO. 3 
DAVIDSON'S ITERATION 
TABLE 111. 
Comparison of the Convergence of Davidson and the Generalized Jacobi-Davidson for H,O at Various 
Geometries in MRCEPA. 
No. of Iterations 
OH Bond Generalized 
Length (Re) Energy (Hattree) Davidson Jacobi-Davidson 
1 .o 
1.5 
2.0 
4.0 
6.0 
8.0 
-76.2590888 
- 76.0731 227 
- 75.9531 400 
- 75.9148751 
- 75.91 47778 
-75.9147748 
8 
10 
12 
12 
12 
12 
The reference space is a5AS containing the 3a,, 4a,, lb,, 2b,, 2b, ,  lb,, and 2b, orbitals. The convergence threshold is 
E - E < 0.1E - 6, where E is the current approximation to the eigenvalue and E is the fully converged eigenvalue. 
The convergence behavior of both methods is 
shown Figure 1. The data for this figure are taken 
from a calculation at 2.0. Re.  In this calculation the 
energy is converged to 14 digits. The figure shows 
that Davidson's method and the GJD method ini- 
tially converge at almost equal rates. However, 
Davidson's method slows down after a few itera- 
tions whereas the GJD method maintains an al- 
most constant convergence rate. 
Dissociation of 0, 
We consider the dissociation of 0, as a second 
test case. The calculations were performed using a 
0 
h -5 
al 
8 
3-10 
- &
M 
-15 
(Ss, 5 p )  - [4s, 2 p ]  double zeta basis.14 The energy 
of the ground state was computed at three geome- 
tries: at the equilibrium bond length Re,  at 1.5 * Re,  
and at 2.0 * Re,  where Re was 2.28 Bohr. The CAS 
contained the 3agr lb,,, lb,,, 3b,,, lb,,, and lb,, 
orbitals. The initial CI vector is computed as de- 
scribed for the H,O test case. The MRCI and 
MRCEPA results are shown in Tables IV and V, 
respectively. 
Again, Davidson's method and the GJD method 
converge at a comparative rate near the equilib- 
rium geometry in the MRCI calculations. How- 
ever, at moderately elongated bonds the GJD 
method converges significantly faster in MRCI as 
well as in MRCEPA calculations. 
I I I I 
Original Davidson 
Generalized 
Jacobi-Davidson 
........ 0." ..... 
3 0 2 2 r4 v1 
Number of iterations 
FIGURE 1. Comparison of the convergence of Davidson and Generalized Jacobi-Davidson. 
JOURNAL OF COMPUTATIONAL CHEMISTRY 271 
VAN DAM ET AL. 
TABLE IV. 
Comparison of the Convergence of Davidson and the Generalized Jacobi-Davidson for 0, at Various 
Geometries in MRSDCI. 
No. of Iterations 
0, Bond 
Length (Re)  Energy (Hartree) Davidson 
Generalized 
Jaco bi-Davidson 
1 .o - 149.8404582 9 9 
1.5 - 149.7729223 13 12 
2.0 - 149.7366070 16 12 
The reference space is a CAS containing the 3a,, lb,,, lb,g, 3b1,, lb,,, and lb,, orbitals. The convergence threshold is 
i - E < 0.1E - 6, where i is the current approximation to the eigenvalue and E is the fully converged eigenvalue. 
TABLE V. 
Comparison of the Convergence of Davidson and the Generalized Jacobi-Davidson for 0, at Various 
Geometries in MRCEPA. 
0, Bond 
Length (Re)  
No. of Iterations 
Generalized 
Energy (Hartree) Davidson Jaco bi-Davidson 
1 .o 
1.5 
2.0 
- 149.8468435 
- 149.7802435 
- 149.7425726 
10 
15 
21 
10 
12 
13 
The reference space is a CAS containing the 3a,, lb,,, lb,,, 3bl,, lb,,, and lb,, orbitals. The convergence threshold is 
€ - E < 0.1E - 6 ,  where € is the current approximation to the eigenvalue and E is the fully converged eigenvalue. 
Conclusions 
We derived a formulation of the modified 
Davidson method due to Sleijpen et al.5 that is 
applicable to realistic MRCI and MRCEPA calcula- 
tions. In calculations on H,O and 0,, the new 
method converges significantly faster than the tra- 
ditional Davidson method for geometries with 
moderately elongated bonds. The improved con- 
vergence properties are due to the more accurate 
calculation of the vacuum part of the update vec- 
tor. Consequently, it is expected that this method 
is effective if the vacuum space is well chosen. 
References 
1. E. R. Davidson, J. Camp. Phys., 17, 87 (1975). 
2. E. R. Davidson, Cam. Phys. Comm., 53, 49 (1989). 
3. C. W. Murray, S. C. Racine, and E. R. Davidson, J. Camp. 
Phys., 103, 382 (1992). 
4. Y. Saad, Numerical Methods for Large Eigenvalue Problems, 
Manchester University Press, Manchester, 1992. 
5. G. L. G. Sleijpen and H. A. van der Vorst, A Generalized 
lacobi-Davidson Iteration Method for Linear Eigenvalue Prob- 
lems, Preprint 856, Dept. Math., Utrecht University, 1994. 
6. E. Kreyszig, Advanced Engineering Mathematics, Wiley, Sin- 
gapore, 1993. 
7. V. R. Saunders and J. H. van Lenthe, Mol. Phys., 48, 923 
(1983). 
8. R. J. Buenker, In Studies in Physical and Theoretical Chemistry, 
R. Carb6, Ed., Vol. 21, Elsevier Scientific Publishing Com- 
pany, Amsterdam, 1982, p. 17. 
9. P. J. A. Ruttink, J. H. v. Lenthe, R. Zwaans, and G. C. 
Groenenboom, I. Chem. Phys., 94, 7212 (1991). 
10. P. E. S. Wormer, F. Visser, and J. Paldus, J. Camp. Phys., 48, 
23 (1982). 
11. G. H. Golub and C. F. van Loan, Matrix Computations, 
North Oxford Academic, Oxford, 1986, p. 363. 
12. C. W. Bauschlicher and P. R. Taylor, J. Chem. Phys., 85,2779 
(1986). 
13. V. R. Saunders and M. F. Guest, ATMOL, 1974. 
14. F. B. v. Duijneveldt, Gaussian Basis Setsfor the Atoms H-Ne 
for use in Molecular Calculations, RJ 945, IBM Research San 
Jose, 1971. 


An Improvement of Davidson&apos;s Iteration Method: Applications to MRCI and MRCEPA Calculations

http://citeseerx.ist.psu.edu/viewdoc/download;jsessionid=8570D5B21805B33B1A3C09940DB051AE?doi=10.1.1.1078.9037&rep=rep1&type=pdf

An Improvement of Davidson's Iteration Method: Applications to MRCI and MRCEPA Calculations

Abstract

Similar works

Full text

Available Versions

CiteSeerX

An Improvement of Davidson&apos;s Iteration Method: Applications to MRCI and MRCEPA Calculations

Abstract

Similar works

Full text

Available Versions

CiteSeerX

An Improvement of Davidson's Iteration Method: Applications to MRCI and MRCEPA Calculations