Computational approaches of relativistic models in quantum chemistry

.


Chapter 1. Introduction 1 QED and relativistic models in Quantum Chemistry
It is now well known, following many experimental and theoretical results, that the use of ab initio relativistic calculations are mandatory if one is to obtain an accurate description of heavy atoms and ions. This is true whether one is considering highly charged ions, inner shells of neutral or quasi neutral atoms or outer shells of very heavy atoms.
From a physics point of view, the natural formalism to treat such a system is Quantum Electrodynamics (QED), the prototype of eld theories. For recent reviews of di erent aspects of QED in few electron ions see, e.g., 16,42,2]. Yet a direct calculation using only QED is impractical for atoms with more than one electron because of the complexity of the calculation. This is due to the slow rate of convergence of the so-called Ladder approximation (1=Z), that in non-relativistic theory amounts to a perturbation expansion using the electron-electron interaction as a perturbation. The only known method to do an accurate calculation is to attempt to treat to all orders the electron-electron interaction, and reserve QED for radiative corrections (interaction of the electron with its own radiation eld, creation of virtual electron-positron pairs). The use of a naive approach however, taking a non-relativistic Hamiltonian and replacing one-electron Schr odinger Hamiltonian by Dirac Hamiltonian fails. This approach does not take into account one of the two main features of relativity: the possibility of particle creation, and leads to severe problems as noted already in Ref. 3] and studied in 52]. This theory, for example, does not preserve charge conservation in intermediate states and leads to divergence already in the second-order of perturbation expansion. The only way to derive a proper relativistic many-electron Hamiltonian is to start from QED. The Hamiltonian of a N electron system can be written formally H = H 0 Ne ? ; 0e + ] + H 1 (N + 1)e ? ; 1e + ] + H 2 (N + 2)e ? ; 2e + ] + (1.1) Keeping only the rst term, the so-called \no-pair" Hamiltonian reads h D (r i ) + X i<j U ij ; (1.2) where (in atomic units) h D (r i ) = c p + mc 2  is the electron-electron interaction of order 1 in = 1=c 1=137, the ne structure constant. This expression is in Coulomb gauge, and is derived directly from QED. Here r ij = jr i ? r j j is the inter-electronic distance, ! ij is the energy of the photon exchanged between the electron i and j, which usually reduces to i ? j where the i are the one-electron energies in the problem under consideration (for exemple diagonal Lagrange multipliers in the case of Dirac-Fock). Note that in (1.6) gradient operators act only on the r ij and not on the following wave functions. The presence of the ! ij in this expression originates from the multi-time nature of the relativistic problem due to the niteness of the speed of light. From this interaction, one can deduce the Breit operator, that contains retardation only to second order in 1=c, in which the ! ij can be eliminated by use of commutation relations between r and the one-particle Dirac Hamiltonian. This operator can then be readily used in the evaluation of correlation, while the higher-order in 1=c in the interaction (1.6) can only be evaluated perturbatively.
Finding bound states of (1.2) is di cult and requires approximations. The different methods of solution are inspired from the non-relativistic problem. The three main categories of methods are the Relativistic Many-Body perturbation theory (RMBPT, see, e.g., 41] for the non-relativistic case), the Relativistic Random Phase Approximation (RRPA, see, e.g., 33]), which has been heavily used for evaluation of photoionization cross-sections, and Multicon guration Dirac-Fock (MCDF). The RMBPT method requires the use of basis sets to sum over intermediate states. The MCDF method is a variational method.
2 Relativistic Many-Body perturbation theory and RRPA In its most general version, the RMBPT method starts from a multidimentional model space and uses Rayleigh-Schr odinger perturbation theory. The concept of model space is mandatory if there are several levels of quasidegenerate energy as in the ground state of Be-like ions (1s 2 2s 2 1 S 0 and 1s 2 2p 2 1 S 0 are very close in energy, leading to very strong intra-shell correlation). In that case one gets would get very bad convergence of the perturbation expansion, because of the near-zero energy denominators, if building the perturbation theory on a single level.
where and P 0 is the projector on P, de ned by We de ne the perturbation potential by We also de ne Q 0 = 1?P 0 as the projector on the orthogonal space Q. We now de ne the wave operator, which build the exact solution of the Hamiltonian equation (2.1) from the 0 = 0 ; (2.5) so that H T = E ; (2.6) with the property P 0 P 0 = P 0 (2.7) The exact eigenenergies can be obtained by the application of the Model-space wave functions on the e ective Hamiltonian H e = P 0 H T P 0 = P 0 H 0 P 0 + P 0 V P 0 = H N 0 + P 0 V P 0 ; (2.8) using Eqs. (2.2) and (2.7), P 2 0 = P 0 and the fact that H N 0 and P 0 commute. This operator, acting on the unperturbed wave functions give the exact eigenenergies: H e 0 = E 0 (2.9) The wave operator obeys the generalized Bloch equation  (1) ; H 0 i P 0 = Q 0 V P 0 (2.12) h (2) ; H 0 i P 0 = Q 0 V (1) P 0 ? (1) P 0 V P 0 (2.13) The RRPA method is based on the solution of the Hamiltonian (1.2) subjected to a time-dependant perturbation (like a classical electromagnetic radiation of known frequency). This time-dependant Dirac-Fock equation is solved over a set of solutions of the unperturbed problem, leading to a set of time-dependant mixing coe cients in the usual fashion of time-dependant perturbation theory. The phases of those coe cients are approximated (leading to the name \Random Phase"), leading to di erential equations very similar to the Dirac-Fock ones. This method include to all orders some classes of correlation contribution that can be easily also evaluated in the framework of RMBPT. It is mostly used for the ground state of atoms and ions to study photoionization. It is more di cult to use for excited states.
This paper is mostly devoted to the MCDF method for atoms and molecules, and to preliminary results for the linear Dirac operator.

The MCDF wave function
We rst start by describing shortly the formalism used to build the Dirac-Fock solutions for a spherically-symmetric system like an isolated atom.
If we de ne the angular momentum operators L = r^p, J = L + 2 , the parity as P, then the total wave function is expressed in term of conguration state functions (CSF) as antisymmetric products of one-electron wave functions so that they are eigenvalues of the parity , the total angular momentum J and its projection M. The label stands for all other values (angular momentum recoupling scheme, seniority numbers, ...) necessary to de ne unambiguously the CSF. For a N-electron system, a CSF is thus a linear combination of Slater determinants: ; (3.1) all of them with the same and M values while the d i 's are determined by the requirement that the CSF is an eigenstate of J 2 .
The total MCDF wave function is constructed as a superposition of CSF's, i.e.
( JM) = NCF X =1 c j JM > ; (3.2) where NCF is the number of con gurations and the c are called the con gurations mixing coe cients. creates important di culties when trying to nd its eigenvalues. The so-called variational collapse is indeed related to this unboundedness property. On the other hand, nite dimensional approximations to this problem may lead to nding spurious solutions: some eigenvalues of the nite dimensional problem do not approach the eigenvalues of the Dirac operator and destroy the monotonicity of the approximated eigenvalues with respect to the basis dimension. These problems seem to be much more acute in molecular than in atomic computations, but they are already present in one-electron systems. In this section we address this di culty for one-electron systems by describing various methods used to deal with this problem. Well-behaved approximation methods should also provide good nonrelativistic limits, that is, variational problems whose eigenvalues and eigenfunctions converge well to those of the corresponding nonrelativistic Schr odinger Hamiltonian.
A way often used to nd good numerical approximations of eigenvalues of an operator A consists in projecting the eigenvalue equation over a well chosen nite dimensional space X N of dimension N, in order to nd an approximation ( N ; x N ) satisfying A N x N = N x N ; (4.3) such that ( N ; x N ) converges to ( ; x) as N ! 1. Then one looks for the eigenvalues of the N N matrix A N and these eigenvalues will converge either to eigenvalues of A or to points in the essential spectrum of A. As N increases, the limit set of the eigenvalues of A N is the spectrum of A.
The di culty with the Dirac operator is that for most physically interesting potentials V , the spectrum of H 0 + V is made of its essential spectrum (?1; ?mc 2 ] mc 2 ; +1) and a discrete set of eigenvalues lying in the gap (?mc 2 ; mc 2 ). Hence, the choice of the nite dimensional space, or equivalently, of the nite basis set, is fundamental if we want to ensure that for some N large, the eigenvalues of (H 0 + V ) N , or at least some of them, will be approximations of the eigenvalues of H 0 + V in the gap (?mc 2 ; mc 2 ). The question of how to choose a good basis set has been addressed in many papers, among which 11,12,22,21,32,35,38], that we will describe with further details in Sections 5 and 6 below. In particular, Section 6 is devoted to the description of numerical techniques based either on discretization or on B-Splines, and shows that with appropriate boundary conditions one can avoid the variational collapse.
When the operator A is bounded from below, it is often possible to characterize its spectrum by variational methods, for instance by looking for critical values of the Rayleigh quotient Q(x) := (Ax; x) (x; x) (4.4) over the domain of A. More concretely, when A is bounded from below, under appropriate assumptions, its ground state energy can be found by minimizing the above Rayleigh quotient. However, this cannot be done directly in the context of the Dirac operator, since it is unbounded from below (and also from above). A large number of works have been devoted to the variational resolution of this problem in view of the Dirac operator. Most of them use the approximation of an e ective Hamiltonian which is bounded from below. The idea hidden behind this kind of techniques is that there is no explicit way of diagonalizing the Dirac Hamiltonian H + V , but this can be done at an abstract level. The diagonalized operator is then approximated via a nite expansion or an iterative procedure. These methods are therefore perturbative and contain an approximation at the operator level. They will be referred to as perturbation theories and e ective Hamiltonian methods and will be described below (see 11,35,21,32] and Section 7 for more details).
Other variational techniques are based on a correspondence between the eigenvalues of A and those of T(A), for some operator function T, like the inverse function Tx = x ?1 (see 29]) or the function Tx = x 2 (see 56,1]). Finally, some authors solve the variational problem in a subspace of the domain in which the operator is bounded from below and 'avoids' the negative continuum. Section y8 will be devoted to these more direct variational approaches, based on either linear or nonlinear constraints.
Before going into the details of the computational methods, let us start with some notations and preliminary considerations. For any with values in C I 4 , if we write = ' , with '; taking values in C I 2 , then the eigenvalue equation 8 is equivalent to the following system: 8 > < > : R = ( ? mc 2 ? V ) ' ; R ' = ( + mc 2 ? V ) ; (4.6) with R=i c (~ :r)= P 3 j=1 i c j @ @x j . Here j , j = 1; 2; 3, are the Pauli matrices.
As long as +mc 2 ?V 6 =0, the system (4.6) can be written as where g = + 2mc 2 58]). For this reason, but also because the principal part of the second order operator in (4.7) is semibounded for not too large potentials V , the partitioned equation has been extensively studied for nding eigenvalues of linear Dirac operators.
To end these preliminary considerations on linear Dirac equations, note that in the case of rotationally invariant potentials, the solutions can be put in the form The dependence on the angular coordinates is contained in the 2-spinors m ( ; '), which are eigenfunctions of the angular momentum operators J, its third component J z (with eigenvalues j(j + 1) and m respectively) and of parity. On the other hand, the radial dependence is contained in the functions f and g which are called the upper and lower radial components of .

Finite basis set approaches
The choice of nite dimensional spaces is essential for the discretization of the operator and the approximation of its eigenvalues. The presence of the negative continuum makes this task di cult in the case of the Dirac operator.
The basic criterium to decide whether a particular space, or a generating basis set, is good, is to check that the approximated eigenvalues found are either negative and lying in the negative continuum or positive. In this case, if they are below the positive continuum, they are approximations of the discrete exact (positive) eigenvalues. Many attempts to construct nite basis sets can be found in the litterature.
In 11], Drake and Goldman introduced the so-called Slater type orbitals (STO): with a particular choice of and which depends on and V . They showed numerical evidence that such a nite basis satis es the above properties in the case of hydrogen-like atoms. Note that the STOs exhibit the same behavior near 0 and at 1 as the exact eigenfunctions. The properties of STO basis sets are made more explicit in 21], where STO basis sets are replaced by orthonormal sets of Laguerre polynomials. The main drawback of this approach is that some of the eigenvalues of the approximated matrix are spurious roots which do not approximate any of the exact eigenvalues.
Another way to construct basis sets with good properties consists in imposing the so-called kinetic balance condition relating the upper and lower components of the functions in the basis set. See for instance 35].
Other types of basis sets proposed in the litterature include those generated by B-splines (see 32]), which have very good properties since, in this approach, the matrices are very sparse: only a nite number (depending on the degree of the splines) of diagonal lines are nonzero. This kind of basis sets has been widely used in atomic and molecular computations (see Sec. 6).
The choice of a good basis set can be quite e ective in some computations, but as it appears clearly in the litterature that we quote, there is very often a risk of nding spurious roots or of variational collapse. In the next subsection we give some more precise exemples of how to use particular basis sets in the context of Dirac operators.

