JOURNAL OF CHEMICAL PHYSICSVOLUME 116, NUMBER 13

1 APRIL 2002

Density matrix variational theory: Application to the potential energy
surfaces and strongly correlated systems
Maho Nakata, Masahiro Ehara, and Hiroshi Nakatsujia)
Department of Synthetic Chemistry and Biological Chemistry, Faculty of Engineering, Kyoto University,
Kyoto 606-8501, Japan

共Received 23 August 2001; accepted 7 January 2002兲
The density matrix variational theory 共DMVT兲 algorithm developed previously 关J. Chem. Phys. 114,
8282 共2001兲兴 was utilized for calculations of the potential energy surfaces of molecules, H4 , H2 O,
NH3 , BH3 , CO, N2 , C2 , and Be2 . The DMVT( PQG), using the P, Q, and G conditions as
subsidiary condition, reproduced the full-CI curves very accurately even up to the dissociation limit.
The method described well the quasidegenerate states and the strongly correlated systems. On theother hand, the DMVT( PQ) was not satisfactory especially in the dissociation limit and its potential
curves were always repulsive. The size consistency of the method was discussed and the G
condition was found to be essential for the correct behavior of the potential curve. Further, we also
examined the Weinhold–Wilson inequalities for the resultant 2-RDM of DMVT( PQG)
calculations. Two linear inequalities were violated when the results were less accurate, suggestingthat this inequality may provide a useful N-representability condition for the DMVT. © 2002
American Institute of Physics. 关DOI: 10.1063/1.1453961兴

I. INTRODUCTION

powerful at that time, so that their methods were applicable
only to extremely small systems from the limitation in thenumber of variational parameters.
Mazziotti6 applied variational principle to Lipkin model
employing positive semidefiniteness of the 4-RDM and thedensity equation. Erdahl and Jin14 considered a merit of using higher 共than 2兲 order RDM in the DMVT for the existence of more effective N-representability condition for themodel system of one-dimensional periodic lattice of electron
pairs. They extended the work of Garrod and Percus forhigher order RDMs, and gave some insights for using higher
order RDM as a basic variable.
In our previous study,15 we could efficiently implement
the DMVT using the semidefinite programming algorithm
共SDPA兲16 –19 and succeeded to calculate the 2-RDM of theground state of different symmetry for many atoms and molecules. We transformed the DMVT to the standard type problem of SDP. We showed that the positive semidefinitenessconditions of the P, Q, 8 and G 9 matrices were very strong
for atoms and molecules, though they are only necessary
conditions of the N-representability.
In the recent work of Mazziotti and Erdahl,20 positive
semidefinite condition of 3- and 4-RDMs were examined forsolving the DMVT coupled with DET. They demonstrated its
performance for a boson model of two-energy-level systemwith N⫽10– 75. Valdemoro et al. also considered the functional reconstruction with respect to the ensemble representability conditions.21
Another promising approach was initiated by one of theauthors.22–25 Since the exact ⌿ is an eigenfunction of the
Hamiltonian that has simple structure composed of only oneand two-body operators, the ⌿ itself should also have asimple structure reflecting the simple structure of the Hamiltonian. Some explicit expressions of the structure of the ex-

The second-order reduced density matrix 共2-RDM兲 completely describes the N-body fermion system since any observable properties of the system can be calculated from the2-RDM.1,2 This fact has motivated us to use 2-RDM as abasic variable of quantum mechanics instead of the wave
function ⌿. If we can determine 2-RDM without using ⌿,we have a closed form of quantum mechanics where thebasic variable is 2-RDM. We refer to such formalism ofquantum mechanics as density matrix theory共DMT兲. There
are two categories in the DMT with respect to the determination of the RDM. One is based on the density equation,3
which is equivalent to the Schrödinger equation in the necessary and sufficient sense. This approach is called density
equation theory 共DET兲. Recently, DET is extensively studied
and developed.4 – 6 They have been summarized in a recent
review paper.7 The other is based on the Ritz variational
principle expressed in terms of 2-RDM. This latter approachis called density matrix variational theory 共DMVT兲. The key
in this approach is how well we can restrict our variable
2-RDM to be N-representable.8
Garrod and Percus9 first formulated the DMVT. Kijewski applied the DMVT to C2⫹ and found that the G condition
was a rather strong condition.10 Garrod et al.11,12 also implemented their method and calculated the ground state of Be
atom very accurately. Erdahl proposed to use the convex
program for solving the DMVT and performed accurate calculation for the He2 molecule.13 Afterwards, the interest forsolving 2-RDM using the DMVT has almost disappeared forabout 20 years. The reasons were probably that there was no
rigorous mathematical and computational algorithm for the
DMVT calculation, and the computer facilities were not soa兲Electronic mail: hiroshi@sbchem.kyoto-u.ac.jp

0021-9606/2002/116(13)/5432/8/$19.00

5432© 2002 American Institute of Physics

Downloaded 16 Apr 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. 116, No. 13, 1 April 2002

Density matrix variational theory

act wave function were given and the theories for the groundand excited states was formulated and applied to a simple
model system.
In this paper, we extensively apply our DMVT to calculations of the potential energy surfaces of molecules. Previously, we applied our DET to calculations of the potential
energy curves of small molecules.26 Though the results wereencouraging around the equilibrium and elongated geometry,the calculation failed to converge at large internuclear distances. Here, special attentions are paid to the performanceof DMVT for describing the electronic state of strongly correlated systems and the multiconfigurational systems. We
also discuss the size-consistency property of the method in
connection with the N-representability condition. We will
also examine the Weinhold–Wilson inequalities27–29 for theobtained 2-RDM and consider their possibilities as another
N-representability conditions in our method.
II. THEORYFirst and second order reduced density matrices 共1-,
2-RDMs兲, ␥ and ⌫, are defined by
共2.1兲

