Second-Order State-Specific Multireference Møller−Plesset

3668

J. Phys. Chem. A 2010, 114, 3668–3682

Second-Order State-Specific Multireference Møller-Plesset Perturbation Theory (SS-MRMPPT) Applied to Geometry Optimization Uttam Sinha Mahapatra† Department of Physics, Taki GoVernment College, Taki, North 24 Parganas-743429, India

Sudip Chattopadhyay* Department of Chemistry, Bengal Engineering and Science UniVersity, Shibpur, Howrah 711103, India

Rajat K Chaudhuri‡ Indian Institute of Astrophysics, Bangalore 560034, India ReceiVed: December 7, 2009; ReVised Manuscript ReceiVed: January 30, 2010

The performance of a numerically oriented gradient scheme for the previously introduced second-order statespecific multireference Møller-Plesset perturbation theory (SS-MRMPPT) has been tested to compute certain geometrical parameters (such as bond lengths and angles). Various examples [H2O, O3, N2H2, C2H4, C2H2F2, 1,3-butadiene (C4H6), cyclobutadiene (C4H4), and 2,6-pyridynium cation (C5NH+ 4 )] have been presented to validate the implementation of the SS-MRMPPT gradient approach. To illustrate the reliability of our findings, comparisons have been made with the previously reported theoretical results. The accuracy of our calculations has further been assessed by comparing with the experimental results whenever available. On the basis of the present work, we arrive at the conclusion that the SS-MRMPPT gradient scheme has substantial potential in computing the geometrical parameters for several rather diverse molecular systems, whether charged or neutral and having the closed-shell ground state or being open-shell radicals or biradicals or strongly perturbed by intruders. It is worthwhile to emphasize that the present work represents the first systematic application of the SS-MRMPPT numerical gradient approach. I. Introduction Energy derivatives are as important as the energies themselves. Since its introduction to quantum chemistry via the pioneering work of Pulay,1,2 the ab initio gradient approach has emerged as a very elegant and reliable avenue as well as a computationally affordable tool to search stationary points (minima and saddle points) in the potential energy surface (PES) of small- to large-sized molecules and hence has been used in the field of investigation of geometries, vibrations, chemical reactions, energy relaxation processes, and dynamics of molecules of arbitrary complexity and generality. It is now a well-accepted fact that for situations such as reaction paths including stretched or dissociating bonds, nonequilibrium (transition state) geometries, or excited states where one cannot logically find a dominant single-reference function, the use of multireference methods is not only useful but also highly desirable. Generally, the key to a successful description of such situations is an accurate and balanced treatment of dynamical and nondynamical electron correlation effects. The development of multireference methods represents important progress in electronic structure theory in the last decades (see refs 3–7 for selected reviews). We now discuss the problem of balancing dynamical and nondynamical correlations in the context of MR-based perturbative theory. * To whom correspondence should [email protected]. † E-mail: [email protected]. ‡ E-mail: [email protected].

be

addressed.

E-mail:

The multireference perturbation theories (MRPT) are among the most common and easily applied MR methods to investigate the above-mentioned cases as they are computationally inexpensive. A considerable effort has been put into developing several MRPT approaches for correlation treatments for MR and open-shell systems. However, there are many problems with implementation of the methods, so that numerous variants of MRPT formulations have appeared. Most of them were built using the framework of the effective Hamiltonian formalism. Unfortunately, effective Hamiltonian based conventional MR methods are often faced with a practical problem, that of intruder states, which can arise regardless of whether a complete or a general model space is used. Hence, the first decisive step for any MR calculation is clearly a proper choice of references. To attenuate this issue, in recent times, the state-specific versions of MR-based methods have made remarkable progress and appear in useful applications for studies of various chemical phenomena and situations. Despite that, in many cases, a single state of the system is the focus of the calculation. The main advantage of the SS approach stems from the fact that the intruder problem is avoided in a natural manner by focusing on a single state at a time. Different formulations of the statespecific MRPT method have been proposed, namely, CASPT28 of Roos and co-workers, MRMPPT of Hirao,9 CIPSI10 of Malrieu and co-workers, SS-MRPT(Mk-MRPT)11 of Mukherjee and co-workers, and works by Davidson,12 Surj´an and coworkers,13 and so on. Hoffmann and co-workers14 have come up with a MRPT approach that is an intruder-free scheme and termed the method the generalized Van Vleck perturbation theory, abbreviated as GVVPT. This approach is distinguishable

10.1021/jp911581f  2010 American Chemical Society Published on Web 02/12/2010

SS-MRMPPT Applied to Geometry Optimization from other MRPT approaches in vogue by the fact that it is subspace-specific (an intermediate Hamiltonian based approach) and it gives an allowance for an interaction of the perturbed states of interest with the unperturbed complementary ones. All of these approaches have been developed with the intention of finding the best possible partitioning scheme because the choice of the zero-order Hamiltonian operator in MRPT is not unique, and a balance between accuracy and efficiency had to be found. Although these methods combine several attractive features, they suffer from the lack of size-extensivity of computing energy even when complete model spaces (CMS) are considered (except CIPSI and SS-MRPT). The size-inextensivity is one of the main sources of errors in the calculations on energies even for medium-sized molecules. As a passing remark, we mention that the MRPT method of Freed and co-workers15 alleviates the loworder intruder state problems simply by imposing a well-defined energy gap between the reference energies of the states corresponding to the model space and the states corresponding to the secondary space. The CAS-based SS-MRPT approach,5,11 which is exactly sizeextensive and explicitly spin-free, seems to be numerically promising and computationally competitive. It thus seems that an effort should be made at further development of the SSMRPT proposed by Mahapatra et al.11 The SS-MRPT method is based on second-order perturbation theory. The first-order wave function is obtained in the SS-MRPT method as an iterative solution to a large set of linear equations. As the SSMRPT method is state-specific in nature, only one of the eigenvalues obtained by diagonalization of the effective Hamiltonian has a physical meaning, and the rest are artificial. We do not present a detailed illustration here but refer to the original papers.11 As that for the other MRPT methods, correct selection of the active space is a very importance issue. All neardegeneracy effects leading to configurations with moderate to large weights must be included at the stage of the zero-order calculation. If this is not done, the first-order wave function will contain large coefficients. If the energy contribution from such a configuration is large, the results are not to be trusted, and a new selection of the active space should be made. The SS-MRPT method allows one to circumvent the intruder state problem in a size-extensive manner using a physically reasonable and mathematically well-defined procedure that ensures the generation of smooth energy curves, even when the curves are in close proximity. It is necessary to note that conventional MRPT is not amenable to calculating globally continuous meaningful energy derivatives due to the perennial intruder state problem. The method is size-consistent with respect to fragment separations using localized orbitals. We should mention here that the use of single-root MR methods is not beneficial when electronic states of different character are close in energy, which may result in an erratic mixing of the different types of wave functions in the reference functions. In this context, we mention the MCQDPT method,16 which is a multiconfiguration basis multiroot PT. It includes MRMPPT,9 which is multiconfiguration but single-state PT, and the conventional quasidegenerate PT (QDPT),17 which is a multistate but single-configuration PT, as special cases. It should be emphasized that neither the multiroot MCQDPT16 nor the single-root MRMPPT9 is rigorously sizeextensive in nature (see ref 18). Here, we want to state that to demonstrate the relative performance, recently proposed various multireference perturbation theories have been compared numerically on representative examples at the second order of approximations by Chaudhuri et al.19 and Hoffmann et al.20 In

J. Phys. Chem. A, Vol. 114, No. 10, 2010 3669 these papers, the authors briefly discuss some of the basic differences among the MRPT methods considered. The SS-MRPT method using Møller-Plesset (MP) as well as Epstein-Nesbet (EN) partitions with respect to the RayleighSchro¨dinger (RS) and a Brillouin-Wigner (BW) perturbation expansions has been tested in a large number of applications, and we have observed that the SS-MRPT approach has been quite successful.5,11,21,22 The method can be viewed as a more economic alternative to MRCI for treatments of electronic state(s) with a change of dominating electronic configurations due to slight distortion of geometry. Very satisfactory performance of the MP-based SS-MRPT (SS-MRMPPT) method within the framework of RS expansion prompted us to extend this formalism to make geometry optimizations feasible for electronic states strongly perturbed by intruders and also for those with pronounced multireference character. The optimized geometry is the geometry which minimizes the strain on a given system. Any perturbation from this geometry will induce the system to change so as to reduce this perturbation unless prevented by external forces. It is now safe to say that no quantum chemical method can survive for long without an associated capability for gradients. Gradients provide structures and the second derivatives to characterize critical points while providing vibrational frequencies. Very recently, Prochnow et al.23 presented a formulation of analytic gradients for the SSMRCC of Mukherjee and co-workers5 along with numerical implementation. To handle larger MR systems for which the application of the SS-MRCC gradient scheme is computationally quite demanding, the development of a gradient scheme for its companion perturbation theory is highly desirable. The development of gradients for the SS-MRMPPT adds greatly to the computational chemist’s ability to study chemical phenomena related to multireference electronic states. The SS-MRMPPT gradient procedure appears to provide the best compromise between accuracy and computational feasibility. Since the SSMRMPPT method works well both in the regimes of quasidegeneracy and also where there are intruders in the traditional effective Hamiltonian formalisms, it should be very useful for computation of geometrical parameter computations via a gradient scheme over a wide range of energy curves. Thus, the SS-MRMPPT gradient approach enables us to describe reaction paths very accurately as long as the target state energy is wellseparated from the energies of the virtual functions. There exist broadly two approaches to compute the derivatives. The simplest is numerically oriented, where the effect of a small external agency on the energy functional of the perturbed correlated state is monitored numerically and the response properties are extracted by finite difference approximations to the various derivatives. The other approach is analytical, which is based on the analytical expressions of the response functions to compute the required derivatives. Both methods have found widespread use. Fully numerical derivatives are inherently less efficient than analytic derivatives; however, the numerical derivatives are more amenable to very efficient coarse-grained parallel computing. Incidentally, the derivation of analytical SS-MRMPPT energy gradients is not at all straightforward, and a computational algorithm for their evaluation has yet to be published. In the present work, we apply the numerical gradient procedure for the SS-MRMPPT method to calculate geometrical parameters such as bond length and bond angles using a simple program implemented within the GAMESS(US)24 suite for computing SS-MRMPPT derivatives numerically and hence for optimizing geometries. Fully numerical first and second derivative codes

3670

J. Phys. Chem. A, Vol. 114, No. 10, 2010

