week ending
29 JULY 2005

PHYSICAL REVIEW LETTERSPRL 95, 050407 (2005)

Analytically Solving the Relativistic Dirac-Coulomb Equation for Atoms and Molecules
Hiroshi Nakatsuji1,2 and Hiroyuki Nakashima1
1

Department of Synthetic Chemistry and Biological Chemistry, Graduate School of Engineering, Kyoto University,
Katsura, Nishikyo-ku, Kyoto 615-8510, Japan
2
Fukui Institute for Fundamental Chemistry, Kyoto University, 34-4 Takano-Nishihiraki-cho, Sakyo-ku, Kyoto 606-8103, Japan(Received 26 April 2005; published 29 July 2005)
A method of solving the Schrödinger equation in an analytically expanded form reported previously isextended to the relativistic case and a general method of exactly solving the Dirac-Coulomb equation has
been proposed for atoms and molecules. Some problems characteristic of the relativistic case are
discussed. Test applications to the hydrogenlike, heliumlike atoms are satisfactory, implying a highpotentiality of the proposed method also for the relativistic case.
DOI: 10.1103/PhysRevLett.95.050407

PACS numbers: 03.65.Pm, 03.65.Ca, 03.65.Ge, 31.15.2p

As noted by Dirac in 1929 [1], the Schrödinger equation(SE) for atoms and molecules,

H

X
i

H E ;
XX
 i =2  ZA =rAi  1=rij
i;A

(1)
(2)

i>j

is a major quantum principle governing chemistry. So,
establishing a general method of exactly solving the SE
is a dream of theoretical chemists. In the SE, the exact
wave function is determined by the Hamiltonian H, sothat the exact may be written as a functional of H applied
to some appropriate function 0 as
 fH 0 :(3)Recently, we have studied the structure of the exact wave
function, i.e., a possible functional form of fH, andproposed the iterative configuration or composites interaction (ICI) method and the simplest extreme coupled cluster(SECC) method as the methods of constructing the exact
wave function [2 – 4]. Generally, when f is expanded in a
Maclaurin series of H, it would include higher powers ofH, but the integrals of the higher powers of H over approximate wave functions do not exist, i.e., diverge because of the singularity of the Coulomb potentials includedin the Hamiltonian [4]. Thus, the singularity difficultycommonly occurs in the method of analytically solving
the SE of atoms and molecules.
We have proposed two methods of solving the singularity problem, i.e., by introducing the inverse Schrödinger
equation (ISE) [3] and the scaled Schrödinger equation(SSE) [4]. Between the two, the latter method was simpler
and more general than the former. Combined with the ICI
method [2], it opened a new field of calculating the exact
wave functions of atoms and molecules in an analyticalexpansion form [4]. We have applied the method successfully to helium atom and hydrogen molecule [4,5] andfurther to three-to-five electron atoms [5]: these calculations have given a basis to confirm the high potentiality ofthe proposed method. A short overview was given as a partof the account of the recent works of one of the authors [6].
0031-9007=05=95(5)=050407(4)$23.00

Dirac equation and Dirac-Coulomb equation. —Whenatoms and molecules include heavier elements, the relativistic effect becomes very important and even dominant. For
a one-electron system, the Dirac equation (DE) describes
exactly the relativistic motion of an electron. It has an
explicit Hamiltonian that obeys the Lorentz transformation. For many-electron systems, however, we do not havesuch relativistic equations that are as compact and accurate
as the DE. The Dirac-Coulomb equation (DCE) and the
Dirac-Coulomb-Breit equation (DCBE) are approximate
and do not exactly obey the Lorentz transformation, butthey are expected to be very accurate for atoms and molecules. For more accurate formulations, we use quantum
electrodynamics, but there, we do not have explicit
Hamiltonians.
The DE for one-electron system and the DCE for N
electron system are written asH^  E ;(4)

where the Hamiltonian H^ is defined by the equations,
H^ 

N
X

H^ N i  I^N i

N
X

H^ N i  I^i1  h^i  I^Ni ;h^i wij ;

i<jvi  c2 I^2c  pi !

(5)c  pi ;vi  c2 I^2

