VOLUME 76, NUMBER 7

PHYSICAL REVIEW LETTERS

12 FEBRUARY 1996

Direct Determination of the Quantum-Mechanical Density Matrix Using the Density Equation
Hiroshi Nakatsuji* and Koji Yasuda†
Department of Synthetic Chemistry and Biological Chemistry, Faculty of Engineering, Kyoto University, Kyoto 606, Japan(Received 9 August 1995)
With the use of the density equation [H. Nakatsuji, Phys. Rev. A 14, 41 (1976)], the second-orderdensity matrices are directly determined without any use of the wave functions. The third- and fourthorder reduced density matrices (RDM’s) are decoupled into lower-order ones using the Green’s function
technique. This method is applied to Be, Ne, H 2O, H3 O1 , NH 3, CH 4, BH42 , NH41 , and CH 3F, andthe results are successfully compared with the full configuration interaction results. The convergenceis fairly good, and the calculated second-order RDM’s almost satisfy the necessary conditions of the
N representability, the P, Q, and G conditions, and the first-order RDM’s are exactly N representable.
These results show that the present method is very promising.
PACS numbers: 31.10.+z, 02.70.Rw, 03.65.Ge, 83.10.–y

It is known that the many-electron wave function
involves more information than we need to know. Therefore the determination of density or density matrices
without using the wave function could be a convenient
alternative to wave mechanics, and along this line thedensity functional theories [1,2] and density matrix approaches [3], etc., have been presented.
The Hohenberg-Kohn theorem [1] ensures the use of
the electron density as a basic variable. Recently, thisdensity functional theory has been much developed andwidely used [1,2]. However, the Hohenberg-Kohn theorem is an existence theorem, and the explicit functional
form is still not known [4].

EG

snd

­

(n
X

n
X

)

Z

(

n
X

)

wsi, jd G 1 sn 1 1d
ysn 1 1d 1
wsi, n 1 1d G sn11d dxn11
i
i.j
i
Z
1
sn12d
1
dxn11 dxn12 ,sn 1 1d sn 1 2d wsn 1 1, n 1 2dG
2ysid 1snd

All of the important physical quantities are obtained
from the second-order density matrix. Since it determines
the energy, we can apply the variational principle in a suitable domain of the density matrices [3]. However, it isstill not completely known what conditions the Pauli principle enforces on the density matrices (N representability condition) [5]. Moreover, since these conditions may
be very complicated, it is impracticable to carry out suchvariational calculations.
On the other hand, a nonvariational approach for the direct determination of the density matrix has been proposedby one of the present authors [6]. He showed that the equation

which is called the density equation, is equivalent to the
Schrödinger equation for each n with n $ 2. That is,
this equation in the domain of the N-representable density
matrices is the necessary and sufficient condition for thecorresponding wave function to satisfy the Schrödinger
equation. The necessity alone was shown by Cho [7]and by Cohen and Frishberg [8]. However, unless theN representability condition is not known, the number of
unknowns exceeds the number of conditions [9], since
the nth-order density equation contains nth-, sn 1 1dth,
and sn 1 2dth-order density matrices. On the other hand,
if we can calculate in some way or even approximately
the sn 1 1dth and sn 1 2dth-order density matrices fromthe nth- and lower-order ones, it becomes possible tosolve the density equation and calculate the density matrix
directly without any use of the wave function. The purposeof this communication is to show that such a method is now
almost established.
Recently, Valdemoro et al. have suggested a decoupling approximation of the higher-order density matrices
0031-9007y 96y76(7)y1039(4)$06.00(1)

in terms of the lower-order ones, based essentially on thefermion’s anticommutation relations [10], and Colmenero
and Valdemoro applied it, using the density equation, tosome four-electron systems [11]. However, as we shallshow later, this method seems to have limitations for more
complex atoms and molecules. We have derived more accurate decoupling approximations of the third- and fourthorder density matrices in terms of the lower-order ones
using the Green’s function technique, and applied them,using the density equation, to some atoms and molecules.
In this communication we briefly outline the method andgive the results of the calculations, comparing them withthe traditional wave function approach. We actually use
the second-order density equation sn ­ 2d.
The density matrix (DM) G snd is related to the Green’s
function G snd as
s2idnn!
3 G snd sx10 02 · · · xn0 02 jx1 01 · · · xn 01 d ,G snd sx10 · · · xn0 jx1 · · · xn d ­© 1996 The American Physical Society