are available in GAMESS(US); therefore, one can optimize molecular geometries with any level of theory. As that of the parent SS-MRMPPT method, the SS-MRMPPT gradient method combines several attractive advantages; (i) it is applicable to a wide class of phenomena and a wide variety of systems, (ii) is much more useful and handy than the corresponding MRCI and MRCC procedures, (iii) is highly stable if reference CAS is appropriately chosen, and (iv) is size-extensive, guaranteeing proper scaling with the number of particles or units in the system. The SS-MRMPT gradient approach utilizes a relaxed coefficient description of the multiconfigurational reference function during the inclusion of dynamical correlation. Allowing for this flexibility in the wave function is often crucial, particularly in cases where the dynamic correlation strongly interacts with the model space contributions (as it happens, for example, in the ground state of the LiF molecule at nonequilibrium distances). However, in the case of the SS-MRMPPT, an unrelaxed scheme, where the coefficients of the model space determinants are held fixed, can also be considered. It may lead to faster convergence and still very accurate calculations. Relatively fewer works exist in the literature which are based on analytical gradients for multireference methods23,25–34 because of the much greater complexity of the formalisms compared to the SR case. There is still no widely applicable and generally available MR-based gradient procedure. Gradients are most easily obtained by using the method of undetermined Lagrangian multipliers.26 In passing, we may say that the application of the single-reference-based gradient theories35 is limited only to the systems where the Hartree-Fock approximation is a good starting point. In this paper, we report only a pilot implementation of a numerical gradient scheme for the SS-MRMPPT and discuss the efficacy of the methodology for the calculations of geometrical parameters for H2O, O3, N2H2, C2H4, C2H2F2, butadiene (C4H6), cyclobutadiene (C4H4), and 2,6-pyridynium cation (C5NH+4 ) systems that require a MR wave function as the zerothorder approximation. Here, we have studied different types of systems including neutral, charged, as well as biradical ones. There have been a number of previous studies of such systems in the literature of electronic structure theories by means of which one can assess the applicability and potentiality of our newly implemented gradient method. Our aim for the present work is not to present a high-precision ab initio study of the geometries of the systems considered here but rather to provide a simple illustration of the effectiveness of the numerical gradient SS-MRMPPT approach. We mention here that the diradicals represent the most clear-cut application of the SFbased formalism and its different variants36 because in these systems, the nondynamical correlation derives from a single HOMO-LUMO pair. SF-based methods are also useful to describe many MR situations. The numerical test of the SS-MRMPT gradient method presented in this paper for various diverse systems clearly demonstrates its ability to generate accurate geometrical parameters, as well as its other advantages, such as the absence of the intruder state problem that often plagues other MR-based PT and CC theories, while at the same time enjoying a very favorable computational cost with the size of the problem. Thus, this work will hopefully rekindle the vanishing interest in the standard numerical gradient schemes. In the present paper, we briefly recall the equations of the SS-MRMPPT method and its underlying ideas in section II. The results of calculations are summarized and discussed in section

Mahapatra et al. III. In section IV, we discuss the problems and future perspectives of the approach. II. Brief Introduction to the SS-MRMPPT Method This section concerns itself with a succinct discussion of the structural features of the SS-MRMPPT method along with the pertinent theoretical issues. For details, we refer to the earlier papers.11,21 In the SS-MRMPPT formalism, the exact wave function ψ is parametrized using the Jeziorski-Monkhorst ansatz37 in a state-specific fashion. In the SS-MRMPPT approach, model space coefficients cµ as well as molecular orbitals (MOs) are determined first through CAS, and then, RayleighSchro¨dinger PT is applied via a multipartitioning MP scheme38 using the CAS wave function as a reference, ψ0 ) ∑µ cµφµ, where φµ is the model space function (configuration state functions, CSF). In an attempt to incorporate compactness of presentation of the working equations and the issues pertaining to the avoidance of intruders, the equation for the first-order cluster amplitudes of SS-MRMPPT is represented as ν*µ

tµl(1)

Hlµ ) + [E0 - Hll]

∑ 〈χlµ|Tν(1)|φµ〉Hµν(c0ν/cµ0) ν

[E0 - Hll]

(1)

where the state energy E0 is defined by ∑ν Hµνcν0 ) E0cν0. Here, the cluster operators in Tµ excite the model function φµ to the virtual functions χµl , and the cµ0 ’s are the combining coefficients of the model functions {φµ}. The cluster amplitude-determining equations (and hence the theory) avoid the intruder problem in a natural fashion as long as E0 is well-separated from virtual functions, which often occurs for low-lying states including the ground state. It should be mentioned that the first term in the eq 1 is strictly first-order in the cluster operator, whereas the last term (coupling term) is not necessarily so. This coupling term is mainly responsible for providing the size-extensive and intruder-free nature of the theory. The working equations of the SU- and SS-based MRPT methods differ in the form of the coupling terms. However, a repeated solution of essentially the same set of equations for each target state is required in the case of the SS-MRMPPT method, in contrast to the conventional SU approach where all states are considered simultaneously. After the solution of the cluster amplitudes, the coefficients and the required energy of the target state are generated by ˜ (2) diagonalizing an effective operator (non-Hermitian) H µν defined in CAS (or CMS)

∑ H˜µν(2)cν ) E(2)cµ

(2)

ν

(2) ˜ µν where H ) Hµν + ∑l Hµνtνl(1)(ν). Since the theory is statespecific, only one root (representing the exact intruder-free energy) is meaningful, while the rest are extraneous, which differentiates it from the state-universal-based method(s). In eq 2, there is a choice; one can use it to compute the energy either an unrelaxed (as an expectation value with respect to the unrelaxed function) or a relaxed mode. The SS-MRMPPT method with fixed reference coefficients is a rigorous secondorder approach. The unrelaxed scheme follows trivially if one fixes the model space coefficients to some preassigned/initial values. In this paper, we have used the relaxed description since it is more general. A typical feature of the present SS-MRMPPT

SS-MRMPPT Applied to Geometry Optimization scheme is the use of the zeroth order coefficients cµ0 to compute the cluster amplitudes and the dressed Hamiltonian but allow the coefficients to relax while computing E since this is obtained via diagonalization. The SS-MRMPPT method applied here is a pseudo-second-order approach as it utilizes a relaxed coefficients description of the multideterminantal reference function. It should be emphasized that the strict separability of the effective Hamiltonian ensures the separability of its eigenvalues provided that CMS and localized orbitals are used. The effective Hamiltonian is computed as a sum of the product of molecular integrals and coupling constants between CSFs divided by the energy difference; therefore, the computer resources required do not depend strongly on the size of the active space and basis set. This presents a stark contrast to the case of MRCI and MRCC. The SS-MRMPPT uses the best traits of the multipartitioning strategy38 as well as those of a rigorously size-extensive formulation. In our MP partitioning scheme, H0 is defined with respect to φ0µ (being the largest closed-shell component of φµ) taken as the vacuum. This leads to a spin-adapted version of SS-MRMPPT, where the cluster operators are defined in terms of the spin-free unitary generators.5 The way that eq 1 has been cast, it becomes evident that the equation is ill-behaved for the very small value of cµ0 . This is not due to an intruder effect (related to the vanishing energy difference in the denominator of the cluster-finding equations leading to convergence difficulties), but it rather due to illconditioning of the working equations. Thus, this is less serious. The variety of systems and CAS that we have explored until now (including the present work) to implement the SSMRMPPT method have all equivocally revealed that the method is free from from any numerical instability of the amplitudedetermining equations even when all of the coefficients (cµ) are included. This may be explained by realizing the fact that if ˜ µνcν/ the reference coefficients tend to 0, the value of 〈χµl |Tν(1)|φµ〉H [〈χµl |Tν(1)|φµ〉Hµνcν0] becomes small and very much akin to the values of the coefficients cµ. As in any reference function (ψ0), if the contribution of some of its component functions (φµ) is very small, then the back coupling from other components in the cluster-finding equation should not be large, at least not larger than the zero-order value, that is, 〈χµl |Tν|φµ〉Hµνcν will be as small as cµ. However, if somehow any sort of numerical instability creeps in due to nearly vanishing values of cµ resulting in an explosion in the magnitude of the amplitudes associated with any nonzero determinant which has been solely created by this reference, one may envisage irregular and incomprehensible behavior of the perturbative corrections. A viable avenue might be to simply delete the cluster operator Tµ corresponding to the model function φµ causing the problem. This, however, might invite questions in the sense that such a deliberate omission eradicates all other determinants which were solely generated from this “special reference” also. However, such a doubt has been clarified recently by Engels-Putzka and Hanrath39 in the context of the N2 system by demonstrating that most of such determinants have a vanishing coefficient and the effect on the energy is negligibly small (see also ref 20). Below, we report SS-MRMPPT gradient results for a molecular systems. III. Results and Discussion In this section, a method for computing SS-MRMPPT energy gradients numerically has been implemented in calculations of the geometrical parameters of the ground state for various interesting systems. A series of calculations concerning the

J. Phys. Chem. A, Vol. 114, No. 10, 2010 3671 ground state of H2O, O3, N2H2, C2H4, C2H2F2, 1,3-butadiene (C4H6), cyclobutadiene (C4H4), and 2,6-pyridynium cation (C5NH4+) systems have been carried out so as to assess the performance of the SS-MRMPPT gradient approach. These systems are usually used to illustrate the efficacy of almost any many-body method. Consequently, such systems are wellstudied, and therefore, comparison of our results is possible with many previous MR-based theoretical calculations to calibrate the SS-MRMPPT gradient method. To illustrate the influence of the multireference treatment, we compare the geometrical parameters determined here with those obtained at the SR level. We emphasize that although comparison of theoretical results involving identical active spaces and bases is physically more appealing, in the present paper, we have cited results of various theories using different schemes to get a feeling about the performance of the present method. Finally, our results have been calibrated against the experimental values since we do not have the exact FCI energies for our comparison. In our numerical analysis, geometries have also been compared with the corresponding CASSCF structures, illustrating the effect of including the correction for dynamical electron correlation. To demonstrate systematically how the nature of the basis set affects the SS-MRMPT geometrical parameters, the various basis sets have been applied to selected cases at the CASSCF and SSMRMPPT levels. It can be seen from the present work that SSMRMPPT geometries show a strong basis set dependence. In all cases, we employ the ground-state CASCSCF orbitals unless stated otherwise. The geometries have been optimized to within a maximum gradient of 10-4. Unless stated otherwise, in most of our calculations, we use CASSCF orbitals optimized with respect to the dominant configuration of the ground state. The calculations using (single-state) CASSCF orbitals are more reliable than those employing RHF MOs since the latter provide a better description in the MR situations. The geometry optimization has been done by interfacing the SS-MRMPPT(RS) method into the GAMESS(US) quantum chemistry software. In this paper, mainly two families (correlation-consistent polarized valence n-tuple zeta, cc-pVnZ, and augmented correlation-consistent polarized valence n-tuple zeta, aug-cc-pVnZ) of basis sets40 have been considered. Spherical harmonics rather than Cartesian (d/f) functions have been used throughout in our present computations. Finally, it should be mentioned that in the case of N2H2, C2H2F2, and 1,3-C4H6, other lowest singlet states apart from the ground state have been considered and briefly discussed. A. Water (H2O). Our first example is the geometry optimization of the ground state of the H2O molecule. We chose the water molecule as our benchmark system for a number of reasons. Water is a prototype system for a variety of spectroscopic and reaction dynamics studies. There are many applications involving the water molecule in which the knowledge of reliable geometry is required. This system has been extensively studied in the past as a critical test to establish the applicability of SR- and MR-based methods.41–46 Thus, the determination of its geometrical parameters continue to be the subject of recent investigations. Due to a reasonable multireference character of its ground state, dynamical as well as nondynamical correlation effects play a significant role when determining the equilibrium geometry of this system. Although the DZP basis is not large enough to provide highly accurate results (as it is not good enough to describe electron correlation and, particularly, the higher excitation effects), to facilitate a comparison with other results,41,42 we employ the

3672

J. Phys. Chem. A, Vol. 114, No. 10, 2010

Mahapatra et al.

TABLE 1: Results of Equilibrium Geometries for the Ground-State H2O Moleculea basis DZP cc-pVDZ aug-ccpVDZ aug-ccpVTZ cc-pVDZ DZP

cc-pvDZ

cc-pvTZ

∞Z

methods

Re

θ

CASSCF(6,6) SS-MRMPPT(6,6) CASSCF(6,6) SS-MRMPPT(6,6) CASSCF(6,6) SS-MRMPPT(6,6) CASSCF(6,6) SS-MRMPPT/CASSCF(6,6) SR-CCSDb 2R-RMR-CCSDb 3R-RMR-CCSDb CISDc CCSDc BDc ODc CISDd CISD[TQ]d CISDTQd SR-SDCIe MP3e CCSD(T)e MRSDCIe MRAQCCe SR-SDCIe MP3e CCSD(T)e MRSDCIe MRAQCCe SR-SDCIe MP3e CCSD(T)e MRSDCIe MRAQCCe experimentf