i i

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

1 2

2

共2.2兲

where a † and a are creation and annihilation operators, respectively. Practical complete N-representability condition is
not known for 2-RDM: we know only some necessary conditions. In the present DMVT, we use P, Q, and G conditions. The P, Q, and G matrices are defined by
i i
P j1 j2 ⫽ 具 ⌿ 兩 a i† a i† a j 2 a j 1 兩 ⌿ 典 ,
1 2
1 2
i i

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

1 2

i i

1

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

1 2

2

lem into the SDP,16 –18 and we employ SDPA19 as a standard
SDP solver. Details were described in Ref. 15.

Weinhold and Wilson,27 Davidson,28 and McRae andDavidson29 derived some other N-representability conditions
that were expressed as linear inequalities using only the diagonal elements of 2-RDM. Among them, the conditions independent from those already used in our present method are
as follows.
Condition VI:jk␥ ii ⫺2⌫ ii jj ⫺2⌫ ik
ik ⫹2⌫ jk ⭓0.

共2.4兲

jk
1⫺ ␥ ii ⫺ ␥ jj ⫺ ␥ kk ⫹2⌫ ii jj ⫹2⌫ ik
ik ⫹2⌫ jk ⭓0.

共2.6兲

where H is the Hamiltonian of the system, P (2) is a set of
2-RDM that satisfy approximate or nearly complete
N-representability condition. We did two types of calculations using the approximate N-representability conditions:
one is with the trivial representability condition plus P and
Q condition, denoted as DMVT( PQ), and the other is withthe trivial condition plus P, Q, and G conditions, denoted as
DMVT( PQG).
For implementing the minimization problem with theselinear and semidefiniteness conditions, we casted this prob-

冉

冊

共2.8兲

Condition VIII: Positive semidefiniteness of the ⍀ matrix

共2.5兲⌫苸P (2)

共2.7兲

Condition VII:

共2.3兲

respectively. We enforce all of these matrices to be positive
semidefinite. We also use seven trivial conditions of 2-RDM,
which are antisymmetric condition, hermiticity, trace condition, number of electrons, number of spins, and expectation
values of S z and S 2 .
In the DMVT, we take 2-RDM as a variational variable,
and minimize the energy within N-representability conditions, namely,
E min⫽ Min Tr H⌫,

FIG. 1. Coordinates for H4 .

B. Additional linear inequalities for the density
matrices

A. Definitions and basic algorithm

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

5433

⍀⫽

␥ 11

12
13
2⌫ 12
2⌫ 13
¯

1t
2⌫ 1t

␥ 11

12
2⌫ 12

␥ 22

23
2⌫ 23
¯

2t
2⌫ 2t

␥ 22

12
2⌫ 12

13
2⌫ 13

␥ 33

3t
2⌫ 3t

␥ 33

]

]

]

¯

2t
2⌫ 2t

3t
2⌫ 3t

␥ tt

]

1t
2⌫ 1t

¯

]

␥ tt

␥ 11

␥ 22

␥ 33

¯

␥ tt

1

.

共2.9兲

These are the representability conditions that may be stronger than and/or may reinforce the P, Q, and G conditions.
Since these conditions are given as linear inequalities, it iseasy to include them into the present DMVT formalism
within the SDP formalism, since SDP is an extension of thelinear programming. In this study, we examine the resultant
2-RDM against these three inequalities, the condition VI,
VII, and VIII, and discuss the possibility of using these conditions as the additional constraints in our DMVT formalism.
C. Size-consistency

The positive semidefiniteness of the G matrix includes anecessary condition for size consistency. In the original nonlinear form, the position representation of the G matrix is
given byDownloaded 16 Apr 2002 to 130.54.33.130. Redistribution subject to AIP license or copyright, see http://ojps.aip.org/jcpo/jcpcr.jsp

5434

J. Chem. Phys., Vol. 116, No. 13, 1 April 2002

Nakata, Ehara, and Nakatsuji

G 共 12兩 1 ⬘ 2 ⬘ 兲 ⫽ 具 共 ␺ † 共 2 兲 ␺ 共 1 兲 ⫺ 具 ␺ † 共 2 兲 ␺ 共 1 兲 典 兲 † 共 ␺ † 共 2 ⬘ 兲 ␺ 共 1 ⬘ 兲 ⫺ 具 ␺ † 共 2 ⬘ 兲 ␺ 共 1 ⬘ 兲 典 兲 典⫽ 具 共 ␺ † 共 1 兲 ␺ 共 2 兲 ⫺ 具 ␺ † 共 2 兲 ␺ 共 1 兲 典 * 兲共 ␺ † 共 2 ⬘ 兲 ␺ 共 1 ⬘ 兲 ⫺ 具 ␺ † 共 2 ⬘ 兲 ␺ 共 1 ⬘ 兲 典 兲 典⫽ 具 ␺ †共 1 兲 ␺ 共 2 兲 ␺ †共 2 ⬘ 兲 ␺ 共 1 ⬘ 兲 典 ⫺ 具 ␺ †共 2 兲 ␺ 共 1 兲 典 *具 ␺ †共 2 ⬘ 兲 ␺ 共 1 ⬘ 兲 典⫺ 具 ␺ † 共 1 兲 ␺ 共 2 兲 典具 ␺ † 共 2 ⬘ 兲 ␺ 共 1 ⬘ 兲 典 ⫹ 具 ␺ † 共 2 兲 ␺ 共 1 兲 典 * 具 ␺ † 共 2 ⬘ 兲 ␺ 共 1 ⬘ 兲 典⫽ 具 ␺ † 共 1 兲 ␺ 共 2 兲 ␺ † 共 2 ⬘ 兲 ␺ 共 1 ⬘ 兲 典 ⫺ 具 ␺ † 共 1 兲 ␺ 共 2 兲 典具 ␺ † 共 2 ⬘ 兲 ␺ 共 1 ⬘ 兲 典 ,where ␺ (i) is a field operator defined by using one-particle
complete basis set 兵 ␺ j 其 ,␺ 共 i 兲⫽ 兺 ␺ j共 i 兲a j .