and wij  1=rij . I^n and I^n are the unit matrices of orders n
and 4n , respectively. h^i is the 4  4 matrix one-electron
Hamiltonian and vi  A ZA =rAi is the nuclear attraction
operator for electron i.  and pi are the ordinary Pauli spin
matrix and the momentum operator of electron i. The
relativistic Hamiltonian H^ is a 4N  4N operator matrix
and the wave function
is a vector composed of 4N
elements. The term E in Eq. (4) is also written as
EI^N : mostly we use the simpler expression. For the
one-electron case (N  1), the DCE reduces to the DE.
In the nonrelativistic limit, the DCE reduces to the SE. So,
the DCE should be very accurate for atoms and molecules.

050407-1

 2005 The American Physical Society

PRL 95, 050407 (2005)

week ending
29 JULY 2005

PHYSICAL REVIEW LETTERS

In this Letter, our aim is to formulate a general method ofsolving the ‘‘exact’’ of the DE and DCE in an analyticalexpansion form.
In the relativistic case, there are no theories dealing withthe structure of the exact wave function, nor the theories
for constructing the exact wave function. So, we first
generalize our theory for the SE to the relativistic case.
The Coulomb singularity problem is common to both thenonrelativistic and relativistic cases and is solved by introducing the scaled equations [4]. Another problem that
often appears in the relativistic field is the variational
instability and collapse. Though some studies using kinetic
and other balances [7,8] were reported, there seems to be
no established method. It was impressive that Pestka andKarwowski closed their important chapter of the
Rychlewski’s book in 2003 [7] by noting ‘‘The Hylleraas
CI approach to solving the Dirac-Coulomb eigenvalue
problem is still in its infancy. However, it certainly offersa perspective of highly accurate benchmark results, at least
for two-electron systems.’’
Formulation of the exact wave function. —Now the firststep of our project is to clarify the structure of the exact
relativistic wave function
and then to give a general
method of constructing it in an analytical form. We follow
the same arguments as in the previous articles [2 – 4]. First,
the variational principle for is written as
h jH^  Ej i  0:

(6)This is a stationary principle for calculating the best possible . The H-square equation is formulated ash jH^  E2 j i  0;

(7)

and is satisfied only with the exact . Then, we have thefollowing theorem.
Theorem: When a includes only one variable matrix C^(diagonal) with the element Cn (n  1 . . . 4N ) and when an
arbitrary variation of satisfies^ H^  E ;  C(8)then this has the structure of the exact wave function.
Proof: When we define H^  E  ,  is a column
vector of the elements n n  1 . . . 4N . When we substitute  given by Eq. (8) into the variational principle,^ H^  Ej i 
Eq. (6), we obtain 0  h jH^  ECN
42
n Cn hjn j i. Since Cn is arbitrary, this meanshjn j2 i  0 for all n, which is equivalent to Eq. (7). So,among with the form of Eq. (8), variationally optimal
wave functions are the exact solutions of the DE and theDCE. Proof ends. This theorem is a generalization of thesimilar theorem given previously for the nonrelativistic
case [2] and the argument is free from the technical arguments like variational collapse, since the H-square equation is valid only for the exact wave function.
Based on the above theorem, we propose the relativistic
ICI method just as in the nonrelativistic case. We define thesimplest ICI (SICI) wave function by the recurrence for-

mula,
n1 

1  C^ n H^  En 

n;

(9)

where n is an iteration number. The variation of n1 , n1  C^ n H^  En  n , becomes equal to that given
by Eq. (8) at the convergence where everything becomesn independent, so that the SICI gives the exact wave
function at convergence. We can also formulate a moregeneral ICI method by dividing the operator H^ in the ND
parts (ICIND) [2], but omit it here for brevity.
The DC Hamiltonian includes the Coulomb potentials
that cause the singularity problem. To avoid it, we introduce the scaled DE and DCE that are given bygH^  EI^N   0;

(10)

where the scaling function g is a positive scalar functionof electron coordinates and can become zero only at the
singular points r0 of the potential V but satisfies
limr!r0 gV  0. By this introduction of the scaled equation, we can avoid the singularity problem as in the nonrelativistic case [4]. Since this idea is very general, we canuse it for singular potentials that appear in other fields ofphysics. Based on the scaled DCE, we introduce the SICI
method that is free from the singularity problem byn1 