Numerical basis sets
This section is devoted to the special case of basis sets whose elements are computed numerically.

Discretisation method
The G oteborg group has developed an e cient technique to obtain basis sets for the Dirac equation 47]. The Dirac equation is discretized and solved on a grid. The atom is placed in a spherical box, large enough not to disturb the bound state wave function considered. The method provides a nite number of orbitals which is complete over the discretized space 48], and resemble lattice gauge eld calculation 57]. The method enables to eliminate spurious states and preserves the Hermiticity of the discretized Hamiltonian. The appearance of spurious states in a discretized method, is traced back to the \fermion doubling", rst encoutered in gauge-eld lattice calculations 34]. On a lattice of dimension (D + 1) (D spatial and one time dimensions), an equation for a massless fermion will describe not one but 2 D ones if no precaution is taken 51].
As an example, let us consider a one-dimensional Dirac equation for a free fermion 0 B @ The derivatives are approximated over the latice points using where h is the space between adjacent lattice sites. Eliminating the large component in (6.1), one gets the following equation in which the left hand side is the kinetic energy operator p 2 x =2m acting on f at the latice point i. Yet this second order derivative does not connect even and odd lattice sites. The highest energy solution over the lattice is the one changing sign at each site so that Using the expression (6.3) acting on this solution gives the same results as if it had no nodes. A high-energy eigenvector thus appears as a spurious state in the low energy part of the spectrum. For low-order derivative two equivalent ways can be used 51,47]. One is to use forward derivatives for f and backward derivatives for g, h : (6.5) The other consists in de ning the large and small components on alternating sites on the lattice.
with h being the separation between g i?2 and g i . In this case the derivative is expressed as with i = 2n ? 1, j = 2n, n = 1; 2; : : : ; N. These methods reduce to the same This six-point formula combined with (6.6) provides a spurious-state-free solution, while using the same lattice for f and g and a forward-backward derivative scheme does not work.
In the spherical case, one needs to use a logarithmic lattice to get a good description of the wave function. The Hermiticity of the Hamiltonian must be preserved by doing the variable change y(r) ! 1 p r y(x) ; x = log(r) : The 1 C A : (6.10) Since the large and small component are de ned on di erent lattices, one needs interpolation formulas to express f(x)= p r and g(x)= p r in the term.
The discretization nally provides a 2N 2N symmetric eigenvalue problem Z(2 + 1) r +3=2 ;(6.14) with = q 2 ? (Z ) 2 , and the equivalent expression for g. To avoid nonlinear terms in the eigensystem (6.11), only the contribution independent of has been kept. This is a good approximation for bound states for which mc 2 .
using the notations of (4.9) (note that in this representation the gap lies between ?2mc 2 and 0), to which suitable boundary conditions are added through From the point of view of the variational principle, is a Lagrange multiplier introduced to ensure that the solutions of the Dirac equation are normalized.
The boundary constraint (6.16) is designed to avoid a hard boundary at the box radius R, following the idea behind the MIT bag model for quark connement, and provides P (R) = Q (R). Forcing P (R) = Q (R) = 0 would amount to introduce an in nite potential at the boundary and possibly leads to the Klein paradox. Other choices of boundary conditions are possible. This particular choice avoids the appearance of spurious solutions. Expanding the radial wave function as P (r) = n X i=1 p i B i;k (r) ; Q (r) = n X i=1 q i B i;k (r) ; (6.17) the variational principle reduces to d(S + S 0 ) dp i = 0 ; d(S + S 0 ) dq i = 0 ; i = 1; 2; : : : ; n : (6.18) This leads to a 2n 2n symmetric, generalized eigenvalue equation Av = Bv ; (6.19) where v = (p 1 ; p 2 ; : : : ; p n ; q 1 ; q 2 ; : : : ; q n ), Diagonalization of (6.20) provides 2n eigenvalues and eigenfunctions, n of which have energies below ?2mc 2 , a few correspond to bound states (typically 5 to 6 for k = 7 to 9) and the rest belongs to the positive energy continuum.