1039

VOLUME 76, NUMBER 7

PHYSICAL REVIEW LETTERS

and the reduced density matrix (RDM) n D is defined byXn
Dsr10 · · · rn0 jr1 · · · rn d ­G snd sx10 · · · xn0 jx1 · · · xn d ,s1 ···sn

where si is a spin variable. Taking lower-order perturbation series of the Green’s functions, we introduce thefollowing decoupling approximations for the third- andfourth-order Green’s functions:

12 FEBRUARY 1996

The UV term in the 3-DM is defined by
UV

1
i1,i2,i3
G s3dj1,j2,j3­√occupied
X

6

!

2

virtualXkk√occupiedX!

2virtualX

k

Usx10 x20 jx1 x2 d 2 Usx10 x20 jx2 x1 d ­ 2G s2d sx10 x20 jx1 x2 d
ÉÉ
G s1d sx10 jx1 dG s1d sx10 jx2 d,
2 s1d 0G sx2 jx1 dG s1d sx20 jx2 d(4)where the second term in the rhs is the second-order DM(2-DM) in the independent-particle approximation [2,6].
The 2P term in the 4-DM is defined by
2P

i1,i2,i3,i4
G s4dj1,j2,j3,j4­

1 i1,i2 i3,i4
Uj1,j2 Uj3,j4 1 · · · ,24where the sum is for all the permuted, antisymmetrized
products of U and consists of 72 terms.
TABLE I.

¥
i1,k
i2,i3
Uj1,j2 Vk,j3 1 · · · .

The U term is defined by Eq. (4), and the V term is thesolution of the linear equation
i1,i2
i1,i2
Vj1,j2 ­ Uj1,j2 1

In Eqs. (2) and (3), we give only the representative graphs,
but we include all the related terms including the exchange
terms. The first two terms of Eq. (2) and the first threeterms of Eq. (3) are the same as those given by Valdemoro
et al., and correspond approximately to the first-orderterms. The last terms of Eqs. (2) and (3), which we call
UV and 2P (two-pair) terms, respectively, are the secondorder terms in electron correlations. These terms are more
explicitly given below, and their physical meanings areclarified.
First we define the collision term U by

≥

≥
¥
s1dk
s1dk
i1,i2
Vj1,k Gj2 2 GHFj2 ,kwhere HF stands for Hartree-Fock. The UV term represents the following process: A particle (or a hole) 2 firsti1,kcollides with 1 sUj1,j2 d, moves some distance, then collides
i2,i3with 3 sVk,j3 d. The sum of all the intermediate states kyields the total amplitude, and the sign reflects whether aparticle or a hole propagates as an intermediate state. The
second term of the linear equation represents the subtraction of some overcounted diagrams. Note that UV term is
not the simple product of the lower-order density matrices.
Since the Green’s function method is quite general,
Valdemoro’s decoupling approximation was derived asa special case. Valdemoro’s procedure, which is quite
interesting, is based on the exact relation between the nRDM’s n D and the n-hole RDM’s (n-HRDM’s) n D, forexample,
j1,j2

i1,i2

1

j1

j2

1

j2

j1

D i1,i2 2 Dj1,j2 ­ 2 sD i1 D i2 2 2 D i1 D i2 d
1

1

i1 i2
i2 i1
2 2 sDj1
Dj2 2 2 Dj1
Dj2 d .

Valdemoro et al. put the RDM (HRDM) term in the lefthand side to be equal to the RDM (HRDM) terms on
the right-hand side. Hence this approximation has thepossibility of neglecting some terms that are common toboth RDM and HRDM. In addition, their approximation
represents the higher-order RDM’s by the simple productsof the lower-order ones, which is generally impossible.

Errors of the various approximations of the 3- and 4-RDM’s for the ground state of Be.
3-RDM

4-RDM

Element aj123 , m123

