Populations of Carbonic Acid Isomers at 210 K from a Fast Two

Oct 6, 2011 - DOI: 10.1021/jp5037469. Sourav Ghoshal and Montu K. Hazra . Autocatalytic Isomerizations of the Two Most Stable Conformers of Carbonic ...
0 downloads 0 Views 1MB Size
Download PDF
ARTICLE pubs.acs.org/JPCA

Populations of Carbonic Acid Isomers at 210 K from a Fast Two-Electron Reduced-Density Matrix Theory Christine A. Schwerdtfeger and David A. Mazziotti* Department of Chemistry and The James Franck Institute, The University of Chicago, Chicago, Illinois 60637, United States ABSTRACT: Parametrization of the 2-electron reduced density matrix (2-RDM) rather than the many-electron wave function yields a new family of electronic-structure methods that are faster and more accurate than traditional coupled electron-pair methods including coupled cluster with single and double excitations. Deriving the parametrization from N-representability conditions generates a 2-RDM that captures significant correlation from triple and higherorder excitations at the cost of double excitations. We apply the parametric 2-RDM method to confirm recent experiments determining the relative thermodynamic populations of the ciscis and cistrans isomers of carbonic acid. In 2010 Bernard et al. showed by infrared spectroscopy that the populations of ciscis and cistrans isomers have a 10:1 ratio at 210 K. By use of the parametric 2-RDM method, we predict a 8:1 ratio at 210 K. Comparable ab initio methods overestimate the stability of the ciscis isomer with 24:1 and 21:1 ratios. These 2-RDM-based methods promise to have significant applications throughout chemistry.

I. INTRODUCTION Calculations of many-electron atoms and molecules generally employ either the many-electron wave function or the oneelectron density as the primary variable. These two classes of methods are known as ab initio wave function methods1 and density functional theory,2 respectively. However, a third class of methods can be defined in which the two-electron reduced density matrix (2-RDM) is the primary variable.3,4 The 2-RDM 2 D can be obtained from the density matrix of a N-electron system by integrating over all electrons except two Z 2 Dð12; 12Þ ¼ Ψð12:::NÞΨð12:::NÞd3:::dN ð1Þ where each number indicates the spatial and spin variables of an electron. Because electrons are indistinguishable with pairwise Coulomb interactions, the energy of a N-electron atom or molecule can be expressed as a linear functional of the 2-RDM without any additional information from the many-electron wave function. As early as fifty years ago, Mayer,5 L€owdin,6 Coleman,7 and others realized that the variational calculation of a molecule’s electronic energy as a functional of the 2-RDM requires significant constraints on the 2-RDM to ensure that it corresponds to a many-electron wave function. The need for such constraints (N-representability conditions)3,4,714 prevented the development of 2-RDM-based electronic structure methods for a long time even after the appearance of methods based on the one-electron density. Several complementary approaches to the direct calculation of the 2-RDM have emerged recently including the parametric variational 2-RDM methods.1521 In these methods the 2-RDM is parametrized in terms of single and double excitations to remain very nearly N-representable. The expression for the r 2011 American Chemical Society

energy is minimized as a functional of the 2-RDM to generate the ground-state energy without additional reference to the many-electron wave function. This 2-RDM-based method is faster and more accurate than comparable electron-pair-based wave function methods including coupled cluster with single and double excitations (CCSD).2224 While the 2-RDM does not explicitly depend on triple excitations, it has an accuracy approaching that of coupled cluster with single, double, and triple excitations (CCSDT).25 Precursors of the parametric 2-RDM methods include Weinhold and Wilson’s parametrization of the 2-RDM26,27 with a double-excitation wave function and the coupled electronpair approximation (CEPA) methods.28,29 Unlike the CEPA methods, however, the parametric 2-RDM methods employ N-representability conditions to guide the parametrization, which results in energy functionals that are both more stable and more accurate. The first calculations with these methods used a parametrization derived by Kollmar15 that yields similar accuracy as CCSD except in the presence of mild multireference correlation where it improves upon CCSD. More recently, a new energy functional, introduced by one of the authors,20,21 yields energies with an accuracy approaching CCSDT at a cost like configuration interaction with single and double excitations. This new parametrization has been applied successfully to a set of test molecules in large basis sets as well as the study of oxywater and its elusiveness.30 In 2010 Bernard et al.31 announced the measurement by infrared spectroscopy of two gas-phase isomers of carbonic acid, the ciscis (CC) and the cistrans (CT) as well as the ciscis dimer. They furthermore determined that the relative thermodynamic populations of CC, CT, and dimer at 210 K are in the Received: June 20, 2011 Revised: August 25, 2011 Published: October 06, 2011 12011