Perturbation theory and e ective Hamiltonians
An alternative way to nd the eigenvalues of the unbounded relativistic operator H consists in looking for a so-called e ective Hamiltonian H e , which is semi-bounded, such that both Hamiltonians have common eigenvalues on an interval above the negative continuous spectrum. Such a Hamiltonian H e cannot usually be found in an explicit way, but can be viewed as the limit of an iterative procedure. This leads to families of Hamiltonians which approach the e ective Hamiltonian and yield approximated eigenvalues for H.
One of the most popular procedure in this direction is due to Foldy and ? + mc 2 ) are bounded from below (resp. above) and have correct nonrelativistic limits. Although this procedure looks very promising, the problem is that is unknown in closed form, and so there is no way of diagonalizing H 0 + V in an explicit way. However, approximations of , and therefore of H FW , can be constructed either by writing a formal series expansion for H FW in the perturbation parameter c ?2 : c ?2k H 2k ; (7.2) and cutting it at level k 0, or by approaching it by an iterative procedure.
In general one identi es the e ective Hamiltonian H e as a solution to a nonlinear equation H e = f(H e ) , which can be solved approximately in an iterative way. By instance, one can produce an equation like the above one by \eliminating" the lower component of the spinor as in (4.7), that is, by partitioning.
Many proposals of e ective Hamiltonians for the Dirac operator can be found in the litterature. Some are Hermitian, some are not, some act on 4 component spinors, others on 2-spinors. A good review about various approaches to this problem and the corresponding di culties has been written by W. Kutzelnigg 37] (see also 36,45,46]). An important di culty arising in this context is that most of the proposed e ective Hamiltonians are quite nice when the potential V is regular, but in the case of the Coulomb potential they contain very singular terms, which are not even well de ned near the nucleus. These serious singularities are avoided by a method used by Chang, P elissier and j( r) 'j 2 1 ? V + + (1 + V ) j'j 2 ! dx (8.2) and the corresponding eigenfunction is the spinor function = ' ?i ( r) ' Note that the idea to build a semibounded energy functional had already been introduced by Bayliss and Peel 1] in another context. It is closely related to previous works of Datta and Deviah 5], and Talman 54], where a particular min-max procedure for the Rayleigh quotient Q V is proposed without proof. We will not give here further details on these theoretical aspects (for tractable numerical applications, see below).
An alternative variational method has been proposed by Dolbeault, Esteban and S er e in 8]. It is based on rigorous results proving that for a very large class of potentials (including all those relevant in atomic models), the ground state of H 0 + V can be found by a minimization problem posed in a class of functions de ned by a nonlinear constraint. The main idea is to eliminate the lower component of the spinor and solve a minimization problem for the upper one. With the notations of the introduction, = ' is an eigenfunction of H 0 +V if and only if (4.7) takes place. The rst equation in (4.7) is an elliptic second order equation for the upper component ', while the second part of (4.7) gives the lower component as a function of ' and the eigenvalue . The dependence of H on = + mc 2 makes this problem nonlinear, since is still to be found, but the di culty of nding the unitary transformation in the Foldy-Wouthuysen approach is now replaced by a much simpler problem.
We may reformulate the question as follows. Let A( ) be the operator de ned by the quadratic form acting on 2-spinors: and consider its lowest eigenvalue, 1 ( ). Because of the monotonicity with respect to , there exists at most one for which 1 ( ) = 0. This is the ground state level.
An algorithm to numerically solve the above problem has been proposed in 10]. The idea consists in discretizing Eq. (8.2) in a nite dimensional space E n of dimension n of 2-spinor functions. The discretized version of (8.4) is A n ( ) x n x n = 0 ; (8.5) where x n 2 E n and A n ( ) is a -dependent n n matrix. If E n is generated by a basis set f' i ; : : : ' n g , the entries of the matrix A n ( ) are the numbers Z IR 3 0 @ ( r) ' i ; ( r) ' j 1 ? V + + (1 ? + V ) (' i ; ' j ) 1 A dx : The ground state energy will then be approached from above by the unique for which the rst eigenvalue of A n ( ) is zero. This method has been tested on a basis of Hermite polynomials (see 10] for some numerical results). More e cient computations have been made recently on radially symmetric con gurations with B-splines basis sets, involving very sparse matrices. Approximations from above of the excited levels can also be computed by requiring successively that the second, third,... eigenvalues of A n ( ) are equal to zero.
Chapter 3. The MCDF method for atoms 9 The Muticon guration Dirac-Fock (MCDF) method The MCDF equations are obtained from (1.2) by a variational principle. The energy functional is written E tot = < JMjH np j JM > < JMjj JM > : (9.1) A Hamiltonian matrix which provides the mixing coe cients by diagonalization is obtained from (9.1) with the help of @ @c E tot = 0 ; (9.2) and a set of integro-di erential equations for the radial wave functions P (r) and Q (r) is obtained from the functional derivatives 8 > < > : P (r) E tot = 0 ; Q (r) E tot = 0 : One assumes the orthogonality condition (restricted Dirac-Fock) Z 1 0 P A (r)P B (r) + Q A (r)Q B (r)] dr = A ; B n A ;n B ; (9.4) in order to make the angular calculations possible. Equation ( ?X P A (r) 1 C A (9.5) where V A is the sum of the nuclear potential and the direct Coulomb potential, while the exchange terms X P A and X Q A include all the two-electron interactions except for the direct Coulomb instantaneous repulsion. The constants A;B are Lagrange parameters used to enforce the orthogonality constraints of (9.4) and thus the summation over B runs only for orbitals with B = A . The exchange terms can be very large if the orbital A has a small e ective occupation (the exchange term is a sum of exchange potentials divided by the e ective occupation of the orbitals). This e ective occupation is the sum where q (A) is the number of electrons in the orbital A in the th con guration.
The numerical MCDF methods are based on a xed-point method, or to be precise on an iteration scheme which provides a self-consistent eld (SCF) state in a way very similar to the method which is used to solve the Hartree-Fock model. Initial wave functions must be chosen, e.g., hydrogenic wave functions, wave functions in a Thomas-Fermi potential or wave functions already optimized with a smaller set of con gurations. One then builds the Hamiltonian matrix (9.2) and obtains the mixing coe cients. Those coe cients and the initial wave functions enter the direct and exchange potential in (9.5), which become normal di erential equations, and are solved numerically for each orbital. A new set of potential terms is then evaluated until all the wave functions are stable to a given accuracy ( 10 ?2 in the rst cycle of diagonalization to 10 ?6 at the last cycle, at the point where the largest variation occurs). A new Hamiltonian matrix is then built and new mixing coe cients are calculated. This process is repeated until convergence is reached. As it is a highly nonlinear process, this can be very tricky, and trial and error on the initial conditions is often required when many con guration and correlation orbitals (i.e. orbitals with very small e ective occupations) are involved. All those calculations are done using direct numerical solutions of the MCDF differential equations (9.5), which has the advantage of providing very accurate results with relatively limited set of con gurations, while MCDF methods using basis set require orders of magnitude more con gurations to achieve similar accuracies.
Explicit expressions for V A , X P A and X Q A can be found in 23,24,6]. All potentials can be expressed in term of the functions Z k i;j (x) = 1 x k Z x 0 dr ij (r) r k ; (9.7) Y k i;j (x) = 1 x k Z x 0 dr ij (r) r k + x k+1 Z 1 x dr ij (r) r k+1 ; (9.8) where ij (r) = P i (r)P j (r)+Q i (r)Q j (r) for the Coulomb part of the interaction, to which are added terms with ij (r) = P i (r)Q j (r) or ij (r) = Q i (r)P j (r) when Breit retardation is included in the self-consistent eld process. These potential terms can be obtained very e ciently numerically by solving a second-order di erential equation (Poisson equation), as a set of two rst-order di erential equations, with the predictor-corrector method prensented in Sec. 10. 10 Numerical solution of the inhomogeneous Dirac-Fock radial equations In order to increase the numerical stability, the direct numerical computation of (9.5) is done by shooting techniques. First one chooses a change of variables to make the method more e cient because bound orbitals exhibit a rapid variation near the origin and exponential decay at large distances. One can choose either t = r 0 log(r) or t = r 0 log(r) + b r : (10.1) The rst choice leads to a pure exponential grid, while the second leads to an exponential grid at short distances and to a linear grid at in nity, and is better suited to represent, e.g., Rydberg states. One then takes a linear grid in the new variable t, t n = nh with h ranging from 0.02 to 0.05. In order to provide the few values needed to start the numerical integration at r = 0, and to have accurate integrals (for evaluation of the norm for exemple) the wave function is represented by its series expansion at the origin, which is of the form 8 > < > : P (r) = r (p 0 + p 1 r + : : : ) ; Q (r) = r (q 0 + q 1 r + : : : ) ; (10.2) where = q 2 ? (Z ) 2 if V N (r) = ?Z=r is a pure Coulomb potential and = j j if V N (r) represents the potential of a nite charge distribution. In this case if > 0, p 0 = p 2 = : : : = 0 and q 1 = q 3 = : : : = 0, and if < 0, p 1 = p 3 = : : : = 0, q 0 = q 2 = : : : = 0.
Predictor-Corrector Methods. In the case of the atomic problem, the use of fancy techniques like adaptative grids is not recommended, as it is much more e cient to tabulate all wave functions over the same grid, particularly if other properties like transition probabilities are calculated as well. One then uses well proven di erential equation solving techniques like predictor-corrector methods and nite di erence schemes. The expansion (10.2) is substituted into the di erential equation (9.5) to obtain the coe cients p i and q i , for i > 0 . These coe cients are used to generate values for the wave function at the few rst n points of the grid, with an arbitrary value of p 0 . Then the value of the function at the next grid point is obtained using the di erential equation solver. At in nity the same procedure is used. An exponential approximation of the wave function is made, and the same di erential equation solver is used downward to some matching point r m , usually chosen close to the classical turning point in the potential V A (r). In the predictor-corrector technique, an approximate value of the function at the mesh point n + 1 is predicted from the known values at the preceding n points. This estimate is inserted in the di erential equation to obtain the derivative that in turn is used to correct the rst estimate, then the nal value may be taken as a linear combination of the predicted and corrected values to increase the accuracy. As an example we consider the ve points Adams' method that has been widely selected because of its stability properties 43]. The predicted, corrected and nal values are given respectively by: p n+1 = y n + (1901y 0 n ? 2774y 0 n?1 + 2616y 0 n?2 ? 1274y 0 n?3 + 251y 0 n?4 )=720 ; c n+1 = y n + (251p 0 n+1 + 646y 0 n ? 264y 0 n?1 + 106y 0 n?2 ? 19y 0 n?3 )=720 ; (10.3) y n+1 = (475c n+ + 27p n+1 )=502 ; where p 0 and y 0 stand for the derivatives with respect to the tabulation variable. The linear combination for the nal value is de ned as to cancel the term of order h 6 , h being the constant interval step of the mesh. In the above equations, y represents either the large or small component of the radial wave function.
Since one starts with a somewhat arbitrary energy and slope at the origin, the components of the wave function obtained by the preceding method are not continuous. A strategy must be devised to obtain the real eigenenergy and slope at the origin from the numerical solution. In the case of an homogeneous equation, one can simply make the large component continuous by multiplying the wave function by the ratio of the inward and outward values of the large component at the matching point and then change the energy until the small component is continuous, using the default in the norm. To rst order the correction to the eigenvalue is = cP(r m ) Q(r ? m ) ? Q(r + m )] R 1 0 P 2 (u) + Q 2 (u)] du ; (10.4) where Q(r m ) are the solutions from each side of the matching point. One then checks that the solution is the desired one by verifying that it has the right number of nodes.
In the inhomogeneous case such a strategy cannot work. In order to obtain a solution which is continuous everywhere, it is possible to proceed in the following way. One uses the well known fact that the solution of an inhomogeneous di erential equation can be written as the sum of a particular solution of the inhomogeneous equation and of the solution of the associated homogeneous equation (in the present case the equation obtained by neglecting the exchange potentials). Thus if P o and P i are respectively the outward and inward solutions for the large component, one obtains, with the same labels for the small  (10.5) where the subscripts I and H stand for the inhomogeneous and homogeneous solutions. The coe cients a and b can be obtained from the di erential equation. Obviously this continuous solution will not be normalized for an arbitrary value of the diagonal parameter A;A of Eq.(9.5). The default in the norm is then used to modify A;A until the proper eigenvalue is found. This method is very accurate but not very e cient since it requires to solve both the inhomogeneous and the homogeneous equations to obtain a continuous solution.
Finite Di erences Methods As seen above, the predictor-corrector method has some disadvantages. In the non-relativistic case the Numerov method associated with tail correction 20] provides directly a continuous approximation (the derivative remains discontinuous until the eigenvalue is found). We consider now alternative methods that easily allow to enforce the continuity of one of the two radial components. Let us de ne the solution at point n+1 as: y n+1 = y n + h(y 0 n + y 0 n+1 ) + n ; (10.6) where n is a di erence correction given, in terms of central di erences, by: n = ?1 12 3 y n+ 1 2 + 1 120 5 y n+ 1 2 ; (10.7) with: 3 y n+ 1 2 = y n+2 ? 3y n+1 + 3y n ? y n?1 ; 5 y n+ 1 2 = y n+3 ? 5y n+2 + 10y n+1 ? 10y n + 5y n?1 ? y n?2 : (10.8) Accurate solutions are required only when self-consistency is reached. Consequently, the di erence correction n can be obtained at each iteration from the wave functions of the previous iteration as it is done for the potential terms. One can then design computationally e cient schemes 7]. We de ne a n = 1 + h 2 r 0 n rn ; u n = P n + h 2 h r 0 n X Q n + r 0 n+1 X Q n+1 i ; b n = ?1 + h 2 r 0 n rn ; v n = P n + h 2 h r 0 n X P n + r 0 n+1 X P n+1 i ; ' n = h 2 n ? V n ] r 0 n ; n = h r 0 n + n ; (10.9) where r 0 stands for dr=dt (to take into account the fact that the tabulation variable t is a function of r) and X P(Q) = X P A (Q A ) + P B6 =A A;B P B (Q B ). All the functions of r are evaluated using wave functions obtained at the previous iteration. Then the system of algebraic equations: a n+1 P n+1 ? n+1 Q n+1 + b n P n ? n Q n = u n ; ' n+1 P n+1 ? b n+1 Q n+1 + ' n P n ? a n Q n = v n ; (10.10) determines P n+1 and Q n+1 if P n and Q n are known. For the outward integration, this system is solved step by step from near the origin to the matching  7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 5 (10.11) and the two column vectors (PQ) and (uv) de ned as: (P Q) = 2 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 4 Q m P m+1 Q m+1 P m+2 : P N?1 Q N?1 P N 3 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 5 : (10.12) As displayed in Eq. (10.11) each row of the matrix M has at most four non-zero elements. To solve this system of equations the matrix M is decomposed into the product of two triangular matrices M = LT in which L is a lower matrix with only three non-zero elements on each row and T an upper matrix with the same property. Introducing an intermediate vector (pq) it is possible to solve L(pq) = (uv) for m , m+1 , : : : , N and then T(PQ) = (pq) for N, N ?1, .., m. The last point of tabulation N is determined by the requirement that P N should be lower than a speci ed small value when assuming Q N = 0. Thus the number of tabulation points of each orbital is determined automatically during the self-consistency process. This elimination process produces, as written here, a large component P that is continuous everywhere. The discontinuity of the small component at the matching point r m can then be used to adjust the eigenvalue A;A . In practice this method works very well for occupied orbitals (i.e. orbitals with e ective occupations at the Dirac-Fock level q (A) o n, n integer larger or equal to 1). Yet it is not su ciently accurate for correlation orbitals and leads to convergence instability. A good strategy 6] is thus to use the accurate predictor-corrector method for the outward integration and the nite di erences method with the tail correction for the inward integration. However the accuracy of the inward integration is increased by computing directly the di erence correction (10.7) from the wave function being computed rather than from the one from the previous iteration.
Diagonal Lagrange multipliers One can use di erential techniques, when the obtention of the eigenenergy AA is di cult. Their evaluation proceeds as follows. One can obtain the rst order variation of the large component P with respect to a change AA of one of the o -diagonal Lagrange multipliers by substituting the development P 0 AA + AA = P 0 AA + AA @P @ AA AA = 0 AA (10.13) (and the equivalent one for the small component Q) into the di erential equation (9.5). De ning p AA = @P @ AA ; q AA = @Q @ AA ; (10.14) leads ?P B (r) 1 C A (10.15) which is very similar to (9.5), with the replacement of X P A (r) (resp. X Q A ) by P B (r) (resp. Q B (r)). This system can be solved in p AA (r) and q AA (r) by the above techniques. With this solutions AA can be calculated in rst order from AA = 1 ? R 1 0 P A (r)P B (r) + Q A (r)Q B (r)] dr 2 R 1 0 p AA (r)P B (r) + q AA (r)Q B (r)] dr : (10.16) Note that such relations could be established to provide the change in the non-diagonal Lagrange multipliers AB as well, if one were to solve for several orbitals of identical symmetry simultaneously.

