JOURNAL OF CHEMICAL PHYSICSVOLUME 114, NUMBER 19

15 MAY 2001

Variational calculations of fermion second-order reduced density matrices
by semidefinite programming algorithm
Maho Nakata, Hiroshi Nakatsuji,a) and Masahiro Ehara
Department of Synthetic Chemistry and Biological Chemistry, Graduate School of Engineering,
Kyoto University, Kyoto 606-8501, Japan

Mitsuhiro Fukuda
Department of Mathematical and Computing Sciences, Tokyo Institute of Technology,
Tokyo 152-8852, Japan

Kazuhide Nakata
Department of Applied Physics, The University of Tokyo, Tokyo 113-8656, Japan

Katsuki Fujisawa
Department of Architecture and Architectural Systems, Kyoto University, Kyoto 606-8501, Japan

共Received 24 July 2000; accepted 9 February 2001兲
The ground-state fermion second-order reduced density matrix 共2-RDM兲 is determined variationallyusing itself as a basic variable. As necessary conditions of the N-representability, we used thepositive semidefiniteness conditions, P, Q, and G conditions that are described in terms of the2-RDM. The variational calculations are performed by using recently developed semidefinite
programming algorithm 共SDPA兲. The calculated energies of various closed- and open-shell atoms
and molecules are excellent, overshooting only slightly the full-CI energies. There was no casewhere convergence was not achieved. The calculated properties also reproduce well the full-CI
results. © 2001 American Institute of Physics. 关DOI: 10.1063/1.1360199兴
I. INTRODUCTION

straightforward consequence of the Ritz variational principle
combined with the fact that the Hamiltonian involves only
one- and two-body operators. The problem here is how well
we can restrict our variable ⌫ (2) to be N-representable.10 The
N-representability condition that is enforced by the Pauliprinciple is not completely known for ⌫ (2) and this is an
obstacle of the DMT in general.
The P, Q, 共Ref. 10兲 and G 共Ref. 11兲 conditions are thewell-known necessary conditions of the N-representability.
They are the semidefiniteness conditions of the matrices derived from ⌫ (2) . Though these three conditions are not complete, they seem to be quite strong to characterize theN-representability of the ground-state 2-RDM. First calculations along this line were performed in a beautiful way by
Garrod et al.12,13 for the ground state of Be, and Mihailović
et al.14 for the nuclear ground state of 15O, 16O, 17O, 18O,
20
Ne, 24Mg, and 28Si. At that time their method was very
heuristic and could hardly be applied to general systems. We
found that this method can be elegantly realized using thesemidefinite programming algorithm 共SDPA兲,15 recently developed in the field of mathematical programming. We calculated the ground-state energies of atoms and molecules
using these three necessary conditions and employing SDPA
as our problem solver.

The ground state of N-body fermion system is completely described by the second-order reduced density matrix
共2-RDM兲 ⌫ (2) because any observable properties of the system can be calculated from the 2-RDM.1,2 This fact led us todesire to use 2-RDM as a basic variable of quantum mechanics instead of the wave function ⌿; if we can determine ⌫ (2)without using ⌿, we have a closed form of quantum mechanics where the basic variable is 2-RDM. We refer to suchformalism of quantum mechanics as density matrix theory
共DMT兲.3
In the nonrelativistic case, the determinative equation for⌿ is the Schrödinger equation 共SE兲. Therefore, to establish
DMT, we have to formulate the equation for the RDM that isequivalent to the SE in the necessary and sufficient sense.3
As such an equation, one of the author derived density equation 共DE兲 共Refs. 4 and 5兲 that has recently been used successfully to calculate the 2-RDMs of atoms and molecules
directly without any use of the wave function.6–9 This approach is called density equation theory 共DET兲 and a review
on the DET in chemical physics has recently been summarized together with some later developments.3
Another equation that is equivalent to the SE but includes 2-RDM alone as a variable is the variational equationof the form,
E g ⭐E 关 ⌫ (2) 兴 ,II. THEORICAL OUTLINE

共1.1兲

First and second order reduced density matrices 共1-,
2-RDMs兲, ␥ and ⌫, respectively, are defined bywhere E g is the exact ground-state energy. This method
called density matrix variational theory 共DMVT兲 is a

␥ ij ⫽ 具 ⌿ 兩 a †i a j 兩 ⌿ 典 ,

a兲

and

Author to whom correspondence should be addressed.

0021-9606/2001/114(19)/8282/11/$18.00

共2.1兲

8282

© 2001 American Institute of Physics

Downloaded 06 Jun 2002 to 130.54.33.130. Redistribution subject to AIP license or copyright, see http://ojps.aip.org/jcpo/jcpcr.jsp

J. Chem. Phys., Vol. 114, No. 19, 15 May 2001

Calculations of second order reduced density matrices

i i⌫ j1 j2 ⫽ 12 具 ⌿ 兩 a i† a i† a j 2 a j 1 兩 ⌿ 典 ,1

1 2

where a † and a denote creation and annihilation operators,
respectively. Note we have simplified ⌫ (2) as ⌫. Throughout
this paper, we assume the elements of 1-RDM and 2-RDM tobe real. Complete N-representability condition is known for
␥,10 but for ⌫, we know only necessary conditions 共theknown complete condition is not practical兲. Some trivial conditions for 2-RDM are:
共1兲 Antisymmetric condition,
i i

i i

i i

1 2

1 2

2 1

⌫ j1 j2 ⫽⫺⌫ j2 j1 ⫽⫺⌫ j1 j2 ;

共2.3兲

j j

1 2

1 2

共2.4兲

2

兺k ⌫ ikjk ⫽ N⫺1 ␥ ij ;共2.5兲共4兲 Number of electrons,兺k ␥ kk ;共2.6兲

Tr N ␣2 ⌫⫽N ␣2 ;

and

共2.7兲

where the operators of N ␣ and N ␣2 are written as

兺i a i†␣ a i ␣ ,N ␣2 ⫽

兺i j a i†␣ a i ␣ a †j ␣ a j ␣ ;i1␴1i2␴2⌫ j ␴ j ␴ ⫽0;

共2.9兲

1 1⬘ 2 2⬘

when ␴ 1 ⫽ ␴ ⬘1 or ␴ 2 ⫽ ␴ 2⬘ and ␴ 1 ⫽ ␴ 2⬘ or ␴ 2 ⫽ ␴ 1⬘ , where ␴denotes spin variable.
共7兲 Expectation value of S 2 ,Tr S 2 ⌫⫽S 共 S⫹1 兲 ,共2.10兲

where the spin-squared operator S 2 is given byS 2 ⫽S z ⫹S z2 ⫹S ⫺ S ⫹
1
2

兺i

⫹

共 a i†␣ a i ␣ ⫺a i†␤ a i ␤ 兲 ⫹兺i j a i†␤ a i ␣ a †j ␣ a j ␤ ;

1
4

冉兺
i

a i†␣ a i ␣ ⫺a i†␤ a i ␤

冊

2

共2.11兲

共8兲 Positive semidefiniteness of P matrix, which is just
2-RDM,

兺 x i i ⌫ ij ij x j j ⭓0,
1 2

1 2

1 2

1 2

共2.14兲

2

These matrices are semidefinite11 and linear to ⌫ as
i i

i

i

i

i

i

i

i

i

Q j1 j2 ⫽ 共 ␦ j1 ␦ j2 ⫺ ␦ j1 ␦ j2 兲 ⫺ 共 ␦ j1 ␥ j2 ⫹ ␦ j2 ␥ j1 兲
1 2

1

2

2

1

1

2

2

i
i
i
i
i i
⫹ 共 ␦ j1 ␥ j2 ⫹ ␦ j2 ␥ j1 兲 ⫺2⌫ j1 j2
2
1
1
2
1 2

共2.12兲

where x i 1 i 2 is an arbitrary geminal.
Note that except for the condition 8, all of these conditions are linear to 2-RDM.
The Q and G matrices are defined by1

共2.15兲

and
i

i

2

i j

1

1 2

共2.16兲

We note that originally the G matrix was written in an
equivalent nonlinear form.16
The Hamiltonian of the system can be written as
i i

i i

1 2

1 2

H j1 j2 ⫽w j1 j2 ⫹

1
i i
i i
共 v 1 ␦ 2 ⫹ v j1 ␦ j2 兲 ,2
1N⫺1 j 1 j 2

共2.17兲

where v and w are 1- and 2-body operators, respectively.
Then, the basic equation of DMVT given by Eq. 共1.1兲 iswritten as the variational minimization of the energy within
our constraints,
共2.18兲

⌫苸 (2) P

where (2) P is the set of 2-RDMs which satisfies the abovenecessary N-representability conditions, namely,
(2)

共2.8兲

共6兲 Spin symmetry,

⫽

1

1 2

E min⫽ min Tr H⌫,

共5兲 Eigenstate of the number of ␣ 共or ␤兲 electrons,

