\documentstyle[pra,aps,twocolumn,psfig]{revtex} \newcommand{\bfa}{\mbox{\boldmath $\alpha$}} \newcommand{\bfs}{\mbox{\boldmath $\sigma$}} \newcommand{\bft}{\mbox{\boldmath $\bigtriangledown$}} \begin{document} \draft \title{Benchmark calculations for lanthanide atoms: calibration of ab initio and density functional methods} \author{ Wenjian Liu$^{1,2}$, and Michael Dolg$^{1}$} \address{$^{1}$ Max-Planck-Institut f\"ur Physik komplexer Systeme, N\"othnitzer Stra{\ss}e 38 D-01187 Dresden, Germany} \address{Email : lwj@mpipks-dresden.mpg.de} % \address{$^{2}$ State Key Laboratory of Rare Earth Materials Chemistry and Applications, College of Chemistry and Molecular Engineering, Peking University, Beijing 100871, P. R. China} % \maketitle \begin{abstract} Relativistic ab initio pseudopotential and fully relativistic density functional benchmark calculations have been carried out for the first to fourth ionization potentials as well as the 'df' charge transfer energies for the whole series of lanthanide atoms. It was found that the two approaches have essentially the same accuracy compared to experimental values. In addition, it is shown that the present (nonrelativistic) density functionals work fairly well within an otherwise relativistic framework even for the rather compact 4f shells, correcting previous opposite statements. \end{abstract} \pacs{ } \narrowtext \section{Introduction} The chemistry of lanthanide elements has received much attention in the past three decades \cite{REhand}. However, the complexity of the open shells of 4f, 5d, 6s and 6p poses a great challenge to theoretical work \cite{dolg3}, e.g., the $^{2S+1}$L$_{J}$ term of the 4f$^n$ subshell may have a spin S as large as 7/2 and an angular momentum L as large as 12. Even more extrem values may result from the coupling of the 4f$^n$ subshell to other partially occupied shells of s, p or d symmetry. Moreover, spin-orbit coupling leads to a large number of energetically adjacent electronic states \cite{ExptRE}. Therefore, the knowledge of the energy levels of free lanthanide atoms and ions is far from being complete. However, a detailed knowledge, at least of the low-lying states, is a necessary prerequisite for understanding the behaviour of lanthanide atoms in molecules or solids. Theoretical first-principles methods are presently at the edge of successfully dealing with such complicated systems containing lanthanides or actinides \cite{dolg3}. It was found that traditional ab initio approaches dealing with relativity at the all-electron Dirac-Coulomb-Breit level and including electron correlation effects by means of coupled-cluster (CC) or configuration interaction (CI) methods need h- or even i-functions in the one-particle basis sets to yield accurate results \cite{Kaldor}. However, such state-of-the-art studies are presently only feasible for atoms by exploiting their spherical symmetry, and to our knowledge, due to its implementation the method is currently applicable only to some special cases, i.e., closed-shell systems, one or two electrons outside a closed shell or one or two holes inside a closed shell. In order to be able to treat all lanthanide atoms and also to be able to deal with molecular systems, compromises have to be made with respect to the treatment of relativity and electron correlation. Two of such approximate schemes are considered in the present work. A very successful approach in relativistic quantum chemistry is the ab initio pseudopotential method \cite{Kutzelnigg}, where the explicit quantum chemical treatment is restricted to the valence electrons and relativistic effects are implicitly accounted for by a proper adjustment of free parameters in the valence model Hamiltonian. Although several sets of such potentials have been published for the lanthanide atoms \cite{dolg1,dolg2,Cundari}, no systematic calibration at the correlated level has been performed for atoms up to now. Another approach, which has gained extensive attention in quantum chemistry during the last decade, is density functional theory (DFT) \cite{DFTreview}. Although in principle the theory based on the works of Hohenberg and Kohn \cite{HK} and Kohn and Sham \cite{KS} is exact, in practice only approximate approaches are at hand. Since most of the density functionals used nowadays, e.g., local density approximation (LDA) and generalized gradient approximation (GGA), take the homogeneous electron gas as input and are therefore expected to work well only for slowly varying charge densities, it is reasonable to doubt their good performance for the rather compact 4f-shells of lanthanide atoms. Indeed, some previous relativistic DFT calculations were not satisfactorily successful in reproducing the term energies of lanthanide atoms, especially when the related states involve occupation changes in the 4f shells \cite{Wang1,Eschrig}. The present authors recently investigated a number of electronic states of Eu and Yb as well as their cations and found out that this failure is, at least partially, an artifact (e.g., basis set error or use of spherically averaged charge densities) \cite{BDFTCA,BDFJCP}. However, in order to establish the reliability of DFT methods for systems containing lanthanides a broader study seems to be needed. Since both relativistic ab initio pseudopotential as well as density functional methods will be the methods of choice for the treatment of systems containing heavy elements in the forseeable future, we decided to investigate their performance for lanthanides in detail. We studied the first to fourth ionization potentials (IP$_{1,2,3,4}$) as well as the 4f$^{n+1}$ 6s$^2$ - 4f$^n$ 5d$^1$ 6s$^2$ excitation energies ('df' charge transfer energies) of the whole series of lanthanide atoms (La to Lu). For the ab initio calculations we applied the energy-consistent quasi-relativistic pseudopotentials (QR-PP) of the groups in Stuttgart and Dresden \cite{dolg2}, whereas for the density functional calculations the newly developed Beijing four-component density functional program package (BDF) \cite{BDFTCA,BDFJCP} was used. We suggest that the DFT results presented here may serve as a benchmark for other DFT calculations using transformed or simplified relativistic Hamiltonians. Our paper is organized as follows. In section II we outline the applied ab initio pseudopotential and density functional methods. In section III we present our results and compare them to the available experimental data to display what accuracy the presently available ab initio and DFT approaches can actually achieve. Finally, in section IV we give our conclusions. \section{The methods} \subsection{QR-PP} The method of quasi-relativistic energy-consistent ab initio pseudopotentials was described in detail elsewhere \cite{dolg1,dolg2} and will be outlined here only briefly. The valence-only model Hamiltonian for an atom or ion with $n$ valence electrons is given as \begin{equation} {\cal H}_{v} = - \frac{1}{2} \sum_{i}^{n} \Delta_i + \sum_{i < j}^{n} \frac{1}{r_{ij}} + V_{av} + V_{so} \qquad . \end{equation} \noindent Here i and j are electron indices. $V_{av}$ denotes a spin-orbit averaged relativistic pseudopotential in a semilocal form \begin{equation} V_{av} = - \sum_i \frac{Q}{r_i} + \sum_i \sum_{l , k} A_{l k} exp(-a_{l k} r_i^2) P_l \qquad , \end{equation} \noindent where $P_l$ is the projection operator onto the Hilbert subspace of angular momentum $l$. The spin-orbit term $V_{so}$ may be written as \begin{equation} V_{so} = \sum_i \sum_{l>0,k} \frac{2}{2 l + 1} B_{l k} exp(-b_{l k} r_i^2) P_l {\bf l}_i {\bf s}_i P_l \qquad . \end{equation} The free parameters $A_{l k}$, $a_{l k}$, $B_{l k}$ and $b_{l k}$ are adjusted to reproduce the valence total energies of a multitude of low-lying electronic states of the neutral atom and its ions. The necessary reference data has been taken from relativistic all-electron calculations. In the present work accurate small-core pseudopotentials for Ce to Yb have been used, e.g., the 1s - 3d shells were included in the pseudopotential core, while the higher shells were treated explicitly. The orbitals were described by medium-sized one-particle basis sets, e.g., the exponents of a (12s10p8d8f) primitive set were optimized for the lowest state of the 4f$^{n+1}$ 6s$^2$ configuration of the neutral atom. The contraction coefficients of a [5s5p4d3f] set were derived from atomic natural orbitals (ANO) of a multi-reference configuration interaction ground state calculation keeping the 4s, 4p and 4d shells frozen. The generalized contraction scheme was applied. A (6g)/[4g] ANO correlation set was then derived in the same way starting from the most important exponents of the f set. Finally, two diffuse functions were added in all symmetries up to g, resulting in (14s12p10d10f8g)/[7s7p6d6f6g] basis sets. For the first to fourth ionization potential of all atoms considered here the basis set errors are less than 0.2 eV at the Hartree-Fock level. The corresponding pseudopotentials errors at the finite difference level are also typically 0.2 eV or less. All scalar-relativistic calculations were carried out with the MOLPRO ab initio program package \cite{MOLPRO}. The atomic orbitals were optimized in state-averaged complete active space multi-configuration self-consistent field calculations (CASSCF). Dynamic correlation was then accounted for by all single and double excitations from the CASSCF reference in averaged coupled-pair functional calculations (ACPF) \cite{ACPF}. The active space in the CASSCF comprised all open-shell orbitals (4f, 5d, 6s), whereas in the ACPF excitations were also allowed from the semi-core orbitals (5s, 5p). No excitations were allowed from the 4s, 4p and 4d shells both in the CASSCF and ACPF, however, the orbitals were optimized for each state. Spin-orbit coupling was taken into account by complete configuration interaction calculations within all open-shell orbitals (COSCI). The corresponding corrections derived from calculations with and without $V_{so}$ were then added to the scalar-relativistic ACPF results. Since spin-orbit contributions were found to amount to at most a few tenths of an electronvolt in the cases considered here, such an additive treatment appears to be justified. Modified versions of the finite-difference programs MCHF \cite{MCHF} and GRASP \cite{GRASP} were used. Due to the use of the state-averaging technique in calculations using MOLPRO and the exploitation of the spherical symmetry in MCHF and GRASP, all ab initio results of this work were obtained with eigenfunctions of the appropriate parity and angular momentum operators. \subsection{BDF} The BDF program package also has been described elsewhere\cite{BDFTCA,BDFJCP}. Briefly, the one-particle Dirac-Kohn-Sham equation (\ref{dks1}) based on the Dirac-Coulomb Hamiltonian under the so-called no-pair approximation are directly solved. % \begin{equation}\label{dks1} (c \bfa \cdot {\bf {p}} + (\beta - 1) c^{2} + V_{ext}({\bf {r}}) + V_{c}({\bf {r}}) + V_{xc}(\rho({\bf {r}})) \varphi_{j}({\bf {r}}) = \epsilon_{j} \varphi_{j}({\bf {r}})\\ \end{equation} % Here $\bf {p}$=-i $\bft$ is the usual momentum operator and c denotes the speed of light, 137.037 au. $\bfa$ and $\beta$ are the Dirac matrices\\ % \begin{equation} \bfa = \left( \begin{array}{cc} 0 & \bfs\\ \bfs & 0 \\ \end{array} \right)\hspace{2mm} \ \ \mbox{and } \ \ \beta= \left( \begin{array}{cc} I & 0\\ 0& -I \\ \end{array} \right) \end{equation} % where $\bfs$ represents the vector of the $2\times2$ Pauli spin matrices ($\sigma_{x}$, $\sigma_{y}$, $\sigma_{z}$) and I is the $2\times2$ unit matrix. The external, Coulomb and exchange-correlation potentials in eq. \ref{dks1} are, respectively, \begin{equation} V_{ext}({\bf {r}}) = - \sum_{A}\frac{Z_{A}}{|{\bf {R}_{A}}-{\bf {r}}|}\\ \end{equation} % \begin{equation} V_{c}({\bf {r}}) =\int \frac {\rho({\bf {r^{'}}})} {|{\bf {r}}-{\bf {r^{'}}}|} d{\bf {r^{'}}} \\ \end{equation} % \begin{equation} V_{xc}(\rho({\bf {r}}))=\frac{\delta E_{xc}(\rho({\bf {r}}))}{\delta \rho} \\ \end{equation} % The charge density reads \\ % \begin{equation} \rho({\bf {r}}) = \sum_{j}^{occ} n_{j}\varphi_{j}^{+}({\bf {r}})\varphi_{j}({\bf {r}})\\ \end{equation} % Since relativistic corrections, e.g., the Breit term, to Coulomb and exchange-correlation potentials have only a very limited influence on valence-electron excitation energies of atoms \cite{Kaldor,Eschrig}, they were not considered in the present calculations. Instead, self-interaction corrections (SIC) to the approximate density functionals are sinificant for the compact shells \cite{Eschrig}. The approximate forms for the exchange-correlation potential $V_{xc}(\rho({\bf {r}}))$ employed in this work are the Perdew-Wang formula \cite{PW92} within the local density approximation (LDA), a SIC term according to Stoll et al. \cite{SIC} and nonlocal exchange corrections (B88) according to Becke \cite{Beck88} as well as nonlocal correlation corrections (P86) according to Perdew \cite{Pd86}. \\ The atoms were treated in the same manner as molecules in the calculations by using the double point D$_{\infty h}^{*}$ group. The jj-coupling scheme was used and Kramer's degeneracy was adopted to carry out moment-polarized calculations for open shells. For the configurations considered here the highest possible moment polarization was always generated. Specifically the 4f shell was occupied as follows: electrons one to three occupy 4f$_{5/2}$ with moment up and electrons four to seven occupy 4f$_{7/2}$ with moment up; then, electrons eight to ten occupy 4f$_{5/2}$ with moment down and finally electrons 11 to 14 occupy 4f$_{7/2}$ with moment down. 5d$_{3/2}$ and 6s$_{1/2}$ were always occupied with moment up when occupied with a single electron. Keeping fixed the highest possible moment polarization, we then used fractional occupation numbers for all moment-polarized subshells with incomplete filling, e.g., for a 4f$^1$ configuration each of the three 4f$_{5/2,m_j}$ spinors with moment up was occupied by 1/3 electrons. A final remark appears to be in order here: although our program works in the jj-coupling scheme we have to account for the fact that the lanthanides are closer to the nonrelativistic LS-coupling scheme. Therefore, instead of filling first 4f$_{5/2}$ and afterwards 4f$_{7/2}$, we used the prescription given above, which also leads to lower total energies. The generalized Gauss-Laguerre quadrature \cite{Yang} and Lebedev quadrature \cite{Lebedev} were employed to calculate the radial and angular integrals, repectively. The numerical accuracy of total energies can be further improved to better than 0.01 eV by the generalized transition-state method \cite{TS}. The frozen-core approximation, i.e., $[$1s$^{2}-$4d$^{10}]$, was employed for all the calculations because relaxation of the 4s, 4p and 4d shells did not change the total energies larger than 0.01 eV \cite{BDFJCP}. Four-component numerical atomic spinors obtained by moment-restricted finite-difference atomic calculations were used for the cores, while the basis sets for the valence orbitals were combinations of the numerical atomic spinors and kinetically balanced double-zeta Slater-type functions (STF). Such basis sets result in errors less than 0.05 eV. \section{Results and discussion} Before we discuss our results in detail we want to emphasize that electron correlation effects turn out to be very important (cf. Fig. \ref{fig1}), e.g., they amount up to about 1 and 4.5 eV for IP$_{1,2}$ and IP$_{3,4}$, respectively. Therefore, uncorrelated Dirac-Hartree-Fock calculations can not be in quantitative agreement with experiment \cite{Matsuoka}. Spin-orbit coupling contributions are typically less than 0.2 eV and 0.5 eV for IP$_{1,2}$ and IP$_{3,4}$, respectively, indicating that for many purposes the use of scalar-relativistic Hamiltonians might be sufficiently accurate. We also want to address a critical point in the DFT studies which is relevant to the fine structure of a multiplet state. The currently existing approximate density functionals lead to unphysical splittings of levels which should be degenerate. These amount to 0.6 eV for the 4f$^{13}$ 6s$^2$ configuration of the Yb atom \cite{BDFJCP}. No remedy is currently at hand to avoid this unpleasant feature. However, in a previous study we found that in cases of a fixed 4f occupation number these splittings are transferable between different states and therefore, energy differences derived for the lowest levels of two states will be only slightly affected by the unphysical splittings due to an error compensation. In cases of a variation of the 4f occupation number the error compensation will certainly be less effective, if present at all. Besides the development of still more accurate density functionals to be used within a single-determinant framework, the extension to multi-determint wavefunctions might partially cure this defect \cite{Savin}. Since in the present work we are only interested in the lowest level of a configuration in order to study the general performance of DFT for ionization potentials and excitation energies, we used the average occupation scheme described above. The calculated ionization potentials IP$_1$ to IP$_4$ are listed in tables 1 to 4, respectively, while the 'df' charge transfer energies are given in table 5. Some additional results from previous studies are also included. The quality of the present work can be judged from the mean absolute errors (m.a.e.) given at the bottom of the tables, which were calculated with respect to the experimental values given by Martin et al. \cite{ExptRE}. Two points have to be taken into account. First, the ab initio and DFT m.a.e. refer to the atoms Ce to Yb and La to Lu, respectively, i.e., they are not defined for the same set of systems. Second, some experimental data is not available, e.g., for Pm, and some other values bear large error bars \cite{ExptRE}. Nevertheless, we think that the m.a.e. given in tables 1 - 5 still allow to roughly judge the performance of the ab initio and DFT approaches presented here. It can be seen from table 1 and 2 that both methods presented here reproduce IP$_1$ and IP$_2$ fairly well, i.e., the m.a.e. are below 0.20 eV for the ab initio and BDF (LDASIC) results. Nonlocal corrections (B88, P86 or BP) do not change the LDASIC results significantly. Forstreuter et al.'s \cite{Eschrig} DFT calculations based on the squared Dirac Hamiltonian (D$^{2}$) using both relativistically\cite{RVxLDA} and SIC \cite{PZSI} corrected local density functionals (RLDASIC) yielded IP$_{1}$ of similar quality to our DFT results. Nevertheless, it is also discernable that the differences for IP$_{1}$ between the BDF (LDASIC) and D$^{2}$ (RLDASIC) results amout to 0.5 eV for lighter elements, although they are in good agreement for heavier elements. However, the BDF results show more systematic errors when compared to experiment. They are are also closer to other DFT results taking relativity into account as a first-order perturbation \cite{Wang1} as well as fully relativistic coupled-cluster calculations\cite{Kaldor} (cf. table 1). The m.a.e. for IP$_{2}$ by D$^{2}$ (RLDASIC) \cite{Eschrig} amounts to 0.50 eV, i.e., it is a factor of 2 larger than the present values. We mention that in contrast to the ab initio method both BDF (cf. table 1 and 2) and D$^{2}$ (RLDASIC) \cite{private} calculations do not reproduce the experimental ground states for Ce/Ce$^{+}$ and Gd/Gd$^{+}$. This might be attributed to the fact that nondynamic correlation effects due to near-degenerate configurations are missing within the single-determinant formulation of DFT. A combination of multireference wavefunctions with DFT might be able to improve the results\cite{Savin}. The rather good performance of both ab initio and DFT approaches for IP$_{1}$ and IP$_{2}$ is mainly related to the fact that for almost all atoms (an exception is Ce) the 4f occupation remains unchanged. Harder tests are the third and fourth ionization potentials where the 4f occupation is changed by one or two electrons. For IP$_{3}$ the m.a.e. of the ab initio data is 0.58 eV, whereas the ones for the BDF results range from 0.45 eV to 0.60 eV at different levels of calculations. A much larger m.a.e. of 2.1 eV as well as larger fluctuations of 1 to 4 eV in the absolute deviations are found for the D$^2$ (RLDASIC) data \cite{Eschrig}. Since in the framework of density functional calculations the squared Dirac Hamiltonian used by these authors should be completely equivalent to the Dirac Hamiltonian used in our work, the differences have to be traced to other sources, except for the small differences between the different density functionals used in these two calculations. It has been found that roughly half of the error stems from the sole use of an incomplete STF basis set \cite{private}, whereas a combination of numerical functions and STFs was used in the present work. In addition, Forstreuter et al.\cite{Eschrig} used spherically averaged charge densities in contrast to the present polarized ones. For the fourth ionization potential the m.a.e. is 0.45 eV for the ab initio results and ranges from 0.42 eV to 0.61 eV for the BDF values. Of course, one has to keep in mind that the error bars in the experimental results amount to 0.7 eV (cf. table 4). Another less robust criterion than the m.a.e to judge the quality of the calculations is the largest relative error (l.r.e.), which is relevant to the amplitude of the fluctuations of the calculated results with respect to the experimental values. From tables 1 to 4 one can see that the largest relative errors are typically about 2 to 6 $\%$ for all the approaches included here, with an exception of 15.3 $\%$ in the D$^2$ (RLDASIC) IP$_3$ results\cite{Eschrig}. It is worthwhile to mention that the present BDF results for Yb, Lu and Pr are very close to those obtained by the fully relativistic all electron coupled-cluster calculations using very large uncontracted basis sets (cf. the footnotes in tables 1-4), whereas the present approach is computationally much cheaper. A real challenge for computational methods appears to be the calculation of 4f$^{n+1}$ 6s$^2$ - 4f$^n$ 5d$^1$ 6s$^2$ excitation energies ('df' charge transfer energies). Ab initio results for 4f$^n$ 5d$^1$ 6s$^2$ are quite difficult to get due to the too large active space resulting from an open d and f shell. Moreover, for Sm, Eu and Yb some states of 4f$^{n+1}$ 6s$^1$ 6p$^1$ are lower in energy than the lowest solutions for 4f$^n$ 5d$^1$ 6s$^2$ and cause root flipping problems in the CI. Convergence also could not be achieved for Pr, Nd and Pm where 4f$^n$ 5d$^2$ 6s$^1$ is nearly degenerate. In fact we were only able to perform the calculations for Ce and Gd with a complete active space: the errors in the excitation energies after correction for spin-orbit coupling are 0.37 eV and 0.86 eV, respectively. In both cases the configuration with the larger 4f occupation number is too high in energy, most likely reflecting the incomplete correlation treatment due to the neglect of higher than g functions in the basis sets. In order to obtain results for the atoms Tb to Tm the occupation of the 5d orbitals in the active space had to be restricted to one in the reference wavefunction. Compared to the lowest experimental levels the m.a.e. of the BDF results range from 0.90 eV to 1.04 eV, whereas 0.43 eV are obtained at the ab initio level. The DFT calculations systematically overestimate the excitation energies and the ab initio calculations underestimate them. We mention that the m.a.e. of nonrelativistic calculations is more than 2 eV \cite{MPrev}. Again, the remaining large deviations of relativistic DFT calculations might be accounted for by the missing nondynamic correlation effects which are expected to be larger for the 'df' charge transfer energies than for the ionization potentials because occupation changes occur to both f and d shells. \section{Conclusions} Benchmark calculations using both ab initio and DFT approaches have been performed for the whole series of lanthanide atoms. The results show that both approaches have essentially the same accuracy when compared to experimental data. Clearly, the current ab initio results might be systematically improved as soon as this is feasible from a computational point of view, e.g., by including higher-order angular momentum basis functions. The presently available approximate (nonrelativisitic) density functionals work fairly well in an otherwise relativistic framework for the rather compact 4f shells, correcting previous opposite statements by other authors. Moreover, a combination of multi-reference wavefunctions with DFT might even be able to further improve the performance of DFT in the open-shell systems studied here. Work along this line is underway in our laboratory. \section{Acknowledement} The authors thank H. Eschrig, J. Forstreuter and M. Richter for valuable discussions. \newpage $ $ \newpage \begin{thebibliography}{999} % \bibitem{REhand} {\it Handbook on the Physics and Chemistry of Rare Earths}, edited by K. A. Gschneidner, Jr. and L. Eyring (Elsevier, Amsterdam, 1978) % \bibitem{dolg3} M. Dolg and H. Stoll, in {\it Handbook on the Physics and Chemistry of Rare Earths}, edited by K. A. Gschneidner, Jr. and L. Eyring (Elsevier, Amsterdam, 1996), vol. 22, p. 607. % \bibitem{ExptRE} W. C. Martin, R. Zalubas, and L. Hagan, {\it Atomic Energy Levels-The Rare Earth Elements} (NSRDS-NBS 60, Washington, DC, 1978). % \bibitem{Kaldor} E. Eliav and U. Kaldor, Phys. Rev. A 52, 291 (1995). % \bibitem{Kutzelnigg} W. Kutzelnigg, Phys. Scr. 36, 416 (1987). % \bibitem{dolg1} M. Dolg, H. Stoll, A. Savin, and H. Preu{\ss}, Theor. Chim. Acta 75, 173 (1989). % \bibitem{dolg2} M. Dolg, H. Stoll, and H. Preu{\ss}, J. Chem. Phys. 90, 1730 (1989). % \bibitem{Cundari} T. R. Cundari and W. J. Stevens, J. Chem. Phys. 98, 5555 (1993). % \bibitem{DFTreview} {\it Modern Density Functional Theory, A Tool for Chemistry}, edited by J. M. Seminario and P. Politzer (Elsevier, Amsterdam, 1995). % \bibitem{HK} P. Hohenberg and W. Kohn, Phys. Rev. 136, B864 (1964). % \bibitem{KS} W. Kohn and L. J. Sham, Phys. Rev. 140, A1133 (1965). % \bibitem{Wang1} S. G. Wang, D. K. Pan, and W. H. E. Schwarz, J. Chem. Phys. 102, 9296 (1995). % \bibitem{Eschrig} J. Forstreuter, L. Steinbeck, M. Richter, and H. Eschrig, Phys. Rev. B 55, 9415 (1997). % \bibitem{BDFTCA} W. Liu, G. Hong, D. Dai, L. Li, and M. Dolg, Theor. Chem. Acc. 96, 75 (1997). % \bibitem{BDFJCP} W. Liu, M. Dolg, and L. Li, J. Chem. Phys. (submitted). % \bibitem{MOLPRO}Program system MOLPRO, P. J. Knowles and H.-J. Werner, Chem. Phys. Lett. 115, 5053 (1985); H.-J. Werner and P. J. Knowles, J. Chem. Phys. 89, 5803 (1988); P. J. Knowles and H.-J. Werner, Chem. Phys. Lett. 145, 514 (1988); P. J. Knowles and H.-J. Werner, Theor. Chim. Acta 84, 95 (1992). % \bibitem{ACPF} R. J. Gdanitz and R. Ahlrichs, Chem. Phys. Lett. 143, 413 (1988). % \bibitem{MCHF} atomic structure code MCHF; C. Froese Fischer, {\it The Hartree-Fock Method for Atoms - A Numerical Approach\/} (Wiley, New York 1976). % \bibitem{GRASP} atomic structure code GRASP; K. G. Dyall, I. P. Grant, C. T. Johnson, F. A. Parpia, and E. P. Plummer, Comput. Phys. Commun. 55, 425 (1989). % \bibitem{PW92} J. P. Perdew and Y. Wang, Phys. Rev. B 45, 13244 (1992). % \bibitem{SIC} H. Stoll, C. M. E. Pavlidou, and H. Preuss, Theor. Chim. Acta 49, 143 (1978); H. Stoll, E. Golka, and H. Preuss, Theor. Chim. Acta 55, 29 (1980). % \bibitem{Beck88} A. D. Becke, Phys. Rev. A 38, 3098 (1988). % \bibitem{Pd86} J. P. Perdew, Phys. Rev. B 33, 8822 (1986); ibid, 34, 7406(E) (1986). % \bibitem{Yang} W. Yang, J. Chem. Phys. 94, 1208 (1991). % \bibitem{Lebedev} V. I. Lebedev, Zh. Vychisl. Mat. Mat. Fiz. 15, 48 (1975); V. I. Lebedev, Zh. Vychisl. Mat. Mat. Fiz. 16, 293 (1976); V. I. Lebedev, Sibirsk. Mat. Zh. 18, 32 (1977). \bibitem{TS} T. Ziegler and A. Rauk, Theor. Chim. Acta 46, 1 (1977). % \bibitem{Matsuoka} H. Tatewaki, M. Sekiya, F. Sasaki, O. Matsuoka, and T. Koka, Phys. Rev. A 51, 197 (1995). % \bibitem{Savin} A Savin, in \begin{it} Recent Developments and Applications of Modern Density Functional Theory \end{it} , edited by J. M. Seminario (Elsevier, Amsterdam, 1996), p.327. % \bibitem{RVxLDA} A. H. MacDonald and S. H. Vosko, J. Phys. C 12, 2977 (1979). % \bibitem{PZSI} J. P. Perdew and A. Zunger, Phys. Rev. B 23, 5048 (1981). % \bibitem{private} H. Eschrig (private communications). % \bibitem{MPrev} R. O. Jones and O. Gunnarsson, Rev. Mod. Phys. 61, 689 (1989). % \bibitem{euo} M. Dolg, H. Stoll, H. Preu{\ss}, Chem. Phys. 148, 219 (1990). % \end{thebibliography} % \newpage \onecolumn \begin{figure}[htb] \hspace{11mm}\psfig{file=fig1.eps,width=130mm,angle=-90}\vspace{-8mm} %\centerline{\psfig{figure=fig1.eps,width=9cm,angle=0}} \vspace{2cm} \small %\caption{Absolute errors in the first to fourth ionization potentials from %scalar relativistic ab initio pseudopotential calculations without (SCF) %and with (ACPF) electron correlation effects. Spin-orbit corrections are %not included.} \label{fig1} \normalsize \end{figure} \mbox{ } \newpage \onecolumn \begin{table} \caption{ The first ionization potential (IP$_1$ in eV) for the the lanthanide atoms from the present fully relativistic density functional calculations (BDF \protect\cite{BDFTCA,BDFJCP}) and ab initio quasirelativistic pseudopotential calculations (QR-PP) \protect\cite{dolg2} in comparison to other theoretical results (D$^2$: squared Dirac equation with relativistically corrected density functionals (RLDASIC) \protect\cite{Eschrig}) and experimental data (Expt. \protect\cite{ExptRE}). LDASIC: local density approximation \protect\cite{PW92} with self-interaction correction \protect\cite{SIC}; B88, nonlocal exchange correction \protect\cite{Beck88}; P86, nonlocal correlation correction \protect\cite{Pd86}; BP, both B88 and P86 \protect\cite{Beck88,Pd86}; ACPF, averaged coupled-pair functional \protect\cite{ACPF} with spin-orbit coupling corrections. The mean absolute error (m.a.e.) and the largest relative error (l.r.e.) are also given.} \begin{center} \begin{tabular}{lcllllllllllllllll} &&&\multicolumn{4}{c}{BDF} && \multicolumn{1}{c}{D$^2$}&&\multicolumn{1}{c}{QR-PP}&& $ $ \\ \cline{4-7}\cline{9-9}\cline{11-11} atom &configurations&&LDASIC & B88 & P86 & BP &&RLDASIC &&ACPF&& Expt.\\ \hline $_{57}$La&$ f^0d^1s^2\rightarrow f^0d^2s^0$&&5.59 & 5.69 & 5.42 & 5.51 && $ $&& &&5.58\\ $_{58}$Ce&$ f^1d^1s^2\rightarrow f^1d^2s^0$$^a$&&5.69 & 5.75 & 5.48 & 5.54 && 5.8&&5.62&&5.54\\ &$ f^2d^0s^2\rightarrow f^2d^0s^1$$^b$&&5.21 & 5.39 & 5.03 & 5.21 && && && \\ $_{59}$Pr&$ f^3d^0s^2\rightarrow f^3d^0s^1$&&5.24 & 5.35 & 5.06 & 5.17 && 5.7&&5.39&&5.46\\ $_{60}$Nd&$ f^4d^0s^2\rightarrow f^4d^0s^1$&&5.29 & 5.40 & 5.09 & 5.21 && 5.8&&5.44&&5.53\\ $_{61}$Pm&$ f^5d^0s^2\rightarrow f^5d^0s^1$&&5.33 & 5.45 & 5.13 & 5.25 && 5.8&&5.48&&5.55\\ $_{62}$Sm&$ f^6d^0s^2\rightarrow f^6d^0s^1$&&5.38 & 5.49 & 5.17 & 5.29 && 5.8&&5.51&&5.64\\ $_{63}$Eu&$ f^7d^0s^2\rightarrow f^7d^0s^1$&&5.68 & 5.80 & 5.45 & 5.58 && 5.9&&5.53&&5.67\\ & &&5.42$^{c}$ & 5.53$^{c}$ & & && && && \\ $_{64}$Gd&$ f^7d^1s^2\rightarrow f^7d^1s^1$$^a$&&5.84 & 5.97 & 5.59 & 5.72 && 6.3 &&6.02&&6.15\\ &$ f^8d^0s^2\rightarrow f^8d^0s^1$$^b$&&6.32 & 6.44 & 5.88 & 5.99 && && && \\ & &&5.87$^{c}$ & 5.99$^{c}$ & & && && && \\ $_{65}$Tb&$ f^9d^0s^2\rightarrow f^9d^0s^1$&&5.65 & 5.77 & 5.44 & 5.56 && 5.8 &&5.58&&5.86\\ $_{66}$Dy&$f^{10}d^0s^2\rightarrow f^{10}d^0s^1$&&5.76 & 5.88 & 5.55 & 5.67 && 5.8 &&5.73&&5.94\\ $_{67}$Ho&$f^{11}d^0s^2\rightarrow f^{11}d^0s^1$&&5.86 & 5.98 & 5.66 & 5.77 && 5.9 &&5.82&&6.02\\ $_{68}$Er&$f^{12}d^0s^2\rightarrow f^{12}d^0s^1$&&5.97 & 6.08 & 5.76 & 5.87 && 6.0 &&5.89&&6.11\\ $_{69}$Tm&$f^{13}d^0s^2\rightarrow f^{13}d^0s^1$&&6.05 & 6.17 & 5.84 & 5.95 && 6.0 &&5.95&&6.18\\ $_{70}$Yb&$f^{14}d^0s^2\rightarrow f^{14}d^0s^1$$^d$&&6.33 & 6.43 & 6.11 & 6.22 && 6.1 &&6.02&&6.25\\ & &&6.06$^{c}$ & 6.18$^{c}$ & & & & & & && \\ $_{71}$Lu&$f^{14}d^1s^2\rightarrow f^{14}d^0s^2$$^e$&&5.24 & 5.30 & 5.12 & 5.18 && 5.7 && &&5.43\\ \multicolumn{2}{c}{m.a.e. (eV)}&&0.17&0.11&0.34&0.23&&0.2&&0.16&& $ $\\ \multicolumn{2}{c}{l.r.e. ($\%$)}&&5.0&3.8&9.1&7.0&&5.0&&5.1&& $ $\\ \end{tabular} \end{center} $^a$ Experimentally measured lowest configurations.\\ $^b$ DFT calculated lowest configurations.\\ $^{c}$ DFT-FORPT : first-order relativistic perturbation theory \protect\cite{Wang1}.\\ $^{d}$ The result of a relativistic all-electron coupled-cluster calculation for Yb with an uncontracted (31s26p21d15f10g6h) basis set is 6.34 eV \protect\cite{Kaldor}.\\ $^{e}$ The result of a relativistic all-electron coupled-cluster calculation for Lu with an uncontracted (34s25p20d15f10g6h) basis set is 5.30 eV \protect\cite{Kaldor}.\\ \end{table} % % \begin{table} \caption{ The second ionization potential (IP$_2$ in eV). For other explanations cf. table 1.} \begin{center} \begin{tabular}{lcrrrrrrrrrrrrrrrllll} &&&\multicolumn{4}{c}{BDF} && \multicolumn{1}{c}{D$^2$}&&\multicolumn{1}{c}{QR-PP}&& $ $ \\ \cline{4-7}\cline{9-9}\cline{11-11} atom &configurations&&LDASIC & B88 & P86 & BP &&RLDASIC &&ACPF&& Expt.\\ \hline $_{57}$La&$ f^0d^2s^0\rightarrow f^0d^1s^0$&&10.85&10.93&10.74&10.81&& $ $&& &&11.06\\ $_{58}$Ce&$ f^1d^2s^0\rightarrow f^2d^0s^0$$^a$&& 9.54& 9.65& 9.50& 9.61&& 8.8&&11.06&&10.85\\ &$ f^2d^0s^1\rightarrow f^2d^0s^0$$^b$&&10.43&10.49&10.26&10.32&& && && \\ $_{59}$Pr&$ f^3d^0s^1\rightarrow f^3d^0s^0$&&10.61&10.68&10.44&10.51&&10.9&&10.57&&10.55\\ $_{60}$Nd&$ f^4d^0s^1\rightarrow f^4d^0s^0$&&10.79&10.86&10.62&10.68&&11.0&&10.73&&10.73\\ $_{61}$Pm&$ f^5d^0s^1\rightarrow f^5d^0s^0$&&10.95&11.02&10.77&10.84&&11.3&&10.87&&10.90\\ $_{62}$Sm&$ f^6d^0s^1\rightarrow f^6d^0s^0$&&11.11&11.17&10.93&10.99&&11.5&&10.98&&11.07\\ $_{63}$Eu&$ f^7d^0s^1\rightarrow f^7d^0s^0$&&11.27&11.35&11.10&11.19&&11.6&&11.11&&11.24\\ $_{64}$Gd&$ f^7d^1s^1\rightarrow f^7d^1s^0$$^a$&&12.26&12.30&12.13&12.17&&12.7&&12.05&&12.09\\ &$ f^8d^0s^1\rightarrow f^8d^0s^0$$^b$&&11.40&11.47&11.20&11.27&& && && \\ $_{65}$Tb&$ f^9d^0s^1\rightarrow f^9d^0s^0$&&11.54&11.61&11.33&11.41&&11.9&&11.14&&11.52\\ $_{66}$Dy&$f^{10}d^0s^1\rightarrow f^{10}d^0s^0$&&11.67&11.75&11.46&11.54&&12.1&&11.41&&11.67\\ $_{67}$Ho&$f^{11}d^0s^1\rightarrow f^{11}d^0s^0$&&11.80&11.88&11.58&11.66&&12.2&&11.57&&11.80\\ $_{68}$Er&$f^{12}d^0s^1\rightarrow f^{12}d^0s^0$&&11.93&12.01&11.70&11.79&&12.3&&11.69&&11.93\\ $_{69}$Tm&$f^{13}d^0s^1\rightarrow f^{13}d^0s^0$&&12.19&12.28&11.96&12.05&&12.5&&11.77&&12.05\\ $_{70}$Yb&$f^{14}d^0s^1\rightarrow f^{14}d^0s^0$$^c$&&12.13&12.23&11.90&12.00&&12.6&&11.73&&12.18\\ $_{71}$Lu&$f^{14}d^0s^2\rightarrow f^{14}d^0s^1$$^d$&&13.86&13.97&13.60&13.71&&14.2&& &&13.90\\ \multicolumn{2}{c}{m.a.e. (eV)}&&0.15&0.19&0.26&0.18&&0.5&& 0.18&& $ $\\ \multicolumn{2}{c}{l.r.e. ($\%$)}&&1.4&1.9&2.9&2.3&&5.1&&3.8&& $ $\\ \end{tabular} \end{center} $^a$ Experimentally measured lowest configurations.\\ $^b$ DFT calculated lowest configurations.\\ $^{c}$ The result of a relativistic all-electron coupled-cluster calculation for Yb with an uncontracted (31s26p21d15f10g6h) basis set is 12.14 eV \protect\cite{Kaldor}.\\ $^{d}$ The result of a relativistic all-electron coupled-cluster calculation for Lu with an uncontracted (34s25p20d15f10g6h) basis set is 14.12 eV \protect\cite{Kaldor}.\\ \end{table} % % \begin{table} \caption{ The third ionization potential (IP$_3$ in eV). For other explanations cf. table 1.} \begin{center} \begin{tabular}{lcrrrrrrrrrrlll} &&&\multicolumn{4}{c}{BDF} && \multicolumn{1}{c}{D$^2$}&&\multicolumn{1}{c}{QR-PP}&& $ $ \\ \cline{4-7}\cline{9-9}\cline{11-11} atom &configurations&&LDASIC & B88 & P86 & BP &&RLDASIC &&ACPF&& Expt.\\ \hline $_{57}$La&$ f^0d^1s^0\rightarrow f^0d^0s^0$&&18.74&18.87&18.54&18.67&&$ $&&&&19.18\\ $_{58}$Ce&$ f^2d^0s^0\rightarrow f^1d^0s^0$&&20.53&20.61&20.26&20.34&&22.6&&19.36&&20.20\\ $_{59}$Pr&$ f^3d^0s^0\rightarrow f^2d^0s^0$&&21.79&21.87&21.53&21.61&&23.9&&21.04&&21.62\\ $_{60}$Nd&$ f^4d^0s^0\rightarrow f^3d^0s^0$&&22.44&22.53&22.19&22.28&&24.9&&21.52&&22.1$\pm$0.3\\ $_{61}$Pm&$ f^5d^0s^0\rightarrow f^4d^0s^0$&&23.39&23.48&23.15&23.23&&25.7&&21.87&&22.3$\pm$0.4\\ $_{62}$Sm&$ f^6d^0s^0\rightarrow f^5d^0s^0$&&24.22&24.31&23.98&24.07&&26.6&&23.14&&23.4$\pm$0.3\\ $_{63}$Eu&$ f^7d^0s^0\rightarrow f^6d^0s^0$&&24.64&24.72&24.41&24.46&&27.3&&24.56&&24.92\\ $_{64}$Gd&$ f^7d^1s^0\rightarrow f^7d^0s^0$&&20.13&20.19&19.92&19.98&&20.9&&20.59&&20.63\\ $_{65}$Tb&$ f^9d^0s^0\rightarrow f^8d^0s^0$&&22.75&22.79&22.37&22.41&&23.7&&21.19&&21.91\\ $_{66}$Dy&$f^{10}d^0s^0\rightarrow f^{ 9}d^0s^0$&&23.62&23.66&23.25&23.30&&24.3&&22.25&&22.8$\pm$0.3\\ $_{67}$Ho&$f^{11}d^0s^0\rightarrow f^{10}d^0s^0$&&23.36&23.42&23.01&23.07&&25.1&&22.04&&22.84\\ $_{68}$Er&$f^{12}d^0s^0\rightarrow f^{11}d^0s^0$&&24.03&24.10&23.69&23.75&&25.6&&21.88&&22.74\\ $_{69}$Tm&$f^{13}d^0s^0\rightarrow f^{12}d^0s^0$&&24.53&24.61&24.21&24.30&&25.9&&22.89&&23.68\\ $_{70}$Yb&$f^{14}d^0s^0\rightarrow f^{13}d^0s^0$&&25.03&25.10&24.70&24.76&&26.3&&24.27&&25.05\\ $_{71}$Lu&$f^{14}d^0s^1\rightarrow f^{14}d^0s^0$$^a$&&21.01&21.10&20.70&20.79&&21.5&& &&20.96\\ \multicolumn{2}{c}{m.a.e. (eV)}&&0.56&0.60&0.45&0.46&&2.1&& 0.58&& $ $\\ \multicolumn{2}{c}{l.r.e. ($\%$)}&&5.7&6.0&4.2&4.4&&15.3&&4.3&& $ $\\ \end{tabular} \end{center} $^{a}$ The result of a relativistic all-electron coupled-cluster calculation for Lu with an uncontracted (34s25p20d15f10g6h) basis set is 20.97 eV \protect\cite{Kaldor}.\\ \end{table} % % \begin{table} \caption{ The fourth ionization potential (IP$_4$ in eV). For other explanations cf. table 1.} \begin{center} \begin{tabular}{lcrrrrrrrll} &&&\multicolumn{4}{c}{BDF} &&\multicolumn{1}{c}{QR-PP}&& \\ \cline{4-7}\cline{9-9} atom &configurations&&LDASIC & B88 & P86 & BP &&ACPF&&Expt.\\ \hline $_{57}$La&$ 5s^25p^6\rightarrow 5s^25p^5 $&&49.18&49.31&48.94&49.07&& &&49.95\\ $_{58}$Ce&$ f^1d^0s^0\rightarrow f^0d^0s^0$&&37.54&37.62&37.22&37.30&&36.05&&36.76\\ $_{59}$Pr&$ f^2d^0s^0\rightarrow f^1d^0s^0$$^a$&&39.10&39.19&38.80&38.88&&38.48&&38.98\\ $_{60}$Nd&$ f^3d^0s^0\rightarrow f^2d^0s^0$&&40.50&40.59&40.20&40.30&&40.26&&40.4$\pm$0.1\\ $_{61}$Pm&$ f^4d^0s^0\rightarrow f^3d^0s^0$&&41.19&41.29&40.91&41.01&&40.81&&41.1$\pm$0.6\\ $_{62}$Sm&$ f^5d^0s^0\rightarrow f^4d^0s^0$&&42.31&42.41&42.04&42.14&&41.26&&41.4$\pm$0.7\\ $_{63}$Eu&$ f^6d^0s^0\rightarrow f^5d^0s^0$&&43.34&43.43&43.07&43.17&&42.73&&42.7$\pm$0.6\\ $_{64}$Gd&$ f^7d^0s^0\rightarrow f^6d^0s^0$&&44.26&44.36&44.00&44.10&&44.86&&44.0$\pm$0.7\\ $_{65}$Tb&$ f^8d^0s^0\rightarrow f^7d^0s^0$&&40.89&40.92&40.47&40.59&&38.96&&39.37\\ $_{66}$Dy&$f^{ 9}d^0s^0\rightarrow f^{ 8}d^0s^0$&&42.03&42.07&41.63&41.67&&40.79&&41.4$\pm$0.4\\ $_{67}$Ho&$f^{10}d^0s^0\rightarrow f^{ 9}d^0s^0$&&43.08&43.13&42.69&42.74&&42.10&&42.5$\pm$0.6\\ $_{68}$Er&$f^{11}d^0s^0\rightarrow f^{10}d^0s^0$&&42.81&42.87&42.43&42.50&&42.06&&42.7$\pm$0.4\\ $_{69}$Tm&$f^{12}d^0s^0\rightarrow f^{11}d^0s^0$&&43.58&43.66&43.23&43.31&&43.32&&42.7$\pm$0.4\\ $_{70}$Yb&$f^{13}d^0s^0\rightarrow f^{12}d^0s^0$&&44.40&44.48&44.05&44.12&&43.11&&43.56\\ $_{71}$Lu&$f^{14}d^0s^0\rightarrow f^{13}d^0s^0$&&45.10&45.17&44.75&44.83&& &&45.25\\ \multicolumn{2}{c}{m.a.e. (eV)}&&0.56&0.61&0.42&0.43&&0.45&& $ $\\ \multicolumn{2}{c}{l.r.e. ($\%$)}&&3.9&3.9&2.8&2.9&&2.0& $ $\\ \end{tabular} \end{center} $^{a}$ The result of a relativistic all-electron coupled-cluster calculation for Pr with an uncontracted (29s23p19d14f10g6h4i) basis set is 38.61 eV \protect\cite{Kaldor}.\\ \end{table} %%%%% \begin{table} \caption{ The 'df' charge transfer energies (eV) defined as $\Delta_{df}$ = $E(f^{n}d^1s^2)$-$E(f^{n+1}d^0s^2)$ (n=0-13 for La-Yb). Density functional results are for an averaged occupation of the open shells. For other explanations cf. table 1.} \begin{center} \begin{tabular}{lrrrrrrrrl} &&\multicolumn{4}{c}{BDF} &&\multicolumn{1}{c}{QR-PP}&& \\ \cline{3-6} \cline{3-6}\cline{8-8} atom&&LDASIC & B88 & P86 & BP &&ACPF&&Expt.\\ \hline $_{57}$La &&-1.34&-1.32&-1.52&-1.49&& &&-1.88\\ $_{58}$Ce && 0.40 &0.46 &0.31& 0.36&& -0.96 &&-0.59\\ $_{59}$Pr && 1.33 &1.34 &1.25 &1.26&& && 0.55\\ $_{60}$Nd && 1.77 &1.86 &1.69 &1.78&& && 0.84\\ $_{61}$Pm && 2.49 &2.52 &2.43 &2.46&& && \\ $_{62}$Sm && 3.13 &3.16 &3.07 &3.10&& && 2.24\\ $_{63}$Eu && 3.74 &3.78 &3.68 &3.72&& && 3.33$^a$\\ $_{64}$Gd && 0.47 &0.45 &0.27 &0.25&& -2.22 &&-1.36\\ $_{65}$Tb && 1.40 &1.39 &1.22 &1.22&& -0.15 && 0.04\\ $_{66}$Dy && 2.23 &2.23 &2.07 &2.07&& 0.60 && 0.94\\ $_{67}$Ho && 1.82 &1.83 &1.69 &1.69&& 0.65 && 1.04\\ $_{68}$Er && 2.43 &2.44 &2.31 &2.32&& 0.36 && 0.89\\ $_{69}$Tm && 2.99 &3.00 &2.89 &2.90&& 1.30 && 1.63\\ $_{70}$Yb && 3.46 &3.50 &3.36 &3.40&& && 2.88\\ \multicolumn{2}{c}{m.a.e. (eV)}&1.02&1.04&0.90&0.93&& 0.43 &&$ $\\ \end{tabular} \end{center} $^{a}$ estimated value \cite{euo}. \end{table} %%%% \newpage \onecolumn \begin{center} {\bf FIGURE CAPTION} \end{center} \setcounter{figure}{0} \begin{figure}[htb] \vspace{2cm} \caption{Absolute errors in the first to fourth ionization potentials from scalar relativistic ab initio pseudopotential calculations without (SCF) and with (ACPF) electron correlation effects. Spin-orbit corrections are not included.} \end{figure} %% \end{document}