O -diagonal Lagrange multipliers
The self-consistent process outlined in Section 3 requires the evaluation of the o -diagonal Lagrange parameters to satisfy the orthonormality constraint (9.4). As in the non-relativistic case, the o -diagonal Lagrange multiplier between closed 1

Solution of the inhomogeneous Dirac-Fock equation over a basis set
It has been found however 30,31] that even the enhanced numerical techniques presented in Sec. 10 would not work for correlation orbitals with very small e ective occupation, particularly when the contribution of the Breit interaction is used in (9.5). This leads to point out that in the numerical MCDF calculations, the projection operators which should be used according to (1.5) are absent, as they have no explicit expression. A new method has been proposed that retains the advantages of the numerical MCDF. The idea is to expand P A , Q A , X P A and X Q A over a nite basis set, e.g., the one based on the B-spline calculated following the method of Sec. 6, using the full MCDF direct potential V A (r). Let us thus assume that one has a complete set of solutions The square of the norm of the solution of (9.5) is then easily obtained as One then can calculate the normalized solution of (9.5) if the o -diagonal Lagrange parameters are known, by solving N ( AA ) = 1 for AA . One can notice the interesting feature of (11.3) that the norm of the solution of the inhomogeneous equation (9.5) has a pole for each eigenenergy of the homogeneous equation. This method has the advantage over purely numerical techniques that by restricting the sums in (11.1) to positive energy eigenstates, one can explicitly implement projection operators, thus solving readily the \no-pair" Hamiltonian (1.2), rather than an ill-de ned equation. More details on this method and on the evaluation of the o -diagonal Lagrange multipliers can be found in 31].