N ␣⫽

i i

G j1 j2 ⫽ 具 ⌿ 兩 a i† a i 2 a †j a j 1 兩 ⌿ 典 .

1 2

共3兲 Trace condition,

Tr N ␣ ⌫⫽N ␣

共2.13兲

1

and

i i

⌫ j1 j2 ⫽⌫ i 1i 2 ;

N⫽

2

1 2

G j1 j2 ⫽ ␦ j2 ␥ j1 ⫺2⌫ j1 i 2 .

共2兲 Hermiticity,
i i

i i

Q j1 j2 ⫽ 具 ⌿ 兩 a i 1 a i 2 a †j a †j 兩 ⌿ 典

共2.2兲

2

8283

P⫽ 兵 ⌫ 兩 P,Q,G matrices are non-negative and theconditions 1 – 7 are satisfied其 .

共2.19兲

Either of the P, Q, and G conditions forms compact convex
set with trace topology,17 and a finite combination of compact convex sets is also compact convex set, therefore thismethod should find a minimum in energy. This method canbe applied to the ground state of any space and spin symmetry.
III. CALCULATION METHOD

The minimization problem with some linear constraints
can be achieved by using semidefinite programming algorithm 共SDPA兲 共Ref. 15兲 as a problem solver. The SDPA has
recently been developed in the field of mathematical programming. In this section, we explain how to apply SDPA to
our problem of solving ⌫ (2) in the constrained variational
method given by Eq. 共2.19兲. The dimensions of the matrices
are n⫻n, if they are not explicitly defined.
A. Simplified problem

First, we introduce a simplified problem which contains
all the essentials, that is
Problem 共a兲: Minimize the total energy of the 2-RDM ⌫subject to the fixed number of electrons and the positive
semidefiniteness of ⌫.
The positive semidefiniteness of ⌫ is the P condition. Notethat this problem gives the exact solution for N⫽2. Problem
共a兲 is written as,

Downloaded 06 Jun 2002 to 130.54.33.130. Redistribution subject to AIP license or copyright, see http://ojps.aip.org/jcpo/jcpcr.jsp

8284

Nakata et al.

J. Chem. Phys., Vol. 114, No. 19, 15 May 2001

Problem 共 a⬘ 兲 :再
再Minimize Tr H⌫
共3.1兲

subject to Tr N⌫⫽N
and ⌫ is positive semidefinite.

Minimize Tr H⌫

Formal expression of the problem15 is,Problem 共 a⬙ 兲 :number of electrons, F2 as the constraint for spin squared
operator, etc. Maximization is altered to minimization by justchanging the sign of F0 . Problem 共2.19兲 is written as

subjected to Tr N⌫⫽N

Minimize F0 "Ysubject to F1 "Y⫽c 1

Tr S 2 ⌫⫽S 共 S⫹1 兲

共3.2兲

where F0 and F1 are constant n⫻n symmetric matrices, andY is n⫻n symmetric variable matrix, c 1 is real constant, and" is an operator such thatF"Y⫽共 F兲 i j 共 Y兲 i j .
兺
i, j

共3.3兲

One can easily confirm that problem (a⬘ ) and problem (a⬙ )are the same when we take Y as ⌫, F0 as the Hamiltonian,and F1 as the number of operator. A generalization of the
problem (a ⬙ ) is called semidefinite programming 共SDP兲.

冦m兺 c ix i兺dual:

Fi x i ⫺F0 ,

YⱰ0,
共3.4兲
where X and Y are n⫻n real symmetric matrices, Fi (1⭐i⭐m) symmetric constraint matrices, c i and x i real constantand variable numbers, respectively, U"V denotes inner product of the matrices, U"V⫽ 兺 i, j Ui, j Vi, j , and XⱰ0 means X to
be positive semidefinite. We assume all the constraint matrices are linearly independent.
Semidefinite programming is usually solved by primaldual interior-point method.18,19 This method is based on the
primal-dual theorem of SDP, which shows an existence ofthe optimal solution and gives a necessary and sufficient condition for the optimal solution 共minimum in primal problem,and maximum in dual problem兲: if there exists (X,Y,x) suchthat they satisfy all the constraints and XⱰ0 and YⱰ0, then
共1兲 SDP has an optimal solution;
共2兲 Necessary and sufficient condition for the optimal
solution (X* ,Y*,x* ) ism

兺 c i x *i ⫽0.

1

1 2

2

2

1

共3.7兲

i i

where a j1 j2 is constant.
1 2

共2兲 Set up the constraint matrix FA such that
i i

i i

1 2

1 2

FA "⌫⫽

subject to Fi "Y⫽c i 共 1⭐i⭐m 兲 ,

i⫽1

1 2

XⱰ0

maximize F0 "Y

X* "Y* ⫽F0 "Y* ⫺i ia j j a i† a i† a j a j ,兺
i i j j共3.8兲共3兲 Then, the constraint is given by the equality,mi⫽1Note that some of the matrices appeared below have four
indices, however, we can reduce them to two indices bymapping indices (i, j) to the composite index k. Imposing
linear constraints for N or S 2 , etc., is straightforward. Constraining the expectation value of the two body operator A to
be c a (Tr A⌫⫽c a ) is done as follows:
共1兲 Explicit expression of A is give by

共 FA 兲 j1 j2 ⫽a j1 j2 .

i⫽1subject to X⫽

and ⌫ (2) Ɒ0, QⱰ0 and GⱰ0.

1 2 1 2The SDPA 共Ref. 15兲 solves the following form of
semidefinite programming and its dual:

SDP

Tr N ␣2 ⌫⫽N ␣2

A⫽

B. Semidefinite programming algorithm „SDPA…

primal: minimize

共3.6兲

Tr N ␣ ⌫⫽N ␣

and Y is positive semidefinite,

共3.5兲

C. Set up of DMVT in SDPA

Our object is to solve the DMVT problem 共2.19兲. It isequivalent to solve the dual of problem 共3.4兲, taking Y as
2-RDM, F0 as the Hamiltonian, F1 as the constraint for the兺 a ij ij Y ij ij ⫽Tr A⌫⫽c a .

i1i2 j 1 j 2

1 2

1 2

1 2

1 2

共3.9兲

For example, we set up the constraint matrix for thenumber of particles N. Explicit expression of N is,N⫽N⫺1

兺i a †i a i ⫽ 2 兺i j a †i a †j a j a i .

共3.10兲

Then, the constraint matrix FN for the number of particle isrepresented by
共 FN 兲 iklj ⫽N⫺1
␦ ik ␦ jl ,2

共3.11兲

and we confirm the following relation:
FN "Y⫽

共 FN 兲 iklj ⌫ iklj
兺
i jkl

N⫺1
␦ ik ␦ jl ⌫ iklj
2

⫽

兺
i jkl

⫽

N⫺1
⌫ij
2 ij ij

⫽

N⫺1
2

兺

2

兺i N⫺1 ␥ ii ⫽ 兺i ␥ ii ⫽N.

共3.12兲

Now we consider how to enforce the 2-RDM to satisfy
the P, Q, and G conditions, simultaneously. We first explain
the case where only P and Q conditions are enforced simultaneously.

Downloaded 06 Jun 2002 to 130.54.33.130. Redistribution subject to AIP license or copyright, see http://ojps.aip.org/jcpo/jcpcr.jsp

J. Chem. Phys., Vol. 114, No. 19, 15 May 2001

Calculations of second order reduced density matrices

We introduce the variable matrix Y in which P and Q
matrices are diagonally arranged,

冉 冊
P

Y⫽

0

0

共3.13兲.

Q

YⱰ0↔PⱰ0and

共3.14兲QⱰ0.

There is a linear relation between ⌫ and Q matrices,
i

i

i

i

Q j1 j2 ⫽ 共 ␦ j1 ␦ j2 ⫺ ␦ j1 ␦ j2 兲 ⫺
1 2

2

1

兺
k k l l1 2 1 2⫽

It is obvious thati i

the second two terms of Eq. 共3.21兲 give,2

1N⫺1 i 1 i 2 k
i k
i
i i
共 ␦ j ⌫ j k ⫹ ␦ j2 ⌫ j1 k 兲 ⫺2⌫ j1 j2 . 共3.15兲
2
1
1
2
1 2
2

兺k

Therefore, we can find a set of linear constraints for eachelement of the Q matrix as
i i

i i

1 2

1 2

1 2

i

i

1

i

i

i

2

i

2

i

i

1

i i

⫺ 共 ␦ j1 ␥ j2 ⫹ ␦ j2 ␥ j1 兲 ⫹Q j1 j2
2

1

1

2

1 2

i i
i i
⫽ ␦ j1 ␦ j2 ⫺ ␦ j1 ␦ j2 .
1
2
2
1

兺
k k l l冉

Using these constraints, the SDP formalism is given by

共3.17兲

1 2

i

1

i

2

i

i i