共2.11兲

j

A static density–density autocorrelation function F(1 兩 1 ⬘ ) 30
corresponds to the G matrix as
F 共 1 兩 1 ⬘ 兲 ⫽ 具 n 共 1 兲 n 共 1 ⬘ 兲 典 ⫺ 具 n 共 1 兲 典具 n 共 1 ⬘ 兲 典⫽ 具 ␺ †共 1 兲 ␺ 共 1 兲 ␺ †共 1 ⬘ 兲 ␺ 共 1 ⬘ 兲 典⫺ 具 ␺ 共 1 兲 ␺ 共 1 兲 典具 ␺ 共 1 ⬘ 兲 ␺ 共 1 ⬘ 兲 典††⫽G 共 11兩 1 ⬘ 1 ⬘ 兲 ,where n(i) is the density operator defined by
n 共 i 兲 ⫽ ␺ †共 i 兲 ␺ 共 i 兲 .
Using the positive semidefiniteness of G(12兩 1 ⬘ 2 ⬘ ),

冕

共2.12兲
共2.13兲

共2.10兲

where x(1) is an arbitrary one-particle function. From thepositive semidefiniteness of F(1 兩 1 ⬘ ), it is shown thatF(1 兩 1 ⬘ ) is everywhere non-negative.
The size consistency requires more strict condition;
when 兩 1⫺1 ⬘ 兩 →⬁, F(1 兩 1 ⬘ ) should asymptotically go tozero, namely,lim兩 1⫺1 ⬘ 兩 →⬁F 共 1 兩 1 ⬘ 兲 ⫽0.

共2.16兲

The positive semidefiniteness of the G matrix guarantees
only the non-negativity of F(1 兩 1 ⬘ ), but does not guarantee
this asymptotical condition. Thus, the DMVT( PQG) includes a necessary condition for the size-consistency, whilein the DMVT( PQ), even F(1 兩 1 ⬘ ) is not necessarily nonnegative.
III. RESULTS AND DISCUSSIONS

x 共 12兲 G 共 12兩 1 ⬘ 2 ⬘ 兲 x 共 1 ⬘ 2 ⬘ 兲 * d ␶ 1 d ␶ 2 d ␶ 1 ⬘ d ␶ 2 ⬘ ⭓0,

共2.14兲
where x(12) is an arbitrary two particle function, F(1 兩 1 ⬘ ) isshown to be also positive semidefinite by integrating the G
matrix with respect to the two particle function x(12) givenby x(12)⫽x(1) ␦ (1⫺2), as

A. H4 system

共2.15兲

First, we applied our DMVT to the potential energy surface of H4 . This system has been frequently used as a benchmark molecule of many methods for the quasidegeneratesituation:31,32 the a g b 2u and a g b 3u configurations become
equivalent for a square geometry and therefore, become degenerate. We used the DZ basis set33,34 for H and defined thepotential energy surface with the coordinates ( ␪ ,R) depicted
in Fig. 1. R gives the size of the molecule and ␪ defines theasymmetry of the structure. We calculated three different
cuts of the potential energy surface that were also tested inRef. 31.
First, we examined the cut of stretching R with ␪ ⫽90°
fixed, namely, the square structure as a function of R. The
results are summarized in Table I. The DMVT( PQG) reproduced the full-CI curve quite accurately. For large R, it gave
almost identical total energy and the errors were within 1
mhartree for R⬎1.0 Å, though the total correlation energies

TABLE I. Total energy and correlation energy in 共%兲 for H4 as a function of
R with ␪ fixed at 90 degree.

TABLE II. Total energy and correlation energy in 共%兲 for H4 as a functionof ␪ with R fixed at the equilibrium value of 0.869 Å.

0⭐
⫽

冕
冕

x 共 12兲 G 共 12兩 1 ⬘ 2 ⬘ 兲 x 共 1 ⬘ 2 ⬘ 兲 * d ␶ 1 d ␶ 2 d ␶ 1 ⬘ d ␶ 2 ⬘x 共 1 兲 ␦ 共 1⫺2 兲 G 共 12兩 1 ⬘ 2 ⬘ 兲 x 共 1 ⬘ 兲 * ␦ 共 1 ⬘ ⫺2 ⬘ 兲 *⫻d ␶ 1 d ␶ 2 d ␶ 1 ⬘ d ␶ 2 ⬘

⫽
⫽

冕
冕

x 共 1 兲 G 共 11兩 1 ⬘ 1 ⬘ 兲 x 共 1 ⬘ 兲 * d ␶ 1 d ␶ 1 ⬘x共 1 兲F共 1兩1 ⬘ 兲x共 1 ⬘ 兲* d ␶ 1 d ␶ 1⬘ ,R(Å)

DMVT( PQ)DMVT( PQG)

Full-CI

Hartree–Fock␪共degrees兲

DMVT( PQ)DMVT( PQG)

Full-CI

Hartree–Fock

0.6
0.8
1.0
1.2
1.4
1.6
1.8
2.0

⫺2.0405(186)⫺2.1485(168)⫺2.1881(177)⫺2.2210(191)⫺2.2407(194)⫺2.2367(183)⫺2.2226(141)⫺2.2055(141)⫺1.9553(104)⫺2.0629(101)⫺2.0693(101)⫺2.0480(100)⫺2.0251(100)⫺2.0087(100)⫺1.9993(100)⫺1.9945(100)⫺1.9511(100)⫺2.0610(100)⫺2.0684(100)⫺2.0474(100)⫺2.0246(100)⫺2.0085(100)⫺1.9992(100)⫺1.9945(100)⫺1.8474(0)⫺1.9330(0)⫺1.9122(0)⫺1.8568(0)⫺1.7939(0)⫺1.7340(0)⫺1.4551(0)⫺1.4818(0)