dx.doi.org/10.1021/jp2057805 | J. Phys. Chem. A 2011, 115, 12011–12016

The Journal of Physical Chemistry A

ARTICLE

ratio 10:1:1. In 2009 Mori et al.32 had discovered the gas-phase CT isomer, but the CC and dimer had previously been undetected. (In 2011 Mori et al.33 also detected the CC isomer by microwave spectroscopy.) Gas-phase carbonic acid may be important in the astrochemistry of comets, planets like Mars, and interstellar space as well as in the Earth’s upper atmosphere,34 and biologically, its aqueous form is an intermediate in the removal of carbon dioxide during respiration.35 Consequently, carbonic acid has been the focus of both experimental3136 and theoretical3742 studies. In this paper we confirm the thermodynamic populations measured by Bernard et al. by computing equilibrium constants at 210 K with electronic energies for the CC and CT isomers from the parametric 2-RDM method. The 2-RDM-based method predicts the CC and CT populations to be in a ratio of 8:1. The data from the 2-RDM methods are compared with ab initio wave function methods including CCSD and CCSD(T), which predict ratios considerably larger than experiment. Section III presents the computational results including absolute and relative energies, relative populations at thermodynamic equilibrium, errors in molecular geometries, natural occupation numbers, and computational timings. Conclusions and salient details of the parametric 2-RDM method with the improved energy functional are given in sections IV and II, respectively.

II. THEORY The energy of any N-electron atom or molecule can be expressed as a functional of the 1- and 2-RDMs as E¼

∑pq 1 K pq 1Dpq

þ

∑ 2 V pqst 2 Dpqst

ð2Þ

pqst

where the matrices 1K and 2V contain the one- and two-electron integrals, 1D and 2D are the 1- and 2-RDMs, and the indices denote the members of a finite set of spin orbitals. The energy cannot be directly minimized as a functional of the 1- and 2-RDMs, however, because the RDMs must be constrained by N-representability conditions to represent N electrons. In the first step of deriving the parametric 2-RDM method21 we parametrize the 1- and 2-RDMs to obtain an energy that is correct through third order of perturbation theory about a meanfield reference. For notational clarity we assume that the orbitals are chosen to be the eigenfunctions of the optimal 1-RDM, known as the natural orbitals, which allows us to neglect the explicit inclusion of single excitations in the energy functional. Practically, the assumption of natural orbitals can be easily relaxed by either: (i) optimizing the orbitals to minimize the energy or (ii) introducing parameters from the second-order part of the 1-RDM, which serve as single excitations, into the functional. In the natural-orbital basis set all second-order parts of the 1- and 2-RDMs can be constructed from the first-order part of the 2-RDM 2T 2

D ¼ T þ Oðλ Þ 2

2

ð3Þ

This third-order parametrization of the energy in fact corresponds to the simplest coupled electron-pair approximation.28 Although the energy expression is size extensive since 2T can be defined to scale linearly with system size, the parametrized 2-RDM is not N-representable, and often the minimum energy is a lower bound to the energy from solving the Schr€odinger equation in the finite orbital basis set, known as full configuration interaction.

This 2-RDM functional for the energy through third order of perturbation theory, however, is not unique since a variety of fourth-order terms can be added. The flexibility in the functional can be understood through the following parametrization of a class of the 2-RDM elements, corresponding to double excitations sffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi 1 2 2 ab 1 ð4Þ Dij ¼ 2 T ab f abcd j2 T cd ij kl j 4 klcd ijkl