and the explicit expression of the element of the constraint
i i k kmatrix (E j1 j2 ) l 1l 2 is given by
1 2

i i k kiij ji ⫹n i ⫹n j ⫹n j ⫹n
共 E j1 j2 兲 l 1l 2 ⫽2 ␦ k1 ␦ k2 ␦ l 1 ␦ l 2 ⫹ ␦ k1 ␦ k2 ␦ l 1 ␦ l 2
1 2 1 2
1
2 1
2
1
2
1
2

N⫺1 i 1 i 2 j 2 k 2 N⫺1 i 2 i 1 j 1 k 2
␦ j ␦k ␦l ␦l ⫹
␦ j ␦k ␦l ␦l
1
1 1
2
2
1 1
2
2
2

N⫺1 i 1 i 2 j 1 k 2 N⫺1 i 2 i 1 j 2 k 2⫺
␦ j ␦k ␦l ␦l ⫺
␦ j ␦k ␦l ␦l
2
1 1
2
1
1 1
2
2
2
共3.19兲
i i
and the constant c j1 j2 in Eq. 共3.16兲 is
1 2
i i

i

i

i

i

1 2

1

2

2

共3.20兲

1

i i

1 2

i i k kk ki i
共 E j j 兲 l l 共 Y兲 l l ⫽c j j兺
k k l l1 2 1 2

兺

冊

2

2

1

1

N⫺1

1

1

2

2

i

1

共3.24兲

2

Y⫽

冉 冊
P

0

0

0

Q

0

0

0

G

共3.25兲

.

We have a linear relation between ⌫ and G,i

2N⫺1 i 1 k
i j
⌫ j k ⫺2⌫ j1 i 2 ,
1
1 2
2

兺k

共3.26兲

element of G matrix as
J j1 j2 "Y⫽0⫽⫺ ␦ j2 ␥ j1 ⫹2⌫ j1 i 2 ⫹G j1 i 2 ,
i i

i

1 2

i

2

1

i j

i j

1 2

1 2

共3.27兲
i i

k k

1 2

1 2

and an explicit expression of the constraint matrix (J j1 j2 ) l 1l 2
is given by
i i

k k

1 2

1 2

i

i

j

i ⫹2n i 2 ⫹2n j 1 ⫹2n j 2 ⫹2n
␦k
␦l
␦l
1
2
1
2

j

共 J j1 j2 兲 l 1l 2 ⫽2 ␦ k1 ␦ k2 ␦ l 1 ␦ l 2 ⫹ ␦ k1
1

⫺

2

1

2

N⫺1 i 2 i 1 i 2 k 2
␦ j ␦k ␦l ␦l ,
2
1 1 2
2

共3.28兲

which is further symmetrized as
i i

i i

j j

1 2

1 2

1 2

J̃ j1 j2 ⫽ 12 共 J j1 j2 ⫹Ji 1i 2 兲 .

共3.29兲

Thus, the DMVT using the P, Q, and G conditions is formulated into SDPA as

1 2

1 2

1 2

1 2

1 2

1 2

1 2

1 2

共3.21兲

subjected to Fi "Y⫽c i
Ẽ j1 j2 "Y⫽ ␦ j1 ␦ j2 ⫺ ␦ j1 ␦ j2
i i

i

1

i

2

i

2

共3.30兲

1

J̃ j1 j2 "Y⫽0.
i i

i
i
j j
i ⫹n i ⫹n j ⫹n j ⫹nk k
共 2 ␦ k1 ␦ k2 ␦ l 1 ␦ l 2 ⫹ ␦ k1 ␦ k2 ␦ l 1 ␦ l 2 兲共 Y兲 l 1l 2
1
2 1
2
1
2
1
2
1 2

i ⫹n,i ⫹n
i i
i i
i i
⫽2 共 Y兲 j1 j2 ⫹ 共 Y兲 j1 ⫹n, j2 ⫹n ⫽2⌫ j1 j2 ⫹Q j1 j2 ,
1 2
1
2
1 2
1 2

i

1 2

as follows. The first two terms of Eq. 共3.21兲 are

k1k2l1l2

共3.23兲

Minimize H"Y

We can confirm Eq. 共3.21兲 holds
E j1 j2 "Y⫽1

1 2

共3.18兲

c j1 j2 ⫽ ␦ j1 ␦ j2 ⫺ ␦ j1 ␦ j2 .

12which is described by a set of linear constraints J j1 j2 for each

1 2

⫹

2

i i

where Ẽ j1 j2 is a symmetric matrix defined by

1 2

2

i

1

1 2

1

i i
i i
j j
Ẽ j1 j2 ⫽ 12 共 E j1 j2 ⫹Ei 1i 2 兲 ,
1 2
1 2
1 2

1

兺k ␦ ij ⌫ ij kk ⫺ 2 兺k ␦ ij ⌫ ij kk

i

2

i i

i

2

N⫺1
2

G j1 j2 ⫽ ␦ j2

Ẽ j1 j2 "Y⫽ ␦ j1 ␦ j2 ⫺ ␦ j1 ␦ j2 .
i2⫽⫺ ␦ j1 ␥ j2 ⫺ ␦ j2 ␥ j1 .

Minimize H"Yi i

1N⫺1 j 2 i 2 j 1 k 2 N⫺1 j 1 i 1 j 2 k 2k k
␦i ␦k ␦l ␦l ⫺
␦ i ␦ k ␦ l ␦ l 共 Y兲 l 1l 2
1
1 1
2
2
1 1
2
1 2
2
2

⫺

⫽⫺

共3.16兲

subject to Fi "Y⫽c i

N⫺1

兺k ␦ ij ⌫ ij kk ⫹ 2 兺k ␦ ij ⌫ ij kk

Combining Eq. 共3.22兲 and Eq. 共3.24兲, we get Eq. 共3.16兲.
Constraining P, Q, and G matrices to be positive
semidefinite is done in essentially the same way as above. In
this case, the variable matrix Y is defined as

E j1 j2 "Y⫽c j1 j2 ⫽2⌫ j1 j2 ⫹ 共 ␦ j1 ␥ j2 ⫹ ␦ j2 ␥ j1 兲
i i

N⫺1
2

i
i
i
i
⫽ ␦ j1 ␥ j2 ⫹ ␦ j2 ␥ j1 ,1
2
2
1

1 2 1 2⫹冊N⫺1 i 1 i 2 j 2 k 2 N⫺1 i 2 i 1 j 1 k 2
k k␦ j ␦k ␦l ␦l ⫹
␦ j ␦ k ␦ l ␦ l 共 Y兲 l 1l 2
1
1 1
2
2
1 1
2
1 2
2
2

and the last term givesN⫺1 i 1 i 2 k
i k
i
共 ␦ j ⌫ j k ⫹ ␦ j2 ⌫ j1 k 兲
2
2
1
1
2

兺k

冉

8285

共3.22兲

1 2

It is convenient to fold our 2-RDM into a compact form,
i i

P j1 j2 → P ij
1 2

共3.31兲

Downloaded 06 Jun 2002 to 130.54.33.130. Redistribution subject to AIP license or copyright, see http://ojps.aip.org/jcpo/jcpcr.jsp

8286

Nakata et al.

J. Chem. Phys., Vol. 114, No. 19, 15 May 2001
TABLE I. Total energy and correlation energy in % in parentheses calculated by the DMVT with P⫹Q andP⫹Q⫹G conditions compared with those obatined by the wave function methods, full CI, and Hartree–Fock.
The basis set is STO-6G except for notice.

System

State

c

1

Be
Be
Bed
Bec
Bed
LiHc
LiH
LiHc
LiH
BeH⫹
BH⫹
BH
CH⫹
CH⫺
CH
NH⫹
NH⫺
NH
OH⫹
OH⫺
OH
HF⫹
HF
BH2
BH2
CH2
CH2
CH2
NH2
NH2
H2O
H2O⫹
FH⫹
2
BH3
CH3
NH3
NH3 共dis兲
H3O⫹

S
S
1
S
3
S
3
S
1 ⫹
⌺
1 ⫹
⌺
3 ⫹
⌺
3 ⫹
⌺
1 ⫹
⌺
2 ⫹
⌺
1 ⫹
⌺
1 ⫹
⌺
3 ⫺
⌺
2
⌸
2
⌸
2
⌸
3 ⫺
⌺
3 ⫺
⌺
1 ⫹
⌺
2
⌸
2
⌸
1 ⫹
⌺
2
A1
2
B1
1
A1
3
B1
3 ⫺
⌺u
2
A1
2
B1
1
A1
2
A1
1
A1
1
A1
2
A2
1
A1
1
A
1
A1
1

ActiveMOa Ele( ␣ ⫹ ␤ ) b
4
5
5
4
5
6
6
6
6
6
6
6
6
6
6
6
6
6
6
6
6
6
6
7
7
7
7
7
7
7
7
7
7
8
8
8
8
8