90.0
89.9
89.5
89.0
88.0
85.0
80.0
70.0

⫺2.1656(170)⫺2.1656(170)⫺2.1656(172)⫺2.1655(174)⫺2.1654(177)⫺2.1673(180)⫺2.1869(185)⫺2.2337(185)⫺2.0711(101)⫺2.0711(101)⫺2.0713(101)⫺2.0718(101)⫺2.0738(101)⫺2.0849(102)⫺2.1120(101)⫺2.1727(101)⫺2.0697(100)⫺2.0697(100)⫺2.0698(100)⫺2.0703(100)⫺2.0721(100)⫺2.0830(100)⫺2.1106(100)⫺2.1721(100)⫺1.9326(0)⫺1.9335(0)⫺1.9372(0)⫺1.9418(0)⫺1.9509(0)⫺1.9777(0)⫺2.0205(0)⫺2.0992(0)Downloaded 16 Apr 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. 116, No. 13, 1 April 2002

Density matrix variational theory

TABLE III. Total energy and correlation energy in 共%兲 for H4 as a functionof ␪ with R fixed at 1.738 Å.

␪共degrees兲

DMVT( PQ)DMVT( PQG)

Full-CI

Hartree–Fock

90.0
89.9
89.5
89.0
88.0
85.0
80.0
70.0

⫺2.2275(174)⫺2.2275(174)⫺2.2276(175)⫺2.2276(175)⫺2.2278(177)⫺2.2289(181)⫺2.2327(188)⫺2.2465(207)⫺2.0016(100)⫺2.0015(100)⫺2.0017(100)⫺2.0019(100)⫺2.0023(100)⫺2.0041(100)⫺2.0087(100)⫺2.0258(100)⫺2.0015(100)⫺2.0015(100)⫺2.0015(100)⫺2.0015(100)⫺2.0018(100)⫺2.0033(100)⫺2.0080(100)⫺2.0255(100)⫺1.6962(0)⫺1.6967(0)⫺1.6988(0)⫺1.7014(0)⫺1.7067(0)⫺1.7231(0)⫺1.7523(0)⫺1.8198(0)were large, for example, 0.51 a.u. for R⫽2.0 Å. The method
was found to give a good description for the quasidegenerate
system. On the other hand, the DMVT( PQ) gave 40– 90%
errors of the correlation energies.
Second, the cut of ␪ ranging from 70.0° to 90.0° withR⫽0.869 Å, which is near equilibrium distance is examinedin Table II. At ␪ ⫽90°, electronic state becomes quasidegenerate. The DMVT( PQG) gave very smooth potential curve
parallel to the full-CI without artificial cusp at ␪ ⫽90°. 31 The
deviations were within 2 mhartree and 2% of the total correlation energy throughout the geometries. For this system, theerrors were constant regardless of the quasidegeneracy.
Last, the cut of ␪ ⫽70° – 90° with R elongated to 1.738
Å, namely, 2⫻R e , is examined in Table III. Surprisingly, theDMVT( PQG) gave almost identical results with the full-CI
ones: the deviations were less than 1 mhartree for all thegeometries. Though the present calculations did not include
polarization functions, the DMVT( PQG) gave very accurate
potential energy surface of H4 .
For this system all the Weinhold–Wilson inequalities
were satisfied for all the potential energy surfaces examined
here. This also supports the high quality of 2-RDM calculated by the DMVT( PQG).

B. Ne and the equilibrium geometry of N2 , CO, C2 ,LiF, and CH4

Next, the DMVT is applied to the ground state of Ne,
N2 , CO, C2 , LiF, and CH4 , which were not calculated in theprevious study.15 In Table IV, we summarized the total energy for these systems. In all calculations, we adopted
STO-6G minimal basis35 and experimental geometries36,37
except for Ne. For Ne, 关 3s2 p 兴 basis set was used. The 1s
orbitals of the second raw atoms were fixed as cores.

5435

