PHYSICAL REVIEW A

VOLUME 56, NUMBER 4

OCTOBER 1997

Direct determination of the quantum-mechanical density matrix using the density equation. II.
Koji Yasuda* and Hiroshi Nakatsuji†
Department of Synthetic Chemistry and Biological Chemistry, Faculty of Engineering, Kyoto University, Kyoto 606, Japan~Received 7 October 1996!
With the use of the density equation @Phys. Rev. A 14, 41 ~1976!#, second-order density matrices of atoms
and molecules are calculated directly without any use of the wave function. The decoupling of the third- andfourth-order density matrices is done using the diagrammatic procedure of the Green’s-function method,
including terms up to second order in electron correlation. In addition to the results given in a previous
communication @Phys. Rev. Lett. 76, 1039 ~1996!#, we give more detailed formulations, discussions, andadditional results for H2O, NH3, CH4, HF, N2, CO, and acetylene using double-z basis sets, and for CH3OH,
CH3NH2, and C2H6 ~staggered and eclipsed! using minimal basis sets. The present density-equation method
gave energies as accurate as, and reduced density matrices ~RDM’s! more accurate than the single and double
excitation configuration interaction method. The present method seems to give better quality results as thesystem becomes large. The convergence was fairly good and the calculated second-order RDM’s almost
satisfied some necessary conditions of the N representability, the so-called P, Q, and G conditions, while thefirst-order RDM’s were exactly N representable. The variety of the molecules calculated and the quality of thecalculated results show that the density-equation method can be a promising alternative to the wave-function
approach in quantum mechanics. @S1050-2947~97!06409-3#
PACS number~s!: 31.10.1z, 02.70.Rw, 03.65.Ge

I. INTRODUCTION