∑

where i, j, k, and l correspond to occupied orbitals, a, b, c, and d correspond to unoccupied orbitals, and the eight-index tensor 15,16,29 Selectf abcd ijkl is sometimes known as the topological factor. abcd ing the elements of f ijkl to be zero yields the coupled electronpair approximation (CEPA0), while selecting the elements to be one yields the configuration interaction method with double excitations (CID). Hence, the topological factor interpolates between CEPA0 and CID with the latter case being the parametrization of the 2-RDM examined by Weinhold and Wilson in 1967.26,27 While the CEPA0 2-RDM is size extensive, it is generally not N-representable, and while the CID 2-RDM is N-representable, it is generally not size extensive, at least not for N > 2. We need a parametrization of the 2-RDM that balances the best of both CEPA0 and CID. Traditional approaches to improving CEPA0 rely on perturbation theory, but in the parametric 2-RDM method the parametrization of the 2-RDM is improved through N-representability conditions. An important set of N-representability constraints, known as 2-positivity conditions,8,1014 corresponds to restricting three matrix forms of the 2-RDM to be positive semidefinite 2

Dg0

ð5Þ

2

Q g0

ð6Þ

2

Gg0

ð7Þ

A matrix is positive semidefinite if and only if it has non-negative eigenvalues. The positive semidefiniteness of 2D—the conventional form of the 2-RDM—corresponds to keeping the probability distribution for finding two particles non-negative while the semidefiniteness of 2Q and 2G corresponds to keeping the probability distributions for finding two holes nonnegative and one particle/one hole nonnegative, respectively. Each of these three matrices can be computed from the others by linear mappings.1013 The 2-positivity conditions imply N-representability conditions known as the CauchySchwarz inequalities such as ij

2 2 2 ab ð2 Dab ij Þ e Dij Dab ij

2 2 2 ab ð2 Q ab ij Þ e Q ij Q ab

ð8Þ ð9Þ

Taking an average of these two inequalities, as described in ref.,21 we obtain the values of the topological factor corresponding to the parametric 2-RDM method derived by Mazziotti (M). Similarly, if we perform this average with four CauchySchwarz inequalities from the 2G g 0 condition, we obtain the topological factor first derived by Kollmar (K) by a related but different argument. Computations show that both the M and K energy functionals are size extensive and nearly N-representable at both equilibrium and nonequilibrium geometries. 12012

dx.doi.org/10.1021/jp2057805 |J. Phys. Chem. A 2011, 115, 12011–12016

The Journal of Physical Chemistry A

ARTICLE

Because the M parametrization of the 2-RDM is derived from N-representability conditions rather than perturbation theory, it contains terms that do not appear in the coupled electron-pair approximations and are not expressible by the topological factor f abcd ijkl in eq 4 without approximation. Importantly, these terms in the M functional separate the family of parametric 2-RDM methods from the coupled electron-pair approximations in terms of both computational stability and accuracy. Practically, the M functional is implementable without approximation as described in eq 62 of a previous study.21 This unabridged implementation scales computationally as r4 with r being the number of orbitals, which causes the M parametrization to remain at the same computational cost as configuration interaction with single and double excitations.