0.9674 0.9621 0.9699 0.9668 0.9685 0.9670 0.9577 0.9574 0.965 0.965 0.965 0.9577 0.9610 0.9609 0.9609 0.9623 0.9670 0.9688 0.9613 0.9618 0.9658 0.9666 0.9669 0.9520 0.9532 0.9578 0.9592 0.9599 0.9490 0.9569 0.9558 0.9569 0.9578 0.957

103.89 103.66 101.81 101.32 102.96 103.23 101.06 102.56 102.0 102.0 102.0 104.63 104.64 104.64 104.1 103.6 103.7

104.5

a

Bond lengths and bond angles are given in Angstroms (Å) and degrees, respectively. b Reference 45. c Reference 44. d Reference 43. e Reference 46. f Reference 48.

often used double-ζ plus polarization (DZP) basis. We also use cc-pVDZ, aug-cc-pVDZ, and aug-cc-pVTZ basis sets40 to demonstrate the effect of basis sets on the quality of the geometrical parameters. In our calculations, we use the C2V symmetry-adapted (6,6) active space [CAS(6,6)], which contains three sigma bonding orbitals and three corresponding antibonding orbitals. In our study, the lowest occupied (1a1) orbitals are frozen. Table 1 compares the results provided by the SS-MRMPPT gradient approach for various basis sets. As the results obtained with standard correlated methods [such as CISD, Brueckner CCD (BD and OD), SRCCSD, R-RMR CCSD, and so forth] are well-cited in the literature,43–47 we tabulate them in the same table for the sake of comparison so that one can get a clear picture about the performance of the SS-MRMPPT gradient approach. It is also worthwhile to compare the SS-MRMPPT results to the popular CISD[TQ] and CISDTQ methods.43 For further calibration of our results, we also compared them with experimental findings.48 Here, we mentioned that the all-electron full-CI equilibrium bond length (Re) for H2O with the cc-pVDZ basis has been determined by Ka´llay et al.49 to 0.96616 Å. The satisfactory quality and reliability of the SS-MRMPPT geometries is clearly evident from the comparative study of the results presented in Table 1. The optimized geometry of the ground state of H2O from SS-MRMPPT calculations is in good agreement with experiment. We now turn our attention to the performance of SR-based methods for the water molecule. Not

surprisingly, the SR-based methods perform well, all being in good agreement with experiment. From our numerical analysis, it is evident that the SS-MRMPPT results are very close to the previous SR-CCSD and RMR-CCSD45 as well as the CISD, BD, and OD43,44 ones. Despite the modest size of the model space, the SS-MRMPPT-predicted geometry is in good agreement with those obtained from the computationally expensive CISD[TQ] and CISDTQ. This is the most remarkable feature of the table. Nonetheless, we must emphasize that the algorithms associated with CC-based methods are computationally very demanding and thus cost-ineffective in comparison to the SSMRMPPT scheme. The SS-MRMPPT gradient method achieving near CCD, OD, and R-MRCCSD quality geometries provides support for its efficacy. B. Ozone(O3). As the next example, we choose ozone, O3, which is one of the most important reactive species in the lower parts of the atmosphere, where a comparison with other singlereference and multireference treatments and with experimental data is available. Ozone is a common example of a molecule in which near-degeneracy effects are very pronounced even in the equilibrium region and represents a particularly challenging PEC (has the double minima), the open structure (C2V) and the equilateral triangle ring structure (D3h).32,45,50–52 The electronic structure of ozone is a mixture of ionic and diradical valence structures, and one can therefore expect some multiconfigurational character for the wave function. Thus, it is difficult to investigate its molecular structure using the conventional SRbased methods. It is now a well-illustrated fact that a highquality and balanced description of the dynamical correlation and near-degeneracy effects is required to yield accurate geometry. The ozone molecule provides a classical example of the fact that neither the SR-based method nor the traditional MR method can satisfactorily describe states having a strong multireference character.45 Therefore, it is again a good example to test the performance of the SS-MRMPPT gradient approach with a small reference space. We optimize its equilibrium geometry using the SS-MRMPPT numerical energy gradient. To get correct molecular geometries of O3, the most crucial step is an appropriate choice of the reference space. In present calculations, we have used C2V symmetry. A minimal active space is the two-electron/two-orbital one, that is, the CAS(2,2) space involving the highest occupied molecular orbital (MO) and the lowest unoccupied MO. The ground state of O3 is known to require a minimum CAS(2,2) reference space for its correct description. We have used DZP and aug-ccpVXZ type (X ) D and T) basis sets.40 Although the DZP40 basis set used for this test study is obviously inadequate for a meaningful comparison with experiment, the examination of ozone with DZP provides an opportunity for further assessment of the SS-MRMPPT method itself with respect to other methods as the study of O3 using DZP basis by various methods has been reported in the literature.32,45,50–52 To establish the accuracy and applicability of our method, it is worth comparing our results with the results of CASPT2. In this context, we want to mention that CASPT2(2,2) and CASPT2(6,6) are inadequate for accurate description of geometrical properties of ozone.50 Table 2 depicts the optimized geometrical structures calculated for ozone using the SS-MRMPPT gradient along with results obtained with other perturbative and nonperturbative methods. A comparison with the CCSD, 2R-RMR-CCSD, MBPT, MRMP, and CASPT2 results provides useful measure for the utility and potentiality of the SS-MRMPPT based gradient methodology. The resulting bond distance via the CASPT2 method using a different CAS space and basis52 is

SS-MRMPPT Applied to Geometry Optimization

J. Phys. Chem. A, Vol. 114, No. 10, 2010 3673

TABLE 2: Results of Equilibrium Geometries for the Ground State of the O3 Molecule, Using the Methods Described in the Texta basis

methods

Re

θ (deg)

DZP aug-ccpVDZ

SS-MRMPPT/CASSCF(2,2) SS-MRMPPT/RHF CASSCF(2,2) SS-MRMPPT(2,2) SS-MRMPPT/RHF CASSCF(2,2) SS-MRMPPT(2,2) MCSCFb MP2b GVVPT2b CISDb SR-CISDc TCSCF-CISDc CISD[TQ]c SR-CCSDc SR-CCSD(T)c SR-CCSDTc 2R-RMR-CCSDc 3R-RMR-CCSDc CASPT2(2,2)d CASPT2(6,6) CASPT2(8,7) CASPT2(12,9) CASPT2(18,12) MBPT2(2)e MBPT2(3)e RASSCFe MRCI 10e CASPT2-De CASPT2-Ne CASPT2(12,9)e POL-CIf MCSCF-CIf MRD-CIf MRMPf expt.g

1.2631 1.2575 1.2445 1.2813 1.2461 1.2424 1.2740 1.2572 1.2791 1.2700 1.2711 1.247 1.271 1.281 1.263 1.287 1.286 1.276 1.277 1.2683 1.2779 1.2882 1.2940 1.2954 1.2920 1.2440 1.2740 1.2710 1.2860 1.2780 1.2775 1.299 1.277 1.29 1.291 1.2717

118.15 118.22 115.41 115.75 118.33 115.83 116.14 115.1 117.3 117.6 116.2 117.7 116.2 116.7 117.4 116.8 116. 116.5 116.7 117.19 116.62 115.97 116.26 116.24 116.7 117.0 116.8 116.6 115.9 116.5 116.70 116.0 116.1 116.0 116.6 116.8

aug-ccpVTZ DZP

cc-pVTZ

ANO(4s3p2d)

a Bond lengths and bond angles are given in Angstroms (Å) and degrees, respectively. b Reference 32. c Reference 45. d Reference 52 [with CASSCF(8,7)]. e Reference 50. f Reference 51. g Reference 53.

also reported in the same table. The best results of CASPT2 for the O3 geometry have been reported by Ljubic´ and Sabljic´52 with the full outer valence active space, CAS(12,9), and extrapolation from Dunning’s weighted core-valence basis sets.40 The results of SS-MRMPPT calculations are compared with those obtained from experiment53 for further assessment. From the results of Table 2, we have found that the SSMRMPPT results are promising with respect to the experimental observations as the difference between them is modest. As illustrated in the table, the geometry obtained by SS-MRMPPT is reasonably close to MRMP9 (which bears close kinship with the SS-MRMPPT method). The results in Table 2 also show that the SS-MRMPT predictions are in close agreement with the results of CASPT2, CASPT2-D (or N), MBPT, GVVPT2 (generalized Van Vleck perturbation theory), MRCI10, and RASSCF (restricted active space SCF) methods. The table clearly depicts the fact that our computed equilibrium geometry is competitive with the prediction of computationally expensive full-blown CCSD, CCSDT, CCSD(T), and R-RMR-CCSD (the “externally corrected” CCSD method, involving a MR-CISD part followed by a corrected SR-CCSD part). Keeping in mind that the active space is small, we see that a SS-MRMPPT/ CAS(2,2) approach describes satisfactorily the equilibrium structure of the ground state of the O3 system. It is to be

Figure 1. Diazene/diimide (N2H2).

emphasized that the SS-MRMPPT/CAS(2,2) approach is immensely cheaper than that of the R-MRCCSD and CASPT2(12,9),52 which reemphasizes the usefulness of the method. C. Diazene/Diimide (N2H2). Diazene/diimide is extensively used for stereospecific reductions of olefins and as a ligand for transition-metal complexes. Knowledge of accurate geometrical parameters for N2H2 is therefore a very important step toward the understanding of the reactions involved in such processes. As the electronic structure of this system has been subjected to exhaustive theoretical investigations (see refs 54–60), it is an excellent probing ground for testing the efficiency of newly proposed approaches in their ability to account for correlation effects. This system can be considered as a “real-life” alternative for the H4 model system as an essentially perfect twoconfiguration reference problem. Very recently, the rigid cis-trans rotation of the molecule has been studied as a test case via the BWMRCC method.56 There have been several considerable theoretical works on the computation of the barrier height associated with the isomerization or formation processes of the system.54,55,57–59 Chaudhuri et al.60 studied the torsional potential energy curve as well as geometries and relative stabilities of local minima of the various isomers of N2H2 using the method employing the improved virtual orbital multireference Møller-Plesset perturbation theory. Since the main thrust of the present attempt is to delineate the efficacy of the numerically oriented SS-MRMPPT gradient method, we restrict our considerations to the determination of the geometrical parameters for the trans, cis, and iso conformers of N2H2 (see Figure 1). In our calculations, we have used CAS(2,2) for the cc-pVTZ basis.40 The geometrical parameters of trans- and cis-N2H2 calculated by the SS-MRMPPT method are collected in Table 3. The corresponding values for iso isomer are summarized in Table 4. For all three isomers, there are several theoretical and a few experimental results reported that may serve for comparison purposes, and thereby, one can easily visualize the capability of the SS-MRMPPT gradient approach to provide geometrical parameters of the N2H2 system. Only experimental geometrical data are available for the trans isomer only. The accuracy of our results has been assessed by comparing with the results of Jensen et al.54 as well as those of Hwang and Mebel.58 The goodness of the calculations via the numerical gradient SS-MRMPPT method may further be assessed by the CCSD(T) results (with extrapolation to the basis set limit and inclusion of inner-shell correlation effects) of Martin and Taylor.57 Also given in these tables are the geometries obtained at MRCI/aug-cc-pVXZ level calculations by Biczysko et al.,59

3674

J. Phys. Chem. A, Vol. 114, No. 10, 2010

Mahapatra et al.

TABLE 3: Results of Equilibrium Geometries for the Ground State of the 1,2-Diimine Molecule (N2H2), Using the Methods Described in the Texta conformers cis

basis cc-pVTZ

aug-cc-pVDZ aug-cc-pVTZ aug-cc-pVTZ cc-pVDZ cc-pVTZ cc-pVQZ best cal. 6-31G* trans

cc-pVTZ