Since all the operators we shall concern ourselves with inquantum mechanics are one- and two-body ones, all the elemental physical quantities can be calculated from thesecond-order density matrix ~2-DM!. Many-electron wave
functions involve more information than we need to know.
Thus a determination of the density or density matrices without using the wave function could be a convenient alternative
to the wave mechanics, and this approach is often calledH(n

EG ~ n ! 5inv~ i ! 1

(
i. j

J

w ~ i, j ! G ~ n ! 1 ~ n11 !

1 21 ~ n11 !~ n12 !

E‘‘wave mechanics without wave.’’ Since the 2-DM determines the energy, we may apply the variational principle in asuitable domain of the density matrices @1#. However, it isstill not completely known what conditions the Pauliprinciple enforces on the density matrices ~theN-representability condition! @2#. Moreover, since these conditions may be very complicated, it would be impractical tocarry out such variational calculations.
One of the authors proposed a nonvariational approach fora direct determination of the density matrix @3#. He showedthat the equationEHnv~ n11 ! 1

J(i w ~ i,n11 ! G ~ n11 !dx n11

w ~ n11,n12 ! G ~ n12 ! dx n11 dx n12 ,~1.1!

*Present address: Graduate School of Human Informatics,
Nagoya University, Chikusa-ku, Nagoya 464-01, Japan.
†
Corresponding author. Also affiliated with the Department of Applied Chemistry, Graduate School of Engineering, The Universityof Tokyo, Hongo, Tokyo 113, Japan, and the Institute for Fundamental Chemistry, 34-4, Takano-Nishihikari-cho, Sakyo-ku, Kyoto
606, Japan.

function to satisfy the Schrödinger equation. The necessity
alone was shown by Cho and by Cohen and Frishberg @4#.
Unfortunately, since the nth-order density equation contains
nth, (n11)th, and (n12)th order DM’s, that is, the numberof the unknowns exceeds the number of the conditions as far
as the N representability condition is not completely known,
the solution of the density equation itself is not unique @5#.
However, the solution within the N-representable space isguaranteed to be exact, and we use the subsidiary
N-representability conditions to eliminate unphysical solutions.
On the other hand, if we can approximate (n11)th- and(n12)th-order DM’s using the nth- and lower-order ones,
we can directly determine the density matrix from the den-1050-2947/97/56~4!/2648~10!/$10.00

2648

which is called the density equation, is equivalent to the
Schrödinger equation in the domain of N-representable
DM’s for each n with n>2. That is, this equation in thedomain of the N-representable space is a necessary and sufficient condition for the corresponding antisymmetric wave

56

© 1997 The American Physical Society

56

DIRECT DETERMINATION OF THE QUANTUM- . . .

sity equation without using the wave function. Actually, theequivalence of the Schrödinger equation plus Pauli principlewith the density equation plus N representability condition@3# implies that the N representability condition has a powerof expressing G (n12) and G (n11) in terms of G (n) . Thus, theabove approximation also works as the N-representability
condition. In principle, the ground-state nth-order DM determines the DM’s of all orders including the wave function,
since it determines the electron density and hence the wave
function ~Hohenberg-Kohn theorem @6#!.
Recently, Valdemoro and co-workers @7,8# reported aninteresting approach for solving the density equation. Theysuggested a decoupling approximation of higher-order reduced density matrices ~RDM’s! in terms of the lower-orderones based essentially on the fermion’s anticommutation relations. They derived the identity containing the n-RDM andn-hole RDM ~HRDM! on the left-hand side ~LHS!, and thelower-order RDM’s and HRDM’s on the right-hand side
~RHS!, and then assumed that the n-RDM ~n-HRDM! on thelhs is equal to the sum of the products of the lower-order
RDM’s ~HRDM’s! on the rhs. We call this approximation
the IPH approximation ~the approximation identifying theparticle and hole parts separately!. With this approximation,
Colmenero and Valdemoro solved the density equation @9#
for some four- and six-electron systems, and calculated theirdensity matrices directly @10#. However, their results seem tohave been limited to small systems, and the applicability togeneral many-electron systems is an open question. Since
their decoupling approximation is accurate for a Be atom
~four-electron system! but less accurate for H2O ~ten-electron
system! @8#, we have to look for a more accurate approximation which is suitable for studying more complex atoms andmolecules.
In a previous paper @11#, we proposed a different approach for solving the density equation: we presented thedecoupling approximation based on the Green’s-function
technique, and applied it using the second-order density
equation @Eq. ~1.1! with n52# to directly determine the2-RDM’s of several atoms and molecules. We found that our
method is very promising, giving as accurate energies as, andmore accurate density matrices than, those of the variational
SDCI ~single and double excitation configuration interaction!
one. We checked some N-representability conditions, andfound that the resultant 2-RDM’s almost satisfied these conditions, while the 1-RDM’s were exactly N representable.
This success opened a possibility of an entirely new method
in quantum mechanics, and gave us some insight into the
N-representability conditions.
In this paper we give more detailed explanations of theapproximations for expressing the 3- and 4-RDM’s in termsof the 1- and 2-RDM’s, which were discussed briefly in theprevious paper @11#, show some new results, and discuss theconvergence properties of our method.

II. APPROXIMATIONS OF 3- AND 4-RDM’S

We use the Green’s-function ~GF! method to derive approximate relations of the higher-order DM’s with the lowerorder ones. The IPH approximation is derived as a specialcase. The DM’s are defined as @12#2649

G ~ n ! ~ x 18 •••x n8 u x 1 •••x n !
51^ f † ~ x 1 ! ••• f † ~ x n ! f ~ x n8 ! ••• f ~ x 18 ! & .
n!
~2.1!

where x i denotes space-spin coordinate of i-th electron,
while f † and f denote creation and annihilation field operators, respectively. The n-RDM’s are defined asn

D ~ r 81 •••r 8n u r 1 •••r n !

5(

s1 . . . sn

G ~ n ! ~ x 81 . . . x 8n u x 1 . . . x n !

~2.2!

where r i and s i denote space and spin coordinates, respectively. The n-particle many-body GF’s are defined as @13#
G ~ n ! ~ x 81 t 81 •••x 8n t 8n u x 1 t 1 •••x n t n !
5 ~ 2i ! n ^ T @ f ~ x 81 t 81 ! ••• f ~ x 8n t 8n ! f † ~ x n t n ! ••• f † ~ x 1 t 1 !# &~2.3!
where T denotes the time ordering operator. The DM’s are
related to the GF’s as
G ~ n ! ~ x 18 •••x n8 u x 1 •••x n !
5~ 2i ! n ~ n !
G ~ x 18 0 2 •••x n8 0 2 u x 1 0 1 •••x n 0 1 ! ,
n!

~2.4!

where 0 1 and 0 2 denote positive and negative infinitesimals, respectively. They determine the time ordering of theoperators in Eq. ~2.3! so as to match the definition of theDM’s given by Eq. ~2.1!.
Later we use matrix notations of these quantities, where
the subscript j k and the superscript i k are associated with theannihilation and creation operators, respectively, and theone-particle basis function is a spin-orbital for the DM and a
spatial orbital for the RDM.
In order to obtain the relations among the GF’s, we firstconsider the lower-order perturbation series, and then include some higher-order effects. We use Feynman’s procedure, which represents the GF as a sum of the topologicallydifferent diagrams @13#. After taking the limit of time variables to get the relations among the DM’s, we rewrite theminto the relations among the RDM’s by summing up the spin
variables. We use the Hartree-Fock ~HF! GF’s as unperturbed zeroth-order GF’s, and the electron correlation istreated as a perturbation. Hence our approximation works
best for closed-shell atoms and molecules, dealing mainlywith the dynamic correlations.

2650

KOJI YASUDA AND HIROSHI NAKATSUJI
A. Approximation of the 4-RDM

There are four kinds of diagrams for the four-particle GF
up to the second order in the perturbation. In a standard
diagrammatic notation, they are represented symbolically as
G~4!
5
~2.5!
in which the coefficient for the kth-order perturbation diagram of G (n) is 6i m2k2n , where m is the number of the HF
GF’s in the diagram, and the relative sign ensures the permutational antisymmetry @13#. With the use of the lowerorder relations which connect G (3) with G (2) and G (1) , andG (2) with G (1) , it is easy to show that the sum of the first,
second, and third terms of Eq. ~2.5! reproduces the IPH approximation for the 4-RDM @8#: the last term missing in the
IPH approximation represents the simultaneous collisions of
two electron pairs, so that we call it two-pair ~2P! term. It isa central term in the independent pair model which is a good
approximation for describing the dynamic correlations in atoms and molecules @14#. The third term in Eq. ~2.5! will beconsidered in the approximation of the 3-RDM.
The 2P term is calculated as follows: for a given 2-DM,
the collision term U is defined as
U ~ x 18 x 28 u x 1 x 2 ! 2U ~ x 18 x 28 u x 2 x 1 !
52G

~2!

~ x 81 x 82 u x 1 x 2 ! 2U56

becomes a unit operator d ss 8 in a spin space. The spin sum
yields a factor of 2 42p for each permuted diagram of
4-RDM, where p is a parity of the permutational group. For
i3,i22
example, the spin sum of U i1,i4
j1,j3 U j2,j4 yields a factor of 2 .
In the present approximation, the 4-RDM includes all fourkinds of diagrams in Eq. ~2.5! and some higher-order terms,
because U is calculated from the 2-DM itself. In the direct
determination of the 2-RDM described in Sec. III, this approximate 4-RDM is used to calculate the rhs of the density
equation ~1.1!, then this equation is used to calculate a new
2-RDM. Through this self-consistency process, the 2-RDM,
U, and the approximate 4-RDM come to include somehigher-order terms.
The lowest-order missing terms in the present approximation of the 4-RDM are represented symbolically as~2.9!

which are third-order and are not represented by the simple
products of the lower-order RDM’s.
B. Approximation for the 3-RDM

U

G ~ 1 ! ~ x 81 u x 1 !

G ~ 1 ! ~ x 81 u x 2 !

G ~ 1 ! ~ x 28 u x 1 !

G ~ 1 ! ~ x 28 u x 2 !

.

There are three kinds of diagrams of the three-particle GF
up to second order in the perturbation, and they are represented symbolically as~2.6!
G ~ 3 !5Except for a coefficient, it is represented as
U ~ x 81 x 82 u x 1 x 2 ! 5~2.7!

in which the bold line is the exact one-particle GF, while therectangle is the exact vertex of G (2) @13#. In the above figure
we replaced the time variables in the external lines with positive and negative infinitesimals. The term U represents thesum of the collisions of two electrons which cannot be written by simple products of the 1-DM’s. Similarly, the 2P term
represents independent collisions of two electron pairs, andis written as the sum of the permuted, antisymmetrized products of U:
2P

G

~4!

1
5 U i1,i2
U i3,i4 1••• .
4! j1,j2 j3,j4~2.10!

With the use of the lower-order relations, we can show thatthe sum of the first and second terms of Eq. ~2.10! reproduces the IPH approximation for the 3-RDM. The last termof Eq. ~2.10! is second order in the perturbation of the electron correlation, and may not be negligible. To approximate
it, including not only the lowest-order term but also somehigher-order terms, we use Dyson’s equation and the property of the unperturbed GF’s.
An exact one-particle GF satisfies the Dyson’s integral
equation~2.11!

~2.8!

To generate all the terms in Eq. ~2.8!, we substitute fouri3,i4
indices ~i 1 , i 2 , i 3 , and i 4 ! of three terms U i1,i2
j1,j2 U j3,j4 ,
i1,i3 i2,i4
i1,i4 i2,i3
U j1,j3 U j2,j4 , and U j1,j4 U j2,j3 in all possible ways ~4!!. Then
we obtain the topologically different diagrams for the 2P
term whose total number is 72.
Since our Hamiltonian is spin independent, and since ourpresent subjects are closed-shell atoms and molecules, G (1)where an open circle means the irreducible self-energy which
is an energy-dependent effective one-body potential due tothe electron interaction @13#. The second term then represents
the multiple scattering by the self-energy.
The HF GF which we use as the zeroth order unperturbed
one is written as @13#

56

DIRECT DETERMINATION OF THE QUANTUM- . . .

With the use of Dyson’s equation ~2.11! and the HF GF’s
property ~2.16!, the last term of Eq. ~2.10! is transformed as

G ~01 ! ~ x 8 t 8 ,xt !
occupied

51i

2651

(k w *k ~ x ! w k~ x 8 ! exp@ 2i e k~ t 8 2t !# ;t 8 <tunoccupied

2i

(k

w*
k ~ x ! w k ~ x 8 ! exp@ 2i e k ~ t 8 2t !# ;t 8 >t,~2.12!

where w k and e k are the HF orbital and orbital energy, respectively. Since this unperturbed GF describes the probability amplitude of each independent electron ~hole! from xt to
x 8 t 8 ~and vice versa!, it satisfies the rule of the probability
amplitude on the subsequent events in the path-integral
theory @15#. Hence it has the propertyEdyG ~01 ! ~ zt 8 ,y0 ! G ~01 ! ~ y0,xt !

~2.17!
where the approximation means that it is valid only when
tt 8 ,0. The case tt 8 .0 will be considered later. The term
expressed by Eq. ~2.17! is called the UV term. The only
unknown in the above figure which we call the V term issimilarly approximated as

5 ~ 2i ! G ~01 ! ~ zt 8 ,xt ! ; t 8 >0>t
~ 1i ! G ~01 ! ~ zt 8 ,xt ! ; t 8 <0<t,~2.13!

in which we integrate all the intermediate states, y0. The
space integration of the two GF’s gives one GF in these timeorderings. From the definition of the HF GF, Eq. ~2.12!, only
the occupied ~unoccupied! orbitals contribute to G (1)0 whent 8 ,t (t 8 .t). Since these two kinds of orbitals are orthonormal to each other, the relative sign in Eq. ~2.13! can betreated with the use of the projector,

S(

unoccupied

P ~ y ! 5ikoccupied2(kD

After manipulating some numerical coefficients, the UV
term given by diagram ~2.17! is written into an ordinary
formula asUVu w k ~ y ! &^ w k ~ y ! u . ~2.14!

tt 8 ,0,

i1,i2
i1,i2 l2
~ 1 ! k2
~ 1 ! k2
V i1,i2
j1,j2 5U j1,j2 1V j1,l2 P k2 ~ G j2 2G 0 j2 !

tt 8 ,0,5U1V P g .

~2.15!

which we represent in a diagrammatic form as;~2.19!

From diagram ~2.18!, the V term is written as

dy G ~01 ! ~ zt 8 ,y0 ! P ~ y ! G ~01 ! ~ y0,xt !
5G ~01 ! ~ zt 8 ,xt ! ;l2 k2,i3
i1,i2,i3
G ~ 3 !j1,j2,j3
5 16 ~ V i1,i2
j1,l2 P k2 U j2,j3 1••• !

5 16 ~ V PU1••• ! .

Then Eq. ~2.13! becomesE~2.18!

~2.16!

where the triangle means the projector P. We note that Eqs.
~2.15! and ~2.16! are valid only when tt 8 ,0: the lhs’s of
Eqs. ~2.15! and ~2.16! vanish identically when tt 8 .0, while
their rhs’s do not. Note that only space integration for y isnecessary, in contrast to the ordinary space-time integrationin the GF theory.
It is well known that the exact GF has an expression similar to Eq. ~2.12! ~the Lehmann representation!. One-particle
functions which contribute to G (1) when t 8 .t (t 8 ,t) correspond to the quasiparticle wave functions of anion ~ionized!
states. However, since these functions are not orthonormal,
the exact G (1) does not have an expression like Eq. ~2.15!.

~2.20!

which is used to calculate the UV term in Eq. ~2.19!. Other
terms denoted as . . . in Eq. ~2.19! are generated by substituting the indices 1, 2, and 3 in all possible ways, and thenby substituting three indices ~i 1 , i 2 , and i 3 ! in all possible
ways. The projector P and the correction to the 1-DM g is
written with the unperturbed first-order density matrix G (1)0as
P lk 52G ~ 1 !0kl 2 d lk ,~2.21!

g ij 5G ~ 1 !ji 2G ~ 1 !0ji ,~2.22!

where d lk is Kronecker’s delta.
By substituting Eq. ~2.20! into Eq. ~2.19!, the UV term
UV (3)
G can be written in an infinite sum as
UV

1
G ~ 3 ! 5 ~ U PU1U P g PU1U P g P g PU1••• !
6~2.23!

2652

KOJI YASUDA AND HIROSHI NAKATSUJI

Since Eqs. ~2.19!–~2.22! do not contain orbital energies,
we can use any idempotent first-order density matrices
whose N eigenvalues are equal to 1, and the rest are equal to
0 as unperturbed ones ~N is the number of electrons!. In thecalculations below, the HF density matrix is used as an unperturbed one. Once we get the UV term of the 3-DM, thespin-sum yields the 3-RDM in the same way as the 4-RDM.
The lowest-order missing term in the present UV approximation is the second-order one, which is attributed to theapproximation used in Eq. ~2.17!: in the last term of Eq.
~2.10!, we omitted the term with two internal times t and t 8
being the same sign, tt 8 .0, because Eqs. ~2.15! and ~2.16!
do not hold in this case. This missing term contributes to theelements G (3)i1,i2,i3
j1,j2,j3 with j 1 – 3 are occupied and i 1 – 3 are unoccupied ~or j 1 – 3 are unoccupied and i 1 – 3 are occupied!. By
applying standard diagrammatic techniques and integrating
the time variables, the missing term is written in an explicit
form as
3 ! i1,i2,i3
5
DG ~j1,j2,j3

1
6D e

unoccupied

1
2
6D e

(k

occupied(k^ ki 2 u V 0 u j 1 j 2 &^i i uV uk j &e k 1 e i2 2 e j1 2 e j2 1 3 0 3^ i 1i 3u V 0u k j 3&
e i1 1 e i3 2 e k 2 e j3

3 ^ ki 2 u V 0 u j 1 j 2 & ,~2.24!

D e 5 e i1 1 e i2 1 e i3 2 e j1 2 e j2 2 e j3 ,where V 0 is the perturbation operator, the difference between
the total Hamiltonian and the Fock operator. We did notinclude this term in the present calculation, because our aim
here is to express the 3-RDM in terms of the 1- and2-RDM’s alone, but this term additionally includes the orbital energies and the perturbation operator. We expect that
this term is small because it is proportional to the three-bodycluster amplitude in the lowest second order. In addition to
the lowest second-order term given by Eq. ~2.24!, we can
include some higher-order terms if we substitute6 ^ ki 2 u V 0 u j 1 j 2 & /( e k 1 e i2 2 e j1 2 e j2 ) in Eq. ~2.24! with U.
In Sec. III we will show the contribution of this term in
actual calculations of atoms and molecules.
In the present formulation, we also neglect the third- andhigher-order terms which cannot be broken into two parts by
cutting an electron propagator, as expressed by~2.25!
The physical meaning of the present approximation of the
3-RDM may be explained as follows. The three-particle collision term in the last term of Eq. ~2.10! is considered toinclude two distinct effects: three particles interact at thesame time ~three-body cluster effect!, or particle ~hole! 2
interacts with particle 1, moves some distance, then interacts
with particle 3. The UV term represents the latter process.
Because of the short-range nature of the electron correlations

56@14#, the three-body cluster effect is relatively small in atoms
and molecules, and is neglected in the present approximation.
To test the accuracy of the approximations for the 4- and3-RDM’s given in this section, we calculated the 4-RDM
from the exact 1-, 2-, and 3-RDM’s, and the 3-RDM fromthe exact 1- and 2-RDM’s for the singlet ground state of Be,
and the results were presented in Ref. @11#. Our approximations were very accurate, and reduced the errors of the IPH
approximations by about two orders of magnitude for both 3and 4-RDM’s.
III. DIRECT DETERMINATION
OF THE DENSITY MATRIXIn Sec. II, we derived a decoupling approximation for the3- and 4-RDM’s whose accuracy is up to about second orderof the perturbation. In contrast to the approximation adopted
by Valdemoro and co-workers, our approximation involves
almost all the second-order perturbation terms, except for theterm given by Eq. ~2.24!. This is necessary for a quantitativedescription of the electronic structure of atoms and molecules. Although our approximation based on the Green’sfunction technique is conceptually simple and systematically
improvable, the validity of the approximation itself must be
checked through numerical calculations. In a previous paper@11# we solved the second-order density equation iteratively
using the decoupling approximation of the 3- and 4-RDM’s
described above.
In this section, we first summarize our calculational
method and then give the calculated results. We used thesecond-order density equation @Eq. ~1.1! with n52# in matrix form, which can be represented asi1,i2
R i1,i2
j1,j2 5ED j1,j2 ,~3.1!

where the density matrix and the energy density matrix are
given by

E i1,i2
j1,j2 51
i1,i2
D i1,i2
j1,j2 5 2 ^ E j1,j2 & ,

~3.2!

1
i1,i2
R i1,i2
j1,j2 5 2 ^ HE j1,j2 & ,~3.3!

( a †i1 s 1 a †i2 s 2 a j2 s 2 a j1 s 1 .

~3.4!

s 1,s 2The original integrodifferential equation was replaced by thematrix equation including the partial trace. The density matrices and the Hamiltonian were represented in matrices
whose one-electron base are the HF orbitals. The generalized
two-electron integrals were used for simplicity.
Equation ~3.1! is equivalent to the two conditions
1 i1,i2
j1,j2
1R i1,i2! 2ED i1,i2~R
j1,j2 50,
2 j1,j2

~3.5!

j1,j2
R i1,i2
j1,j2 2R i1,i2 50~3.6!

because of the Hermiticity of the 2-RDM,
j1,j2
D i1,i2
j1,j2 5D i1,i2 .

~3.7!

56

DIRECT DETERMINATION OF THE QUANTUM- . . .

In the present calculation, we used the Hermite part of thedensity equation, Eq. ~3.5!, imposing the Hermiticity, Eq.
~3.7!, and did not use the two-particle Brillouin condition@Eq. ~3.6!#.
We approximate the 3-RDM from 1- and 2-RDM’s, andthen 4-RDM from 1-, 2-, and 3-RDM’s, using the decoupling
approximations given in Sec. II. We then put these 3- and4-RDM’s to Eq. ~3.5!, and solve it iteratively. Since thepresent decoupling approximations of the 3- and 4-RDM’s
are correct essentially up to the second order in the perturbation of the electron correlation, we refer to the present
method as DEPT2 ~density equation with the 3- and4-RDM’s correct to second-order in correlation perturbation!, or simply DE2. In this notation, the solution based on
the IPH approximation ~original Valdemoro’s approximation! may be referred to as DEPT1 or DE1, since the IPH
approximation of the 3- and 4-RDM’s corresponds to thefirst-order approximation.
Since the 1-RDM was represented by the 2-RDM, and the3- and 4-RDM’s by the sums of the products of the 1- and2-RDM’s in the present approximation, the density equationin matrix form becomes a multidimensional nonlinear equation. This equation was solved by Newton’s method, which
is often used for this kind of problem. The iterative procedure is as follows. We use the HF 2-RDM as an initial guess,normalize the trial 2-RDM, then calculate the energy and the3- and 4-RDM’s from the 1- and 2-RDM’s, and substitutethem into Eq. ~3.5!. Then we calculate a correction of the2-RDM ~and hence a new 2-RDM! by Newton’s method.
This procedure is repeated until the convergence is obtained.
In some cases, interpolation techniques are useful for an efficient convergence: the previous and present 2-RDM’s areaveraged with a given weight. We used the symmetry propi2,i1
erty (D i1,i2
j1,j2 5D j2,j1 ) of the 2-RDM in addition to the Hermiticity to reduce the number of the variables.
We calculate the coefficient matrix of the linear equationin the Newton procedure as follows. From its definition thiscoefficient matrix is determined using the linearlized density
equation. The density equation with the decoupling approximations of the 3- and 4-RDM’s contains linear and higherorder terms of the 2-RDM. The coefficient matrix elements
derived from the linear term are calculated in every iteration.
The rest are calculated using the HF RDM’s at the first iteration. For a difficult convergence case, they are calculated
at each iteration using the IPH approximations for both 3and 4-RDM’s.
In the previous communication @11# we have applied theDE2 method to several atoms and singly bonded molecules;
Be, Ne, H2O, NH3, H3O1, CH4, BH42 , NH41 , and CH3F,
and found that our method shows a good convergence, giving energies as accurate as and density matrices more accurate than the SDCI one. In this paper we apply the DE2
method to more complex systems ~better basis sets and larger
molecules!. New molecules include singly bonded ones
~CH3OH,CH3NH2 ,C2H6!, and multiply bonded ones
~N2 ,CO,C2H2!. Next we give a comparison of the accuraciesof the various decoupling approximations from DE1 to DE2:
the results of the density-equation method are compared withthe full-CI results.
Comparative wave-function calculations were done at the
HF, SDCI @16#, full-CI @17# ~in case of relatively small basis2653

sets!, and fourth-order Mo” ller-Plesset perturbation method~MP4, for larger cases!. Experimental molecular geometries@18# were used and some 1s core orbitals of C, N, O, and F
were frozen.
Tables I and II show the summary of the present results.
In Table I we compare the energies and the properties calculated by the present DE2 method with those of the wavefunction approaches. The results for the double-z basis @19#
and the minimal Slater-type-orbital STO-6G basis @20# weregrouped separately. The DE2 method includes the UV and2P terms in the 3- and 4-RDM’s, respectively, but does notinclude the term given by Eq. ~2.24!.
It is known that the Hartree-Fock theory can be formulated as an equation for the first-order density matrix, and theuncorrelated, Hartree-Fock density matrix can be calculated
without using the wave function @3,4,12,21#. However, few
attempts have been successful for a direct determination of
the correlated density matrices. The variational method imposing approximate N-representability conditions on thesecond-order density matrix does not work well, except forsome special systems with a very small number of electrons@1,22#. With the use of the nonvariational density-equation
method and accurate decoupling approximations, we havebeen able to determine the correlated energy and the correlated density matrices of general atoms and molecules quantitatively without any use of the wave function. The use of
the density equation and the density matrix instead of the the
Schrödinger equation and the wave function is the central
idea in the present density-equation approach.
As seen in Table I, the energies of the DE2 method are as
accurate as and even more accurate than the SDCI results.
This is true even for acetylene with double-z basis distributing ten electrons within 20 molecular orbitals in the active
space. This is the largest calculation ever made using thedensity equation. The full-CI calculations were difficult forthis class of systems, so that the comparison was made with
the MP4 results. Because of the limitation of our current
program, we did not apply this method to larger molecules
and to larger active spaces, but this limitation should beovercome. The calculated dipole moment of the CO molecule by the DE2 method was 0.0344 a.u., which agrees wellwith the SDTQ ~single, double, triple, and quadruple!-CI
value, 0.0445 a.u., and the experiment, 0.044 a.u. @23#. The
calculated rotational barrier of ethane was 3.31 kcal/mol,which also agrees well with the MP4 value, 3.34 kcal/mol,and the experiment, 2.9 kcal/mol @24#.
Even though the energy errors shown in Table I are sometimes larger than those of the SDCI method, the dipole orquadrupole moments calculated from the present density matrices, are more accurate than those of the SDCI method.
This result is understood from the facts that the densityequation method is a nonvariational method while the SDCI
method is a variational one, and that the density-equation
method is a method directly solving the densities themselves.
In Table II we compare the calculated energies and density matrices obtained by the various decoupling approximations of the 3- and 4-RDM’s from the DE1 to the DE2 level.
The effects of the term given by Eq. ~2.24! are further investigated. It gives more detailed information than that given inKOJI YASUDA AND HIROSHI NAKATSUJI

2654

56

TABLE I. Energies and properties calculated by the density-equation method and the wave-function method.
Density equation
Molecule

DE2a

Wave function
Hartree-Fock

SDCI

Full-CIb

Energy in a.u. ~correlation energy error in percentage!
dipole or quadrupole moment in a.u.d

Active Number of
space, electronsc
H2O
539, 10NH34310, 8CH44312, 8
HF
4312, 8N25311, 10CO
5311, 10C2H25315, 10276.148 27(6.48)
0.9801
256.298 88(4.24)
0.9098240.295 82(3.73)
02100.218 76(2.85)
0.79202109.079 09(12.27) e
2.046
2112.872 93(12.99) e
0.0344276.980 06(6.37) e
4.642

Double-z basis ~Ref. @19#!
276.009 84(100)1.0142
256.175 99(100)0.9226240.185 46(100)02100.021 79(100)0.82152108.878 26(100)1.808
2112.685 05(100)20.1652276.799 07(100)5.315

276.150 01(5.30)
0.9826256.297 17(5.58)
0.9133240.294 05(5.27)
02100.216 22(4.10)
0.79312109.082 19(10.92) e
1.993
2112.873 82(12.58) e
0.0586276.975 59(8.65) e
4.817

276.157 87(0)0.9762
256.304 33(0)0.9081240.300 09(0)02100.224 53(0)0.79192109.107 19e~0!
2.054f2112.900 98e~0!
0.0445f
276.992 30e~0!
4.762gCH3OH735, 14
CH3NH2736, 14
C2H6 ~staggered!
737, 14
C2H6 ~eclipsed!
737, 142114.711 44(5.24)
0.6024295.067 81(3.76)
0.5634279.209 24(0.24) e
0.4279279.204 05(0.20) e
0.4843

STO-6G basis ~Ref. @20#!
2114.589 82(100)0.6548294.938 49(100)0.5754279.062 33(100)0.4589279.056 86(100)0.51592114.710 82(5.72)
0.6058295.065 05(5.81)
0.5633279.204 48(3.47) e
0.4307279.199 25(3.46) e
0.48762114.718 16(0)0.5991295.072 86(0)0.5619279.209 59e~0!
0.4284f
279.204 35e~0!
0.4851fa

The UV and 2 P terms in the 3- and 4-RDM’s, respectively, are included.
Full-CI only for small active space calculations.
cNumber of electrons in the active space.
d
Quadrupole moment is given in case of zero dipole moment.
eFourth-order Mo” ller-Plesset values.
f
SDTQ-CI values.
g
SDT-CI value.
bRef. @11#. We used the double-z STO basis @25# expanded by
six Gaussian-type orbitals for Be, and the minimal STO-6G
basis @20# for molecules.
We see that unless both of the UV and 2P terms areincluded in the 3- and 4-RDM’s ~DE2!, the results are notsatisfactory. When both terms are included, the convergenceof the iterative calculation is remarkably improved, giving
satisfactory results for both energy and density matrices incomparison with the full-CI results. Original Valdemoro’s
approximation ~IPH, IPH!: DE1 and other combinations~1UV, IPH! and ~IPH, 12P! did not give good results, anddid not even converge in some cases. The combinationsof ~IPH, IPH! and ~IPH, 12P! give only small fractions of
the correlation energies for molecules, as well as too large acorrelation energy, twice as large as the exact value for Be,
and the density errors are as large as those of the HartreeFock method. The combination ~1UV, IPH!, which contains
the first, second, and the third terms of Eq. ~2.5! for the
4-RDM, gives better results than the above ones. It actually

gives the best energies among the density-equation method
for NH41 and CO, but it does not give the best densities: thedensity errors in the ~1UV, 12P!: DE2 approximation aresmallest in both density-equation and wave-function methods.
Strictly speaking, the DE2 method must be called, say,DE1.9, because it does not include all the second-orderterms, and the term given by Eq. ~2.24! is missing. We thenconsider the effect of this term. When this term is further
included, and hence the method is strictly DE2, the energy
becomes slightly better than that of the DE1.9 method~which is written in Table II as DE2!. For the CO moleculethe missing term improves the energy considerably, and reduces the correlation energy error up to half of that of DE1.9.
However this may be a special case because the effect of thesame term is small for the triply bonded molecules N2 andC2H2. In all cases, the density matrices calculated by including the missing term are similar or a bit worse than those
calculated by the DE1.9 method.

DIRECT DETERMINATION OF THE QUANTUM- . . .

56

2655

TABLE II. Comparison of the density-equation method using the various decoupling approximations and the wave-function method.
Numbers in square brackets indicate powers of 10.
Density-equation method
3

D
4
D

DE1
IPH
IPH

1UV
IPH

IPH
12P

SystemActive spaceelectronsa
Be
232
4
H2O
532
10
CH4
434
8
NH41
434
8
N2
533
10
CO
533
10
C2H2
535
10
a

Wave function method

DE2
1UV
12P

Hartree-Fock

SDCI

Full-CI

1UV1Eq. (2.24)
12P

EnergyCorrelation energy error in percentage
2
D errorb
214.597 75
2106.4
1.165@ 22 #not conv.

214.581 85
5.93
8.768@ 24 #
275.722 69
12.6
9.132@ 22 #
240.185 37
6.38
3.818@ 22 #256.485 93
22.96
4.956@ 22 #not conv.

214.598 50
2111.6
1.195@ 22 #not conv.

214.582 69
0.014
4.823@ 25 #275.727 51
3.00
7.122@ 23 #
240.124 26
240.111 36 240.187 26
82.4
98.5
4.02
2.948@ 21 #3.713@ 21 # 9.757@ 23 #256.423 54
256.404 41 256.479 90
72.2
95.3
4.31
2.912@ 21 #3.394@ 21 # 1.164@ 22 #not conv.
not conv. 2108.673 19
17.00
7.601@ 22 #not conv. 2112.445 29 not conv. 2112.419 26
1.93
16.77
3.959@ 21 #1.202@ 21 #not conv.
276.752 13
not conv.
276.757 17
14.30710.98
1.852@ 21 #6.493@ 22 #214.582 69
24.60@ 23 #4.784@ 25 #275.727 70
2.63
1.254@ 22 #
240.187 53
3.68
1.007@ 22 #256.480 46
3.63
1.338@ 22 #2108.676 16
15.13
7.106@ 22 #2112.431 25
8.15
3.398@ 21 #276.758 28
10.34
6.353@ 22 #214.568 53 214.582 69 214.582 69
100
0.028
0
4.631@ 22 #1.141@ 24 #0275.678 84 275.728 29 275.729 02
100
1.44
0
3.154@ 21 #1.030@ 22 #0240.110 15 240.187 72 240.190 49
100
3.45
0
3.702@ 21 #2.819@ 22 #0
256.400 47 256.480 87 256.483 47
100
3.14
0
3.671@ 21 #2.757@ 22 #02108.541 77 2108.687 46 2108.700 11
100
7.99
0
7.159@ 21 #1.282@ 21 #02112.303 32 2112.430 00 2112.442 61
100
9.05
0
7.454@ 21 #1.618@ 21 #0276.603 02 276.760 57 276.776 19
100
9.02
0
8.864@ 21 #1.711@ 21 #0

Number of electrons in the active space.
Square norm of the difference between the calculated 2-RDM and the full-CI one.

b

It is remarkable that the correlation energy errors of thepresent density-equation method are improved as the system
becomes larger ~in size and in basis set!. This is seen bycomparing the results given in Table I with those of Table II,
the former being for larger systems. Though the energy errors of the density-equation method given in Table II arealmost always worse than those of the SDCI method, thosein Table I are opposite except for H2O, N2, and CO. This iscertainly an encouraging result for the density-equation
method.
Recently, Colmenero and Valdemoro applied the density
equation, using their approximations for 3- and 4-RDM’s, tofour-electron systems like Be, B1, etc. and a six-electron
system Betta @10#. They further used corrections adopting the
3-RDM spin-symmetry, and keeping some diagonal elementsof the 3-RDM non-negative. The error in the correlation energy for Be was 0.27%, which is comparable to our present
result, showing that their additional corrections worked well.
In our calculations, some diagonal elements of the approximate 3-RDM were slightly negative, though they must
be non-negative for the N-representable one. Hence we further tried to keep all the diagonal elements of the 3-RDM
and 3-HRDM non-negative, and at the same time to keepthem spin symmetry adapted. The spin symmetry of theRDM corresponds to the permutational antisymmetry of the
explicit ~unintegrated! variables of the DM. Our correctionused here is similar to the one used by Valdemoro and coworkers @8#: if a negative diagonal element appeared in the
3-RDM, it was corrected to be zero, and some nondiagonal
elements were modified to keep the spin symmetry. Similarcorrections were applied to the 3-HRDM because the3-RDM and the 3-HRDM must satisfy an exact relation @8#.
However, these corrections either had no effect or resulted ina divergence at least from our simple convergence technique.
Hence we did not use them in the calculations given in
Tables I and II. Although these corrections keep the spinsymmetry of the 3-RDM and hence the antisymmetry of the
explicit variables of the 3-DM, it may not ensure the antisymmetry of the integrated variables of the 3-DM.
We next show the convergence property of our method.
Figure 1 shows the energy, the errors in the density equation,and the 2-RDM for each iteration of the calculation of
CH 3 F. The error in the density equation means the square
norm of the difference between the rhs and lhs of the density
equation, and the error in the 2-RDM means the square normof the difference between the calculated 2-RDM and the exact one. We see that our iterative method gave a nice convergence which requires about ten iterations, in contrast toKOJI YASUDA AND HIROSHI NAKATSUJI

2656

FIG. 1. Energy and errors of the density equation and 2-RDM
for each iteration of the calculation of CH3F.

the previous one @10# which required longer iterations.
We examine the N representability of the calculated
2-RDM. We calculated some necessary requirements, the P,
Q, and G conditions @1,2#, which are the non-negativities of
the RDM, HRDM, and the g matrix defined as
g i j,kl 523 2 D iljk 1 d ki 1 D lj 2 1 D ij 1 D kl .

~3.8!

The G condition is equivalent to the non-negativity of thenorm,^ B † B & 5x *i j x kl g i j,kl >0,~3.9!

in which x i j are arbitrary scalars and B is the operator defined as
B5( x i j ~ a i†s a j s 2 1 D ij ! .

i j, s~3.10!

Table III shows the lowest eigenvalues of the three matri-

56

FIG. 2. Minimum eigenvalues of the 2-RDM and 2-HRDM, andg matrix and trace errors of the n-RDM’s for each iteration of thecalculation of CH3F.

ces P, Q, and g for several molecules. Though they must be
non-negative for the N-representable 2-RDM’s, the lowest
values are slightly negative but the absolute values are small.
Note that the sums of the eigenvalues of the P and Q matrices are quite large, which are N C 2 and M C 2 , respectively,
where N being the number of the electrons and M the number of the holes. Note also that the sum of the eigenvalues of
the g matrix in the Hartree-Fock approximation isN(N12M )/2, which is also quite large. Figure 2 showsthese lowest eigenvalues for each iteration of the calculationof CH3F. The deviations from zero are large in the first few
iterations, but become very small after convergence, that is,the calculated 2-RDM’s almost satisfy these three conditions. A small deviation from the exact N representability is
due to an inaccuracy of our decoupling approximation of the
RDM’s. As the approximation becomes more accurate, thecalculated 2-RDM’s should satisfy the N-representability
condition more precisely. Figure 2 also shows the trace er-

TABLE III. Minimum eigenvalues of the 2-RDM, 2-HRDM, and g matrix, and the range of the eigenvalues of 1-RDM calculated by theDE2 method. Numbers or square brackets indicate powers of 10.
Molecule

Activespace

Number ofelectronsa

2-RDM

2-HRDMg matrix

1-RDM

H2O
NH3
CH4
CH3OH
N2
CO
C2H2

539
4310
4312
735
5311
5311
5315

10
8
8
14
10
10
10

23.34@ 24 #24.59@ 24 #25.71@ 24 #25.39@ 24 #21.02@ 23 #27.83@ 24 #21.52@ 23 #21.85@ 24 #21.46@ 24 #21.59@ 24 #24.35@ 24 #25.31@ 24 #23.37@ 24 #25.94@ 24 #23.28@ 24 #22.70@ 24 #23.22@ 24 #25.39@ 24 #23.70@ 23 #26.40@ 23 #23.10@ 23 #1.94@ 24 # ;1.9998
3.00@ 24 # ;1.9838
2.35@ 24 # ;1.9815
1.78@ 22 # ;1.9974
4.30@ 24 # ;1.9866
5.51@ 24 # ;1.9884
3.40@ 25 # ;1.9796

aNumber of electrons in the active space.

56

DIRECT DETERMINATION OF THE QUANTUM- . . .

rors of the n-RDM from the exact trace, N C n , which also
became very small after the convergence.
We finally check the N representability of the 1-RDM’s.
The necessary and sufficient N-representability condition for
the 1-RDM is that all the eigenvalues are in the range of zero
to two @2#. Table III shows the range of the eigenvalues of
the 1-RDM’s for several molecules. All the eigenvalues of
the calculated 1-RDM’s are in the range of zero to two,showing that the calculated 1-RDM’s are exactly N representable. This was true for all atoms and molecules calculated by the present DE2 method.
IV. CONCLUSION

The second-order density matrices of some atoms andmolecules were directly calculated without any use of the
wave function. We solved the density equation using thedecoupling approximations for the third- and fourth-orderdensity matrices. The decoupling approximations were derived based on the Green’s-function technique and using the
Feynman diagram. Previous approximations by Valdemoro
and co-workers @7,8# were derived as a special case. The
present approximation involves up to the second-order and@1# C. Garrod, M. V. Mihailovic, and M. Rosina, J. Math. Phys.
~N.Y.! 16, 868 ~1975!.
@2# A. J. Coleman, Rev. Mod. Phys. 35, 668 ~1963!; C. Garrod and
J. Percus, J. Math. Phys. ~N.Y.! 5, 1756 ~1964!; H. Kummer, J.
Math. Phys. ~N.Y.! 8, 2063 ~1967!; T. L. Gilbert, Phys. Rev. B
12, 2111 ~1975!.
@3# H. Nakatsuji, Phys. Rev. A 14, 41 ~1976!.
@4# S. Cho, Sci. Rep. Gumma Univ. 11, 1 ~1962!; L. Cohen and C.
Frishberg, Phys. Rev. A 13, 927 ~1976!.
@5# E. Davidson ~private communication! to one of the authors~H.N.! at Tomakomai in 1976; J. E. Harriman, Phys. Rev. A
19, 1893 ~1979!.
@6# P. Hohenberg and W. Kohn, Phys. Rev. B 864, 136 ~1965!; W.
Kohn and L. J. Sham, Phys. Rev. A 1133, 140 ~1965!.
@7# C. Valdemoro, Phys. Rev. A 45, 4462 ~1992!.
@8# F. Colmenero, C. Perez del Valle, and C. Valdemoro, Phys.
Rev. A 47, 971 ~1993!; F. Colmenero, and C. Valdemoro, ibid.
47, 979 ~1993!.
@9# Valdemoro called the density equation a ‘‘contracted Schrödinger equation.’’ However, this name is not appropriate since
it seems to represent only the necessary part of the density
equation. We think the necessity alone is trivial, and what is
important is the necessity and sufficiency relation of the density equation with the Schrödinger equation.
@10# F. Colmenero and C. Valdemoro, Int. J. Quantum Chem. 51,
369 ~1994!; for the recent development of the approach, see C.
Valdemoro, L. M. Tel and E. Pérez-Romero, Adv. Quantum
Chem. 28, 33 ~1997!.
@11# H. Nakatsuji and K. Yasuda, Phys. Rev. Lett. 76, 1039 ~1996!.

2657

some higher-order perturbation terms in electron correlations, and is more accurate than the previous one @7,8#,
which neglects some important second-order terms. The calculated 2-RDM’s were more accurate than the SDCI results,
and the calculated energies were comparable with or better
than the SDCI ones: the SDCI method is a variational
method which minimize the energy, while the densityequation method is a nonvariational method directly determining the 2-RDM. The density-equation method is found to
give better quality results as the system becomes larger. The
convergence was fairly good and the calculated 2-RDM’s
almost satisfied some necessary conditions of the
N-representability, P, Q, and G conditions, while the1-RDM’s were exactly N representable. These results showthat the density-equation method is very promising.

ACKNOWLEDGMENTS

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 financial support
from the Japan Society for the Promotion of Science.

@12# R. G. Parr and W. Yang, Density-Functional Theory of Atoms
and Molecules ~Oxford University Press, New York, 1989!.
@13# A. A. Abrikosov, L. P. Gor’kov, and E. Dzyaloshinskii, Methods of Quantum Field Theory in Statistical Physics ~PrenticeHall, Englewood, Cliffs, NJ, 1963!.
@14# O. Sinanoglu, J. Chem. Phys. 36, 706 ~1962!.
@15# R. P. Feynman and A. R. Hibbs, Quantum Mechanics andPath Integrals ~McGraw-Hill, New York, 1965!
@16# M. Dupuis, A. Farazdel, S. P. Karna, and S. A. Maluendes,
HONDO8.1, IBM Corporation.
@17# P. J. Knowles and N. C. Handy, Chem. Phys. Lett. 111, 315~1984!.
@18# L. E. Sutton, D. G. Jenkin, and A. D. Mitchell, Tables ofInteratomic Distances ~Chemical Society, London, 1958!.
@19# T. H. Dunning, J. Chem. Phys. 53, 2823 ~1970!.
@20# W. J. Hehre, R. F. Stewart, and J. A. Pople, J. Chem. Phys. 51,
2657 ~1969!.
@21# R. McWeeny, Rev. Mod. Phys. 32, 335 ~1960!.
@22# M. Rosina, B. Golli, and R. M. Erdahl, in Density Matrices
and Density Functionals, edited E. Erdahl and V. H. Smith, Jr.
~Reidel, Dordrecht, 1987!, p. 231.
@23# A. Szabo and N. S. Ostlund, Modern Quantum Chemistry:
Introduction to Advanced Electronic Structure Theory ~Macmillan, New York, 1982!.
@24# S. H. Pine, Organic Chemistry ~McGraw-Hill, New York
1987!.
@25# E. Clementi and C. Roetti, At. Data Nucl. Data Tables 14, 428~1974!. The exponents were reoptimized by T. Koga and S.
Watanabe ~private communication!.