1  C^ n gH^  En 

n:

(11)

The SICI includes only one variable matrix C^ n at each
iteration and its convergence may be slow [2]. A better
performance is expected by introducing the free ICI
method as in the nonrelativistic case [4]. In the free ICI,
we take all the independent functions fi gn from the righthand side of Eq. (11) and give independent coefficients toall such functions i , i.e.,
n1 

all
X

ci i :

(12)i

From the variational point of view, this free ICI should
converge faster than the SICI. Note that i is a column
vector of 4N dimension and ci is a diagonal matrix with 4N
elements. In the free ICI, each iteration does not depend on
the results of the former iterations and so no accumulationof error occurs. In some choice of g, the right-hand side ofEq. (11) may include diverging functions but they are
eliminated in the free ICI process since such functions
are inappropriate for describing the wave function.
Variational calculations. —Now, we have an analytical
function having the structure of the exact wave function in
the form of free ICI [Eq. (12)]. The next step is the variational calculation of the variables fci g and the corresponding energy, which we call the diagonalization step. Here,an obstacle that often appears in the relativistic field is theso-called variational collapse. Since the lowest electronic
state of the DE and the DCE is not the lowest state of theseequations, the so-called Ritz-type property does not hold.
However, a method to recover Ritz-type property was

050407-2

week ending
29 JULY 2005

PHYSICAL REVIEW LETTERSPRL 95, 050407 (2005)

proposed by Hill and Krauthauser [9] for one-electron DC
calculations with basis set expansion method. They introduced the inverse Hamiltonian and wrote the DE asH^ 1  E1 :(13)

As in our case of the inverse Schrödinger equation [3], it iseasy to show the equivalence between the original equationand the inverse one. For the inverse DE, the electronic
ground state is the highest solution against the complete
vacuum and therefore we have the Ritz-like variational~ is aprinciple E~1  h ~jH^ 1 j ~i=h ~j ~i  E1
0 , where
variational trial function for and E0 is the true energy ofthe ground electronic state. Now, how do we explicitly
write down the inverse Hamiltonian? A clever trick [9] is^
to take our variational function in the form of ~  H’
with ’ representing a free variation. Then, the above variational equation is rewritten as^E~ 1  h’jHj’i=h’j
H^ 2 j’i  E1
0 ;

(14)

and the variance is formulated as2  h ~jH^ 1  E1 2 j ~i=h ~j ~i^ h’j’i  2E1 h’jHj’i
 E2 h’jH^ 2 j’i=h’jH^ 2 j’iwhich we can calculate. For the exact wave function, 2must be zero: this is the inverse H-square theorem [3].
At the diagonalization step of the relativistic free ICI, we
utilize the above Ritz-like variational equation to ensure
the bound property. Though it was introduced originallyfor the one-electron DE [9], we utilize it for solving the
many-electron DCE. We know from the theory of ICI that
our free ICI function ’n approaches the exact wave function, and therefore it is easy to show that our variational
function ~n also approaches the exact wave function. SinceH^ must not be singular to define its inverse, we often haveto shift it by a constant W, i.e., H^ 0  H^  W I^N . Further,since Eq. (14) includes the integral of H^ 2 that is more
^ we have to eliminate a larger
strongly singular than H,
number of diverging functions from our ICI basis than in^ This may make thethe ordinary case that involves only H.