aug-cc-pVTZ aug-cc-PVQZ aug-cc-pVTZ cc-pVQZ cc-pVDZ cc-pVTZ cc-pVQZ best cal. 6-31G* cc-pVDZ

cc-pVTZ

∞Z

methods

RNN

CASSCF(2,2) SS-MRMPPT(2,2) CASSCF(4,4) SS-MRMPPT(4,4) MRCI(12,10)b MRCI(12,10)b IVO-CASCI/CAS(8,8)c CCSD(T)d CCSD(T)d CCSD(T)d CCSD(T)d G2M(MP2)e CASSCFf CASSCF(2,2) SS-MRMPPT(2,2) CASSCF(4,4) SS-MRMPPT(4,4) MRCI(12,10)b MRCI(12,10)b IVO-CASCI(8,8)c IVO-CASCI(8,8)c CCSD(T)d CCSD(T)d CCSD(T)d CCSD(T)d G2M(MP2)e CASSCFf SR-SDCIg MP3g CCSDg MRSDCIg MRAQCCg SR-SDCIg MP3g CCSDg MRSDCIg MRAQCCg SR-SDCIg MP3g CCSDg CCSD(T)g MRSDCIg MRAQCCg expt.g

1.2370 1.2501 1.2368 1.2484 1.2597 1.2558 1.2388 1.2592 1.2512 1.2481 1.2456 1.2609 1.2660 1.2362 1.2531 1.2579 1.2579 1.2577 1.2563 1.2393 1.2481 1.2643 1.2536 1.2494 1.2468 1.2670 1.2660 1.2435 1.2488 1.2546 1.2655 1.2662 1.2274 1.2510 1.2386 1.2544 1.2556 1.2233 1.2308 1.2354 1.2458 1.2472 1.2487 (1.252, 1.247, 1.2457)

a

e

Bond lengths and bond angles are given in Angstroms (Å) and degrees, respectively. Reference 58. f Reference 54. g Reference 46.

which include the Davidson size-consistency and zero-point energy corrections for further assessment of our findings. Table 3 also contains results of the various single-reference and multireference theories after Shepard et al.46 The tables clearly display that the geometry calculations reported here converge fairly satisfactorily at the SS-MRMPPT/aug-cc-pVDZ level. Tables 3 and 4 display that the geometries of trans, cis, and iso isomers of N2H2 calculated at the SS-MRMPPT level with CAS(2,2) agree well with the best theoretical estimates. In the case of the trans isomer, a good performance can also be seen for the SS-MRMPPT method [even for the CAS(2,2) reference] with respect to the experimental observations. At this point, it should be noted that the CAS(2,2) used in the present work is much smaller than what is used in benchmark calculations on multireference methods reported in Table 3, and thus, SSMRMPT is a good scheme to predict the geometry of various local minima of the PES of the N2H2 system. To visualize the effect of the dimension of the CAS, we have computed the geometries for three conformations implementing the CAS(4,4)

∠HNN

RNH 1.0141 1.0320 1.0148 1.0250 1.0587 1.0481 1.0178 1.0501 1.0360 1.0343 1.0331 1.0369 1.0450 1.0120 1.0271 1.0107 1.0251 1.0426 1.0418 1.0097 1.0152 1.0447 1.0310 1.0294 1.0281 1.0320 1.0450 1.0348 1.0385 1.0439 1.0446 1.0453 1.0186 1.0227 1.050 1.0306 1.0315 1.0169 1.0217 1.0239 1.0270 1.0283 1.0300 (1.028,1.029, 1.0286) b

Reference 59.

c

Reference 60.

112.47 111.71 112.43 111.45 111.5 111.6 112.3 111.48 111.64 111.79 111.88 111.9 104.8 107.48 105.75 106.70 105.74 105.9 106.0 107.3 106.4 104.97 105.73 106.05 106.17 105.1 104.8

106.3 d

Reference 57.

and tabulated the results in Table 3. Clearly, the results show no significant improvement vis-a-vis the CAS(2,2). For this system, as a byproduct of these calculations, we have also presented the trans-cis and trans-iso energy gap, that is, the isomerization energy, (∆E) without considering the zeropoint energy correction to assess the wide applicability of our SS-MRMPPT gradient approach, although computation of ∆E is not the primary focus of the present work. In Table 5, we have listed the calculated energy gap along with the best calculated previously reported theoretical results. From the literature database, it has been found that the CASSCF calculations fail drastically in predictions of relative energies of trans-iso-N2H2 (see ref 54). As shown in the table, the cis-trans energy gap provided by the SS-MRMPPT(2,2)/(4,4) is satisfactorily close to those of IVO-MRMP, MRCI, and CCSD(T) values. On the other hand, for the trans-iso gap, the SS-MRMPPT(4,4) result is in close proximity with those of CCSD(T), MRCI, and MRCI+Q in comparison to the SSMRMPPT(2,2) scheme. To be modest, we may say that the

SS-MRMPPT Applied to Geometry Optimization

J. Phys. Chem. A, Vol. 114, No. 10, 2010 3675

TABLE 4: Results of Equilibrium Geometries for the Ground State of the iso-1,2-Diimine Molecule (N2H2), Using the Methods Described in the Texta basis cc-pVTZ

aug-cc-pVDZ aug-cc-pVTZ aug-cc-pVTZ cc-pVDZ cc-pVTZ cc-pVQZ best cal. 6-31G*

methods

RNN

RNH

∠HNN

CASSCF(2,2) SS-MRMPPT(2,2) CASSCF(4,4) SS-MRMPPT(4,4) MRCI(12,10)b MRCI(12,10)b CASSCF(12,12)c IVO-CASCI(12,12)c CCSD(T)d CCSD(T)d CCSD(T)d CCSD(T)d G2M(MP2)e CASSCFf

1.210 1.216 1.211 1.214 1.2245 1.2166 1.2166 1.2129 1.2280 1.2214 1.2194 1.2172 1.2220 1.2240

1.013 1.036 1.013 1.032 1.0569 10.523 1.0523 1.0362 1.0528 1.0370 1.0351 1.0342 1.2299 1.0380

122.08 123.85 122.82 123.5 124.1 124.2 124.2 123.5 124.46 123.70 123.49 123.49 124.1 124.4

a Bond lengths and bond angles are given in Angstroms (Å) and degrees, respectively. b Reference 59. c Reference 60. d Reference 57. e Reference 58. f Reference 54. g Reference 46.

TABLE 5: Calculated Relative Energy ∆E (in kcal/mol) for N2H2 ∆E

methods

∆Etrans-cis

cis-trans

CASSCF(2,2)/cc-pVTZ SS-MRMPPT(2,2)/cc-pVTZ SS-MRMPPT(4,4)/cc-pVTZ IVO-MRMPa CCSD(T)/cc-pV∞Zb CASSCFb G2M(MP2)/6-31G*b MRCI/aug-cc-pVQZb MRCI+Q/aug-cc-pVQZb CASSCF(2,2)/cc-pVTZ SS-MRMPPT(2,2)/cc-pVTZ SS-MRMPPT(4,4)/cc-pVTZ IVO-MRMPa CCSD(T)/cc-pV∞Zb CASSCFb G2M(MP2)/6-31G*b MRCI/aug-cc-pVQZb MRCI+Q/aug-cc-pVQZb

5.64 5.31 5.91 5.87 5.21 6.68 4.81 5.03 5.05 27.90 26.95 23.32 24.37 24.12 34.59 24.65 24.02 24.04

iso-trans

a

Reference 60. b Reference 59.

performance of the SS-MRMPPT gradient approach is very satisfactory in the case of N2H2. D. Ethylene (C2H4). Next, we consider the ethylene molecule, which plays a fundamental role in the understanding of photoisomerization processes.61 Theoretical calculations on ethylene have a very long tradition.60–64 We have considered the planar C2H4 (D2h) molecule. The ground state (1Ag) of ethylene is a closed-shell singlet state. At the equilibrium geometry, the two carbon p orbitals perpendicular to the molecular plane form bonding π and antibonding π* orbitals. In the ground state, it possesses MR character. The π2 f π*2 excitation is quite important and has a coefficient of -0.213 in the 2RCISD/DZP calculation at the computed equilibrium geometry.63 To get correct geometrical parameters, a high-level treatment for dynamical and nondynamical correlation effect is important for this system. Therefore, this system can be employed to test the performance of the SS-MRMPPT gradient approach. It has been reported65 that σ-π correlation is indeed very important, and it is difficult to recover it when it is absent in the zero-order wave function. In our calculations, the carbon (1s)-like orbitals have been constrained to be doubly occupied. For the this molecule, we have used two D2h symmetry-adapted

active spaces, namely, CAS(2,2) and CAS(4,4), which also allow us to investigate the dependence of the equilibrium bond distance and angles on the reference space size. We used a cc-pVTZ basis and froze the 1s core orbitals of carbon. The optimized geometrical parameters for the C2H4 molecule of our calculations are summarized in Table 6. To get a clear illustration about the potentiality of the SS-MRMPPT gradient method, we compare the results against those obtained with various standard correlated methods.60,61,63,64 As a further test of the robustness of the formalism, we have also tabulated the experimental results.66 We found that the geometrical parameters are only weakly affected by basis set expansion in the case of CCSD and IVO-CASCI methods. It is also worth noting that the main problem with the results of IVO-CASCI, CCSD, MRSDCI, and SF-based methods is that ∠CHH is severely underestimated. It is clear from the computed geometries that they are fairly reproducible between the SS-MRMPPT and computationally expensive CCSD(T)63 levels of theory. Inspection of the data in the table suggests that our calculated geometry is considerably closer to experiment. This appears to validate our computed geometry by any standard. However, a close observation of our results in Table 6 clearly indicates that in order to obtain a reliable value for the ground-state geometry of the C2H4 molecule, it is essential to employ the CAS(2,2). Keeping in mind the experimental results, in this case, the SSMRMPPTs with CAS(2,2) and CAS(2,2) approaches seem to perform almost identically (for different basis sets), the former one being more economical. E. 1,2-Difluoroethene (C2H2F2). Computation of the groundstate geometry of 1,2-difluoroethene is our next example that benefits from a multireference description. The determination of the structure and properties of this system has attracted considerable attention from environmentalists because of its significant role as an air pollutant. This is one of the reasons for regular studies on this type of molecule. Several workers have devoted considerable effort to cis- and trans-dihaloethylenes in order to gain a better understanding of their electronic structure. Despite the obviously greater repulsion of the C-F bond dipoles in the cis configuration than that in the trans configuration, the cis isomer has a lower energy (known as the cis effect (in contrast to chemical intuition)). The anomalous stability of the cis isomer of 1,2-difluoroethylene has been attributed to the formation of mesomeric resonance structure(s). A number of theoretical calculations have been done to explore this effect.67–70 Investigation of the origin of the cis effect is not the aim of this study. The reasons behind the consideration of this molecule is that there are a lot of theoretical studies available which can be used to calibrate our theoretical results. Various theoretical investigations have illustrated that a correct interpretation of the cis effect depends on the accuracy of the geometrical parameters obtained from molecular orbital calculations. Consequently, the computed geometries can be strongly dependent on the calculation level and basis sets used. In our applications, the geometry of 1,2-difluoroethylene has been calculated using augmented correlation consistent basis sets of spd quality (aug-ccpVDZ) and CAS(4,4). In our calculation, the 1s orbitals for each C have been excluded from the correlation treatment. Selected geometrical parameters of 1,2-difluoroethene obtained by the SS-MRMPPT gradient scheme are given in Table 7. To assess the quality of our computed geometrical parameters, the optimized geometrical parameters at various levels of theory are summarized in the same table. The experimental values69 are also given in the same

3676

J. Phys. Chem. A, Vol. 114, No. 10, 2010

Mahapatra et al.