Chapter 4. Numerical relativistic methods for molecules
Most of molecular methods that include relativistic corrections are based on the expansion of the molecular orbitals in terms of basis sets (most of the time taken to be Gaussian functions). We shall not review these methods here but refer the interested reader to a book to be published soon 49]. Let us just point out that the sometimes observed lack of convergence to upper bounds in the total energy (the so-called variational collapse) is not unambigously related to the Dirac negative energy continuum. Indeed this attractive explanation is unfortunately unable to explain the appearance of spurious solutions. Both the existence of spurious solutions and the lack of convergence to expected levels can be traced back to originate from poor basis sets and bad nite matrix representations of the operators (in particular the kinetic energy). For an extensive discussion see 15]. Numerical methods succesfully used are brie y sketched in the next two paragraphs.
12 Fully numerical two-dimensional method = (r 1 + r 2 )=R ; = (r 1 ? r 2 )=R ; (12.2) where r 1 and r 2 are the distances between the electron and each of the nucleus, R is the inter-nuclear distance. The third variable ' is the azimutal angle around the axis through the nuclei.
As usual for molecular calculations, the variational collapse is avoided by de ning the small component in terms of the large one 35]. Starting from the Dirac equation in a local potential V one possibility is to use: S = c :p L = h 2c 2 + E ? V i : (12.3) After this substitution, the large component is given as solution of a second order di erential equation that can be solved using well known relaxation methods 55].
For e ciency, the distribution of integration points must be chosen as to accumulate points where the functions are rapidly varying. It was found that the transformation, = arccosh( ) ; = arccosh( ) ; (12.4) which yields a quadratic distribution of points near the nuclei, is some kind of optimum to reduce the number of points needed to achieve a given accuracy. Then the derivatives of the Laplace operator are approximated by n-point nite di erences. In so doing, the di erential equations are replaced by a set of linear equations that can be written in a matrix form as (A ? ES) X = B ; (12.5) where the matrix A, that represents the direct part of the Fock operator, is diagonal dominant but has non-diagonal elements arising from the discretization of the Laplace operator. Here E is the energy eigenvalue, S is the overlap diagonal matrix and B a vector due to the exchange part of the Fock operator whose values change during iterations. Then the relaxation method can be viewed as an iterative method to nd the x i component of X such that (A ? ES) x i = b i ; (12.6) each iteration n being associated with a linear combination of the initial and nal estimate of x i at iteration n ? 1, i.e.
x initial n+1 i = (1 ? !)x initialn i + !x finaln i : (12.7) It was found that with overrelaxation (i.e. ! > 1), the method may be slow in convergence but it is quite stable. Applications of the method outlined above may be found in 53] and in references therein.