inverse variation method rather slowly converging.
Nevertheless, the bound property of this variational
method is useful in actual calculations where we do notknow the true energy E0 . In the ICI formulation, theexcited states are calculated as the higher solutions: anappropriate choice of 0 is necessary for obtaining areasonable convergence speed.
Another method of avoiding the variational collapse is toensure kinetic balance, etc. on the basis functions [7,8,10].
From the DE, we obtain a relation connecting small
and large components as S  cp=E  V  c2  L .
A strict use of this relation is called atomic balance and an
approximate one kinetic balance. We note here that the balancing similar to the atomic balance is automatically donein the ICI formalism. In ICI, n1 is generated by applying
the DC operator to n [Eq. (11)] and this is essentially abalancing that is theoretically valid even for many-electron
systems. We call this balancing ICI balance. In an earlier
stage of the ICI calculation, the number of the basis functions is small and the ICI balancing would be insufficient.
So, practically, we have to do an additional balancing
before diagonalization. Detailed examinations of such balancing methods will be given in the forthcoming paper.
Hydrogen isoelectronic atoms.—The first application of
the present theory is to hydrogen isoelectronic atoms forwhich the exact analytical wave functions are known. Wetake H (Z  1) and Fe25 (Z  26): the relativistic effectis larger for the latter. For the ground state of these oneelectron atoms, the DE is written with the radial part alone,
separated from the angular part, as [11]
c2  Z=rc d=dr  1  k=r
E ;c d=dr  1  k=r
c2  Z=rwhere k   j  1=2 with j being the quantum numberof the total angular momentum and c  137:03596 [12].
The relativistic exact wave function has a mild singularity
at the position of the nucleus [11].
In the free ICI, we took the starting wave function asLS
0  0  exp$r with $  1:5Z instead of $  Z.
The scaling function was taken to be g  r99=100 instead ofg  r in order to describe, ad hoc, the mild singularity. The
free ICI iteration process is shown in Table I. The varia-

TABLE I. Relativistic free ICI energy (a.u.) of hydrogenlike atoms (Z  1 and 26).anb

Mc

Energy (H)Energy(Fe25 )L (Fe25 )S (Fe25 )

0
1
2
3
4
5

2
6
12
20
30
42

0:375 033 696
0:493 381 703
0:499 822 155
0:500 004 337
0:500 006 642
0:500 006 657
0:500 006 657

268:305 491
336:985 794
340:983 740
341:096 421
341:097 830
341:097 839
341:097 839

3:3  101
9:8  102
1:6  102
1:6  103
1:2  104
7:2  106
0.0

3:3  101
7:1  102
1:1  102
1:2  103
9:0  105
5:3  106
0.0

a

The bottom row shows the exact values.
Iteration number.
cNumber of the independent functions for expanding . Half and half for the large and small components.
b050407-3

PRL 95, 050407 (2005)

week ending
29 JULY 2005

PHYSICAL REVIEW LETTERS

tional method was an ordinary one. For both H and Fe25 ,the energy at fifth iteration agrees
with the exact value to
R ICInine figures.  defined by
j n  exact j2 d& 1=2 approaches zero for both large and small components asshown for heavier Fe25 : the ICI wave function itself
also converges to the exact one. We have performed
many different calculations, but as far as we do ICI, wehave never experienced the variational collapse for the oneelectron case. It was unnecessary to resort to the inverse
Hamiltonian as Hill et al. considered originally [9].
Helium isoelectronic atoms.—The next application is toheliumlike two-electron atoms, whose DCE is the matrix
equation of dimension 16. We calculated He (Z  2) andTh88 (Z  90). The wave function is composed of ll ,ls ,sl , andss , each being 4-vector. They satisfyll 1;2   ll 2;1, ls 1;2   sl 2;1, and ss 1;2  ss 2;1: ls and sl are not independent. The initial
function was taken as 0ll  As1=2;1=2  s1=2;1=2 , 0ls ssands1=2;1=2  p1=2;1=2  s1=2;1=2  p1=2;1=2 ,0 Ap1=2;1=2  p1=2;1=2  with A being an antisymmetrizer,
and a function   r*1 expZr was used with * 
k2  Z=c2 1=2 that describes a mild singularity at thenucleus. We performed the free ICI using the g function,g  1  r1  r2  r12 , and therefore, we dropped out
from our expansion bases the functions that lead to diverging integrals. We used the inverse Hamiltonian method at
the diagonalization step to ensure the bound property. We
used W  c to avoid the singularity originating from thels part.
Table II shows the result. When we used the inverse
Hamiltonian at the diagonalization step, the ground stateenergy monotonically converged from above (actually, theinverse energy converged from below) toward the literature
value [8,12] for both Z  2 and 90. For Z  2, the energyof Ref. [12] is more accurate than that of Ref. [8]: theformer included r12 explicitly, but the latter did not. For
Z  90, the value cited [8] is due to the calculation notincluding r12 explicitly. Since we used the inverse
Hamiltonian, the present energy given in Table II should
be an upper bound of the true energy. When we used theregular Hamiltonian, instead of the inverse Hamiltonian, atthe diagonalization step, the variational collapse occurred
for both Z  2 and 90.
Conclusion.—We have proposed in this Letter a methodof solving the DCE for atoms and molecules in an analytical expansion form. The relativistic ICI automatically generates the wave function having the exact structure. The
exact ICI balancing is done in the iteration process. For the
one-electron atoms, we experienced no variational problem as far as we used ICI. For the two-electron atoms, we
could avoid the variational collapse by introducing theinverse Hamiltonian. Thus, we propose the relativistic
free ICI method based on Eqs. (11) and (12) and the inverse
variational method having the Ritz-like property as a general method of calculating the exact solutions of the DCE
of many-electron atoms and molecules. This Ritz-like