4(2⫹2)
4(2⫹2)
4(2⫹2)
4(3⫹1)
4(3⫹1)
4(2⫹2)
4(2⫹2)
4(3⫹1)
4(3⫹1)
4(2⫹2)
5(3⫹2)
6(3⫹3)
6(3⫹3)
8(5⫹3)
7(4⫹3)
7(4⫹3)
9(5⫹4)
8(5⫹3)
8(5⫹3)
10(5⫹5)9(5⫹4)
9(5⫹4)
10(5⫹5)7(4⫹3)
7(4⫹3)
8(4⫹4)
8(5⫹3)
8(5⫹3)
9(5⫹4)
9(5⫹4)
10(5⫹5)9(5⫹4)
10(5⫹5)8(4⫹4)
9(5⫹4)
10(5⫹5)10(5⫹5)10(5⫹5)DM( P⫹Q)

DM( P⫹Q⫹G)

Full CI

Hartree–Fock⫺14.5934(176)⫺14.5579(103)⫺14.6064(200)⫺13.3168(120)⫺14.3346(177)⫺8.0034(139)⫺7.9731(104)⫺7.8997(167)⫺7.8554(191)⫺14.8452(106)⫺24.8169(151)⫺25.1234(211)⫺37.9618(227)⫺37.9834(148)⫺38.2472(240)⫺54.4510(248)⫺54.5292(161)⫺54.8280(144)⫺74.7805(138)⫺74.8127(100)⫺75.1164(158)⫺99.1376(153)⫺99.5258(100)⫺25.7549(235)⫺25.7317(233)⫺38.9301(294)⫺38.9043(214)⫺38.8836(187)⫺55.4134(244)⫺55.4856(243)⫺75.7953(232)⫺75.4912(262)⫺99.8894(244)⫺26.4681(258)⫺39.6375(290)⫺56.2061(334)⫺56.1808(326)⫺75.9422(276)⫺14.5827(100)⫺14.5561(100)⫺14.5895(100)⫺13.3146(100)⫺14.3241(100)⫺7.9924(100)⫺7.9724(100)⫺7.8939(98)⫺7.8552(97)⫺14.8439(100)⫺24.8015(100)⫺25.0630(106)⫺37.8896(107)⫺37.9714(99)⫺38.1917(111)⫺54.3957(111)⫺54.5150(99)⫺54.8160(100)⫺74.7719(100)⫺74.8112(95)⫺75.1013(99)⫺99.1278(100)⫺99.5229(89)⫺25.7089(115)⫺25.6837(113)⫺38.8228(119)⫺38.8566(107)⫺38.8358(103)⫺55.3570(111)⫺55.4195(108)⫺75.7310(104)⫺75.4218(106)⫺99.8305(103)⫺26.3932(120)⫺39.5283(117)⫺56.0617(115)⫺56.0394(115)⫺75.8636(103)⫺14.5827(100)⫺14.5561(100)⫺14.5895(100)⫺13.3146(100)⫺14.3241(100)⫺7.9922(100)⫺7.9723(100)⫺7.8940(100)⫺7.8552(100)⫺14.8438(100)⫺24.8015(100)⫺25.0593(100)⫺37.8853(100)⫺37.9718(100)⫺38.1871(100)⫺54.3914(100)⫺54.5151(100)⫺54.8161(100)⫺74.7720(100)⫺74.8127(100)⫺75.1014(100)⫺99.1279(100)⫺99.5258(100)⫺25.7031(100)⫺25.6783(100)⫺38.8110(100)⫺38.8534(100)⫺38.8342(100)⫺55.3525(100)⫺55.4157(100)⫺75.7290(100)⫺75.4192(100)⫺99.8294(100)⫺26.3827(100)⫺39.5178(100)⫺56.0516(100)⫺56.0293(100)⫺75.8621(100)⫺14.5685(0)⫺14.5034(0)⫺14.5725(0)⫺13.3036(0)⫺14.3105(0)⫺7.9635(0)⫺7.9519(0)⫺7.8854(0)⫺7.8549(0)⫺14.8226(0)⫺24.7712(0)⫺25.0015(0)⫺37.8251(0)⫺37.9477(0)⫺38.1443(0)⫺54.3510(0)⫺54.4920(0)⫺54.7887(0)⫺74.7491(0)⫺74.7851(0)⫺75.0756(0)⫺99.1096(0)⫺99.4998(0)⫺25.6649(0)⫺25.6383(0)⫺38.7497(0)⫺38.8089(0)⫺38.7772(0)⫺55.3101(0)⫺55.3670(0)⫺75.6789(0)⫺75.3748(0)⫺99.7879(0)⫺26.3287(0)⫺39.4547(0)⫺55.9855(0)⫺55.9622(0)⫺75.8166(0)aNumber of active MOs.
Number of electrons with the number of ␣ and ␤ electrons in parentheses.
cBasis set is double-␨.
dBasis set is triple-␨.
bby renumbering i⫽i 1 ⫹ 兵 关 i 2 (i 2 ⫺1) 兴 /2其 if i 1 ⬎i 2 and discarding P when i 1 ⭐i 2 . This helps to cut down unnecessary
variables and to automatically assume that 2-RDM has antisymmetric property. Similarly, the Q matrix and other linear
constraints are also folded. Note that the G matrix does nothave such a symmetry property, so that we use all the elements.
The present method involves very large number of linear
constraints and may not be efficient; a merit is that the SDPA
program is used without any modification. However, if wemake a problem-specific SDP solver, it would be much more
efficient than the present one, and such study is now inprogress.
The DMVT formulated above has been applied to theground states of different space and spin symmetries of neutral and charged species of 16 different atoms and molecules.
They are Be( 1 S), Be( 3 S), LiH( 1 ⌺ ⫹ ), LiH( 3 ⌺ ⫹ ), BeH⫹

BH⫹, BH, CH⫹, CH, CH⫺, NH⫹, NH⫺, NH, OH⫹, OH⫺,
OH, HF⫹, HF, BH2( 2 A 1 ), BH2( 2 B 1 ), CH2( 1 A 1 ), CH2( 3 B 1 ),⫹2
2⫹
linear CH2( 3 ⌺ ⫺u ), NH2( A 1 ), NH2( B 1 ), H2O, H2O , FH2 ,
BH3 , CH3 , NH3 , NH3共dis兲 共‘‘dis’’ stands for distorted in thesense that one bond length is shortened by 0.9 time, anotherone is lengthened by 1.1 time兲 and H3O⫹.
We used three different basis sets, double and triple-␨s-type GTOs and STO-6G, for Be, and double-␨ s-type
GTOs by Huzinaga20 and Dunning21 and STO-6G for LiH.
For all the other molecules, we used the STO-6G basis set.22
The geometries we used are the experimental ones.23,24
IV. RESULTS AND DISCUSSION

We show in Table I the total energy of the system calculated by the present method and in parentheses the calculated correlation energy in percentage relative to the

Downloaded 06 Jun 2002 to 130.54.33.130. Redistribution subject to AIP license or copyright, see http://ojps.aip.org/jcpo/jcpcr.jsp

J. Chem. Phys., Vol. 114, No. 19, 15 May 2001

Calculations of second order reduced density matrices

TABLE II. Dipole moments calulated by the DMVT with P⫹Q and P⫹Q⫹G conditions compared with those obtained by the wave function
method. The basis set is STO-6G except for notice.
Molecule
LiHa
LiH
LiHa
LiH
BeH⫹
BH⫹
BH
CH⫹
CH⫺
CH
NH⫹
NH⫺
NH
OH⫹
OH⫺
OH
HF⫹
HF
BH2
CH2
CH2
NH2
NH2
H2O
H2O⫹
FH⫹
2
NH3
NH3 共dis兲H3O⫹
a

TABLE III. Virial coefficients calulated by the DMVT with the P⫹Q andP⫹Q⫹G conditions compared with those obtained by the wave function
methods. The basis set is STO-6G except for notice.

State DM( P⫹Q) DM( P⫹Q⫹G) Full CI Hartree–Fock

System

⌺⫹
⌺⫹
3 ⫹
⌺
3 ⫹
⌺
1 ⫹
⌺
2 ⫹
⌺
1 ⫹
⌺
1 ⫹
⌺
3 ⫺
⌺
2
⌸
2
⌸
2
⌸
3 ⫺
⌺
3 ⫺
⌺
1 ⫹
⌺
2
⌸
2
⌸
1 ⫹
⌺
2
A1
1
A1
3
B1
2
A1
2
B1
1
A1
2
A1
1
A1
1
A1
1
A
1
A1

Bea
Be
Beb
Bea
Beb
LiHa
LiH
LiHa
LiH
BeH⫹
BH⫹
BH
CH⫹
CH⫺
CH
NH⫹
NH⫺
NH
OH⫹
OH⫺
OH
HF⫹
HF
BH2
BH2
CH2
CH2
CH2
NH2
NH2
H2O
H2O⫹
FH⫹
2
BH3
CH3
NH3
NH3 共dis兲
H3O⫹