TABLE 6: Results of Equilibrium Geometries for the Ground State of the Planar Ethylene Molecule, Using the Methods Described in the Texta basis cc-pVTZ

aug-cc-pVDZ aug-cc-pVDZ aug-cc-pVTZ DZP

aug-cc-pVDZ cc-pvTZ

∞Z

cc-pVDZ cc-pVTZ best estimated

e

∠HCH

methods

R(CdC)

R(CH)

CASSCF(2,2) SS-MRMPPT/CAS(2,2) CASSCF(4,4) SS-MRMPPT/CAS(4,4) CASSCF(4,4) SS-MRMPPT CCSDb IVO-CASCI/CAS(8,8)b CCSDb IVO-CASCI/CAS(8,8)b VOO-CCDc SF-CISc SF-B3LYPc SF-CIS(D)c SF-ODc MRDC1c MR-CISD/SA-3-CAS(2,2)d SR-SDCIe MP3e CCSD(T) MRSDCIe MRAQCCe SR-SDCIe MP3e CCSD(T)e MRSDCIe MRAQCCe CCSD(T)f CCSD(T)f CCSD(T)f experimentg

1.3355 1.3356 1.3359 1.3361 1.335 1.336 1.348 1.323 1.327 1.314 1.358 1.352 1.339 1.345 1.351 1.328 1.350 1.3198 1.3250 e 1.3330 1.3385 1.3386 1.3168 1.3227 1.3306 1.3348 1.3350 1.35162 1.33713 1.3307 1.334,1.330,1.336

1.0728 1.0815 1.0728 1.0825 1.072 1.082 1.094 1.083 1.078 1.076 1.102 1.077 1.088 1.084 1.088 1.10 1.090 1.0734 1.0758 1.0788 1.0834 1.0834 1.0770 1.0766 1.0797 1.0824 1.0825 1.09841 1.083 18 1.0809 1.081,1.076

a Bond lengths and bond angles are given in Angstroms (Å) and degrees, respectively. Reference 46. f Reference 63. g Reference 66.

b

Reference 60.

c

121.54 121.48 121.35 121.42 121.35 121.41 116.9 116.8 117.1 116.6 116.84 117.1 116.8 116.9 117.1 117. 117.2

121.492 121.453 121.4 121.32,120.95 Reference 64.

d

Reference 61.

TABLE 7: Results of Equilibrium Geometries for the Ground State of the 1,2-Difluoroethylene Molecule, Using the Methods Described in the Texta conformers cis

basis

methods

R(CdC)

R(C-H)

R(C-F)

∠CCH

∠CCF

aug-cc-pVDZ

CASSCF(4,4) SS-MRMPPT(4,4) MP2b MP2b MP2c MP2c B3LYPc B3PW91c B3LYPc B3PW91c MP2d CCSD(T)d experimentd experimentd CASSCF(4,4) SS-MRMPPT(4,4) MP2b MP2b MP2c B3LYPc B3PW91c B3LYPc B3PW91c MP2d CCSD(T)d experimentd experimentd

1.3467 1.3466 1.341 1.339 1.330 1.332 1.327 1.327 1.328 1.328 1.341 1.336 1.324 1.331 1.3463 1.3455 1.340 1.339 1.332 1.327 1.327 1.327 1.327 1.340 1.336 1.329 1.316

1.0755 1.0861 1.089 1.091 1.082 1.082 1.084 1.084 1.083 1.084 1.089 1.083 1.089 1.084 1.0761 1.0862 1.089 1.092 1.082 1.085 1.084 1.084 1.084 1.090 1.085 1.080 1.080

1.3287 1.3557 1.356 1.343 1.343 1.356 1.343 1.383 1.348 1.341 1.356 1.342 1.335 1.335 1.3348 1.3629 1.363 1.346 1.363 1.342 1.383 1.355 1.348 1.363 1.346 1.334 1.352

122.93 123.71 122.8 121.7 123.2 123.6 122.8 122.6 123.4 123.1 122.9 122.4 124.0 121.6 125.75 126.75 125.0 124.2 126.4 124.9 122.6 126.1 125.8 125.8 125.1 129.3 126.3

122.20 121.31 122.3 123.1 122.0 122.0 122.7 122.9 122.5 122.67 122.2 122.4 122.1 123.7 119.48 118.50 119.5 120.5 119.3 120.5 122.9 119.9 120.0 119.5 119.8 119.3 119.2

cc-pVDZ 6-31G* 6-31+G* 6-31G* 6-31+G* aug-cc-pVDZ cc-pVTZ trans

aug-cc-pVDZ cc-pVDZ 6-31+G* 6-31G(d) 6-31+G* aug-cc-pVDZ cc-pVTZ

a

Bond lengths and bond angles are given in Angstroms (Å) and degrees, respectively. b Reference 70. c Reference 68. d Reference 69.

SS-MRMPPT Applied to Geometry Optimization

J. Phys. Chem. A, Vol. 114, No. 10, 2010 3677

TABLE 8: Calculated Relative Energy, ∆Etrans-cis (in kcal/mol) for 1,2-Difluoroethene

a

basis

methods

∆Etrans-cis

aug-cpVDZ

CASSCF SS-MRMPPT MP2a IVO-MRMPa a H3rd v CCSDa CCSD(T)a

0.24 0.53 1.07 2.38 1.19 0.88 0.88

Reference 70. b Reference 68.

table. From previous theoretical studies,67–69 it has been observed that the incorporation of polarization functions in the basis set not only leads to the lowering of the ∠CCH bond angle but also that of the C-F bond length. On the other hand, the inclusion of correlation effects explicitly increases both the CdC and C-H bond lengths. Consequently, to get correct agreement between the computed geometrical parameters and the experimental observation, a balanced and simultaneous inclusion of these effects is very important. In our present study, we observed that there is no significant change in C-H and CdC bond lengths between the cis and trans conformers of the compound. On the other hand, the situation in the C-F bond length is different. From the inspection of the values in Table 7, we find that the C-F bond length is smaller in the case of the cis conformer in comparison to that of its trans one in all of the levels of theory considered here. This is due to the fact that in the cis form, both of the highly electronegative fluorine atoms are in the same plane, and the repulsion between them leads to strengthening of the C-F bond. Our computed values show that all of the bond angles change if one moves from the cis to trans forms. This has also been observed in previous studies.68–70 It is found that the performances of MP2/aug-ccpVDZ and SSMRMPPT/aug-ccpVDZ are almost identical. However, unfortunately, like other single-reference methods, the standard MP2 method gives a very poor performance when quasidegenerate boundary orbitals are encountered. It is important to note that the SS-MRMPPT method works very well in such a case. Satisfactory agreement between the SS-MRMMPT and other theoretical results demonstrates the efficacy of our gradient scheme. The favorable agreement of our predicted geometries with experiment suggests that the computed bond angles and bond lengths for the ground state of C2H2F2 should be reliable. To get a better feeling regarding the performance of our SSMRMPPT method, in addition to performing geometry optimizations for the ground state of the cis and trans isomers of C2H2F2 at equilibrium, we have also reported the energy gap between them (∆Etrans-cis). The results are summarized in Table 8. Here, we also reported the CASSCF results along with the SS-MRMPPT one to illustrate the influence of dynamical correlation (incorporated via the SS-MRMPPT method) on the computed energy gap. Kanakaraju et al.68 showed that the HF and MP2 methods with the 6-31G* basis set predict the trans form of the molecule as the minimum energy conformer. On the other hand, it is observed that the same theories predict the correct stability order of conformers with the addition of diffuse function(s) to the 6-31G* basis set. Thus, the nature of the basis has a strong influence to predict the correct stable conformers, which is compatible with the known fact that the dynamic correlation increases the barrier. Very recently, Chaudhuri et al.70 demonstrated that the stability of the cis isomer with respect to its trans structure, and ∆Etrans-cis increases with the size of the basis set for various methods [such as HF, MP2, IVO-

Figure 2. 1,3-Butadiene (C4H6).

MRMP, H v3rd, CCSD, and CCSD(T)]. In the present paper, we have not tried to delineate the influence of basis sets on the barrier height. Perusal of the presented data shows that the energy barrier spans a wide range of 0.24-2.38 kcal/mol depending on the basis set and method used, implying that the problem has not been settled yet, despite the fact that the electronic structure of the system has been subjected to exhaustive theoretical investigations. Table 8 clearly indicates that cis-1,2-difluoroethene is energetically more stable than the corresponding trans isomer (termed the cis effect) as the ∆Etrans-cis value is positive and the predicted ∆Etrans-cis from the SS-MRMPPT method agrees with the predictions of other advanced theoretical results. Our results compare rather well with the results obtained from CCSD(T) employing the augcc-pVDZ basis set,70 but CCSD(T) takes more time than other perturbative schemes. Very recent results70 obtained with IVOMRMP and Hv methods yield slightly elevated values, ranging from 2.38 to 1.19 kcal/mol. Thus, it is difficult to draw an unambiguous conclusion concerning the goodness of the SSMRMPPT value of the barrier height. The analysis of our computed geometries (including the cis-trans energy gap) of C2H2F2 along with experimental and other computed results reveals in a convincing way the efficacy of our SS-MRMPPT gradient approach. F. 1,3-Butadiene (C4H6). In our next set of test calculations, we have studied the ground state of the cis and trans conformers of the 1,3-butadiene (C4H6) system (see Figure 2). The trans form has the lowest energy, and the cis form is a rotational transition state. As the smallest conjugated polyene, the electronic structure of the ground and excited electronic states of 1,3-butadiene has received a great deal of experimental71 and theoretical72–76 attention for many decades. Conformational studies related to the internal rotation of 1,3-butadiene is one subject that has been widely investigated.73–76 The fully optimized geometries of stationary points cis, trans, and iso on the potential energy curves have been computed via the SSMRMPPT gradient approach. These computed results have been compared with experiment and with other previously published correlated methods. From previous theoretical study, it has been found that electron correlation has a pronounced effect on the computation of the geometry and properties of 1,3-butadiene. It thus represents an excellent probing ground for testing the efficiency of the numerically oriented gradient-based SSMRMPPT method and various MR-based approaches in their ability to account for correlation effects. We have used CAS(4,4) and the cc-pVDZ basis set in our calculations, and we froze 1s orbitals on carbon in all correlated calculations. In addition, calculations have been performed using the basis set 6-311G** as in refs 74 and 76, so that one can compare our results with the findings of Gong et al., Xiao,74 and Chattopadhyay et al.76 Table 9 displays our computed ground-state optimized geometrical parameters of the cis and

3678

J. Phys. Chem. A, Vol. 114, No. 10, 2010

Mahapatra et al.

TABLE 9: Results of Equilibrium Geometries for the Ground State of the 1,3-Butadiene Molecule, Using the Methods Described in the Texta conformers trans

basis

methods

R(Cb-Cb)

R(CbdCb)

∠Ca-Cb-Cb

cc-pVDZ 6-311G**

SS-MRMPPT(4,4) SS-MRMPPT(4,4) MP2b CCSDc CASPT2d B3LYPd MP2d experimente SS-MRMPPT(4,4) SS-MRMPPT(4,4) MP2b CCSDc CASPT2d B3LYPd MP2d experimente

1.4641 1.4535 1.460 1.469 1.454 1.456 1.456 1.476 1.4762 1.4717 1.472 1.483 1.468 1.470 1.470 1.472

1.3547 1.3486 1.345 1.344 1.348 1.339 1.343 1.337 1.3554 1.3474 1.344 1.343 1.351 1.339 1.342 1.349

122.78 123.67 123.6 123.7 123.6 124.3 123.7 122.9 123.15 125.26 123.8 126.4 126.7 127.3 126.5 124.4

6-31G(d) 6-31G(d) 6-31G(d) cis

cc-pVDZ 6-311G** 6-31G(d) 6-31G(d) 6-31G(d)