TABLE II. Relativistic free ICI energy (a.u.) of the heliumlike
atoms (Z  2 and 90) using the inverse Hamiltonian at the
diagonalization step.
na

Mb

Energy (He)

Energy (Th88 )

0
1
2
3
4
5
6
Ref. [8]Ref. [12]

3
12
38
105
224
482
852

2:75 008 563
2:88 771 973
2:90 307 277
2:90 375 350
2:90 384 265
2:90 385 116
2:90 385 405
2:90263 6
2:90385 7

9166:575 433
9166:809 415
9166:858 084
9166:881 423
9166:894 218
9166:903 050
9166:907 227
9166:927 2aIteration number.
Number of the independent functions for expanding . Moredetails of M (Mll , Mls , Mss ) are 31; 1; 1, 124; 4; 4, 3812;
13; 13, 10529; 35; 41, 22466; 69; 89, 482145; 139; 198, and852257; 223; 342.
bproperty is useful in actual calculations where we do notknow the true energy E0 .
Our method of calculating the analytical exact wave
functions has been confirmed to have a high potentiality
for solving the SE of atoms and molecules [4 –6] andextended here to solve the relativistic DCE. We have still
many things to do in future to develop this new methodology as a really useful method in quantum chemistry.
This study has been supported financially by a Grant for
Creative Scientific Research from the Ministry of
Education, Science, Culture, and Sports of Japan.

[1] P. A. M. Dirac, Proc. R. Soc. A 123, 714 (1929).
[2] H. Nakatsuji, J. Chem. Phys. 113, 2949 (2000);
H. Nakatsuji and E. R. Davidson, ibid. 115, 2000 (2001);
H. Nakatsuji, ibid. 115, 2465 (2001); H. Nakatsuji, ibid.
116, 1811 (2002); H. Nakatsuji and M. Ehara, ibid. 117, 9(2002); H. Nakatsuji and M. Ehara, ibid. 122, 194108(2005).
[3] H. Nakatsuji, Phys. Rev. A 65, 052122 (2002).
[4] H. Nakatsuji, Phys. Rev. Lett. 93, 030403 (2004).
[5] H. Nakatsuji, H. Nakashima, and Y. Kurokawa (to be
published).
[6] H. Nakatsuji, Bull. Chem. Soc. Jpn. (to be published).
[7] G. Pestka and J. Karwowski, in Explicitly Correlated
Wave Functions in Chemistry and Physics — Theory
and Applications, edited by J. Rychlewski (Kluwer,
Dordrecht, 2003), pp. 331–346.
[8] A. Kolakowska, J. D. Talman, and K. Aashamar, Phys.
Rev. A 53, 168 (1996).
[9] R. N. Hill and C. Krauthauser, Phys. Rev. Lett. 72, 2151
(1994).
[10] G. Pestka, Phys. Scr. 69, 203 (2004).
[11] J. J. Sakurai, Advanced Quantum Mechanics (AddisonWesley, Reading, MA, 1967).
[12] A. Kolakowska, J. Phys. B 30, 2773 (1997).

050407-4