Numerical integrations with linear combinations of atomic orbitals
A widely used approximation in molecular calculations is to expand the molecular orbitals as a linear combination of atomic orbitals. If these atomic orbitals are chosen as the numerical solutions of some kind of Dirac-Fock atomic calculations, then small basis sets are su cient to achieve good accuracy. The main disavantage of this choice is that all multi-dimensional integrals have to be calculated numerically. This is compensated by two advantages: rst the kinetic energy contribution can be computed by a single integral using the atomic Dirac equations (thus avoiding numerical di erenciation), second, by including only positive energy atomic wave functions, no \variational collapse" will occur.
In this method, the molecular wave functions are expanded in terms of symmetry molecular orbitals as: = X c ; (13.1) while the symmetry molecular orbitals are taken to be linear combinations of atomic orbitals ' : = X i d i ' i : (13.2) The coe cients d i are given by the symmetry of the molecular orbital and are obtained from the irreducible representations of the double point groups.
Computing all necessary integrals (overlaps, matrix elements of the Dirac operator, the Coulomb interaction, etc...) the Dirac-Fock equations are reduced to a generalized matrix eigenvalue problem that determines both the eigenvalues and the c coe cients of equation ( where the weight function !(r i ) can be considered as an integration weight corresponding to a local volume per point. This function is also constrained to force the error momenta to vanish on the grid points following the work of Haselgrove 27]. Furthermore the set of the sampling points r i ] must be chosen to preserve the symmetries of the system under con guration (this is accomplished by taking a set of sampling points that includes all points Rr i , R standing for operations of the symmetry group). A full description of the DVM can be found in the references given above.