1

1

1.64451.73720.6225
1.58971.32030.04950.28330.68930.18260.6016
0.89370.13590.47300.99880.06200.44970.96000.54200.00370.24350.08380.51700.64330.59930.87181.03680.69030.66601.4162

1.61641.75230.62581.5906
1.31880.02230.29350.67640.19250.48780.87290.13110.49950.9741
0.06370.47380.99930.53830.03280.50570.08570.54070.6816
0.64600.98571.04290.6901
0.6634
1.4286

1.61921.75190.62581.59071.31960.02230.29940.69050.19290.50440.88040.13210.4996
0.97420.06200.47450.99990.54200.03440.5293
0.09340.55090.68960.6487
0.99201.04370.6922
0.67671.42892.0764
1.93390.62611.5915
1.29870.01970.38060.72530.16690.44060.87890.1431
0.52330.98750.07250.5166
1.0786
0.52280.04660.62240.1006
0.55800.72000.69271.0724
1.05600.69350.69371.4320Basis set is double-␨.

Hartree–Fock 共0%兲 and full-CI 共100%兲 results. Two types of
SDP relaxation calculations are performed. One uses the P
and Q conditions together with the seven conditions given
by Eqs. 共2.3兲–共2.10兲; it is referred to as DM( P⫹Q). The
other uses the G condition additionally and it is denoted asDM( P⫹Q⫹G).
We first examine the results of DM( P⫹Q) calculations.
We see that the results for OH⫺ and HF are excellent, but
this is not a good news but simply due to the too restrictive
variational space: 10 electrons are distributed into six orbitals and therefore in this case P⫹Q condition gives the complete N-representability condition 共2 hole system兲.25 Similarly, the extent of overshooting is relatively small because
the variational space is too restrictive. When the variation isreasonably free, the DM( P⫹Q) energy overshoots too muchthe full-CI energy up to 334% of the full-CI correlation energy for NH3. This result shows that the P⫹Q condition
together with the above seven conditions is still too far fromthe complete N-representability condition.
When we impose further the G condition, we obtain the
results shown under DM( P⫹Q⫹G). They are much improved in comparison with the results of DM( P⫹Q). The
calculated correlation energy percentages range within
100%–110% for atoms and diatomic molecules, while theyrange in 110%–120% for triatomic molecules. This meansthat the G condition is a nice restrictive condition for the

8287

a

State DM( P⫹Q) DM( P⫹Q⫹G) Full CI Hartree–Fock
1

S
S
1
S
3
S
3
S
1 ⫹
⌺
1 ⫹
⌺
3 ⫹
⌺
3 ⫹
⌺
1 ⫹
⌺
2 ⫹
⌺
1 ⫹
⌺
1 ⫹
⌺
3 ⫺
⌺
2
⌸
2
⌸
2
⌸
3 ⫺
⌺
3 ⫺
⌺
1 ⫹
⌺
2
⌸
2
⌸
1 ⫹
⌺
2
A1
2
B1
1
A1
3
B1
3 ⫺
⌺u
2
A1
2
B1
1
A1
2
A1
1
A1
1
A1
2
A2
1
A1
1
A
1
A1
1

1.99751.96212.00171.74591.97742.00381.98321.99081.95792.0036
1.99181.95742.0044
1.9393
1.97812.0164
1.95961.99412.01991.9672
1.99672.02182.00011.9722
1.96991.98491.98891.9882
1.99501.9956
1.99662.01762.0183
1.9835
1.9941
1.9981
1.9973
1.9972

1.99891.96142.00061.74611.9766
1.99771.9826
1.9875
1.95772.0031
1.99191.95652.00401.9392
1.97782.0156
1.95971.99392.0194
1.96711.99652.02122.00011.97271.9702
1.98401.9884
1.98761.9943
1.9955
1.99682.01522.01591.98361.99391.9985
1.99761.9941

1.99891.96142.00061.74611.9766
1.99771.9826
1.9875
1.95772.0031
1.99181.95662.00391.9393
1.97772.0154
1.95971.99382.0194
1.9672
1.99652.02122.00011.97311.97051.98401.9886
1.98781.9942
1.9955
1.99682.01512.01581.9843
1.9944
1.9984
1.9974
1.9941

1.9994
1.95582.00001.7464
1.97591.99291.98371.9826
1.95742.00411.9931
1.95502.00251.9396
1.97732.01441.96011.99392.0183
1.96781.99652.0201
1.99991.97381.9712
1.98411.9886
1.98771.99401.9956
1.99672.01382.0143
1.9853
1.99481.9985
1.99761.9934

Basis set is double-␨.
Basis set is triple-␨.

b

N-representability. We investigated distorted ammonia to examine whether the spatial symmetry affects the
N-representability condition, however, this calculation showsthat there is no effect by such a small distortion; the accuracies of the two calculations are almost the same.
The SDP variational method should give, in principle, alower bound for energy, however, compared to the full-CI
results, the breakdown where the calculated SDPA energy ishigher than the full-CI energy occurs for LiH( 3 ⌺, STO-6G),
LiH( 3 ⌺, double- ␨ ), CH⫺, NH⫺, OH, OH⫺, and HF, though
the violations are within 1 mhartree. It seems that thesebreakdowns are related to the numerical errors in the SDPA
procedure, which we discuss later.
In Table II, we show the 共nonzero兲 dipole moments ofthe molecules calculated here. The dipole moment obtainedat the level of DM( P⫹Q) is not so good. In particular, those
for CH, NH⫹, CH2 , H2O, and H2O⫹ are worse than the
Hartree–Fock results. At the DM( P⫹Q⫹G) level, how-

Downloaded 06 Jun 2002 to 130.54.33.130. Redistribution subject to AIP license or copyright, see http://ojps.aip.org/jcpo/jcpcr.jsp

8288

Nakata et al.

J. Chem. Phys., Vol. 114, No. 19, 15 May 2001
TABLE IV. Number of the constraints and the numerical errors of the DM( P⫹Q) calculations. The basis set
is STO-6G except for notice.

System
Bec
Be
Bed
Bec
Bed
LiHc
LiH
LiHc
LiH
BeH⫹
BH⫹
BH
CH⫹
CH⫺
CH
NH⫹
NH⫺
NH
OH⫹
OH⫺
OH
HF⫹
HF
BH2
BH2
CH2
CH2
CH2
NH2
NH2
H2O
H2O⫹
FH⫹
2
BH3
CH3
NH3
NH3 共dis兲
H3O⫹

State
1

S
S
1
S
3
S
3
S
1 ⫹
⌺
1 ⫹
⌺
3 ⫹
⌺
3 ⫹
⌺
1 ⫹
⌺
2 ⫹
⌺
1 ⫹
⌺
1 ⫹
⌺
3 ⫺
⌺
2
⌸
2
⌸
2
⌸
3 ⫺
⌺
3 ⫺
⌺
1 ⫹
⌺
2
⌸
2
⌸
1 ⫹
⌺
2
A1
2
B1
1
A1
3
B1
3 ⫺
⌺u
2
A1
2
B1
1
A1
2
A1
1
A1
1
A1
2
A2
1
A1
1
A
1
A1
1

ActiveaMOsEle( ␣ ⫹ ␤ ) b
constraints

No. ofconstraints

Primal feasible
error

Dual feasible
error

Gap

4
5
5
4
5
6
6
6
6
6
6
6
6
6
6
6
6
6
6
6
6
6
6
7
7
7
7
7
7
7
7
7
7
8
8
8
8
8

4(2⫹2)
4(2⫹2)
4(2⫹2)
4(3⫹1)
4(3⫹1)
4(2⫹2)
4(2⫹2)
4(3⫹1)
4(3⫹1)
4(2⫹2)
5(3⫹2)
6(3⫹3)
6(3⫹3)
8(5⫹3)
7(4⫹3)
7(4⫹3)
9(5⫹4)
8(5⫹3)
8(5⫹3)
10(5⫹5)9(5⫹4)
9(5⫹4)
10(5⫹5)7(4⫹3)
7(4⫹3)
8(4⫹4)
8(5⫹3)
8(5⫹3)
9(5⫹4)
9(5⫹4)
10(5⫹5)9(5⫹4)
10(5⫹5)8(4⫹4)
9(5⫹4)
10(5⫹5)10(5⫹5)10(5⫹5)

183
440
440
183
440
911
911
911
911
911
911
911
911
911
911
911
911
911
911
911
911
911
911
1692
1692
1692
1692
1692
1692
1692
1692
1692
1692
2897
2897
2897
2897
2897

