Computer Physics Communications ELSEVIER
Computer Physics Communications 120 (1999) 155161 www.elsevier.nl/ locate/cpc
Lanczos algorithm and energyweighted sum rules for linear response C.W. Johnson a'l, G.F. Bertsch b, W.D.
Hazelton b
a Department of Physics and Astronomy, Louisiana State University, Baton Rouge, LA 708034001, USA b Department of Physics and Institute Nuclear Theory, Box 351560, University of Washington, Seattle, WA 98915, USA
Received 14 December 1998
Abstract
An iterative algorithm is presented for solving the RPA equations of linear response. The method optimally computes energyweighted moments of the strength function, allowing one to match the computational effort to the intrinsic accuracy of the basic meanfield approximation, avoiding the problem of solving very large matrices. For local interactions, the computational effort for the method scales with the number of particles Np as O(N3). © 1999 Elsevier Science B.V. All rights reserved. PACS: 21.60.Jz; 31.15.Ew Keywords: Lanczos algorithm; Meanfield approximation; Linear response
1. I n t r o d u c t i o n
In a number of branches of physics, mean field theory gives a remarkably effective approximation to the ground state. Similarly, for the response of the system to small perturbations, timedependent mean field theory is a useful extension. This is the experience in nuclear physics [ 1,2], atomic and molecular physics [ 3  8 ] and condensed matter physics [9,10]. The meanfield approximation is inherently limited; thus, rather than squander the large computational resources required for systems of interest, one is motivated to look for algorithms that best match the computational effort. We use convergence of energyweighted moments of the response as a criterion of efficiency. In this paper we describe the algorithm, prove that it preserves, after a welldefined number of iterations, certain energyweighted sum rules of the response, and then give a simple example wherein we illustrate the convergence graphically. We take our inspiration from the Lanczos algorithm [ 11 ], which is best known in manybody physics for extracting lowlying eigenstates of very large Hamiltonian matrices [ 12,13]. When dealing with large spaces, the computational question often comes down to the number of times the Hamiltonian operates on a state vector. I Corresponding author; email:
[email protected] 00104655/99/$  see front matter (~) 1999 Elsevier Science B.V. All rights reserved. PII S00104655 (99) 002489
C W. Johnson et al./Computer Physics Communications 120 (1999) 155161
156
Depending on the Hamiltonian and the starting vector, the Lanczos algorithm is able to extract an accurate ground state vector in a basis of 1056 states with a few hundred Hamiltonian operations. The algorithm may be viewed as a numerically stable technique [ 14] to compute moments of the Hamiltonian with respect to some initial state !/0, that is, /zk  (g'01 f/k i~0) " For large k, /Zk is dominated by the extremal eigenvalues [15], which are thus available for recovery. The Lanczos algorithm has also been applied to many other topics in atomic, molecular, solid state, and nuclear physics, including computation of the Smatrix [16], timeevolution of wave packets [ 17], level densities [ 18], and the continuedfraction expansion of the resolvent or Green's function [19]. Particularly relevant to us is the application to strength functions. The strength function S for an operator Q on a state i is defined as
S ( E )  ~'~ 8 ( E  E: + Ei) I(fl 0 1012 • f
(1)
A powerful technique to calculate the strength function, successfully applied to the nuclear shell model [ 13,20], uses the Lanczos algorithm with a starting vector I~0) = (~ Ii)/(ilO21i) ~/2. The Lanczos algorithm implicitly computes the moments
Mk = /
d E ( E  El)kS(E)
of the strength function. After a few tens of iterations one can accurately reconstruct the distribution of the exact strength function. In many cases, however, the matrix elements of the operator Q are sensitive to correlations in the ground state, and then the size of the wave function basis in the straightforward Hamiltonian approach becomes problematic. In this situation, the timedependent meanfield theory offers a reasonable compromise. The small amplitude theory, the RPA or linear response, can be cast into a matrix form in a particlehole basis. However, the RPA matrix is not symmetric as required by the Lanczos algorithm. The matrix equation is commonly written as B
A
= to
'
(2)
where A and B are particlehole Hamiltonian matrices, to is the eigenfrequency, and x and y are the vectors of positive and negativefrequency particlehole amplitudes, respectively. An important property of the RPA equation is that eigenvectors come in conjugate pairs: in Eq. (2) (y, x) is also an eigenvector with eigenfrequency  t o . For the linear response, the matrix element between the RPA ground state 10) and an excited state Ito) may be expressed as
(tolO_'O) = (q,q) " ( xy)
,
(3)
where q is the vector of particlehole matrix elements and the vector (x,y) is normalized as 1 = x . x  y . y . There are a number of ways to introduce a Lanczostype algorithm for the RPA matrix. The method we describe here has the advantages that it preserves the form equation (2) of the RPA matrix and it produces strength functions that respect sum rules. The algorithm is similar to that of Olsen, Jensen, and Jorgenson [21] (who apply it to the response of ethylene and neon), but we can show explicitly the conservation of sum rules derived below. Starting from Eq. (2) we seek a new basis of vectors [Zi) := (Xi, Yi) where the matrices of column vectors U : ( X 1 , X 2 , X 3 . . . . ) and V := (Y], Y2, Y3 . . . . ) transform the RPA matrix as
v:)
(o

V
Bl
AI
,I
(4)
C.W. Johnson et al./Computer Physics Communications 120 (1999) 155161
157
where the transformed matrices A ~ and B t are now tridiagonal,
al 0
At
=
e2 a2
a2 e3
,
Bt =
bl 0
d2 b2 b2 d3
".o
(5) "..
The Lanczos basis vectors and matrix elements are generated iteratively as follows. Suppose we have the vectors IZl) ..... IZn) already computed, together with the transformed matrix up to enl, dnl, an1 and bnI. The iteration starts by applying the RPA matrix in Eq. (2) to the vector IZ,),
Xt
(6)
 B X n  AYn / " The diagonal elements en and dn are now easily computed,