Valdemoro
error

1UV
error

Full-CIvalue

Element aj1234 , m1234

Valdemoro
error

12P
error

Full-CIvalue

2,3,3,2,3,3
2,2,3,2,2,3
2,2,3,1,1,3
1,1,3,1,1,3
1,3,3,1,3,3
3,3,2,2,3,3
3,2,2,2,2,3

3.356 3 1024
23.174 3 1024
22.233 3 1024
22.015 3 1024
1.752 3 1024
21.678 3 1024
1.587 3 1024

7.250 3 1027
24.751 3 1027
22.552 3 1027
22.814 3 1027
4.459 3 1027
23.625 3 1027
2.376 3 1027

1.758 3 1024
2.029 3 1024
2.235 3 1024
3.187 3 1024
3.364 3 1024
28.789 3 1025
21.014 3 1024

3,3,2,2,3,3,2,2
3,3,2,2,3,3,1,1
3,3,2,2,2,3,2,3
3,3,1,1,3,3,1,1
4,4,2,2,3,3,1,1
4,4,3,3,2,2,1,1
4,3,2,2,3,3,1,1
3,3,2,2,1,3,1,3
4,4,1,1,3,3,1,1

27.928 3 1025
25.581 3 1025
3.964 3 1025
23.922 3 1025
23.740 3 1025
23.652 3 1025
2.973 3 1025
2.791 3 1025
22.628 3 1025

2.931 3 1027
1.169 3 1027
21.466 3 1027
1.779 3 1027
1.869 3 1028
21.731 3 1027
23.329 3 1028
25.844 3 1028
6.203 3 1028

3.933 3 1025
5.591 3 1025
21.967 3 1025
7.948 3 1025
3.742 3 1025
3.657 3 1025
22.977 3 1025
22.796 3 1025
1.631 3 1026a

All the 3- and 4-RDM elements with the Valdemoro errors larger than 1.5 3 1024 and 2.5 3 1025 , respectively, are listed.

1040

VOLUME 76, NUMBER 7

PHYSICAL REVIEW LETTERS

12 FEBRUARY 1996

TABLE II. Calculated energies and density errors in the density-equation method and the wave function approach.
Density-equation method
Valdemoro
SystemActive space
Number ofelectrons

Be

232
4

Ne

336
6

H2 O

532
10

not conv.

NH3

533
10

not conv.

H 3 O1533
10

not conv.

CH4

434
8

BH2
4

434
8

NH1
4

434
8

CH3 F

734
11

240.12426
82.4
2.948 3 1021
226.80563
99.5
3.198 3 1021
256.42354
72.2
2.912 3 1021
not conv.

present

Wave function approachaHartree-Fock

SDCI

Full-CI

214.58269
0.028
1.141 3 1024
2128.57726
1.33
6.943 3 1023
275.72829
1.44
1.030 3 1022
256.05344
2.53
2.039 3 1022
276.09706
2.18
1.793 3 1022
240.18772
3.45
2.819 3 1022
226.87534
2.58
2.831 3 1022
256.48087
3.14
2.757 3 1022
2138.52491
4.91
6.118 3 1022

214.58269
002128.57835
00275.72902
00256.05513
00
276.09833
00240.19049
00226.87719
00256.48347
002138.53004
00Energycorrelation energy error(%) b2
D error c
214.59775
2106.4
1.165 3 1022
not conv.

214.58269
0.014
4.823 3 1025
128.57862
20.3311.373 3 1022
276.72751
3.00
7.122 3 1023
256.05293
3.30
6.231 3 1023
276.09629
3.55
9.977 3 1023
240.18726
4.02
9.757 3 1023
226.87450
3.75
1.072 3 1022
256.47990
4.31
1.164 3 1022
2139.52339
6.37
4.514 3 1022

214.56853
100
4.631 3 1022
2128.49637
100
1.947 3 1021
275.67884
100
3.154 3 1021
255.98860100
3.473 3 1021
276.03996
100
3.067 3 1021
240.11015
100
3.702 3 1021
226.80526
100
3.813 3 1021
256.40047
100
3.671 3 1021
2138.42548
100
5.614 3 1021a