1.87⫻10⫺14
7.87⫻10⫺14
5.49⫻10⫺14
4.12⫻10⫺11
7.86⫻10⫺12
7.24⫻10⫺14
8.52⫻10⫺14
3.80⫻10⫺12
3.41⫻10⫺11
1.16⫻10⫺13
6.69⫻10⫺14
4.43⫻10⫺14
3.62⫻10⫺14
5.49⫻10⫺12
5.91⫻10⫺14
3.15⫻10⫺14
8.57⫻10⫺12
9.79⫻10⫺12
6.77⫻10⫺12
7.50⫻10⫺12
5.54⫻10⫺12
2.36⫻10⫺12
6.82⫻10⫺12
6.43⫻10⫺14
5.66⫻10⫺14
3.25⫻10⫺14
3.98⫻10⫺14
4.12⫻10⫺14
3.63⫻10⫺14
3.99⫻10⫺14
3.33⫻10⫺14
3.31⫻10⫺14
3.46⫻10⫺14
6.57⫻10⫺14
3.87⫻10⫺14
4.07⫻10⫺14
4.38⫻10⫺14
4.13⫻10⫺14

9.24⫻10⫺11
4.35⫻10⫺11
1.41⫻10⫺10
1.17⫻10⫺7
1.12⫻10⫺7
9.24⫻10⫺9
1.08⫻10⫺9
1.75⫻10⫺6
4.57⫻10⫺8
3.02⫻10⫺11
1.86⫻10⫺10
8.29⫻10⫺10
1.30⫻10⫺9
8.97⫻10⫺8
2.23⫻10⫺99.98⫻10⫺10
7.06⫻10⫺8
4.97⫻10⫺7
8.16⫻10⫺8
5.25⫻10⫺7
1.35⫻10⫺7
8.42⫻10⫺8
1.12⫻10⫺7
4.94⫻10⫺10
8.16⫻10⫺12
3.26⫻10⫺11
5.94⫻10⫺11
6.59⫻10⫺11
1.21⫻10⫺9
2.46⫻10⫺11
2.88⫻10⫺10
3.02⫻10⫺11
3.19⫻10⫺10
7.62⫻10⫺10
1.12⫻10⫺11
5.16⫻10⫺11
5.87⫻10⫺10
5.29⫻10⫺10

1.17⫻10⫺98.03⫻10⫺13
1.89⫻10⫺11
3.87⫻10⫺8
9.58⫻10⫺9
3.46⫻10⫺9
4.76⫻10⫺11
3.64⫻10⫺7
1.89⫻10⫺8
5.73⫻10⫺13
3.26⫻10⫺12
2.41⫻10⫺11
9.62⫻10⫺12
6.05⫻10⫺92.80⫻10⫺11
1.49⫻10⫺11
3.70⫻10⫺9
1.85⫻10⫺8
2.49⫻10⫺97.32⫻10⫺8
3.48⫻10⫺92.96⫻10⫺9
6.65⫻10⫺9
4.65⫻10⫺11
3.20⫻10⫺12
3.19⫻10⫺10
5.75⫻10⫺10
5.88⫻10⫺11
1.22⫻10⫺10
3.91⫻10⫺10
1.29⫻10⫺10
2.64⫻10⫺10
1.12⫻10⫺10
3.37⫻10⫺11
8.68⫻10⫺10
1.02⫻10⫺10
1.62⫻10⫺10
6.77⫻10⫺11aNumber of active MOs.
Number of electrons with the number of ␣ and ␤ electrons in parentheses.
cBasis set is double-␨.
dBasis set is triple-␨.
bever, the dipole moments are drastically improved and all the
results well reproduce the full-CI ones, except for NH3 andNH3共dis兲 for which even Hartree–Fock calculations give
good results and the deviations are very small.
In Table III, we show the virial coefficient 具 V 典 / 具 T 典 ,where 具 V 典 and 具 T 典 denote average potential and kinetic energies, respectively, which must be two for completely variational wave function. When we use DM( P⫹Q⫹G) approximation, the calculated virial is almost completely identicalwith the full-CI result.
Next, we discuss the numerical accuracy of the SDP
method. In Tables IV and V, we summarize the number ofthe constraints and the numerical errors of the DM( P⫹Q)
and DM( P⫹Q⫹G) calculations. The primal feasible error
is defined by

再 冏冋 兺
m

max

X⫺

i⫽1

Fi x i ⫹F0册冏: p,q⫽1,2,...,n

pq

冎

共4.1兲

and the dual feasible error is defined by
max兵 兩 Fi "Y⫺c i 兩 : i⫽1,2,...,m 其 .

共4.2兲

The gap denotes the difference between the primal and dual
functions defined by

冏兺
m

i⫽1

冏

c i x i ⫺F0 "Y .

共4.3兲

These three quantities give criteria of the accuracy of the
SDPA. In the SDPA, our object is the minimization of the

Downloaded 06 Jun 2002 to 130.54.33.130. Redistribution subject to AIP license or copyright, see http://ojps.aip.org/jcpo/jcpcr.jsp

J. Chem. Phys., Vol. 114, No. 19, 15 May 2001

Calculations of second order reduced density matrices

8289

TABLE V. Number of the constraints and the numerical errors of the DM( P⫹Q⫹G) calculations. The basis
set is STO-6G except for notice.

System
Bec
Be
Bed
Bec
Bed
LiHc
LiH
LiHc
LiH
BeH⫹
BH⫹
BH
CH⫹
CH⫺
CH
NH⫹
NH⫺
NH
OH⫹
OH⫺
OH
HF⫹
HF
BH2
BH2
CH2
CH2
CH2
NH2
NH2
H2O
H2O⫹
FH⫹
2
BH3
CH3
NH3
NH3 共dis兲
H3O⫹

State
1

S
S
1
S
3
S
3
S
1 ⫹
⌺
1 ⫹
⌺
3 ⫹
⌺
3 ⫹
⌺
1 ⫹
⌺
2 ⫹
⌺
1 ⫹
⌺
1 ⫹
⌺
3 ⫺
⌺
2
⌸
2
⌸
2
⌸
3 ⫺
⌺
3 ⫺
⌺
1 ⫹
⌺
2
⌸
2
⌸
1 ⫹
⌺
2
A1
2
B1
1
A1
3
B1
3 ⫺
⌺u
2
A1
2
B1
1
A1
2
A1
1
A1
1
A1
2
A2
1
A1
1
A
1
A1
1

ActiveaMOsEle( ␣ ⫹ ␤ ) b

No. ofconstraints

Primal feasible
error

Dual feasible
error

Gap

4
5
5
4
5
6
6
6
6
6
6
6
6
6
6
6
6
6
6
6
6
6
6
7
7
7
7
7
7
7
7
7
7
8
8
8
8
8

4(2⫹2)
4(2⫹2)
4(2⫹2)
4(3⫹1)
4(3⫹1)
4(2⫹2)
4(2⫹2)
4(3⫹1)
4(3⫹1)
4(2⫹2)
5(3⫹2)
6(3⫹3)
6(3⫹3)
8(5⫹3)
7(4⫹3)
7(4⫹3)
9(5⫹4)
8(5⫹3)
8(5⫹3)
10(5⫹5)9(5⫹4)
9(5⫹4)
10(5⫹5)7(4⫹3)
7(4⫹3)
8(4⫹4)
8(5⫹3)
8(5⫹3)
9(5⫹4)
9(5⫹4)
10(5⫹5)9(5⫹4)
10(5⫹5)8(4⫹4)
9(5⫹4)
10(5⫹5)10(5⫹5)10(5⫹5)

983
2365
2365
983
2365
4871
4871
4871
4871
4871
4871
4871
4871
4871
4871
4871
4871
4871
4871
4871
4871
4871
4871
8993
8993
8993
8993
8993
8993
8993
8993
8993
8993
15 31315 31315 31315 31315 313

6.55⫻10⫺11
7.02⫻10⫺12
7.92⫻10⫺11
1.12⫻10⫺10
1.41⫻10⫺11
5.72⫻10⫺12
5.77⫻10⫺12
3.53⫻10⫺11
5.79⫻10⫺11
2.42⫻10⫺11
7.84⫻10⫺12
2.80⫻10⫺11
3.75⫻10⫺12
1.43⫻10⫺11
6.06⫻10⫺12
4.15⫻10⫺12
6.89⫻10⫺12
8.69⫻10⫺11
7.98⫻10⫺12
2.54⫻10⫺11
1.01⫻10⫺11
1.09⫻10⫺11
1.39⫻10⫺11
2.50⫻10⫺12
9.48⫻10⫺12
4.14⫻10⫺12
1.57⫻10⫺10
4.72⫻10⫺12
7.65⫻10⫺12
2.36⫻10⫺12
1.54⫻10⫺12
1.12⫻10⫺11
3.78⫻10⫺12
5.02⫻10⫺11
4.26⫻10⫺12
1.65⫻10⫺12
3.29⫻10⫺12
1.59⫻10⫺12