III. RESULTS Relative electronic energies at 0 K with zero-point vibrational energies and relative thermodynamic populations at 210 K are shown for CC, CT, and their transition state (TS) in Figure 1 and Table I. Computations of the electronic energies of the CC and CT isomers of carbonic acid and their transition state (TS) were performed by the parametric 2-RDM method with the M energy functional, implemented in a Fortran 2003 code developed by the authors, and the wave function methods CCSD and CCSD with a perturbative treatment of triple excitations [CCSD(T)],2224,43,44 implemented in the GAMESS electronic structure package.45 For all electronic structure methods we employed the correlation consistent polarized valence double-ζ (cc-pVDZ), triple-ζ (cc-pVTZ), and quadruple-ζ (cc-pVQZ) basis sets46 with extrapolation to the complete basis-set limit (BL).47,48 Details of the extrapolation are given in ref 30. Molecular geometries from all methods were optimized except in the cc-pVQZ basis set where geometries were taken from the cc-pVTZ basis set. Thermodynamic equilibrium constants were calculated with vibrational energies from the normal-modes approximation and rotational energies from the standard classical approximation.49 In the basis-set limit (BL) the parametric 2-RDM method predicts the CC isomer to be 1.2 kcal/mol lower than the CT isomer while the CCSD and CCSD(T) methods predict 1.6 and 1.6 kcal/mol isomerization energies, respectively. For all three methods the isomerization energies decrease as a function of basis-set size. While the variations in energies with method are small, these variations significantly affect the predicted relative thermodynamic populations at 210 K. The parametric 2-RDM method predicts an 8:1 ratio for the relative population of CC to CT, which is in close agreement with the 10:1 ratio determined from the 2010 spectroscopic studies of Bernard et al.31 In contrast, both CCSD and CCSD(T) predict the CT isomer to be much scarcer than the CC isomer with relative population ratios of 24:1 and 21:1, respectively. The inclusion of perturbative triples does not significantly affect either the isomerization energy or the relative population. The activation energy from all methods is computed to be in the range of 9.0 to 9.5 kcal/mol, which is sufficiently small for isomerization to occur at 210 K. Figure 1 summarizes the relative energies at 0 K from the 2-RDM method in the basis-set limit. Mean-field (HartreeFock) electronic energies and correlation energies at 0 K are shown for CC, CT, and their TS in Table II. At all basis sets including the BL the correlation energies from the parametric 2-RDM method with the M functional are

Figure 1. The potential energy curve (kcal/mol), connecting the CC and CT through their TS, is shown from the parametric 2-RDM method in the extrapolated basis set limit (BL). The isomerization energy from CC to CT is 1.2 kcal/mol, which after conversion to free energy predicts a 8:1 ratio of populations in agreement with the 10:1 ratio from Bernard et al.’s experiment.31 Illustration by Kasra Naftchi-Ardebili, Mazziotti Group, The University of Chicago, 2011. Used with permission.

much closer to those from CCSD(T) than CCSD. The 2-RDM parametrization does not explicitly contain triple excitations, and yet its parametrization by single and double excitations implicitly accounts for many of the triple excitations’ degrees of freedom without the cost of their explicit inclusion. Both the CCSD and 2-RDM calculations yield 1-RDMs whose eigenvalues are the occupation numbers of the natural orbitals (eigenfunctions of the 1-RDM). Comparison of the naturalorbital occupation numbers from the CCSD and 2-RDM methods shows that the 2-RDM method captures greater correlation of the orbitals. Occupation numbers are not available within GAMESS from CCSD(T), which simply computes a perturbative energy correlation for CCSD. For CC, CT, and TS the Nth natural occupation numbers are 0.974, 0.974, and 0.972 from the CCSD method and 0.969, 0.968, and 0.966 from the 2-RDM method. The occupation numbers from the 2-RDM method are lower than those from CCSD, reflecting a small increase in correlation. Both methods show the TS to be slightly more correlated than the isomers of carbonic acid. The optimized molecular geometries from the 2-RDM and CCSD methods in the cc-pVTZ basis set are assessed relative to the geometries from CCSD(T) in Table III. For each molecule we compute the root-mean-squared deviations of the 2-RDM and CCSD geometries relative to the geometry of CCSD(T) with the geometry deviations being reported separately for bond lengths and angles in Ångstroms and degrees, respectively. Results indicate that the molecular geometries of the parametric 2-RDM method with the M functional are consistently closer than the geometries of CCSD to those of CCSD(T). The largest deviation of 2-RDM from CCSD(T) is seen in the bond angle of the CC isomer where the root-mean-squared deviation of CCSD is about four times that of the 2-RDM method. Computational timings for energy calculations at single geometries (single-point calculations) from the parametric 2-RDM, CCSD, and CCSD(T) methods are reported in Table IV as 12013

dx.doi.org/10.1021/jp2057805 |J. Phys. Chem. A 2011, 115, 12011–12016