e

a Bond lengths and bond angles are given in Angstroms (Å) and degrees, respectively. Reference 71.

trans conformers of 1,3-butadiene along with experimental results. In addition to the MP2 and CCSD results, a comparison of the SS-MRMPPT result with the data obtained with the CASPT2 and B3LYP methods has also been presented in the same table. The results, which are presented in Table 9, are rather encouraging. Comparing the bond lengths and angles listed in Table 9, one sees that the geometries of the cis and trans isomer of 1,3-butadiene agree well with experiment. The calculated results for the cis and trans conformers are in close proximity with other high-level theoretical calculations. From this comparison, it is seen that our state-specific calculations show the geometrical parameters to be very close to the CASPT2 ones. The close agreement between SR-based and MR-based methods can be attributed to the reasonable SR nature of the stable conformers of the system considered. We now focus our attention on the value of the barrier height between the cis and trans isomers of 1,3-butadiene. As we have already mentioned, the focus of this study is not to provide the benchmark calculations for the barrier height between the cis and trans isomers by using SS-MRMPPT. In our calculations, we have observed that there is a substantial energy difference between the cis and the trans conformers for the cc-pVDZ basis [3.23 and 3.60 kcal/mol for CASSCF(4,4) and SS-MRMPPT(4,4), respectively, in the case of 6-311G** basis; these are 3.23 kcal/mol with CASSCF and 3.85 kcal/mol with SSMRMPPT]. Considering that the energy difference is so small, it is instructive to perform calculations involving complete geometry optimization at a higher level of theory to observe whether more definitive predictions can be made to validate our results. It should be mentioned that there is a large scatter of the theoretical estimates of the energy barrier height in the literature depending on the method employed, implying that the problem is not settled as of yet, which should be rectified by using state of the art methodology and computations. G. Cyclobutadiene (C4H4). In order to further explore its range of applicability and level of performance, we apply the present gradient method to the correlation problem of cyclobutadiene. Geometry optimization of rectangular cyclobutadiene (CB) is a rather challenging example.14,28,77–82 Cyclobutadiene has been a challenge for both theoretical28,78–82 and experimental83 chemists for many decades, not due to its seemingly simple electronic structure but rather due to its exhibiting of some very interesting features. Neutral cyclobutadiene is a very unstable,

b

Reference 74.

c

Reference 76.

d

Reference 75.

nonaromatic species that appears to have a rectangular structure (D2h) composed of alternating single and double bonds. The earlier experimental observations indicate clearly that the lifetime of cyclobutadiene is certainly so short that isolation appears hopeless. As per a recent study, cyclobutadienes can be obtained as monomers by the use of the matrix isolation technique.83 Thus, the compound is amenable to study. The ground state of cyclobutadiene is found to be a closedshell singlet, implying that it is well-described by singlereference methods. However, various theoretical studies demonstrate the fact that the MR character of the rectangular conformation in its ground state has noticeable effect on the equilibrium geometry. CB has, therefore, aroused a lot of interest of theoreticians to assess the performance of various SR- and MR-based methods in their ability to describe its electronic structure. We use the cc-pVDZ basis set and freeze 1s orbitals on carbon in our calculations, where the reference function is of the CAS(4,4) type. In addition, we also use the aug-cc-pVDZ basis. In order to demonstrate the usefulness of the SS-MRMPPT gradient method, the repercussions of our findings with regard to recent theoretical results have been discussed. The optimized geometries as obtained in the SS-MRMPPT numerical gradient calculations with the cc-pVDZ and aug-cc-pVDZ basis sets are given in Table 10, where it is compared to the geometries obtained by other theoretical methods. There are numerous theoretically measured structures of CB, and all have not been commented upon here. It is worthwhile to note that experimental geometrical data for this system are not available (as per the current state of our knowledge) due to its elusive nature. Table 10 also gives information on the dependence of the geometries on the basis set size. Again, the findings given in the table indicate the shortcomings of the SR-based treatments and serve as an argument for the use of the present MR-based approach. The results presented in Table 10 demonstrate that SR-based approaches are not able to provide reliable description of the equilibrium geometry for the cyclobutadiene system. The results obtained in this paper by the SS-MRMPPT gradient method are consistent with the results of previous theoretical calculations.28,78–82 The SS-MRMPPT geometries are very close to the results obtained by the computationally expensive MR-BWCC and Mk-MRCC methods (using the cc-pVDZ basis set) by means of which one can judge the goodness of the SS-MRMPPT

SS-MRMPPT Applied to Geometry Optimization

J. Phys. Chem. A, Vol. 114, No. 10, 2010 3679

TABLE 10: Optimal Geometries and Energies of the Rectangular Cyclobutadiene Using CAS(4,4)a basis cc-pVDZ aug-cc-pVDZ cc-pVDZ

cc-pVDZ

cc-pVTZ(mixed) cc-pVTZ 3s2p1d/2s 3s2p1d/2s cc-pVDZ cc-pvTZ aug-cc-pVDZ cc-pVDZ

method

RC-C

RCdC

RC-H

HCCa

CASSCF(4,4) SS-MRMPPT(4,4) CASSCF(4,4) SS-MRMPPT(4,4) SCFc MP2c MP4c CCSDc CCSD(T)c MR-BWCCSD(n.c.)d MR-BWCCSD(a.c.)d MR-BWCCSD(i.c.)d MR-MkCCSDd MR-BWCCSD(T)(n.c.)d MR-BWCCSD(T)(a.c.)d MR-BWCCSD(T)(i.c.)d MR-MkCCSD(T)d EOM-SF-CCSDe CCSD(T)e EOM-SF-CCSDe MR-CCSDe MRAQCC/SA-4-CASSCFf MRAQCC/SA-4-CASSCFf MRAQCC/SA-4-CASSCFf GVVPT2g

1.5513 1.5674 1.5571 1.5756 1.5686 1.5764 1.5837 1.5795 1.5831 1.577 1.570 1.570 1.573 1.583 1.575 1.574 1.579 1.551 1.566 1.566 1.570 1.573 1.562 1.574 1.574

1.3584 1.3753 1.3514 1.3636 1.3230 1.3545 1.3598 1.3528 1.3644 1.357 1.364 1.365 1.361 1.363 1.363 1.372 1.366 1.344 1.343 1.370 1.367 1.367 1.349 1.368 1.372

1.0786 1.0892 1.07622 1.0872 1.0795 1.0906 1.0930 1.0913 1.0953 1.093 1.093 1.092 1.093 1.095 1.095 1.095 1.095 1.081 1.074 1.091 1.103 1.093 1.077 1.091 1.095

134.87 134.88 134.90 134.8 134.95 135.03 134.95 134.95 134.94 134.86 134.86 134.86 134.87 134.89 134.87 134.88 134.88 134.9 134.9 134.8 134.7 134.9 134.9 135.0 134.88

a Bond lengths and bond angles are given in Angstroms (Å) and degrees, respectively. b Exterior angle between the C-H and the longer C-C bonds. c Reference 79. d Reference 81. e References 28 and 78. f Reference 80. g Reference 14.

Figure 3. 2,6-Pyridynium cation (C5NH+ 4 ).

gradient method. We note that CCSD(T)/cc-pVTZ calculations of Levchenko and Krylov78 and the MR-AQCC/SA-4-CASSCF values are in good agreement with our findings for the bond lengths and angle. As that of the other theoretical observations,81 we observed that there is a decreasing trend in the CdC and C-H bond lengths as well as the ∠HCC with the increase in size of the basis set, whereas the opposite trend is displayed in the case of the single C-C bond length. We believe that our estimate of the rectangular structure characterized by the bond lengths and bond angle is reliable too. The distance together with the angle generated by the SS-MRMPPT gradient approach is also a good indicator of the usefulness of the method for the description of the multireference character. H. 2,6-Pyridynium Cation (C5NH4+). Finally, for the sake of diversity, we also consider a charged species, such as the 2,6-pyridynium cation, C5NH4+ (that possesses pronounced diradical character), which has been extensively studied in the past. Several authors have performed density functional theory and wave-function-based calculations to explore its electronic structure.23,84–86 Among all six possible isomers, here, we focus

only on the 2,6 isomer (see Figure 3). Formally, it can be generated by replacing one C-H group in benzyne with the N-H+ group. It represents a consistently challenging subject for electronic structure calculations, having a singlet ground state that is poorly described by SR methods, thus representing a good testing ground for MR-type approaches. Consequently, the prediction of equilibrium geometrical parameters of this system is a natural choice to illustrate the efficacy of our SSMRMPPT gradient method. Very recently, this system has been investigated by Prochnow et al.23 using various basis sets via the Mk-MRCC/SS-MRCC (SS-MRMPPT was originally developed by Mukherjee and co-workers87 as a perturbative approximant to the rigorously size-extensive SS-MRCC theory) analytical gradient approach in order to show the potentiality of the newly developed gradient scheme. Mk-MRCC data of geometrical parameters published by Prochnow et al.23 on this system have been used as a benchmark for testing and calibrating the SS-MRMPPT gradient method. For an easier comparison, we use the same basis as that in ref 23. As a basis set, the correlation-consistent-core-valence polarized set, ccpCVDZ, of Woon and Dunning has been used.88 It should be mentioned that the experimental structure of the C5NH+ 4 system has yet to be offered. In this paper, we apply a relatively small CAS as the reference, with two active electrons on two active orbitals (a1 (HOMO) and b1 (LUMO)), denoted as CAS(2,2). These two closed-shell determinants are used as references, with the orbitals taken from CASSCF calculations. The optimized structural parameters of the ground state of C5NH4+ obtained at the SS-MRMPPT level of theory with the cc-pCVDZ basis are presented in Table 11. It also summarizes some of the previously reported geometrical parameters calculated by SR, MR, as well as DFT methods. Comparative studies shown in this table illustrate the efficacy of the SS-MRMPPT gradient approach to handle difficult systems like C5NH4+ characterized by a significant quasidegeneracy of the reference state. The HF-SCF methods provide the incorrect closed-shell

3680

J. Phys. Chem. A, Vol. 114, No. 10, 2010

Mahapatra et al.

TABLE 11: Structural Parameters of the Ground State of the 2,6-Pyridynium Cation (C5NH4+)a parameters

SCFb

TCSCFb

CCSDb

CCSD(T)b

Mk-MRCCSDb

B3LYP

BPW91c

CASCF(8,8)c

SS-MRMPPT

N1-C2 C2-C3 C3-C4 C2-C6 N1-H C3-H C4-H ∠C2-N1-C6 ∠N1-C2-C3 ∠C2-C3-C4 ∠C3-C4-C5 ∠C2-C3-H

1.2993 1.3592 1.4290 1.4959 1.0170 1.0763 1.0844 70.29 164.11 103.90 113.69 128.65

1.3297 1.3587 1.3979 2.2299 1.0060 1.0796 1.0820 114.09 126.75 116.67 119.07 120.47

1.3370 1.3738 1.4127 2.1633 1.0215 1.0929 1.0958 107.99 130.81 117.30 115.79 119.78

1.3502 1.3814 1.4174 2.2387 1.0236 1.0949 1.0973 112.12 127.65 117.72 117.14 119.84

1.3418 1.3734 1.4127 2.1983 1.0214 1.0929 1.0949 111.26 128.43 117.43 117.02 120.01

1.3291 1.3602 1.4129 2.1015 1.0242 1.0933 1.0955 104.47 133.59 116.88 114.59 120.54

1.347 1.371

1.332 1.375

2.242

2.235

1.3415 1.3784 1.4103 2.2053 1.0243 1.0939 1.0955 110.56 128.66 117.95 116.21 119.38