1.87⫻10⫺7
4.93⫻10⫺7
4.15⫻10⫺7
7.44⫻10⫺7
1.65⫻10⫺7
2.33⫻10⫺5
7.55⫻10⫺8
6.35⫻10⫺7
6.56⫻10⫺7
1.88⫻10⫺7
7.50⫻10⫺8
1.43⫻10⫺5
2.41⫻10⫺7
3.73⫻10⫺7
8.54⫻10⫺6
3.64⫻10⫺6
3.39⫻10⫺6
1.88⫻10⫺5
4.05⫻10⫺7
5.65⫻10⫺6
2.78⫻10⫺6
1.40⫻10⫺6
7.37⫻10⫺6
8.84⫻10⫺8
1.23⫻10⫺6
3.54⫻10⫺6
4.05⫻10⫺5
3.98⫻10⫺7
4.30⫻10⫺6
3.22⫻10⫺8
1.05⫻10⫺7
6.86⫻10⫺6
5.14⫻10⫺7
5.34⫻10⫺7
6.16⫻10⫺7
4.62⫻10⫺7
1.42⫻10⫺6
2.30⫻10⫺7

2.87⫻10⫺6
1.42⫻10⫺6
2.48⫻10⫺6
3.08⫻10⫺6
2.94⫻10⫺7
4.50⫻10⫺6
1.69⫻10⫺6
2.41⫻10⫺5
2.58⫻10⫺6
1.93⫻10⫺6
7.82⫻10⫺7
1.61⫻10⫺8
9.68⫻10⫺7
2.13⫻10⫺5
1.91⫻10⫺7
4.98⫻10⫺7
1.37⫻10⫺52.90⫻10⫺6
1.90⫻10⫺6
2.72⫻10⫺5
6.15⫻10⫺6
6.07⫻10⫺6
3.99⫻10⫺58.00⫻10⫺7
1.59⫻10⫺6
2.69⫻10⫺7
1.28⫻10⫺7
1.58⫻10⫺6
1.30⫻10⫺6
2.16⫻10⫺6
3.63⫻10⫺7
7.63⫻10⫺7
6.07⫻10⫺7
2.88⫻10⫺7
9.01⫻10⫺8
4.39⫻10⫺7
4.14⫻10⫺6
2.24⫻10⫺6aNumber of active MOs.
Number of electrons with the number of ␣ and ␤ electrons in parentheses.
cBasis set is double-␨.
dBasis set is triple-␨.
bdual form of the problem, so that the dual feasible error is an
important quantity, indicating the numerical accuracy of thecalculation.
For DM( P⫹Q), the dual feasible error is in the range of10⫺7 – 10⫺12, while for DM( P⫹Q⫹G), it ranges
10⫺5 – 10⫺8 . As the number of the constraints increases drastically in the P⫹Q⫹G calculations, the numerical accuracy
becomes much worse in the DM( P⫹Q⫹G) results. The gap
value shows the same tendency. The worst five are
HF(3.99⫻10⫺5 ), OH⫺ (2.72⫻10⫺5 ), LiH共double-␨, 3 ⌺;2.41⫻10⫺5 兲, CH⫺ (2.31⫻10⫺5 ), and NH⫺ (1.37⫻10⫺5 ).
We notice that they have the DM( P⫹Q⫹G) energies
higher than the full-CI ones, though these values must be
lower than the full-CI values. There seems to be some relation between the gap value and the numerical accuracy in the
SDPA technique. Another reason is certainly the too small

variational freedom in the calculations of HF and OH⫺; actually in these cases P and Q conditions are already sufficient; the number of holes is 2, so the 2-hole system with Q
condition is just like performing variational calculation forthe 2-electron system with the P condition. Therefore enforcing P, Q, and G conditions is essentially the same as
enforcing P and Q conditions.
The primal feasible values are very small(10⫺12 – 10⫺14) for DM( P⫹Q) calculations and also small(10⫺10 – 10⫺12) for DM( P⫹Q⫹G) calculations. We do notfind any relationship between the accuracies of the present
calculations and the primal feasible errors. So the accuracyof the present calculation seems to be related only to that ofthe primal problem.
In Table VI, we show the occupation numbers 共eigenval-

Downloaded 06 Jun 2002 to 130.54.33.130. Redistribution subject to AIP license or copyright, see http://ojps.aip.org/jcpo/jcpcr.jsp

8290

Nakata et al.

J. Chem. Phys., Vol. 114, No. 19, 15 May 2001
TABLE VI. Occupation number calculated by the DMVT with P⫹Q and P⫹Q⫹G conditions compared withthose obatined by the wave function methods, full CI, and Hartree–Fock, for Be, H2O, and CH2.
System, state, basis

DM( P⫹Q)

DM( P⫹Q⫹G)

Full CI

Hartree–Fock

Be, S, STO-6G

0.036 641⫻4,0.036 642⫻20.890 243⫻2
0.999 832⫻20.035 909⫻60.892 275⫻20.999 997⫻20.035 901⫻60.892 298⫻20.999 998⫻2

0⫻6
1⫻4

Be, 1 S, triple-␨

0.000 202⫻20.001 041⫻2
0.006 149⫻2
0.993 649⫻20.998 959⫻20.000 064⫻20.000 595⫻20.004 153⫻2
0.995 837⫻20.999 352⫻20.000 055⫻20.000 590⫻20.004 119⫻2
0.995 879⫻20.999 357⫻20⫻60.000 004
0.000 187
0.000 573
0.000 645
0.000 702
0.001 511
0.998 534
0.999 252
0.999 293
0.999 299

0.000 000
0.000 001
0.000 011
0.000 013
0.000 707
0.000 712
0.999 280
0.999 286
0.999 990
0.999 998

0.000 000
0.000 000
0.000 007
0.000 009
0.000 707
0.000 711
0.999 284
0.999 287
0.999 995
1.000 000

0⫻60.029 795⫻20.031 585⫻20.970 205⫻2
0.970 570⫻20.998 733⫻20.999 114⫻2
0.999 998⫻20.013 850⫻20.014 766⫻20.986 433⫻20.987 475⫻2
0.998 702⫻20.998 776⫻20.999 999⫻2

0.013 304⫻20.013 509⫻20.986 732⫻20.988 323⫻2
0.998 973⫻20.999 161⫻20.999 999⫻20⫻40.037 796⫻40.314 251⫻20.685 804⫻20.962 204⫻4
0.999 945⫻20.014 336⫻20.016 294⫻2
0.069 501⫻20.929 480⫻20.984 098⫻2
0.986 297⫻2
0.999 993⫻20.010 854⫻20.012 979⫻2
0.050 589⫻2
0.947 809⫻20.987 470⫻2
0.990 307⫻20.999 993⫻2

0⫻6

1

Be, 3 S, triple-␨

H2O, 1 A 1 , STO-6G

CH2 , 1 A 1 , STO-6G

ues of 1-RDM兲 for Be( 1 S, STO-6G), Be( 1 S, triple- ␨ ),Be( 3 S, triple- ␨ ), H2O( 1 A 1 , STO-6G), and CH2( 1 A 1 ,STO-6G). For Be( 1 S, STO-6G), the occupation numbers ofthe 2p orbitals should be sixfold degenerate. Although we
did not impose such constraints, this degeneracy accurately
holds in both DM( P⫹Q) and DM( P⫹Q⫹G) calculations.
For singlet states, both DM( P⫹Q) and DM( P⫹Q⫹G) calculations reproduced the degeneracy of the twofold occupation without constraints. Generally, the occupation numbersof the DM calculations are much more distributed over all
the natural orbitals than those of the full CI. Although such atrend is reduced for the DM( P⫹Q⫹G) calculation, it contradicts our expectation; the occupation numbers are ex-

1⫻4

1⫻4

1⫻10

1⫻8

pected to be less distributed in the calculations with lesssufficient N-representability conditions. An extreme case
was CH2 , this tendency is very amplified and the accidental
degeneracy of occupation are found in the DM( P⫹Q) calculation.
In Table VII, the root-mean-square 共rms兲 deviation d ofthe 2-RDM from the full CI,d⫽冑兺 兵i1i2 j 1 j 2

i i

i i

1 2

1 2

共 ⌫ calculated兲 j1 j2 ⫺ 共 ⌫ full CI兲 j1 j2 其 2

共4.4兲

is presented for the systems examined in Table VI. The deviations of the 2-RDM are quite small in DM( P⫹Q⫹G)

TABLE VII. RMS deviations of the 2-RDMs calculated by DMVT with P⫹Q and P⫹Q⫹G conditions fromthose by full CI for Be, H2O, and CH2.
System, state, basis

DM( P⫹Q)

DM( P⫹Q⫹G)

Full CI

Hartree–Fock

Be, 1 S, STO-6G
Be, 1 S, triple-␨Be, 3 S, triple-␨
H2O, 1 A 1 , STO-6G
CH2 , 1 A 1 , STO-6G

0.049 208
0.049 615
0.029 401
0.467 084
1.604 712

0.000 162
0.003 567
0.000 715
0.029 694
0.153 503

0
0
0
0
0

0.526 569
0.100 331
0.039 123
0.266 154
0.484 788