The Journal of Physical Chemistry A

ARTICLE

Table I. Reaction Energies and Barriers (ΔE) with Vibrational Zero-Point Energies (kcal/mol), Relative to the CC Isomer, Are Given in the cc-pVDZ (DZ), cc-pVTZ (TZ), and cc-pVQZ (QZ) Basis Sets as Well as the Extrapolated BLa CT (relative to CC) ΔE

TS (relative to CC) ΔG°

ΔE

population BL

ΔG°

method

DZ

TZ

QZ

BL

BL

DZ

TZ

QZ

BL

BL

CCSD

2.00

1.64

1.59

1.59

1.33

1:24

9.92

9.45

9.45

9.49

9.20

CCSD(T) M

1.97 1.45

1.58 1.13

1.52 1.13

1.56 1.19

1.26 0.89

1:21 1:8

9.91 9.47

9.41 9.35

9.40 9.39

9.44 9.47

9.16 9.20

a

The free energy G° (kcal/mol) of the CT isomer and the TS and the population of the CT isomer, relative to the CC isomer, are shown in the basis set limit.

Table II. Correlation Energies [Hartrees (H)] of the CC, the TS, CT Are Presented from the M Parametric 2-RDM, CCSD, and CCSD(T) Methods in the cc-pVDZ (DZ), cc-pVTZ (TZ), and cc-pVQZ (QZ) Basis Setsa energy (H) molecule basis set CC

TS

CT

HF

Table IV. Computational (Wall) Timings for Energy Calculations at Single Geometries (Single-Point Calculations) from the Parametric 2-RDM, CCSD, and CCSD(T) Methods Are Reported as Functions of Basis-Set Sizea

correlation energy (H) CCSD

2-RDM

relative timings

CCSD(T)

molecule

basis set

CCSD

2-RDM

CCSD(T)

CC

DZ

1.0

0.2

1.2

TZ QZ

24 1100

4.3 92

36 1300

DZ

263.680028 0.686028 0.698601 0.705513

TZ QZ

263.763835 0.857617 0.877897 0.892326 263.784796 0.916791 0.939825 0.955582

BL

263.791788 0.966667 0.992062 1.008622

DZ

1.0

1.0

1.2

DZ

263.662103 0.686480 0.700156 0.705997

TZ

22

38

31

TZ

263.746592 0.858143 0.878952 0.892924

QZ

1100

980

1200

QZ

263.767632 0.917236 0.940738 0.956105

DZ

1.0

0.6

1.3

BL

263.768299 0.973377 0.999173 1.015408

TZ

27

16

38

DZ

263.676101 0.686448 0.699238 0.705958

QZ

1200

270

1400

TZ QZ

263.760311 0.858196 0.878641 0.892994 263.781399 0.917331 0.940447 0.956211

BL

263.788444 0.966721 0.992538 1.009139

TS

CT

a

For each molecules timings are relative to those from CCSD in the cc-pVDZ basis set. All methods use Intel’s Math Kernel Library for optimized basic linear algebra subroutines (BLAS).

a

These energies (BL) are extrapolated to the complete-basis-set limit using a two part extrapolation scheme. The Hartree Fock energy is extrapolated by fitting to an exponential function47 while the correlation energy is fitted to an inverse power function.48

Table III. Root-Mean-Squared (rms) Deviations of the 2-RDM and CCSD Geometries Relative to the Geometry from CCSD(T) Were Computed for Each Moleculea rms deviation of geometric parameters molecule

parameter type

CCSD

2-RDM

CC TS

bond length

0.0029 0.0022

0.0010 0.0006

0.0023

0.0006

CT CC

bond angle

0.2493

0.0622

TS

0.1052

0.0611

CT

0.1692

0.0253

a

These deviations are reported separately for bond lengths and angles in Ångstroms and degrees, respectively.

functions of basis set size. Theoretically, the parametric 2-RDM and CCSD methods scale as kN2(r  N)4 in floating-point operations where the prefactor k is smaller for the 2-RDM method than for CCSD. The perturbative correction for triple