Present approximation corresponds to Valdemoro 1UV for 3-RDM and Valdemoro 12P for 4-RDM.
Error percentage in the correlation energy.
cThis error means the square norm of the differences between the calculated 2-RDM and the full-CI one.
b

Note that we have shown previously [6] that, by putting
the independent particle model approximation∂µ
1
1 i2 i1
i1,i2
i1 i2
Dj1,j2 ­Dj2 2
Dj1
Dj1 Dj2
2
2
into the density equation, we get the HF equation as aresult.
Based on the above decoupling approximation of the3- and 4-RDM’s, we solved the second-order density
equation [Eq. (1) with n ­ 2] iteratively. The HF orbitals
were used as one-electron basis and the HF 2-RDM wasused as an initial guess. The iterative procedure is asfollows: We first symmetrize and normalize a trial 2RDM, then calculate the energy and the 3- and 4-RDM’sfrom the 2-RDM, and substitute them in the density
equation. Since only the exact density matrices satisfy
it, its right-hand side will differ from the left-hand sidein general. From this residue we calculate a new 2RDM by the Newton method. The above procedure isrepeated until convergence is obtained. In some cases,
interpolation techniques are necessary for an efficientconvergence: The previous and the present 2-RDM’s areaveraged with some given weight.
We applied the present density-equation method toseveral atoms and molecules: Be, Ne, H2 O, H3 O1 , NH3 ,CH4 , BH 42 , NH41 , and CH3 F. The comparative wave
function approach was used at the HF, configurationinteraction SDCI, and full CI levels. The following
basis sets were used: For atoms, the double-zeta Slatertype orbitals (STO’s) [12] were used and expanded by
six Gaussian-type orbitals (GTO’s) [13]. For molecules,
the minimal STO-6G basis [13] and the experimental
molecular geometries [14] were used. The 1s and 2s
orbitals of Ne and the 1s orbitals of CH4 , BH 42 , NH41 ,and CH3 F were frozen as cores.
First, we examine the accuracy of the present decoupling approximations as given by Eqs. (2) and (3). We
calculated the 3-RDM from the exact 1- and 2-RDM’s

1041

VOLUME 76, NUMBER 7

PHYSICAL REVIEW LETTERS

and the 4-RDM from the exact 1-, 2-, and 3-RDM’s. The
exact RDM’s were calculated from the full CI wave function. Table I shows the results for Be. In going from
Valdemoro’s approximation to the present one, the errors
are reduced from the order of 1024 to 1027 for the 3RDM, and from the order of 1025 to 1027 for the 4-RDM.
Thus our approximations significantly improve the Valdemoro’s ones.
We next compare the results of the density-equation
method at different levels of approximations. Table IIshows the results of the energies and the 2-RDM’s forthe calculated systems. In comparison with the exact fullCI results, the Valdemoro’s approximation shows poor
results, while our approximations for the 3- and 4-RDM’s
improve over it remarkably. Some calculations with the
Valdemoro’s approximation were difficult to converge,
while our approximation led to convergence within a few
iterations.
Recently, Colmenero and Valdemoro have applied thedensity equation, using the Valdemoro’s approximations
for the 3- and 4-RDM’s, to four-electron systems like Be,
B1 , etc. [11]. They further used corrections making the
RDM’s spin-symmetry adapted and keeping the diagonal
elements nonnegative. The error in the correlation energy
for Be was 0.27%, which is comparable to our 1UV 1
2P term result. Their special corrections seem to work
well for the four-electron atoms. Note, however, that for
four-electron systems the 4-RDM is essentially equal to
the wave function itself.
On the other hand, our method including the UV and 2P
terms works well even for many-electron systems having
more than four electrons. This is in marked contrast to
the Valdemoro’s approximation which either gives very
poor results or does not even give a convergence in thecalculations. For atoms, the error in the correlation energy
is only 0.014% for Be and 20.33% for Ne (the minus
sign means overshooting). For molecules, it is 3–4 % forhydrides and 6% for CH3 F. Although these energy errors
are larger than those of the SDCI method except for theatoms, the errors of the 2-RDM’s measured by the square
norms of the deviations from the exact ones are smaller byan order of magnitude than those of the SDCI method. Weinterpret this as follows: The SDCI method is a variational
method, while the present density-equation method is anonvariational method directly calculating the 2-RDM’s.
We finally checked the N representability of the calculated 2-RDM’s using some necessary conditions, P, Q,
and G conditions [5], which are the nonnegatives of theRDM, the HRDM, and the G matrix. We found that our
2-RDM’s almost satisfy these N representability conditions. We also found that all the eigenvalues of the 1RDM’s are in the range of zero to 2, so that the calculated
1-RDM’s are exactly N representable [15].
The present study is summarized as follows. Withthe use of the density equation the second-order density
matrices were directly determined without any use of the1042