Generally, the results of DMVT( PQG) calculations
were satisfactory except for C2 and CH4 . The DMVT( PQ)
calculations overshoot the energy of these molecules, especially for C2 , by 802%. The DMVT( PQG) recovered it up
to 117%. The deviation is still not small, however, the convergence to the exact value is encouraging since the groundstate of C2 is known to be quasidegenerate even at the equilibrium geometry. We also obtained remarkable improvementfor other systems by requiring the G condition.
The Weinhold–Wilson inequalities VI, VII, and VIII
were examined for these systems. For CO, LiF, and Ne, all of
the inequalities were satisfied. For CH4 and C2 , the inequality type VI and VII were violated, and for N2 , the inequality
type VI was violated but others were satisfied, though theviolations were very small as ⫺0.000 297–⫺0.002 063
type VI),⫺0.000 127– 0.001 803(CH4 , type VII),(CH4 ,⫺0.004 169–⫺0.019 446(C2 , type VI兲, ⫺0.000 169–
⫺0.030 852(C2 , type VII兲 and ⫺0.002 83–⫺0.00551(N2 , type VI), respectively. These violations were parallelto the errors of the DMVT( PQG) calculations: the deviations CH4 (124%), C2 (117%), and N2 (108%) were larger
than those of the other systems. This implies that the inequalities VI and VII may be adopted as one of the additional
N-representability conditions for the DMVT.

C. Potential curves of H2 O, NH3 , and BH3

Double dissociation of H2 O and triple dissociation of
NH3 and BH3 are interesting examples, since four and six
electrons are correlated in the bond dissociation processes.
We calculated the potential curves for the symmetric stretching mode of these systems at several points within R⫽0.5– 5.0 Å, and the results were shown in Figs. 2– 4. We
used STO-6G basis set and kept 1s orbitals of O, N, and B tobe frozen. Spectroscopic constants of equilibrium distance(r e ), harmonic frequency ( ␻ e ), and dissociation energy(D e ) were summarized in Table V. The potential energy
curve was fit with the 6th extended Rydberg function forsome points near the equilibrium geometry and the ␻ e wascalculated by the Dunhum method.38 The H–O–H and
H–N–H angles were fixed at the experimental values andonly the H–O and H–N bonds were symmetrically stretched:␻ e was defined for this coordinate and therefore different
from that of the normal mode analysis.
For H2 O and NH3 , DMVT( PQG) simulated the full-CI
curves very accurately even up to the dissociation limit andthe two curves almost overlapped. H2 O and NH3 dissociate

TABLE IV. Total energy and correlation energy in 共%兲 Ne, CO, N2 , LiF, C2 , and CH4 at equilibrium geometry.

a

SystemState

MOaAct. ElebDMVT( PQ)DMVT( PQG)

Full-CI

Hartree–Fock

Ne
CO
N2
LiF
C2
CH4

1

8共9兲
8共10兲
8共10兲
8共10兲
8共10兲
8共10兲

8
12
12
10
10
8

⫺129.2430(705)⫺113.1163(584)⫺109.4466(571)⫺106.7727(568)⫺77.3387(802)⫺40.4335(403)⫺128.6292(105)⫺112.4544(108)⫺108.7123(108)⫺106.4448(102)⫺75.4793(117)⫺40.2100(124)⫺128.6245(100)⫺112.4426(100)⫺108.7002(100)⫺106.4435(100)⫺75.4340(100)⫺40.1905(100)⫺128.5224(0)⫺112.3033(0)⫺108.5418(0)⫺106.3731(0)⫺75.1626(0)⫺40.1102(0)

S
⌺
1 ⫹
⌺g
1
⌺
1 ⫹
⌺g
1
A1
1

Number of active MOs, with the number of total MOs in parentheses.
Number of active electrons.

b

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

5436

J. Chem. Phys., Vol. 116, No. 13, 1 April 2002

Nakata, Ehara, and Nakatsuji

FIG. 2. Potential curve for the double dissociation of H2 O.

FIG. 4. Potential curve for the triple dissociation of BH3 .

into O( 3 P)⫹H( 1 S)⫹H( 1 S) and N( 4 S)⫹H( 1 S)⫹H( 1 S)
⫹H( 1 S), respectively, and at the dissociation limit, the electronic state becomes a multiconfigurational state.
DMVT( PQG) accurately described these multiconfigurational states, namely, static electron correlations: the deviations from the full-CI were less than 0.5 mhartree in the
dissociation limit. On the other hand, while DMVT( PQ)
reproduced the curves in the short bond region, it failed at
the large internuclear distances. DMVT( PQ) curve did notbound. For BH3 , even the DMVT( PQG) curve slightly deviates from the full-CI curve for R B–H⬎2.0 Å. The dissociation limit of BH3 is heavily quasidegenerate: the electronic
state is represented by several configurations including quadruple excitations. DMVT( PQ) curve for BH3 has a humpat around 1.5 Å, and, the potential curve is repulsive in nature.
Since DMVT( PQG) calculations gave accurate potential curves, their spectroscopic constants were also accurate.
For these systems, the deviations from the full-CI were
within 0.003 Å and 30 cm⫺1 , for r e and ␻ e , respectively.
The dissociation energies (D e ) were estimated slightly larger
by 0.06, 0.27, and 0.26 eV, for H2 O, NH3 and BH3 , respectively. This is because DMVT( PQG) calculations overshootthe full-CI energy around the equilibrium geometries rather
than deviate in the dissociation limit. Since the dissociationsof these systems are homolytic and include multiple bonds,
the Hartree–Fock description of the dissociation limit was ofcourse very crude.
We examined the Weinhold–Wilson inequalities for
BH3 , since the deviation from the full-CI was large for thismolecule. Actually, the violations of the conditions VI and
VII at the equilibrium distance ranged ⫺0.000 019–
⫺0.003 018 and ⫺0.000 391– ⫺0.002 608, respectively, andthose of the conditions VI and VII were ⫺0.004 804–
⫺0.013 975, and ⫺0.001 448– ⫺0.001 956 at the dissociation limit (R⫽5 Å).
D. Potential curves of CO, C2 , N2 , and Be2

Next, we apply the DMVT to the potential energy curvesof CO, C2 , N2 , and Be2 , since their electronic states are
TABLE V. Spectroscopic constants of H2 O, NH3 , BH3 , C2 , N2 , and CO.
Methodr e (Å)␻ e (cm⫺1 )D e (eV)H2 O

Hartree–FockFull-CI
DMVT( PQG)

1.824
1.895
1.894

3952
3253
3276

18.471
6.162
6.227

NH3

Hartree–FockFull-CI
DMVT( PQG)

1.025
1.057
1.057

3750
3324
3291

33.00810.68610.956BH3Hartree–FockFull-CI
DMVT( PQG)

1.154
1.178
1.181

3115
2883
2854

31.28414.28014.537

C2

Hartree–FockFull-CI
DMVT( PQG)

1.233
1.257
1.299

2207
2035
1679

16.876
6.790
7.212

N2

Hartree–FockFull-CI
DMVT( PQG)

1.129
1.210
1.199

2715
2061
1980

31.211
6.220
6.622

CO

Hartree–FockFull-CI
DMVT( PQG)

1.146
1.193
1.201

2461
2063
1990

12.692
9.328
9.540

SystemFIG. 3. Potential curve for the triple dissociation of NH3 .

Downloaded 16 Apr 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. 116, No. 13, 1 April 2002

Density matrix variational theory

5437

FIG. 5. Potential curve of CO.

FIG. 7. Potential curve of N2 .

very characteristic. In the potential curve of CO, theHartree–Fock configuration is dominant at around the equilibrium distance, but its weight decreases as the distance increases and finally becomes zero at the dissociation limit. C2
has unoccupied p ␴ MO, therefore, the ground state is always
quasidegenerate even in the equilibrium geometry. N2 includes triple-bond dissociation, therefore its potential curveis highly quasidegenerate at large internuclear distance. Be2
has no bonding interaction. Potential curves of these molecules were calculated for R⫽0.5– 5.0 Å. Minimal STO-6G
basis set was used and the 1s orbitals were kept as frozen.
The potential curves were shown in Figs. 5– 8 and the spectroscopic constants were given in Table V.
As in other systems, DMVT( PQG) curves almost overlapped with the full-CI curves, while DMVT( PQ) curves
were calculated as repulsive. The deviations of
DMVT( PQG) from the full-CI increases in the order of N2 ,CO, and C2 . Though it is true that the description of thequasidegeneracy of C2 is difficult, there is another factor in
the accuracy. Since we used minimal basis set, the calculations of N2 and CO were for 16 spin orbitals with 12 electrons, namely, 4 hole spin orbitals, while those of C2 are for
16 spin orbitals with 10 electrons; 6 hole spin orbitals. We

think this also affected the accuracy of the results. The potential curve of Be2 was repulsive, since van der Waals interaction was not described by the present basis set. The
DMVT( PQ) gave better description than other systems.
The DMVT( PQG) results for the spectroscopic constants of these diatomic molecules were less accurate thanthose for H2 O and NH3 . The deviations were ⬃0.01 Å and⬃80 cm⫺1 for r e and ␻ e , respectively, for CO and N2 . For
C2 , the errors were as large as 0.04 Å and 350 cm⫺1 . These
results reflect the quality of the DMVT around the equilibrium geometry.
We also calculated the Weinhold–Wilson inequalities forC2 and CO. As expected, large violations occurred for theinequalities VI and VII. For C2 , the violations were calculated as ⫺0.001 810– ⫺0.027 667 and ⫺0.000 900–
⫺0.011 494 for conditions VI and VII, respectively, at R
⫽1.5 Å, and ⫺0.003 047– ⫺0.002 500 only for condition
VI at R⫽5.0 Å: the violations at R⫽1.5 Å were larger thanthose at R⫽5.0 Å. The ground state of C2 is quasidegenerate even at the equilibrium geometry and this is the reason of
the crude spectroscopic constants for C2 by the
DMVT( PQG).
For CO, the violations ranged⫺0.008 185– ⫺0.008 434 for condition VI at the dissociation

FIG. 6. Potential curve of C2 .

FIG. 8. Potential curve of Be2 .

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

5438

J. Chem. Phys., Vol. 116, No. 13, 1 April 2002

Nakata, Ehara, and Nakatsuji

TABLE VI. Examination of size consistency for Ne and H2 O, NH3 , BH3 ,Be2 , CO, and C2 .
SystemMethodE molaE atomb⌬E

H2 O

DMVT( PQG)
Full-CI

⫺75.4589⫺75.4588⫺75.4589⫺75.4588

0.00000.0000NH3DMVT( PQG)
Full-CI

⫺55.6622⫺55.6622⫺55.6622⫺55.6622

0.00000.0000BH3DMVT( PQG)
Full-CI

⫺25.8680⫺25.8482⫺25.8482⫺25.8482

0.01980.0000Be2DMVT( PQG)
Full-CI

⫺29.1655⫺29.1654⫺29.1654⫺29.1654

0.00010.0000CODMVT( PQG)
Full-CI⫺112.1153⫺112.1095⫺112.1095⫺112.1095

0.00580.0000C2DMVT( PQG)
Full-CI

⫺75.2187⫺75.1854⫺75.1854⫺75.1854

0.03330.0000N2DMVT( PQG)
Full-CI⫺108.4982⫺108.4982⫺108.4982⫺108.4982

0.00000.0000aEnergy of molecule at the dissociation limit.
Sum of the energies of the isolated atoms.

blimit and ⫺0.000 546 for condition VI at R⫽1.3 Å. There
were no errors of conditions VII and VIII for CO.
E. Size consistency

As we discussed in Sec II C, the G condition is very
important for the size-consistent property of the method.
Here, we examine the size consistency of the results. In
Table VI, the total energy of the molecule in the dissociation
limit and the sum of the total energies of the isolated atoms
are compared for H2 O, NH3 , BH3 , N2 , C2 , CO, and Be2 .
For H2 O, NH3 , and N2 , the total energies calculated by
DMVT( PQG) agree within numerical accuracy, whichshows the size consistency holds for these systems. In thesesystems, Weinhold–Wilson inequalities were also satisfied
and the calculations were quite accurate. As seen from thepotential curves, DMVT( PQ) calculations gave miserable
results from the standpoint of the size consistency. The G
condition is apparently indispensable for the size-consistent
property. For other systems, BH3 , C2 , CO, and Be2 , the size
consistency of DMVT( PQG) was not satisfactory. This isbecause the G condition is not a sufficient condition for thesize consistency. Note that the Weinhold–Wilson inequalities
VI and VII were not satisfied for these systems.

F. Artificially correlation enhanced system

It is interesting to see the performance of the present
method for the strongly correlated system. We here introduced the model Hamiltonian in which the electron correlations are controlled by a parameter. The Hamiltonian is partitioned into F, Fock operator and the rest, V:
共3.1兲

H⫽F⫹␭V,

where ␭ is a real parameter that controls the strength of theelectron correlations and ␭⫽1 corresponds to the original
Hamiltonian. We adopted Be and H2 O and changed ␭ form
0.1– 10 000 and the results were shown in Tables VII and
VIII, respectively.
For Be system, the DMVT( PQG) reproduced the exact
correlation energy quite accurately and the deviations wereeven tempered for the variation of ␭ within 1%. On the other
hand, the DMVT( PQ) gave random errors for the variationof ␭. For H2 O, the errors of the DMVT( PQG) became
large, but, were within 15% relative to the total electron correlations. The DMVT( PQG) calculations converged evenfor the heavily correlated systems (␭⫽10 000), though theabsolute errors were not small.
We did not see the Weinhold–Wilson violations in the
Be system for all ␭. The violations for H2 O were not sosimple. For ␭⫽1.0, the violations of the condition VI occurred as ⫺0.000 54– ⫺0.003 68, but no violations occurredfor ␭⫽2.0, which has the largest correlation energy error inpercent. For ␭⫽10 000, we got large violations as⫺0.008 71– ⫺0.001 036.

IV. CONCLUSION

The DMVT was applied to the calculations of the potential energy surfaces of the atoms and small molecules, Ne,
H4 , H2 O, NH3 , BH3 , CO, N2 , C2 , and Be2 . This is the firststudy in which the bond dissociation was properly describedby the DMVT. In the previous DET study of potential
curves,26 the results were good up to R⬃2R e , but at large
distances, the calculations failed to converge. Generally, theDMVT( PQG) calculation reproduced the full-CI curves
very accurately and they sometimes overlapped even in the
dissociation limit, though the potential curves for BH3 andC2 were less accurate than others. The quasidegenerate states
were well described by the DMVT( PQG) calculations. On

TABLE VII. Total energy for the model Hamiltonian, H⫽F⫹␭V of Be and the correlation energy in 共%兲.
␭

DMVT( PQ)DMVT( PQG)

Full-CI

Hartree–Fock

0.10
0.50
1.00
2.00
3.00
4.00
5.00
10.00
10 000.0⫺10.5324(249)⫺12.3365(226)⫺14.6064(200)⫺19.2097(151)⫺24.1197(115)⫺29.6662(139)⫺37.1084(212)⫺76.5823(165)⫺81 838.44(102)⫺10.5322(100)⫺12.3317(100)⫺14.5895(100)⫺19.1600(100)⫺24.0491(100)⫺29.2180(101)⫺34.7026(100)⫺68.0991(100)⫺81 294.14(100)⫺10.5322(100)⫺12.3317(100)⫺14.5895(100)⫺19.1596(100)⫺24.0469(100)⫺29.2115(100)⫺34.6922(100)⫺68.0414(100)⫺81 290.74(100)⫺10.5321(0)⫺12.3278(0)⫺14.5725(0)⫺19.0619(0)⫺23.5513(0)⫺28.0407(0)⫺32.5301(0)⫺54.9771(0)⫺44 904.03(0)Downloaded 16 Apr 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. 116, No. 13, 1 April 2002

Density matrix variational theory

5439

TABLE VIII. Total energy for the model Hamiltonian, H⫽F⫹␭V of H2 O and the correlation energy in 共%兲.
␭
0.10
0.50
1.00
1.50
2.00
3.00
4.00
5.00
10.00
10 000.0

DMVT( PQ)DMVT( PQG)

Full-CI

Hartree–Fock⫺41.1669(234)⫺56.5293(233)⫺75.7953(232)⫺95.1938(235)⫺115.0471(222)⫺155.4737(137)⫺196.0891(124)⫺236.8429(120)⫺441.1663(116)⫺417 093.18(124)⫺41.1664(100)⫺56.5158(102)⫺75.7310(104)⫺95.0064(106)⫺114.5560(115)⫺154.6553(101)⫺195.0768(101)⫺235.6115(101)⫺438.6811(102)⫺411 089.69(102)⫺41.1664(100)⫺56.5155(100)⫺75.7290(100)⫺94.9978(100)⫺114.4863(100)⫺154.6348(100)⫺195.0347(100)⫺235.5285(100)⫺438.3447(100)⫺410 645.30(100)⫺41.1661(0)⫺56.5051(0)⫺75.6789(0)⫺94.8526(0)⫺114.0264(0)⫺152.3740(0)⫺190.7215(0)⫺229.0690(0)⫺420.8068(0)⫺383 512.82(0)the other hand, the curves by DMVT( PQ) were always repulsive, which showed the potential importance of the G
condition.
We examined the size consistency of the present method.
The G condition is found to be related to the size consistencyof the method and shown to be essential to the behavior of
the potential curves of DMVT( PQG), especially in the dissociation limit.
We also examined the Weinhold–Wilson inequalities forthe 2-RDM of DMVT( PQG) calculations where the results
were less accurate, and found that the inequalities VI and VII
were violated. We think these inequalities may be new candidates for the N-representability condition of the DMVT.
Since these are linear conditions, it would be easily includedin the conditions of the DMVT relaxed with the SDP. Suchstudy is now in progress.
ACKNOWLEDGMENTS

This study has been supported by the Special Grant for
Scientific Research from the Japanese Ministry of Education,
Science, Culture, and Sports. The authors give warm thanks
to Mr. Mitsuhiro Fukuda, Dr. Katsuki Fujisawa, Dr. Kazuhide Nakata, and Professor Masakazu Kojima, for useful discussions and encouragements.
K. Husimi, Proc. Phys. Math. Soc. Jpn. 22, 264 共1940兲.
P.-O. Löwdin, Phys. Rev. 99, 1474 共1955兲.
3
H. Nakatsuji, Phys. Rev. A 14, 41 共1976兲.
4
C. Valdemoro, Phys. Rev. A 45, 4462 共1992兲; F. Colmenero, C. Pérez delValle, and C. Valdemoro, ibid. 47, 971 共1993兲; F. Colmenero and C.
Valdemoro, ibid. 47, 979 共1993兲.
5
H. Nakatsuji and K. Yasuda, Phys. Rev. Lett. 76, 1039 共1996兲; K. Yasuda
and H. Nakatsuji, Phys. Rev. A 56, 2648 共1997兲.
6
D. A. Mazziotti, Phys. Rev. A 57, 4219 共1998兲; Chem. Phys. Lett. 289,
419 共1998兲.
7
H. Nakatsuji, in Many-Electron Densities and Reduced Density Matrices,
edited by J. Cioslowski 共Kluwer Academic, New York, 2000兲.
8A. J. Coleman, Rev. Mod. Phys. 35, 668 共1963兲.
9
C. Garrod, and J. K. Percus, J. Math. Phys. 5, 1756 共1964兲.
10
L. J. Kijewski and J. K. Percus, Phys. Rev. A 2, 1659 共1970兲; L. J.
Kijewski, ibid. 6, 1659 共1972兲; L. J. Kijewski, ibid. 9, 2263 共1974兲.
11
C. Garrod, M. V. Mihailović, and M. Rosina, J. Math. Phys. 16, 868
共1975兲.
12
C. Garrod and M. A. Fusco, Int. J. Quantum Chem. 10, 495 共1976兲.
1
2

R. M. Erdahl, Rep. Math. Phys. 15, 147 共1979兲.
R. M. Erdahl, and B. Jin, J. Mol. Struct.: THEOCHEM 527, 207 共2000兲;
R. M. Erdahl and B. Jin, in Many-Electron Densities and Reduced Density
Matrices, edited by J. Cioslowski 共Kluwer Academic, New York, 2000兲.
15
M. Nakata, H. Nakatsuji, M. Ehara, M. Fukuda, K. Nakata, and K.
Fujisawa, J. Chem. Phys. 114, 8282 共2001兲.
16
Y. Nesterov and A. S. Nemirovskii, Interior Point Polynomial Method inConvex Programming: Theory and Applications 共SIAM, Philadelphia,
1993兲.
17
L. Vandenberghe and S. Boyd, SIAM Rev. 38, 4 9 共1996兲.
18
M. Kojima, Semidefinite Programming and Interior-Point Methods,
http://www.is.titech.ac.jp/˜ kojima/wabun.html, 1996 共in Japanese兲.
19
K. Fujisawa, M. Kojima, and K. Nakata, SDPA 共SemiDefinite Programming Algorithm兲 User’s Manual Version 5.00, August 1999.
http://is-mj.archi.kyoto-u.ac.jp/˜ fujisawa/software.html
20
D. A. Mazziotti and R. M. Erdahl, Phys. Rev. A 63, 042113 共2001兲.
21
C. Valdemoro, L. M. Tel, and E. Perez-Romero, Phys. Rev. A 61, 032507
共2000兲.
22
H. Nakatsuji, J. Chem. Phys. 113, 2949 共2000兲.
23
H. Nakatsuji and E. R. Davidson, J. Chem. Phys. 115, 2000 共2001兲.
24
H. Nakatsuji, J. Chem. Phys. 115, 2465 共2001兲.
25
H. Nakatsuji J. Chem. Phys. 116, 1811 共2002兲.
26
M. Ehara, M. Nakata, H. Kou, K. Yasuda, and H. Nakatsuji, Chem. Phys.
Lett. 305, 483 共1999兲.
27
F. Weinhold and E. B. Wilson, Jr., J. Chem. Phys. 47, 2298 共1967兲.
28
E. R. Davidson, J. Math. Phys. 10, 725 共1969兲.
29
W. B. McRae and E. R. Davidson, J. Math. Phys. 13, 1527 共1972兲.
30
J.-P. Hansen and I. R. McDonald, Theory of Simple Liquids 共Academic,
Harcourt Brace, New York, 1990兲.
31
T. V. Voorhis and M. Head-Gordon, J. Chem. Phys. 113, 8873 共2001兲.
32
J. Paldus, P. Piecuch, L. Pylypow, and B. Jezirovski, Phys. Rev. A 47,
2738 共1993兲.
33
S. Huzinaga, J. Chem. Phys. 42, 1293 共1970兲; T. H. Dunning, Jr., ibid. 53,
2823 共1970兲; T. H. Dunning, J. and P. J. Hay, In Method of Electronic
Structure Theory, edited by H. F. Schaefer III 共Plenum, New York, 1977兲,Vol. 2.
34
Basis sets were obtained from the Extensible Computational Chemistry
Environment Basis Set Database, Version 4/22/01, as developed and distributed by the Molecular Science Computing Facility, Environmental
and Molecular Sciences Laboratory, the Pacific Northwest Laboratory,
P.O. Box 999, Richland, WA 99352. http://www.emsl.pnl.gov:2080/forms/
basisform.html
35
W. J. Hehre, R. F. Stewart, and J. A. Pople, J. Chem. Phys. 51, 2657
共1969兲.
36
J. H. Callomon, E. Horita, K. Kuchitsu, W. J. Lafferty, A. G. Maki, and C.
S. Pote, Landolt–Börnstein 共Springer-Verlag, Berlin, 1976兲.
37
K. P. Huber and G. Herzberg, Molecular Spectra and Molecular Structure
IV, Electronic Constants of Diatomic Molecules 共Van Nostrand Reinhold,
New York, 1979兲.
38
J. L. Dunhum, Phys. Rev. 41, 713, 721 共1932兲; H. M. Hulburt and J. O.
Hirschfelder, J. Chem. Phys. 9, 61 共1941兲.
13
14

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