excitations in CCSD(T) has a higher computational scaling. Practically, the timings reported in Table IV reflect differences in implementation, which include (i) the use of Abelian pointgroup symmetry by the parametric 2-RDM method in contrast to the coupled-cluster implementation in GAMESS and (ii) the use of disk storage for the larger coupled-cluster calculations in contrast to the parametric 2-RDM methods which perform integral storage in memory. The implementation of the parametric 2-RDM methods, still preliminary, has also not been fully optimized. With these important caveats, nonetheless, the following conclusions can be made: (i) the current implementation of the parametric 2-RDM method is competitive with the coupled cluster methods in GAMESS; (ii) the 2-RDM calculations in the cc-pVQZ basis set are generally 23 orders of magnitude more expensive than calculations in the cc-pVDZ basis set; (iii) the 2-RDM method is fastest for the molecule with the highest point-group symmetry.

IV. DISCUSSION The idea of computing the 2-RDM directly without the manyelectron wave function was significantly promoted by Charles Coulson at an international conference in Boulder, Colorado in 1959.50 Despite the 2-RDM’s promise of both accuracy and computational efficiency, he and others recognized that significant progress toward practical 2-RDM methods would require 12014

dx.doi.org/10.1021/jp2057805 |J. Phys. Chem. A 2011, 115, 12011–12016

The Journal of Physical Chemistry A the development of nontrivial constraints on the 2-RDM—Nrepresentability conditions—to ensure that it corresponds to an N-electron wave function.10 By parametrizing the 2-RDM through N-representability conditions, the parametric 2-RDM method enables a direct calculation of the 2-RDM that is faster and more accurate than traditional wave function methods based on single and double excitations. In its application to the two isomers of carbonic acid, CC and CT, the 2-RDM method predicts in the limit of an extrapolated basis set that the CT isomer is 1.2 kcal/mol higher in energy than the CC isomer. At 210 K the 2-RDM method predicts the thermodynamic populations of the CC and CT isomers to be in an 8:1 ratio, which agrees with the 10:1 ratio determined experimentally in 2010 by Bernard et al.31 The wave function methods CCSD and CCSD(T) yield higher isomerization energies of 1.6 and 1.6 kcal/mol and relative populations of 24:1 and 21:1, respectively. While the parametric 2-RDM methods have connections to earlier methods like the coupled electron-pair approximation, they also possess significant differences that result in better stability and accuracy. The simplest coupled electron pair approximation and its variations are based solely on arguments from perturbation theory, especially the effect of higher order diagrams on lower order diagrams. In contrast, the parametric 2-RDM methods develop the parametrization of the 2-RDM according to not only perturbation theory but also N-representability constraints in the form of CauchySchwarz inequalities.21 The N-representability conditions enable the degrees of freedom in the parametrized 2-RDM to model the effects of single and double excitations as well as higher excitations. Both the parametrizations of Kollmar (K)15,16 and Mazziotti (M)20,21 improve upon CCSD for moderately correlated systems as occur at transition states and during bond dissociations. Furthermore, the M energy functional approaches the accuracy of CCSDT at both equilibrium and nonequilibrium geometries. Importantly, the M functional adds fourth-order terms to the energy functional, suggested by the N-representability conditions, that have never appeared in coupled electron pair approximations. These additional terms can be added at a cost of r4 with r being the number of orbitals,21 and hence, the M parametric 2-RDM method retains the same overall computational scaling N2(r  N)4 of configuration interaction with single and double excitations. Both the K and M functionals, like the N-representability conditions, also possess particlehole symmetry, meaning that the particles and holes can be exchanged without altering the results of the functionals while the coupled pair methods typically break particle-hole symmetry. In comparison with the best ab initio electron-pair wave function methods, the parametric 2-RDM methods have a number of advantages. First, the computations in this paper as well as earlier benchmark calculations20,21,30 demonstrate that the parametric 2-RDM methods, particularly the M parametrization, generate significantly more accurate energies, geometries, and properties at a reduced computational cost (the benchmark calculations included comparisons with CCSD(T) and full configuration interaction). Second, the optimization of the energy can be performed variationally, which has several benefits including the validity of the HellmannFeynman theorem51 for the efficient calculation of derivatives for optimizing molecular geometries. Finally, the reduced computational cost of the parametric 2-RDM methods relative to coupled cluster with single and double excitations can also facilitate efficient parallelization of the code. Fewer intermediate matrices means that the