Downloaded 06 Jun 2002 to 130.54.33.130. Redistribution subject to AIP license or copyright, see http://ojps.aip.org/jcpo/jcpcr.jsp

J. Chem. Phys., Vol. 114, No. 19, 15 May 2001

Calculations of second order reduced density matrices

TABLE VIII. Comparison of largest eigenvalues P, Q, G-matrices andsmallest G-matrix calculated by DMVT with P⫹Q and P⫹Q⫹G conditions compared with those obtained by the wave function methods, full CI,
and Hartree–Fock for Be, H2O, and CH2.
System, state,basis

DM( P⫹Q)

DM( P⫹Q⫹G)

Full CI

Hartree–Fock

Be, 1 S, STO-6G
Largest P
1.005 284
Largest Q
2.264 502
Largest G
3.675 970
Smallest G⫺0.007 288

1.000 907
2.259 861
3.682 340
0.000 000

1.000 874
2.259 834
3.682 407
0.000 00012
40

Be, 1 S, triple-␨Largest P
1.004 751
Largest Q
2.008 541
Largest G
3.978 475
Smallest G⫺0.002 482

1.000 672
2.002 209
3.985 686
0.000 000

1.000 612
2.002 072
3.985 726
0.000 00012
4
0

Be, 3 S, triple-␨Largest P
1.002 019
Largest Q
2.005 387
Largest G
3.994 883
Smallest G⫺0.001 439

0.999 991
2.000 230
3.997 833
0.000 000

0.999 995
2.000 233
3.997 849
0.000 00012
40

H2O, 1 A 1 , STO-6G
Largest P
1.034 378
Largest G
9.805 708
Smallest G⫺0.188 367

1.010 041
9.900 236
0.000 000

1.008 652
9.904 3120.000 0001
100

CH2 , 1 A 1 , STO-6G
Largest P
1.090 030
Largest Q
2.227 026
Largest G
6.942 144
Smallest G⫺0.208 966

1.024 402
2.059 546
7.679 238
0.000 001

1.019 529
2.042 417
7.746 013
0.000 0001280calculation especially for small systems, where DM( P⫹Q⫹G) give the identical total energy and virial coefficient tofull CI. However, DM( P⫹Q) calculations gave worse results and even worse than Hartree–Fock for H2O and CH2.
In Table VIII, we compare the largest eigenvalues of P,
Q, and G-matrices and smallest eigenvalues of G-matrix, forthe same systems. Largest eigenvalues of P and Q-matrices
become smaller as the calculation quality becomes better,
while those of G-matrix become larger. In DM( P⫹Q) calculations, smallest eigenvalues of G-matrix are negative. As
we expected, smallest eigenvalue of G-matrix becomes
smaller when electron correlation gets larger. We did notshow the smallest eigenvalues of P and Q-matrices since inany case, they are almost zero 共absolute values are smaller
than 10⫺6 兲. The deviation of these values that are large for
the CH2 关largest eigenvalue of G-matrix for DM( P⫹Q⫹G)兴 calculation is 7.679 238 compared to the full CI’s one
7.746 013, while the SDPA errors are small 共primal and dual
feasibilities are 4.14⫻10⫺12 and 3.54⫻10⫺6 , respectively,
and gap is 2.69⫻10⫺7 兲. Therefore, we conclude that the error originates from the insufficiency of the N-representability
conditions rather than that of the SDPA.
The trace of the Q matrix is normalized to (r⫺N)⫻(r⫺N⫹1), where r is number of MO 共or rank of 1-RDM兲 and
N is the number of the electrons. This condition is satisfied
when we impose the constraint for the number of the electrons.
Lastly, we note that we find essentially no problem in

8291

finding the minimum and this should be the case for othersystems. This is certainly a merit of the present method.
V. CONCLUSION

The DMVT is developed systematically by using SDPA
as a problem solver. This technique is very stable and there
were no example where we could not get a convergence. In
addition to several trivial conditions, the P⫹Q condition is
insufficient, while the P⫹Q⫹G condition gives satisfactory
results, the extent of overshooting the full-CI energy being
small for the systems presently examined. The dipole moment and the virial coefficient calculated by the DM( P⫹Q⫹G) method are also very close to the full-CI values. Thismethod is applicable to the ground state of any spin- andspace symmetry of closed and open-shell systems.
In this DMVT approach, the calculated energy is a lower
bound of the exact energy. The errors of the present DM( P⫹Q⫹G) method are permissible in both energy and properties. Though most quantum chemical method available
give the ground-state energy higher than the full-CI one, thepresent method giving lower energy is equally permissible as
an approximate quantum chemical method, if it is stable andfeasible in cost performance. For the second requirement, thepresent stage of the theory is an infant stage, but muchprogress is expected in future.
ACKNOWLEDGMENT

This study has been supported by the Grand-in-Aid for
Scientific Research from the Japanese Ministry of Education,
Science, Culture, and Sports.
K. Husimi, Proc. Phys. Math. Soc. Jpn. 22, 264 共1940兲.
P.-O. Löwdin, Phys. Rev. 99, 1474 共1955兲.
3
H. Nakatsuji, in Many-Electron Densities and Reduced Density Matrices,
edited by J. Cioslowski 共Kluwer Academic, New York, 2000兲.
4H. Nakatsuji, Phys. Rev. A 14, 41 共1976兲.
5
H. Nakatsuji, Theor. Chem. Acc. 102, 97 共1999兲.
6
H. Nakatsuji and K. Yasuda, Phys. Rev. Lett. 76, 1039 共1996兲.
7
K. Yasuda and H. Nakatsuji, Phys. Rev. A 56, 2648 共1997兲.
8
M. Ehara, M. Nakata, H. Kou, K. Yasuda, and H. Nakatsuji, Chem. Phys.
Lett. 305, 483 共1999兲.
9
M. Nakata, M. Ehara, K. Yasuda, and H. Nakatsuji, J. Chem. Phys. 112,8772 共2000兲.
10A. J. Coleman, Rev. Mod. Phys. 35, 668 共1963兲.
11
C. Garrod and J. Percus, J. Math. Phys. 5, 1756 共1964兲.
12
C. Garrod and M. A. Fusco, Int. J. Quantum Chem. x, 495 共1976兲.
13
C. Garrod, M. V. Mihailović, and M. Rosina, J. Math. Phys. 16, 868
共1975兲.
14
M. V. Mihailović and M. Rosina, Nucl. Phys. A 237, 229 共1975兲.
15
K. Fujisawa, M. Kojima, K. Nakata, SDPA 共SemiDefinite Programming
Algorithm兲 User’s Manual Version 5.00, August 1999, ftp://
ftp.is.titech.ac.jp/pub/OpRes/software/SDPA/5.00/.
16
R. M. Erdahl and C. Garrod, in Proceedings of the Conference on Density
Matrices, Kingston, June 1974, edited by A. J. Coleman and R. M. Erdahl.
17
H. Kummer, J. Math. Phys. 8, 2063 共1967兲.
18
Yu. Nesterov and A. S. Nemirovskii, Interior Point Polynomial Method inConvex Programming: Theory and Applications 共SIAM, Philadelphia,
1993兲.
19
M. Kojima, Semidefinite Programming and Interior-Point Methods,
http://www.is.titech.ac.jp/˜kojima/wabun.html, 1996
共in Japanese兲.
20
S. Huzinaga, J. Chem. Phys. 42, 1293 共1965兲.
21
T. H. Dunning, Jr., J. Chem. Phys. 53, 2823 共1970兲; T. H. Dunning, Jr.
and P. J. Hay, in Method of Electronic Structure Theory, edited by H. F.
Schaefer III 共Plenum, New York, 1977兲, Vol. 2; Actually, basis sets are1
2

Downloaded 06 Jun 2002 to 130.54.33.130. Redistribution subject to AIP license or copyright, see http://ojps.aip.org/jcpo/jcpcr.jsp

8292

J. Chem. Phys., Vol. 114, No. 19, 15 May 2001

taken from EMSL Gaussian Basis Set Order Form, http://
www.emsl.pnl.gov:2080/forms/basisform.html.
22
W. J. Hehre, R. F. Stweart, and J. A. Pople, J. Chem. Phys. 51, 2657
共1969兲.
23
K. P. Huber and G. Herzberg, Molecular Spectra and Molecular Structure

Nakata et al.
IV, Electronic Constants of Diatomic Molecules 共Van Nostrand Reinhold,
New York, 1979兲.
24
J. H. Callomon, E. Horita, K. Kuchitsu, W. J. Lafferty, A. G. Maki, andC.S. Pote, Landolt–Börnstein 共Springer-Verlag, Berlin, 1976兲.
25
A. J. Coleman, Rep. Math. Phys. 4, 113 共1973兲.

Downloaded 06 Jun 2002 to 130.54.33.130. Redistribution subject to AIP license or copyright, see http://ojps.aip.org/jcpo/jcpcr.jsp