12 FEBRUARY 1996

wave function. The 3- and 4-RDM’s were decoupled into
lower-order ones using the Green’s function technique.
The results for the atoms and molecules, Be, Ne, H2 O,
H3 O1 , NH3 , CH4 , BH42 , NH1
4 , and CH3 F were quite
accurate and promising. The convergence was fairly
good, and the calculated 2-RDM’s almost satisfy somenecessary conditions of N representability, while the 1RDM’s are exactly N representable. These results show
this method to be very promising. We are now exploring
both more accurate and more efficient methods. Moredetails of the theory and the calculations will be publishedin near future.
This study was supported in part by a Grant-in-aid
for Scientific Research from the Japanese Ministry of
Education, Science and Culture. K. Y. acknowledges thefellowship of the Japan Society for the Promotion of
Science.

*To whom correspondence should be addressed, and alsoaffiliated with the Institute for Fundamental Chemistry,
34-4, Takano-Nishihikari-cho, Sayko-ku, Kyoto 606,
Japan.
†
Research Fellow of the Japan Society for the Promotionof Science for Young Scientists.
[1] P. Hohenberg and W. Kohn, Phys. Rev. B 864, 136(1965); W. Kohn and L. J. Sham, Phys. Rev. A 1133, 140(1965).
[2] R. G. Parr and W. Yang, Density-Functional Theory ofAtoms and Molecules (Oxford University Press, New
York, 1989).
[3] C. Garrod, M. V. Mihailovic, and M. Rosina, J. Math.
Phys. 16, 868 (1975).
[4] H. Nakatsuji and R. G. Parr, J. Chem. Phys. 63, 112(1975); H. Nakatsuji, J. Chem. Phys. 67, 1312 (1977).
[5] A. J. Coleman, Rev. Mod. Phys. 35, 668 (1963); C. Garrod
and J. Percus, J. Math. Phys. 5, 1756 (1964); H. Kummer,
J. Math. Phys. 8, 2063 (1967).
[6] H. Nakatsuji, Phys. Rev. A 14, 41 (1976).
[7] S. Cho, Sci. Rep. Gumma Univ. 11, 1 (1962).
[8] L. Cohen and C. Frishberg, Phys. Rev. A 13, 927 (1976).
[9] E. Davidson, private communication to one of the authors(H. N.) at Tomakomai in 1976; J. E. Harriman, Phys. Rev.
A 19, 1839 (1979).
[10] C. Valdemoro, Phys. Rev. A 45, 4462 (1992); F. Colmenero, C. Perez del Valle, and C. Valdemoro, Phys. Rev. A
47, 971 (1993); F. Colmenero and C. Valdemoro, Phys.
Rev. A 47, 979 (1993).
[11] F. Colmenero and C. Valdemoro, Int. J. Quant. Chem. 51,
369 (1994).
[12] E. Clementi and C. Roetti, Atom Data Nucl. Data Tab. 14,428 (1974); the exponents were reoptimized by T. Koga
and S. Watanabe (private communication).
[13] W. J. Hehre, R. F. Stewart, and J. A. Pople, J. Chem. Phys.
51, 2657 (1969).
[14] L. E. Sutton, D. G. Jenkin, and A. D. Mitchell, Tables ofInteratomic Distances (The Chemical Society, London,
1958).
[15] T. L. Gilbert, Phys. Rev. B 12, 2111 (1975).