ARTICLE

2-RDM-based methods require less memory communication, which is typically a significant bottleneck in the treatment of electron correlation in electronic structure computations. With the improvements on earlier electron-pair theories and these advantages, the parametric 2-RDM methods provide a direct route to the calculation of 2-RDMs with a balance of accuracy and efficiency, as seen in the present work and previous benchmark calculations, that promises significant impact in computing energies and properties of many-electron systems in both chemistry and physics.

’ AUTHOR INFORMATION Corresponding Author

*E-mail: [email protected].

’ ACKNOWLEDGMENT C.A.S. and D.A.M. thank Kasra Naftchi-Ardebili for his illustration of the carbonic acid isomerization in Figure 1. D.A.M. gratefully acknowledges the NSF, ARO, Microsoft Corporation, Dreyfus Foundation, and David-Lucile Packard Foundation for support. ’ REFERENCES (1) Jensen, F. Introduction to Computational Chemistry, 2nd ed.; Wiley: West Sussex, 2007. (2) Parr, R. G.; Yang, W. Density Functional Theory of Atoms and Molecules; Oxford University Press: New York, 1994. (3) Reduced-Density-Matrix Mechanics: With Application to ManyElectron Atoms and Molecules; Mazziotti, D. A., Ed.; Wiley: New York, 2007; Vol. 134. (4) Coleman, A. J.; Yukalov, V. I. Reduced Density Matrices: Coulson’s Challenge; Springer-Verlag: New York, 2000. (5) Mayer, J. E. Phys. Rev. 1955, 100, 1579. (6) L€ owdin, P. Phys. Rev. 1955, 97, 1474. (7) Coleman, A. J. Rev. Mod. Phys. 1963, 35, 668. (8) Garrod, C.; Percus, J. J. Math. Phys. 1964, 5, 1756. (9) Harriman, J. E. Phys. Rev. A 1978, 17, 1257. (10) Mazziotti, D. A.; Erdahl, R. M. Phys. Rev. A 2001, 63, 042113. (11) Nakata, M.; Nakatsuji, H.; Ehara, M.; Fukuda, M.; Nakata, K.; Fujisawa, K. J. Chem. Phys. 2001, 114, 8282. (12) Zhao, Z.; Braams, B. J.; Fukuda, H.; Overton, M. L.; Percus, J. K. J. Chem. Phys. 2004, 120, 2095. (13) Mazziotti, D. A. Phys. Rev. Lett. 2004, 93, 213001. (14) Shenvi, N.; Izmaylov, A. Phys. Rev. Lett. 2010, 105, 213003. (15) Kollmar, C. J. Chem. Phys. 2006, 125, 084108. (16) DePrince, A. E., III; Mazziotti, D. A. Phys. Rev. A 2007, 76, 042501. (17) DePrince, A. E., III; Kamarchik, E.; Mazziotti, D. A. J. Chem. Phys. 2008, 128, 234103. (18) DePrince, A. E., III; Mazziotti, D. A. J. Chem. Phys. 2010, 132, 034110. (19) DePrince, A. E., III; Mazziotti, D. A. J. Chem. Phys. 2010, 133, 034112. (20) Mazziotti, D. A. Phys. Rev. Lett. 2008, 101, 253002. (21) Mazziotti, D. A. Phys. Rev. A. 2010, 81, 062515. (22) Cizek, J.; Paldus, J. Phys. Scr. 1980, 21, 251. (23) Purvis, G. D.; Bartlett, R. J. J. Chem. Phys. 1982, 76, 1910. (24) Piecuch, P.; Kucharski, S. A.; Kowalski, K.; Musial, M. Comput. Phys. Commun. 2002, 149, 71. (25) Hoffmann, M. R.; Schaefer, H. F. Adv. Quantum Chem. 1986, 18, 207. (26) Weinhold, F.; Wilson, E. B. J. Chem. Phys. 1967, 46, 275. (27) Weinhold, F.; Wilson, E. B. J. Chem. Phys. 1967, 47, 2298. 12015