a All of the calculations are carried out with cc-pCVDZ (128 CGTOs). Bond lengths and bond angles are given in Angstroms (Å) and degrees, respectively. SS-MRMPPT geometry optimizations are carried with a (2 × 2) CAS constructed from a1(HOMO) and b1(LUMO) orbitals (two-configurations SCF). b Reference 23. c Reference 84.

bicyclic structure (as is evident from the C2-C6 distance, ∼1.5 Å).23,86 It stems from an artifact associated with this level of theory, it being a single-determinant picture. However, with the inclusion of nondynamical correlation, the bicyclic structures inevitably open to monocyclic ones, as is evident from the TCSCF caculations by Prochno et al.23 Akin to the CCSD, CCSD(T), and Mk-MRCCSD methods,23 the SS-MRMPPT gradient too provides the monocyclic structure with a C2-C6 distance of 2.2053 Å for this system. Various sophisticated theoretical studies on this system reveal that the system bears a monocyclic structure in the ground state. We now focus on the geometrical parameters (the C2-C6 length and the ∠C2-N1-C6) generated using several methods reported in the table. This is needed since the extent of the MR nature of the system is best diagnosed by these geometrical parameters. The smallest ∠C2-N1-C6 obtained from CCSD is 107.99°, and CCSD(T), Mk-MRCCSD, and SS-MRMPPT favor wide values of 112.12, 111.26, and 110.56°, respectively. This is in complete agreement with the value of the C2-C6 length where the CCSD treatment yields a distance of 2.163 Å, and the CCSD(T), MkMRCCSD, and the SS-MRMPPT methods give values of 2.204, 2.1983, and 2.205 Å, respectively. The table also points out the fact that the incorporation of the triple excitation influences the performance of the SRCC method significantly. These observations call for a genuine MR description for such systems. The most significant observation from the data in Table 11 is that the SS-MRMPPT gradient method performs very close to Mk-MRCCSD values as compared to the proximity of values obtained from the CCSD and CCSD(T) results. Overall, the SSMRMPPT gradient exhibits very satisfactory performance even at CAS(2,2) for a very complicated system like C5NH4+ by predicting a correct equilibrium geometry (open-shell monocyclic) and geometrical parameters at a fraction of the computational cost of SR- or MR-based CCSD methods. Thus, the SS-MRMPPT level of calculations may be expected to be of considerable utility in investigations of aryne biradicals. The foregoing description demonstrated the viability of the SS-MRMPPT gradient method by several rather diverse applications. Our results presented in the various tables reveal that the inclusion or not of polarization functions and the proper incorporation or not of electronic correlation in the ab initio computations are the two main issues explaining the total variation in the results for the geometrical parameters. Thus, to obtain predicted geometries in good agreement with the experimental ones, a democratic and synchronous incorporation of these effects is imperative. From the observations, we have noted that the geometrical parameters obtained by the numerical-gradient-based SS-MRM-

PPT method compare quite well with the experimental values. It produces satisfactory results (despite our use of the conceptually minimal choice of the reference space) that are superior or very close to those obtained with other highly sophisticated methods. The SS-MRMPPT approach represents a very promising approach to the many-electron correlation problem since it enables a high-level treatment of both dynamical and nondynamical correlation effects in an intruder-free and size-effective manner for a wide class of closed- and open-shell systems, enjoying at the same time a very favorable scaling of its computational cost with the size of the problem. Thus, the present method has the potential to generate accurate bond lengths and bond angles for different molecular systems. The present work suggests that the SS-MRMPPT gradient approach can be viewed as an efficient and useful companion of the complete Mk-MRCC gradient method.23 Finally, though convenient, the numerical gradient approach is not free from objection(s) as it often gives rise to problems of numerical stability. In the numerical gradient approach, much higher accuracy of the computed energies is needed for getting stable results than in customary energy calculations. As a result, many more iterations are needed in all iterative procedures during the course of computation, which makes the scheme numerically expensive in comparison to the corresponding analytical one. Although the above-presented results are very encouraging, much work remains to be done in order to achieve its routine exploitation as it is not able to predict satisfactorily the geometrical parameters in the situations of highly stretched or distorted geometries. These situations illustrate that the SSMRMPPT numerical gradient method has its limitations. To obviate this issue, there is much need for improved numerical algorithms, both in terms of efficiency and flexibility, as well as their improved convergence characteristics. The development of an analytical gradient approach may ameliorate the situations. It is worth mentioning that the numerical gradient approaches are inherently less efficient than the corresponding analytic ones. It is hoped that the present work will provide some impetus for the development of a more general analytical gradient approach. IV. Concluding Remarks It is now well-known that the key to a successful description of geometrical parameters (such as bond length and angle) is the accurate and balanced treatment of dynamical and nondynamical correlation effects. Balancing dynamical and nondynamical correlation effects in molecular systems in a proper manner is one of the main advantages of the size-consistent second-order SS-MRMPPT method. In this paper, we have

SS-MRMPPT Applied to Geometry Optimization explored the implementation of a numerically oriented gradientbased SS-MRMPPT procedure that could be applied to at least some of the most frequently used multireference systems [such as H2O, O3, N2H2, C2H4, C2H2F2, butadiene (C4H6), cyclobutadiene (C4H4), and the 2,6-pyridynium cation (C5NH4+))] by reporting test calculations for the geometrical parameters such as bond lengths and angles. The numerical first derivative SSMRMPPT procedure calculates energy gradients at the multireference and single-state second-order perturbation theory levels by numerical differentiation. Comparison of our results on the various systems with literature theoretical results continues to support the usefulness of the SS-MRMPPT gradient method. We have seen that even for the states having a strong multireference character, as in the case of ozone, the SSMRMPPT method s relying on a minimal active space s provides a highly satisfactory description of geometries that is much more effective than that provided by the CCSD/CCSD(T) methods or, in fact, the MR-based methods exploiting large reference spaces. Our results have been compared to those obtained from single-reference-based calculations, and thereby, we have illustrated the utility of the multireference-based method. Our numerical analysis depicts that the SS-MRMPPT gradient approach (with the conceptually minimal choice of the reference space) works satisfactorily. The results obtained from the SS-MRMPPT gradient scheme agree within chemical accuracy with the experimental estimates whenever available. Thus, one can say that the numerical-gradient-based SSMRMPPT is a method applicable for satisfactory calculation of geometrical parameters of various systems. The present applications of the SS-MRMPPT gradient procedure lead to good hope for its wide applicability in the near future. However, we do not claim that the SS-MRMPPT numerical gradient scheme is a very elegant method for general use to provide geometrical parameters. Finally, the SS-MRMPPT gradient method has enough room for further development and wide applications. Acknowledgment. The authors gratefully acknowledge the Department of Science and Technology, India, for financial support of this research under Grant No. SR/S1/PC-32/2005. References and Notes (1) Pulay, P.; Meyer, W. J. Mol. Spectrosc. 1971, 40, 59. (2) Yamaguchim, Y.; Osamura, Y.; Goddard, J. D.; Schaefer, H. F., III. A New Dimension in Quantum Chemistry: Analytic DeriVatiVe Methods in Ab-initio Molecular Electronic Structure Theory; Oxford University Press: New York, 1994. (3) Freed, K. F. In Many-Body Methods in Quantum Chemistry; Keldor, U., Ed.; Springer-Verlag, Berlin, Germany, 1988. (4) Crawford, T. D.; Schaefer, H. F., III. ReV. Comput. Chem. 2000, 14, 33. (5) Pahari, D.; Chattopadhyay, S.; Das, S.; Mukherjee, D.; Mahapatra, U. S. In Theory and Applications of Computational Chemistry: The First 40 Years; Dytkstra, C. E., Frenking, G., Kim, K. S., Scuseria, G. E., Eds.; Elsevier: Amsterdam, The Netherlands, 2005; p 581. (6) Piecuch, P.; Woch, M.; Lodriguito, M.; Gour J. R. In Progress in Theoretical Chemistry and Physics; Recent Advances in the Theory of Chemical and Physical Systems; Julien, J.-P., Maruani, J., Mayou, D., Wilson, S., Delgado-Barrio, G., Eds; Springer: Berlin, Germany, 2006; Vol. 15, pp 45-106. (7) Krylov, A. I. Acc. Chem. Res. 2006, 39, 83. (8) Andersson, K.; Malmqvist, P. Å; Roos, B. O.; Sadlej, A. J.; Wolinski, K. J. Chem. Phys. 1990, 94, 5483. (9) Hirao, K. Chem. Phys. Lett. 1993, 59, 201. (10) (a) Evangelisti, S.; Daudey, J. P.; Malrieu, J.-P. Chem. Phys. 1983, 75, 91. (b) Cimiraglia, R.; Persico, M. J. Comput. Chem. 1987, 8, 39. (11) Mahapatra, U. S.; Datta, B.; Mukherjee, D. J. Phys. Chem. 1999, 103, 1822. (12) Kozlowski, P. M.; Davidson, E. R. J. Chem. Phys. 1994, 100, 3672.