en = Xt • Xn  Yt " Yn, dn = Xt • Yn  Yt " Xn.
(7)
We next project out IZ~), the component of IZt) that is orthogonal to the space IZI) ..... IZn). This can be done conveniently by using the matrix elements in (5) that have already been calculated,
,Z~)= f Xtt) (XtenXn+dnYnanlXnl +bn,Yn,) \ Ytt = Yt dnXn + enYn bnlXn1 + anlYnI
"
(8)
The norm of the vector IZt) is then computed as
.IV" = X't • Xtt  Y~t" Ytt.
(9)
The norm can be negative, and the definition of the new vector IZn+l) depends on the sign. In fact, because we are actually doing blockLanczos, implicitly operating not only on the vector (X, Y) but also its RPA conjugate (Y, X) simultaneously, there is a degree of freedom, corresponding to a hyperbolic rotation, in choosing the new vector. The simplest choice for the vectors and corresponding RPA matrix elements is
IZn+l)
=~ 1 (
Xtt ) yt t ,
an+l = V~,
bn+l = 0 ,
an+] = 0
b n. + , = x /  ~ ,
.A/" > 0
(10)
and
Izo+,)
1(
Ytt ) . xt,
.
.
.A/" < 0
(11)
This completes the iteration cycle. Finally, in analogy with the application to strength functions in the nuclear shell model, we start with the vector given by
With such a starting vector, we now prove that the algorithm manifestly preserves the energyweighted sum rules, Mk = Z for k odd.
to,k (toul 0 10): ,
(13)
C.W. Johnson et al./Computer Physics Communications 120 (1999) 155161
158
/ \ Represent the ~,th eigenvector by I x ~ ] . The RPA matrix is expanded in the eigenvectors \ Y~ J
B
=
09~
® (x~'
y~
y~)
®(y~,
x~
x~)
as
.
(14)
Unity (1) has a similar representation but with no 09. Therefore, the kth power of the RPA matrix is
(" B
k 09~
=
{() x~ y~
®(x~,
y~)(l)
~
() Y~ x~
®(y~,
x~)
}
.
(15)
P
Therefore, for odd k, A ½(q'q)
B
= ~1
=~
09ku(q,q).
x~ y~
®(x~,
09~ (09~1010)= = Mk,
w~ (q. (x + y ) ) 2 = ~ P
y~).
x2 (16)
P
by Eq. (2). In the algorithm described above, the nth iterate respects the oddk sum rules for k < 2n  1. If .A/" in Eq. (9) is zero or very small, there is the potential danger of numerical instability due to Eqs. (10), ( 1 1 ) blowing up; in the case of computing the linear response, however, one should stop and iterate no further because one has found all the states which Q can excite. We now illustrate the method and the rapid convergence due to the sum rules with a very simple model, a collective particlehole interaction fragmented by singleparticle energies. We consider states i = 1 ..... N with matrix elements Aij = Ei~ij "4 Kqiqj and Bij = K q i q j . Here e represents the energy spacing of the particlehole configurations, t¢ is the strength of the collective coupling to the field Q, and the components of the vector qi ~ i ( N  i) x r, where the r are Gaussian distributed random amplitudes, and normalize Iql 2  1. The factor i(Ni) weights the collective response towards the middle of the excitation spectrum. The parameter K should be positive for a repulsive collective interaction such as the Coulomb that generates plasmons. In Fig. 1 we show the strength function for such an RPA matrix in a space of 500 states, with parameter values given in the caption. The parameters were chosen to obtain moderate collectivity, with a strong but broadlyfragmented collective excitation distributed over the spectrum. Fig. 1 also displays the n = 3, 10, and 50 approximants to the strength function, where n is the number of Lanczos vectors IZi), or, equivalently, the number of multiplications with the RPA matrix. One sees that with a handful of states, one state closely approximates the collective excitation and the others distribute themselves over the remaining spectrum. A better way to see the convergence of the strength function is to plot its integral, 1(09) = ~ 0(09 09~)(09.lO10) 2. This is shown in Fig. 2 for n = 3 and 10. After 50 iterations the integral of the strength function is virtually indistinguishable from the exact solution. We mention that the algorithm does not explicitly preserve the total strength M0. If there were no correlations in the ground state, that is, if the vectors Y/all vanished, then the total strength would be Iql 2  1. The nontrivial deviations from 1 in our examples are related to the amount of correlations in the ground state. This is illustrated in Fig. 3, which is the integrated strength function for a model identical to that in Figs. 1, 2 except that the collective interaction is attractive rather than repulsive. Here the total strength is about 3.7, i.e. quite different from 1. Fig. 3 illustrates how with n = 3 and 10 the total approximate strength converges rapidly to the exact value. (In the repulsive model of Fig. 2 the total strength had already converged by n = 3.) Although we cannot prove this rapid convergence in all cases, it seems likely in light of the strong constraints imposed by the oddk sum rules. We anticipate that the algorithm will be particularly useful in problems which require a singleparticle dimensionality of the order of tens or hundreds of thousands, but which allow a sparse matrix approximation
C.W. Johnson et aL /Computer Physics Communications 120 (1999) 155161
0.40  
T
T
N=3
n
1
159
0.20
N=IO
0.30
0.20
0.10
0.10
s: 0.00
0.00 =
=
0.015
0.010
0.05
0.005
0.000
0.0
20.0
40.0
0.0
20.0
0.00
40.0
excitation energy Fig. 1. Strength function for the model described in the text with x = 10 and e = 0.1 (in arbitrary units) for 500 states, and the Lanczos approximants for 5, 10, and 50 Lanczos vectors. The scales for the abscissae are different because the strength is fragmented over a different number of states.
0.8 .... N=3      N=10 exact (N=500) 

/
0.6
== 0.4 = = ._=
,
I i
i
i t"
0.2
0"00.0
...... i; 10.0
20.0 30.0 excitation energy
40.0
50.0
Fig. 2. Integrated strength function for the model described in Fig. 1. For 50 Lanczos vectors the integrated strength is virtually indistinguishable on this graph from the full calculation.
C. W. Johnson et al. / Computer Physics Communications 120 (1999) 155161
160
4.0
'
T
,
•
,
•
,
:................. r . . . . . . . . . 1
1
'  . . . .
............
3.0
== "~ 2.0
.... ............  
N=3 N=5 N=10 exact
1.0
i
0"%.0
10.0
i
i
20.0 30.0 excitation energy
I
40.0
50.0
Fig. 3. The same as Fig. 2, except with the collective interaction is attractive, K =  1 0 . Notice that the total strength is not constrained, as
described in the text, but has converged by the 10th iteration.
for the Hamiltonian, such as the local density approximation. This applies to molecular and condensed matter physics modeled with the KohnSham equations, and to nuclear physics for excitations in deformed nuclei [ 2]. With the LDA Hamiltonian, an efficient particlehole representation can be constructed from the orbital representation of holes and the coordinatespace representation of particles [22]. The computational difficulty for the basic matrixvector multiplication then scales as the number of particles Np and the dimensionality of the singleparticle space M as MN2p ,,~ N3p. Only a fixed number of these operations, of the order of ten, are needed to obtain the strength function to the accuracy of the fundamental mean field approximation. Thus the overall scaling of the method is O(N3).
Acknowledgements This study arose in the program at the Institute for Nuclear Theory, "Numerical methods for strongly interacting quantum systems", and we wish to thank J. Carlson and R. Wiringa for providing that forum. G.B. also thanks K. Yabana for many discussions. Financial support was provided by the INT under Department of Energy Grant FG0690ER40561 and by Department of Energy Grant DEFG0296ER40985.
References [ 1] E Ring, P. Schuck, The Nuclear ManyBody Problem (Springer, New York, 1980). [2] E E Bortignon, R.A. Broglia, Nucl. Phys. A 371 (1981) 405; G.E Bertsch, E E Bortignon, R.A. Broglia, Rev. Mod. Phys. 55 (1983) 287; E Ring, L. M. Robledo, J.L. Egido, M. Faber, Nucl. Phys. A 41 (1984) 261; J.L. Egido, H.A. Weidenmiiller, Phys. Rev. C 39 (1989) 2398. [3] K. Yabana, G.E Bertsch, Phys. Rev. B 54 (1996) 4484. [4] A. Zangwill, E Soven, Phys. Rev. A 21 (1980) 1561. [5] A. Rubio et al., Phys. Rev. Lett. 77 (1996) 247. [6] C. Yannouleas et al., J. Phys. B 27 (1994) L642. [7] Y. Luo et al., J. Phys. Chem. 98 (1994) 7782.
C.W. Johnson et al./Computer Physics Communications 120 (1999) 155161
[8] [91 [ 10] [ 11 ] [12] [13] [14] [15]
[16] [17] [18] [19]
[20]
[21] [22]
161
C. Jamorski et al, J. Chem. Phys. 104 (1996) 5134. X. Blase et al., Phys. Rev. B 52 (1995) R225. A.A. Quong, A.G. Eguiluz, Phys. Rev. Lett. 70 (1993) 3955. J.K. Cullum, R.A. Willoughby, Lanczos Algorithms for Large Symmetric Eigenvalue Computations (Birkhauser, Basel, 1985). R. Haydock, Solid State Phys. 35 (1980) 215. R.R. Whitehead et al., Adv. Nucl. Phys. 9 (1977) 123. C.C. Paige, J. Inst. Math. Applic. 10 (1972) 373. R.R. Whitehead, A. Watt, J. Phys. G 4 (1978) 835; R.R. Whitehead, in: Theory and Applications of Moment Methods in ManyFermion Systems, B.J. Dalton et al., eds. (Plenum Press, New York, 1980) p. 235. W. Yang, W.H. Miller, J. Chem. Phys. 91 (1989) 3504. T.J. Park, J.C. Light, J. Chem. Phys. 85 (1986) 5870. J.S. Dehesa, A. Zarzo, Europhys. Lett. 8 (1989) 589. R. Haydock, J. Phys. A 7 (1974) 2120; H.D. Meyer, S. Pal, J. Chem. Phys. 91 (1989) 6195; J. Engel, W.C. Haxton, P. Vogel, Phys. Rev. C 46 (1992) R2153. E. Caurier, A. Poves, A.P. Zuker, Phys. Lett. B 252 13 (1990); W.C. Haxton, C.W. Johnson, Phys. Rev. Lett. 65 (1990) 1325; E. Caurier, A. Poves, A.P. Zuker, Phys. Rev. Lett. 74 (1995) 1517. J. Olsen et al., J. Comp. Phys. 74 (1988) 265. J.R. Chelikowsky et al., Phys. Rev. B 50 (1994) 11355.