dx.doi.org/10.1021/jp2057805 |J. Phys. Chem. A 2011, 115, 12011–12016

The Journal of Physical Chemistry A

ARTICLE

(28) Ahlrichs, R. Comput. Phys. Commun. 1979, 17, 31. (29) Ahlrichs, R.; Scharf, P.; Ehrhardt, C. J. Chem. Phys. 1985, 82, 890. (30) Schwerdtfeger, C. A.; DePrince, A. E., III; Mazziotti, D. A. J. Chem. Phys. 2011, 134, 174102. (31) Bernard, J.; Seidl, M.; Kohl, I.; Liedl, K. R.; Mayer, E.; Galvez, O.; Grothe, H.; Loerting, T. Angew. Chem., Int. Ed. 2010, 49, 1. (32) Mori, T.; Suma, K.; Sumiyoshi, Y.; Endo, Y. J. Chem. Phys. 2009, 130, 204308. (33) Mori, T.; Suma, K.; Sumiyoshi, Y.; Endo, Y. J. Chem. Phys. 2011, 134, 044319. (34) Hage, W.; Liedl, K.; Hallbrucker, A.; Mayer, E. Science 1998, 279, 1332–1335. (35) Adamczyk, K.; Premont-Schwarz, M.; Pines, D.; Pines, E.; Nibbering, E. T. J. Science 2009, 326, 1690–1694. (36) Terlouw, J. K.; Lebrilla, C. B.; Schwarz, H. Angew. Chem., Int. Ed. 1987, 26, 354. (37) Nguynen, M. T.; Ha, T. K. J. Am. Chem. Soc. 1984, 106, 599. (38) Janoschek, R.; Csizmadia, I. G. J. Mol. Struct. 1993, 300, 637. (39) Wight, C. S.; Boldyrev, A. I. J. Phys. Chem. 1995, 99, 12525. (40) Loerting, T.; Tautermann, C.; Kroemer, R. T.; Kohl, I.; Hallbrucker, A.; Mayer, E.; Leidl, K. R. Angew. Chem., Int. Ed. 2000, 39, 891. (41) Kumar, P. P.; Kalinichev, A. G.; Kirkpatrick, J. J. Chem. Phys. 2007, 126, 204315. (42) Murillo, J.; David, J.; Restrepo, A. Phys. Chem. Chem. Phys. 2010, 12, 10963. (43) Raghavachari, K.; Trucks, G. W.; Pople, J. A.; Head-Gordon, M. Chem. Phys. Lett. 1989, 157, 479. (44) Wloch, M.; Gour, J. R.; Kowalski, K.; Piecuch, P. J. Chem. Phys. 2005, 122, 214107. (45) Schmidt, M. W.; Baldridge, K. K.; Boatz, J. A.; Elbert, S. T.; Gordon, M. S.; Jensen, J. H.; Koseki, S.; Matsunaga, N.; Nguyen, K. A.; Su, S.; Windus, T. L.; Dupuis, M.; Montomery, J. A. J. Comput. Chem. 1993, 14, 1347. (46) Dunning, T. H., Jr J. Chem. Phys. 1989, 90, 1007. (47) Klopper, W.; Helgaker, T. Theor. Chem. Acc. 1998, 99, 265. (48) Valeev, E. F.; Allen, W. D.; Hernandez, R.; Scherrill, C. D.; Schaefer, H. F., III. J. Chem. Phys. 2003, 118, 8594. (49) McQuarrie, D. A. Statistical Mechanics; Harper and Row: New York, 1976. (50) Coulson, C. A. Rev. Mod. Phys. 1960, 32, 170. (51) Feynman, R. P. Phys. Rev. 1939, 56, 340.

12016

dx.doi.org/10.1021/jp2057805 |J. Phys. Chem. A 2011, 115, 12011–12016