J. Phys. Chem. A, Vol. 114, No. 10, 2010 3681 ´ .; Surja´n, P. R. J. Chem. Phys. 2003, (13) (a) Rolik, Z.; Szabados, A ´ .; Ko¨halmi, D Ann. 119, 1922. (b) Surja´n, P. R.; Rolik, Z.; Szabados, A Phys. (Leipzig) 2004, 13, 223. (14) (a) Hoffmann, M. R. Chem. Phys. Lett. 1992, 195, 127. (b) Khait, Y. G.; Hoffmann, M. R. J. Chem. Phys. 1998, 108, 8317. (c) Jiang, W.; Khait, Y. G.; Hoffmann, M. R. J. Phys. Chem. A 2009, 113, 4374. (15) Sheppard, M. G.; Freed, K. J. Chem. Phys. 1981, 75, 4507. (16) Nakano, H. J. Chem. Phys. 1993, 99, 7983. (17) (a) Brandow, B. H. ReV. Mod. Phys. 1967, 39, 771. (b) Lindgren, I. Int. J. Quantum Chem. 1978, S12, 33. (c) Hose, G.; Kaldor, U. J. Phys. B 1981, 12, 3827. (18) (a) Rintelman, J. M.; Adamovic, I.; Varganov, S.; Gordon, M. S. J. Chem. Phys. 2005, 122, 044105. (b) Azizi, Z.; Roos, B. O.; Veryazova, V. Phys. Chem. Chem. Phys. 2006, 8, 2727. (19) Chaudhuri, R. K.; Freed, K. F.; Hose, G.; Piecuch, P.; Kowalski, K.; Wloch, M.; Chattopadhyay, S.; Mukherjee, D.; Rolik, Z.; Szabados, A.; Toth, G.; Surjan, P. R. J. Chem. Phys. 2005, 122, 134105. (20) Hoffmann, M. R.; Datta, D.; Das, S.; Mukherjee, D.; Szabados, ´ .; Rolik, Z.; Surja´n, P. R. J. Chem. Phys. 2009, 131, 204104. A (21) (a) Mahapatra, U. S.; Chattopadhyay, S.; Chaudhuri, R. K. J. Chem. Phys. 2008, 129, 024108. (b) Mahapatra, U. S.; Chattopadhyay, S.; Chaudhuri, R. K. J. Chem. Phys. 2009, 13, 014101. (22) Evangelista, F. A.; Simmonett, A. C.; Schaefer, H. F., III; Mukherjee, D.; Allen, W. D. Phys. Chem. Chem. Phys. 2009, 11, 4728. (23) Prochnow, E.; Evangelista, F. A.; Schaefer, H. F., III; Allen, W. D.; Gauss, J. J. Chem. Phys. 2009, 131, 064109. (24) Schmidt, M. W.; Baldridge, K. K.; Boatz, J. A.; Elbert, S. T.; Gordon, M. S.; Jensen, J. H.; Koseki, S.; Mastunaga, N.; Nguyen, K. A.; Su, S.; Windus, T. L.; Dupuis, M.; Montgomery, J. A. J. Comput. Chem. 1993, 14, 1347. (25) Shepard, R. The analytic gradient method for configuration interaction wave functions. In Modern Electronic Structure Theory; Advanced Series in Physical Chemistry; Yarkony, D. R., Ed.;World Scientific: Singapore, 1995; Vol. 2, pp 345-345. (26) Szalay, P. G. Int. J. Quantum Chem. 1995, 55, 151. (27) (a) Ajitha, D.; Pal, S. J. Chem. Phys. 1999, 111, 3832. (b) Shamasundar, K. R.; Pal, S. J. Chem. Phys. 2001, 115, 1979. (28) Levchenko, S. V.; Wang, T.; Krylov, A. I. J. Chem. Phys. 2005, 122, 224106. (29) Stevens, J. E.; Chaudhuri, R. K.; Freed, K. F. J. Chem. Phys. 1996, 105, 8754. (30) Chaudhuri, R. K.; Freed, K. F. J. Chem. Phys. 2007, 126, 114103. (31) Nakano, H.; Hirao, K.; Gordon, M. S. J. Chem. Phys. 1998, 108, 5660. (32) Dudley, T. J.; Khait, Y. G.; Hoffmann, M. R. J. Chem. Phys. 2003, 119, 651. (33) Celani, P.; Werner, H. J. J. Chem. Phys. 2003, 119, 5044. (34) (a) Demel, O.; Pittner, J. J. Chem. Phys. 2006, 124, 144112. (b) Pittner, J.; Sˇmydke, J. J. Chem. Phys. 2007, 127, 114103. (35) (a) Handy, N. C.; Schaefer, H. F., III. J. Chem. Phys. 1984, 81, 5031. (b) Bartlett, R. J.; Stanton, J. F.; Watts, J. D. In AdVances in Molecular Vibrations and Collision Dynamics; Bowman, J., Ed.; Elsevier: New York, 1991; Vol. 1, p 139. (c) Gauss, J.; Cremer, D. AdV. Quantum Chem. 1992, 23, 205. (d) Nakajima, T.; Nakatsuji, H. Chem. Phys. Lett. 1997, 280, 79. (36) Krylov, A. I. Chem. Phys. Lett. 2001, 338, 375. (b) Krylov, A. I. Chem. Phys. Lett. 2001, 350, 522. (37) Jeziorski, B.; Monkhorst, H. J. Phys. ReV. 1981, 24, 1668. (38) Malrieu, J. -P.; Heully, J.-L.; Zaitevskii, A. Theor. Chim. Acta 1995, 90, 167. (39) (a) Hanrath, M. J. Chem. Phys. 2005, 123, 84102. (b) Engels-Putzka, A.; Hanrath, M. Mol. Phys. 2009, 107, 143. (40) (a) Feller, D. J. Comput. Chem. 1996, 17, 1571. (b) Schuchardt, K. L.; Didier, B. T.; Elsethagen, T.; Sun, L.; Gurumoorthi, V.; Chase, J.; Li, J.; Windus, T. L. J. Chem. Inf. Model. 2007, 47, 1045. see: www.emsl.pnl.gov/ forms/basisform.html. (41) Bauschlicher, C. W.; Taylor, P. R. J. Chem. Phys. 1986, 85, 2779. (42) Szalay, P. G.; Bartlett, R. J. J. Chem. Phys. 1995, 103, 3600. (43) King, R. A.; Sherrill, C. D.; Schaefer, H. F., III. Spectrochim. Acta, Part A 1997, 53, 1163. (44) Sherrill, C. D.; Krylov, A. I.; Byrd, E. F. C.; Head-Gordon, M. J. Chem. Phys. 1998, 109, 4171. (45) Li, X.; Paldus, J. J. Chem. Phys. 1999, 110, 2844. (46) Shepard, R.; Kedziora, G. S.; Lischka, H.; Shavitt, I.; Mu¨ller, T.; Szalay, P. G.; Ka´llay, M.; Seth, M. Chem. Phys. 2008, 349, 37. (47) (a) Varandas, A. J. C. J. Chem. Phys. 1996, 105, 3524. (b) Piecuch, P.; Woch, M.; Varandas, A. J. C. Theor. Chem. Acc. 2008, 120, 59. (48) Cameron, M. R.; Kable, S. H.; Bacskay, G. B. J. Chem. Phys. 1995, 103, 4476. (49) Ka´llay, M.; Gauss, J.; Szalay, P. G. J. Chem. Phys. 2003, 119, 2991. (50) Borowski, P.; Anderson, K.; Malmqvist, P. A.; Roos, B. O. J. Chem. Phys. 1992, 97, 5568. Roos, B. O. In Theory and Applications of

3682

J. Phys. Chem. A, Vol. 114, No. 10, 2010

Computational Chemistry: The First 40 Years; Dytkstra, C. E., Frenking, G., Kim, K. S., Scuseria, G. E., Eds.; Elsevier: Amsterdam, The Netherlands, 2005; p 725. (51) Tsuneda, T.; Nakano, H.; Hirao, K. J. Chem. Phys. 1995, 103, 6520. (52) (a) Ljubic´, I.; Sabljic´, A. J. Phys. Chem. A 2002, 106, 4745. (b) Ljubic´, I.; Sabljic´, A. Chem. Phys. Lett. 2004, 214, 385. (53) Tanaka, T.; Morino, Y. J. Mol. Spectrosc. 1970, 33, 538. (54) Jensen, H. J. A.; Jorgensen, P.; Helgaker, T. J. Chem. Phys. 1986, 85, 3917. (55) (a) Whitelegg, D.; Wooley, R. G. J. Mol. Struct.: THEOCHEM 1990, 209, 23. (b) Kim, K.; Shavitt, I.; Del Bene, J. E. J. Chem. Phys. 1992, 96, 7573. (56) Mach, P.; Ma´sˇsik, J.; Urban, J.; Hubacˇ, I. Mol. Phys. 1998, 94, 173. (57) Martin, J. M. L.; Taylor, P. R. Mol. Phys. 1999, 96, 681. (58) Hwang, D. Y.; Mebel, A. M. J. Phys. Chem. A 2003, 107, 2865. (59) Biczysko, M.; Poveda, L. A.; Varandas, A. J. C. Chem. Phys. Lett. 2006, 424, 46. (60) Chaudhuri, R. K.; Freed, K. F.; Chattopadhyay, S.; Mahapatra, U. S. J. Chem. Phys. 2008, 128, 144304. (61) Barbatti, M.; Paier, J.; Lischka, H. J. Chem. Phys. 2004, 121, 11614. (62) (a) Freed, K. F. Chem. Phys. Lett. 1968, 2, 255. (b) Chaudhuri, R. K.; Mudholkar, A.; Freed, K. F.; Martin, C. H.; Sun, H. J. Chem. Phys. 1997, 106, 9252. (63) Martin Jan, M. L.; Lee, T. J.; Taylor, P. R.; Franois, J.-P. J. Chem. Phys. 1995, 103, 2589. (64) (a) Krylov, A. I.; Sherrill, C. D.; Byrd, E. F. C.; Head-Gordon, M. J. Chem. Phys. 1998, 109, 10669. (b) Shao, Y.; Head-Gordon, M.; Krylov, A. I. J. Chem. Phys. 2003, 118, 4807. (65) Davidson, E. R. J. Phys. Chem. 1996, 100, 6161. (66) (a) Herzberg, G. In Electronic Spectra and Electronic Structure of Polyatomic Molecules; Von Nostrand Reinhold: New York, 1966. (b) Duncan, J. L. Mol. Phys. 1974, 28, 1177. (c) Kuchitsu, K. J. Chem. Phys. 1966, 44, 906. (d) Clabo, D. A.; Allen, W. D.; Remington, R. B.; Yamaguchi, Y.; Schaefer, H. F. Chem. Phys. 1988, 123, 187. (67) (a) Craig, N. C.; Abiog, O. P.; Hu, B.; Stone, S. C.; Lafferty, W. J.; Xu, L. J. Phys. Chem. 1996, 110, 5310. (b) Dixon, D. A.; Smart, B. E.; Fukunaga, T. Chem. Phys. Lett. 1986, 125, 447. (c) Yamamoto, T.; Kaneno, D.; Tomoda, S. Chem. Lett. 2005, 34, 1190.

Mahapatra et al. (68) Kanakaraju, R.; Senthilkumar, K.; Kolandaivel, P. J. Mol. Struct.: THEOCHEM 2002, 589, 95. (69) Bosco, J.; da Silva, P.; Ramos, M. N. J. Braz. Chem. Soc. 2004, 15, 43. (70) Chaudhuri, R. K.; Hammond, J. R.; Freed, K. F.; Chattopadhyay, S.; Mahapatra, U. S. J. Chem. Phys. 2008, 129, 064101. (71) (a) Aston, J. C.; Szasz, C.; Wooley, H. W.; Brickwedde, F. G. J. Chem. Phys. 1946, 14, 67. (b) Skaarup, S.; Boggs, J. E.; Skancke, P. N. Tetrahedron 1976, 32, 1179. (72) (a) Brink, M.; Jonson, H.; Ottosson, C.-H. J. Phys. Chem. A 1998, 102, 6513. (b) Orlandi, G.; Zerbetto, F.; Zgierski, M. Z. Chem. ReV. 1991, 91, 867. (73) (a) Skaarup, S.; Boggs, J. E.; Skancke, P. N. Tetrahedron 1976, 32, 1179. (b) Alberts, I. L.; Schaefer, H. F., III. Chem. Phys. Lett. 1989, 161, 375. (c) Bock, C. W.; Panchenko, Y. N. J. Mol. Struct. 1989, 187, 69. (b) Wiberg, K. B.; Rosenberg, R. E. J. Am. Chem. Soc. 1990, 112, 1509. (74) Gong, X.; Xiao, H. Int. J. Quantum Chem. 1998, 69, 659. (75) Page, C.; Olivucci, M. J. Comput. Chem. 2003, 24, 298. (76) Chattopadhyay, S.; Chaudhuri, R.; Mahapatra, U. S. J. Chem. Phys. 2008, 129, 244108. (77) Balkova´, A.; Bartlett, R. J. J. Chem. Phys. 1994, 101, 8972. (78) Levchenko, S. V.; Krylov, A. I. J. Chem. Phys. 2004, 120, 175. (79) Sancho-Garcı´a, J. C.; Pittner, J.; Cˇa´rsky, P.; Hubacˇ, I. J. Chem. Phys. 2000, 112, 8785. (80) Eckert-Maksic´, M.; Vazdar, M.; Barbatti, M.; Lischka, H.; Maksic´, Z. B. J. Chem. Phys. 2006, 125, 064310. (81) Bhaskaran-Nair, K.; Demel, O.; Pittner, J. J. Chem. Phys. 2008, 129, 184105. (82) Li, X.; Paldus, J. J. Chem. Phys. 2009, 131, 114103. (83) (a) Maier, G. Angew. Chem., Int. Ed. Engl. 1974, 27, 309. (b) Eckert-Maksic´, M.; Maksic´, Z. B. In The Chemistry of Cyclobutanes; Rappoport, Z., Liebman, J. F., Ed.; Wiley: Chichester, U.K., 2005; p 16. (84) Debbert, S. L.; Cramer, C. J. Int. J. Mass Spectrom. 2000, 201, 1. (85) Schreiner, P. R.; Prall, M.; Lutz, V. Angew. Chem., Int. Ed. 2003, 42, 5757. (86) Li, X.; Paldus, J. J. Chem. Phys. 2008, 129, 174101. (87) Mahapatra, U. S.; Datta, B.; Mukherjee, D. J. Chem. Phys. 1999, 110, 6171. (88) Woon, D. E.; Dunning, T. H., Jr. J. Chem. Phys. 1995, 103, 4572.

JP911581F