15354
J. Phys. Chem. 1995, 99, 15354-15368
Property Evaluation Using the Two-Reference State-Universal Coupled-Cluster Method Piotr Piecuch and Josef Paldus*3’ Department of Applied Mathematics, University of Waterloo, Waterloo, Ontario, Canada N2L 3G1 Received: February 14, 1995; In Final Form: April 28, 1995@
The recently developed and implemented orthogonally spin-adapted state-universal (SU) coupled-cluster (CC) theory using a model space spanned by two closed-shell configurations and involving singly and doubly excited clusters (SU CCSD) is applied to calculate static properties of a few typical quasidegenerate systems, for which the range of quasidegeneracy can be continuously varied by changing their geometries. Electrostatic multipole moments and polarizabilities are calculated for the two lowest totally symmetric singlet states of the so-called H4 model consisting of two interacting hydrogen molecules in various geometrical arrangements and for methylene at equilibrium geometry. In both cases, the double-5 plus polarization basis set is employed and the properties are evaluated by the finite-field method. We discuss the role of orbital relaxation in the SU CC property calculations and compare our results with available single-reference (SR) CC, many-body perturbation theory, and configuration interaction (CI) data, including the full CI results providing the exact solution for the given models. The studied systems enable us to examine several important aspects that are encountered in property calculations when using the SU CC approach. In particular, the strongly degenerate region of the H4 model provides us with several physically interesting situations, involving broken-symmetry solutions and a wrong sign or a wrong order of magnitude of the multipole moments at the Hartree-Fock (HF) or even SR CC level of approximation. Our results indicate that SU CCSD provides accurate values for various electrostatic properties in both degenerate and nondegenerate regimes, regardless of whether the relaxed or nonrelaxed orbitals are employed. At the same time, it gives very good property values for excited states. Finally, even when HF or SR CCSD results are qualitatively wrong due to the symmetry breaking, SU CCSD is capable of correcting this behavior.
I. Introduction In recent years, coupled-cluster(CC) theory’-4 has often been used to describe quantum-mechanical many-body systems. The primary advantage of this approach, particularly when compared to its configuration interaction (CI) counterpart, is its size extensivity and relatively low computational cost of including highly excited configurations, which arise in CC theory as products of lower-order excitations, such as singly and doubly excited clusters (CCSD method). At present, applications of CC theory range from the electronic structure of atoms and molecule^^-^ and nuclear physics9 to condensed matter theory, the theory of electron gas, and quantum fluid dynamics,1° to mention only a few examples (cf., also, ref 11). Many important chemical problems involving small and medium size systems can be nowadays successfully addressed thanks to the availability of general purpose CC codes, such as GAUSSIAN 9212 or ACES II.I3 In fact, a number of research groups, including our Waterloo group, have developed over the years their own original general purpose programs allowing for a large variety of CC calculations. However, in spite of this remarkable progress in the development and coding of various CC techniques, the vast majority of the present day applications focuses on closed-shell (CS) nondegenerate ground states, where one can rely on single-reference (SR) CC methods. CC calculations for general open-shell (OS) or highly quasidegenerate CS situations, such as multiple bond breaking phenomena and excited states, where the standard SR CC methods cannot be employed, are far less numerous and far less routine. In spite of the great progress made in the formulation of required To whom all correspondence should be addressed. E-mail: @ theochem.uwaterloo.ca. Also at: Department of Chemistry and (GWC)?-Waterloo Campus. University of Waterloo, Waterloo, Ontario, Canada N2L 3GI. Abstract published in Aduance ACS Absrracrs. October 1. 1995,
pal$
OS or multireference (MR) CC theories (see, e.g., refs 8, 11, and 14-20, and references cited therein), neither our understanding of the cluster structure of general OS wave functions nor the computational efficiency of the existing procedures can be regarded as satisfactory. In fact, a number of physically important OS or highly degenerate CS cases cannot be properly handled within the framework of the CC formalism at this time, despite numerous efforts by a number of research groups (cf. refs 8, 11, and 14-20). The success of the SR CC theory is largely due to various extensions of this formalism to calculations of atomic and molecular electronic properties and, in general, analytic energy derivatives, which allow for vibrational frequency and property evaluation as well as for geometry optimization and search of the potential energy In particular, a number of different algorithms have been suggested for the SR CC calculations of static properties, such as multipole moments and (hyper)polarizabilities. They include direct evaluation of expectation v a l ~ e s , ~ variational . * ~ . ~ ~ and bivariational (biorthogonal) appro ache^,^^-^^ a stationary Lagrangian technique,3Afinite field (FF) method^,^'.^* and a linear-response a p p r ~ a c h . * ~ , ~ ~ Similar procedures have been proposed for the determination of dynamic properties of nondegenerate CS systems,21.29.40-44 but so far the generalization of all these methods to OS or MR CC theories, that are applicable to OS and highly degenerate CS electronic states, has attracted little a t t e n t i ~ n . ~ In ” ~fact, ~.~~ very little is known about the applicability of genuine MR CC methods to molecular properties. For example, Pal and Ghose recently discussed a response-type MR CC approach within the framework of the so-called Fock-space or valence-universal (VU) f o r m a l i ~ m , but ~ ~ .they ~ ~ have restricted themselves to formal aspects of the theory. A number of static properties have been recently computed for OS singlet and triplet states of watefi7 and ozoneA6using the so-called two-determinantal CC
0022-3654/95/2099-15354$09.00/0 0 1995 American Chemical Society
Coupled-Cluster Properties approach,48which represents a special case of a more general Hilbert-space or state-universal (SU) MR CC method.14 It should be noted, however, that the OS singlet and triplet states are in principle SR one-dimensional (though multiconfigurational) problems, which can be efficiently handled using one of the SR-like spin-adapted state-specific (SS) CC theories, such as the unitary group based SS CC approach of ref 49. Clearly, for general OS situations, none of the above approaches is applicable and a genuine MR formalism is called for. In spite of several advances in the development and implementation of the SU CC theory'6,50-55 and a number of as well as actua155,61-65 molecular applications to mode150~52~55-60 systems (cf., also, ref 66 for the fvst application of the linearized SU CC theory'), the applicability of the MR SU CC theory to molecular properties of highly quasidegenerate systems, for which a MR description is essential, has never been thoroughly investigated. We have thus decided to employ the SU CC theory for the determination of electrostatic properties of a few typical quasidegenerate systems, for which the range of quasidegeneracy can be continuously varied by changing their geometries. The results we report include calculations of multipole moments and polarizabilities for the two lowest-lying totally symmetric singlet states of the double-g plus polarization (DZP) basis H4 modeF7 consisting of two slightly stretched interacting hydrogen molecules in various geometrical arrangements,68as well as static property calculations for the two lowest totally symmetric singlets of the CH;! molecule, as described by the DZP basis of ref 69 (cf., also, ref 70). Both systems enable us to address several important problems regarding the usefulness of the SU CC approach in property calculations. The quasidegeneracy of the two lowest-lying totally symmetric singlets (which are overall second and fourth electronic states) of methylene is well-known. In this way, we can explore the ability of the SU CC method to describe properties of rather highly excited electronic states. On the other hand, the H4 model provides us with a number of physically interesting situations, including degenerate and nondegenerate regimes, broken-symmetry solutions, and a wrong sign or a wrong order of magnitude of the multipole moments at the Hartree-Fock or even SR CC level of approximationin the strongly degenerate region (see section I11 for details). Again, we may use this system to investigate the applicability of the SU CC method to properties of excited states and the shape of property surfaces when we break a single chemical bond. Both methylene and H4 DZP systems are small enough to be analyzed using a variety of quantum-chemical approaches, including the full CI (FCI) method providing exact solutions for these models. The SU CC results will thus be compared with the available SR CC, many-body perturbation theory (MBFT), and limited and full CI data. Clearly, investigation of properties using the SU CC theory and a comparison of the results with those obtained with other approaches may shed some light on the performance of the SU CC method in general, since the properties are often more sensitive to the quality of the electronic wave function than are the corresponding energies. In actual calculations, we use the most complete version of the recently formulated and i m ~ l e m e n t e dorthogonally ~~ spinadapted (OSA)71-73SU CC theory involving single and double excitations (SU CCSD) from the two CS configuration model space. The advantages and shortcomings of the spin-orbital us OSA version of the CC (in particular, SU CC) theory were discussed earlier in a number of places (cf., e.g., refs 51 and 71 -73). We thus only mention that the OSA formulation has the advantage of using the minimum number of cluster amplitudes and an obvious formal appeal of exploiting the spin
J. Phys. Chem., Vol. 99, No. 42, 1995 15355 symmetry of the Hamiltonian. It also offers a direct connection with the corresponding spin-adapted CI approaches. Although the OSA CC theories are more complex than their spin-orbital counterparts and require a higher level of sophistication in the corresponding mathematical considerations, they offer substantial computational memory and time savings and, in consequence, an opportunity to study relatively large molecular systems or small systems described by large molecular orbital (MO) bases (cf., e.g., ref 64). Recall that the OSA SU CCSD approach of ref 55 includes previously neglected'6,50.51,56-62,65 cubic and quartic coupling terms. Numerical evidence indicates that these terms could safely be ignored,55 but, on the other hand, the cost of their computation is very small so that their inclusion does not represent any significant problem. In the present paper, properties are calculated as numerical energy derivatives using the FF method.74 In order to study the role of orbital relaxation in the presence of electric field, we performed both relaxed and nonrelaxed SU CCSD (as well as SR CCSD) calculations. The paper is organized as follows. An outline of the OSA SU CCSD method of ref 55 and the nonrelaxed and relaxed SU CCSD (as well as SR CCSD) models for static property evaluation are presented in section II. Computational details and studied systems are described in section III. The uncorrelated results and the problems arising in property calculations for the studied systems are discussed in section IV,while the correlated results of property calculations using the OSA SU CCSD method are given in section V. In particular, section VB discusses an interesting case of symmetry breaking by the SU CC theory in the strongly degenerate limit of the H4 model. The last section VI summarizes the main conclusions concerning the performance of the SU CC method and discusses future research projects. 11. State-Universal Multireference Coupled-Cluster Theory and Static Property Evaluation
The SU CC theory was described in a number of places (cf. refs 8, 14, 16, and 50-62) so that we restrict ourselves to bare essentials that are relevant for property calculations. We then focus on the complete two-reference OSA SU CCSD theory of ref 55 and discuss relaxed and nonrelaxed SU CC models. A. Basic Formalism. All genuine MR CC formalisms involve two basic concepts, namely, that of the multidimensional reference or model space 4,providing a reasonable zero-order approximation of the exact manifold & which is spanned by M quasidegenerate eigenstates IW,), v = 1, ..., M , of the electronic Hamiltonian H,and that of the wave operator U, which transforms the model space & into the exact space Al (cf. refs 8, 11, and 14-20, and references cited therein). To define the reference configurations I@,),p = 1, ...,M , spanning &, one partitions the employed MO set into the three disjoint subsets of core, active, and virtual orbitals. The core orbitals are occupied and the virtual ones unoccupied in the references. The configurations I@,) differ in the occupancies of active orbitals. All possible distributions of active electrons among the active orbitals result in a complete model space, which is invariant with respect to unitary transformations of active orbitals.I4 In the presence of spin and spatial symmetries, & factors into the direct sum of invariant subspaces, each of which is complete for a given symmetry specie^.^.'^,^^,^^,^^ When defining the wave operator U , we rely on the Bloch f o r m a l i ~ m , ~ , ' ~ .in' ~which . ' ~ . ~U ~ is an inverse of a bijective map P = P &obtained as a restriction of the projection operator P3& & to where 9% represents the entire N-electron Hilbert space, assuming simultaneously that
-
J. Phys. Chem., Vol. 99, No. 42, 1995 15357
Coupled-Cluster Properties
involving two active orbitals k and 1 that belong to different symmetry species of the spatial symmetry group (or one of its subgroups) of the system (cf., also, refs 16, 50, 51, and 5658), while the core consists of a number of doubly occupied orbitals,
Here, XK,m(XK,m)designate the creation (annihilation) operators associated with the spin-orbital basis IK,m) = I K ) 8 I1/2,m),m = &l/2, ( K , 1,..., label spatial MOs employed) and 10) is the true vacuum (zero particle) state. Within the OSA formalism, the M = 2 model space
4=SP~~l@I)~l@2)1
+ p [ q q+ (-1)5qE3, i = 0, 1
$$(i) = 1/2(2i
(e))*
H(A>=H+AW
(20)
is complete for the totally symmetric singlet eigenstates of H (see, e.g., refs 16, 55, and 57), in contrast to the spin-orbital formulation, in which case the states of different multiplicity are involved and M = 4. We can thus express the mono- and biexcited components of the cluster operators P)in terms of the OSA particle-particle-hole-hole coupled single- and double-excitation operators (cf. refs 16 and 76; see also refs 71-73)
SQ, = 2 - L / 2 q
energy expectation values of the SR CC theory, involve at most quadratic terms. The role of different terms of the SU CCSD formalism and various ways of handling the effective Hamiltonian and bilinear coupling terms, were studied in considerable detail in refs 5558. While the inclusion of bilinear terms [such as the direct and coupling terms] is crucial for the elimination of the singular behavior of the linearized version of the SU CCSD f o r m a l i ~ m the ? ~ effect ~ ~ ~ of ~ ~cubic ~ and quartic terms, especially of the cubic and quartic coupling terms, is very Computation of the latter terms is, however, relatively inexpensive, so that in the following we consider the complete tworeference SU CCSD formalism including all higher-order cubic and quartic terms. B. Relaxed and Nonrelaxed SU CC Formalism for Property Evaluation. To calculate static properties, we consider the perturbed Hamiltonian
(21) (22)
(26)
where H H(0) represents a “zero-order” Hamiltonian describing a given N-electron system in the absence of the perturbation, W is the property in which we are interested, and 1 is the perturbation parameter, which has usually a simple physical interpretation. For example, when calculating the dipole moment, p = cUx, py, pd = bx,,px2, pxg),or diagonal components of the dipole polarizability, a,,, ( i = 1-3), W is one of the components of the dipole moment operator and 2 is the negative of the electric field (for the off-diagonal components a,,,i j , it would be necessary to consider a multiple perturbation in all three components of the electric field). When W is the quadrupole moment operator with Cartesian components &, (i, j = 1 - 3), 1 corresponds to the negative of the electric field gradient, etc. We assume throughout this paper that W is a spin-independent one-electron operator, so that
*
respectively, where
w = wA$?q
(27)
I12
zq=
c xK,mxa,m
(23)
m=-1/2
are the standard orbital unitary group generators77and a,/3, ... (e,u,...) designate the occupied (unoccupied) orbitals in lop). Using Einstein’s summation convention, we can express the oneand two-body cluster operators in the formS1355.71-73
q)=
@)
a @)$ a
(24)
where wi = ( ~ 1 ~ 1are 1 ) the corresponding property matrix elements. We next consider the field-dependent perturbed problem,
which we tackle using the SU CC theory. Here IW,(O)) = IY,) are the eigenstates of the M-dimensional unperturbed (fieldfree) problem (cf. section IIA),
HIY,) = E,,IY,),
where @)tiand 0’ = 0, 1) designate the corresponding OSA mono- and bixecited cluster amplitudes. As mentioned earlier, the completeness of &,eq 20, requires that the terms involving only active orbital indictes k and 1 must be excluded. Explicit equations of the OSA SU CCSD theory, based on the model space of eq 20, were given in ref 55. The left-hand side direct terms Ap’(G!), eq 14, where GI = s“, or S$(i)), are exactly the same as those appearing in the OSA SR CCSD theory of ref 73, except that now we must write them for two reference configurations, while excluding equations involving all-active excitations. These terms are at most quartic in cluster amplitudes, as are the coupling term A,@’(G:),eq 15, (cf. ref 55) and the off-diagonal elements of the effective Hamiltonian matrix, while the diagonal elements I$, being essentially the
v = 1, ...,M
(29)
and E,(O) = E, are the corresponding eigenvalues. The energies E&) can be expanded in a Taylor series,
where the corrections
represent various static properties of the system for the electronic state IY,). Specifically, when W is the multipole moment operator, the “first-order” correction E‘,“ is simply the value of the corre-
15358 J. Phys. Chem., Vol. 99, No. 42, 1995
Piecuch and Paldus
TABLE 1: Comparison of Basic Elements (Model Spaces, Cluster and Wave Operators, Wave Function Ansatz, and Effective Hamiltonians) of Relaxed and Nonrelaxed SU CC Methods [SU CC(R) and SU CC(N) Approaches, Respectivelyp SU CC(N) SU CC(R)
pp'(A)=
@'tI(A) 'p'GI(A) I
q=l
q= I
'I=I
The functional dependence of the quantity X on the parameter is designated by X(A). If a given quantity X must be calculated in the absence of the field [cf. the SU CC(N) model], we indicate this fact by writing X ( 0 ) (see the text for the mathematical meaning of various symbols). sponding multipole moment in the quantum state described by the wave function IYJ. Using the Hellmann-Feynman theorem for the exact unperturbed wave function IY,,), we can simply write
Of course, in most practical applications, IY!,) is only an approximate solution of the Schrodinger equation (eq 29), so that the Hellmann-Feynman theorem and eq 32 no longer hold. However, the more accurate the approximate wave function, the smaller the difference between the correction I$" and the expectation value, eq 32. Higher energy derivatives E ': (eq 31 with k > 1) describe higher-order properties. For example, when W is the dipole moment operator, the corrections E':', and describe diagonal components of (hyper)polarizabilities a,/?, and y for the electronic state lYv). In particular, for the dipole field when W = p,, (i = 1 - 3), we have that
e',
where we indicated a given electronic state by adding the argument v to a,,,. Clearly, the zero-order correction d:' represents the energy of the system in the absence of the field, i.e. E ': = E,,. Generally, we can employ two different SU CC approaches to solve the perturbed problem of eq 28 and to evaluate the corrections Depending on whether the orbitals defining the reference configurations I CDp) and excitation operators @!GI are determined in the presence or in the absence of the perturbation, we obtain the relaxed or nonrelaxed SU CC models. The relaxed SU CC method, in which MOs are allowed to relax with changing 1,will be designated by an acronym SU CC(R). In the nonrelaxed SU CC theory, hereafter designated as SU CC(N), the perturbed quasidegenerate problem, eq 28, is solved using the perturbed Hamiltonian H(A), eq 26, but the MOs used to construct the references I@'p) and excitations @)GI are those of the unperturbed (A = 0) problem. Clearly, in the limit when all possible excitations for a given atomic orbital
g').
(AO) basis set are included, the relaxed and nonrelaxed approaches become equivalent and correspond to the full CC or FCI model. However, for truncated SU CC approaches, such as our two-reference OSA SU CCSD method described in section IIA, the relaxed and nonrelaxed models will differ. This situation reminds us of the distinction between the relaxed and nonrelaxed SR CC methods [designated, respectively, as SR CC(R) and SR CC(N)] or between relaxed and nonrelaxed CI models78 [denoted by CI(R) and CI(N), respectively]. Recall that, in the SR CC theory, the orbital relaxation effects play a minor role,38 since the orbital relaxation is effectively included in the nonrelaxed theory through the presence of nonlinear higher-order terms involving TI clusters (for the analogous role of monoexcited configurations in limited CI property calculations, cf. ref 78). This feature of the SR CC method can be useful for the development of efficient analytical methods for property evaluation, such as the linear-response CC theory of M o n k h ~ r s t . ~It' ,will ~ ~ thus be interesting to examine the role of the orbital relaxation effects in the context of the MR SU CC formalism. The essential equations that distinguish the relaxed and nonrelaxed SU CC methods are summarized in Table 1. The main purpose of this table is to highlight the quantities which explicitly depend on the parameter 1. We could use this table to derive analytic expressions for the individual corrections E',"' (more precisely, corrections to the effective Hamiltonian matrix) in terms of the corrections to cluster operators PI [derivatives of PI(,?) with respect to 1 at A = 01. This would parallel the SR linear-response CC f ~ r m a l i s m or ~ lthe . ~ ~response VU CC t h e ~ r y . ~Alternatively, ~,~~ we can determine the using numerical differentiation techindividual corrections niques. In this case, the SU CC equations are solved for a number of values of I. near 0, and the resulting functions &,(A) are differentiated numerically (the so-called FF method). The computational procedures used to calculate the derivatives of &(A) are described in the next section. In this study, all property calculations are carried out using our two-reference OSA SU CCSD theory.s5 We must note, however, that this theory has certain restrictions and is not applicable to all types of perturbations or fields (unless we
e:'
J. Phys. Chem., Vol. 99,No. 42, 1995 15359
Coupled-Cluster Properties decide to use it in an ad hoc nonrigorous manner). It may happen that the applied field (or the perturbation W) will break the spatial symmetry of the unperturbed (field-independent) Hamiltonian H and lower its spatial symmetry. If this lowering is such that the active orbitals k and 1 no longer span different irreducible representations of the spatial symmetry group of the perturbed Hamiltonian H(L), the OSA SU CCSD theory described in section IIA cannot be applied. In this case, we would have to consider a three-reference theory obtained by enriching the model space of eq 20 with the OS singlet configuration (cf. ref 55) (34) For example, we can use the two-reference OSA SU CCSD method to calculate the z-component of the dipole and the polarizability component a, for the two lowest-lying totally symmetric singlets of methylene [ G Usymmetry; here the z-axis is the principal C2 axis and (y, z) is the molecular plane], since the relevant pzoperator preserves the C2usymmetry of the fieldindependent Hamiltonian H. However, we cannot calculate ax, since the px operator is not CzUtotally symmetric [it transforms according to the B I ( C ~ irreducible ~) representation]. The px operator lowers the CzUsymmetry of the Hamiltonian to C, [generated by the a(xz) plane bisecting the molecule], and, as a result, the al(C2,) and b1(C2~) symmetries of the active orbitals (see the next section for details) lower to a’(C2) so that the active orbitals are no longer of different symmetry. We can calculate, however, the a, component, despite the fact that the cly operator is not C2u totally symmetric. In this case, the Hamiltonian symmetry also lowers to C, [generated in this case by the a(yz) plane], but the symmetry species of active orbitals remain distinct [the corresponding irreducible representations being a’(C,) and a”(C,)]. We can also use the two-reference OSA SU CCSD theory to determine all three components of the quadrupole moment QxAz (i = 1 - 3) for CH2 in one of its totally symmetric singlet states, since the corresponding quadrupole moment operators are CzUtotally symmetric, but we must be aware that this may not be possible in general.
III. Model Description and Computational Details
In order to examine the performance of the SU CCSD method in property calculations, we considered selected components of the dipole and quadrupole moments and polarizability tensor for the two lowest-lying totally symmetric singlets of the H4 model system at a few representative geometries and for the two lowest totally symmetric singlets of CH2 at the equilibrium. Both systems exhibit large configurational quasidegeneracy involving two CS-type configurations and, thus, are ideally suited for our study. Throughout the present paper, we use the ground-state restricted Hartree-Fock (RHF) orbitals to generate pertinent configuration spaces so that one of the two reference configurations is always the ground-state RHF wave function. The H4 model consists of two interacting hydrogen molecules arranged in an isosceles trapezoidal configuration, with all nearest-neighbor intemuclear separations fixed and equal to a (ref 68; see Figure 1). As in previous studies using this model system,5°,55-57,68 the hydrogen molecules are slightly stretched and the H-H distance a is chosen to be 2.0 au. The geometry of the H4 model is described by a single parameter a E [0, 0.51 (cf. refs 57, 67, and 68) so that the a = 0 limit represents the square conformation and a = 0.5 corresponds to the linear arrangement of hydrogen atoms (cf. Figure 1). The equilibrium geometry of CH2 is taken from refs 69 and 70.
I L
S=an
Figure 1. Geometry of the H4 model.
In general, molecular properties are tensor quantities whose components depend on the choice of the molecular coordinate system (see ref 79; see, also, ref 80, and references cited therein). In particular, when the origin of the molecular coordinate frame is shifted to another position, the components of the quadrupole moment tensor describing neutral species change by terms which are proportional to dipole moment components79(in this case, dipole moment components remain unchanged). In this paper, we assume that both the H4 system and the CH2 molecule lie in the (y, z) plane, with the positive z-axis bisecting both systems, so that the x-axis is perpendicular to the molecular plane. For the H4 model, the origin of the molecular frame is placed at the center of the H(2)-H(3) bond (see Figure l), so that the H(2) and H(3) atoms are situated along the y-axis, independently of the value of a. We do not use the center of mass as the origin, since the center of mass shifts its position from the center of the H(2)-H(3) bond for a = 0.5 to the center of the square formed by all four hydrogen atoms for a = 0. For methylene, we follow the prescription of refs 69 and 70 and place the C atom at the origin of the coordinate system. This enables us to directly compare our property results (e.g. quadrupole moments) with results obtained by other authors (cf. ref 81) who used the same convention. In any case, the choice of coordinate system is not of primary importance for the conclusions drawn in this paper, since we always compare our CC results with exact or very accurate results obtained using the FCI or CISDTQ (CI with singles, doubles, triples, and quadruples) methods. Both systems are described using the A 0 basis set of the DZP quality. The basis set for H4 is taken from ref 67. For CH2 we use the same basis set as that in ref 69, except that we employ the six Cartesian components of the carbon d orbitals rather than the five spherical components. In particular, the lowest core orbital (-C 1s orbital) is kept frozen. Neither the presence of the 3s combination of 3d functions in our basis set nor the freezing of the core are expected to play a significant role in the CC (SR, SU CC) property calculations for CH2 (cf., e.g., refs 39 and 81). Although in this case the use of much larger bases would be recommended to obtain “experimental” accuracy (cf., e.g., refs 64 and 82), the DZP basis set is large enough to enable us to draw valuable conclusions regarding the applicability of the SU CCSD method in property calculations. The DZP H4 model is simple enough to allow numerous computations using several different methods, including FCI providing the exact solution. For this reason, this model, as well as its simpler minimum basis set (MBS) counterpart, were often applied in the past to test novel quantum-chemical approaches (cf., e.g., refs 50, 55-57, 59, 67, and 68, and references cited therein). The same features characterize the CH2 molecule, thanks to the availability of the FCI results in the DZP basis, both for the energy70 and for various proper tie^?^ at equilibrium geometry. Again, a number of many-body treatments, including the SR and MR CI and multiconfigura-
Piecuch and Paldus
15360 J. Phys. Chem., Vol. 99, No. 42, 1995 tional self-consistent-field (MC SCF) methods,69s70.83.84 as well as the MR MBPT based on the complete active space (CAS) MC SCF the generalized valence bond (GVB) method,89and the SR MBFTKC theory (cf., e.g., refs 81, 90, and 91), have been tested on this methylene DZP model. We should also mention here our most recent results for this model obtained with the OSA SU CCSD and SS CCSD methods developed in our g r o ~ p . Of ~ ~particular . ~ ~ importance to us are the results of refs 69 and 81, where various limited CI, MC SCF, and MBPT results for the properties of CH2 are collected. The spatial symmetry of the H4 and CH2 systems is Czij. Independently of the value of the parameter a, a E [0, 0.51, the ordering of the MOs for the H4 model, corresponding to increasing orbital energy, is as follows:
(35) where the superscripts 2 or 0 indicate the orbital occupancies in the ground-state RHF configuration. For methylene at the equilibrium, the energy ordering of the MOs is
In both cases, the highest occupied MO (HOMO) and the lowest unoccupied MO (LUMO) (indicated by bold font) transform according to different irreducible representations of the Czi group so that the two-reference OSA SU CCSD formalism of section 11, employing HOMO and LUMO as active orbitals, applies. We thus use it to calculate a number of static properties of H4 and CH2, such as the z-component of the dipole moment, p:, the zz diagonal component of the polarizability tensor, a::, or any nonvanishing component of the quadrupole moment tensor, Qr2,8(i = 1-3), for the two lowest IAl(Cz,) states (designated as 1'Al and 2'Al states). The corresponding perturbations W do not change the C?lrspatial symmetry of the Hamiltonian so that there is no need to include the OS singlet configuration, eq 34, in the model space (cf. the discussion at the end of section IIB). This configuration corresponds to a separate 'Bl(C2,) (for CH2) or 'B2(C2,) (for H4) subproblem, which is not the subject of this study. In the following, we consider the traceless quadrupole moment^,'^ (37) with the trace
calculations for stronger and stronger fields A, we used the "analytical" continuation of the relevant CC solutions; namely, the resulting amplitudes for A = 0 were used as an initial guess for A = fO.OO1 au, the converged amplitudes for A = f O . O O 1 au served as the starting approximation for A = f0.002 au, etc. We also used the CzL spatial symmetry to speed up computations. This was possible, since the applied dipolar or quadrupolar fields preserve this symmetry (cf. the discussion above). The CC energies were always converged to a very high accuracy, such as 14-15 decimal places for the energy or 1314 figures for cluster amplitudes. In this way, we were able to obtain highly accurate property values, with up to 10 stable decimal figures for the first-order properties and up to 6 stable decimal figures for the second-order properties. In all cases, the nonrelaxed SR CCSD values of dipole and quadrupole moments were verified analytically by employing our newly developed linear-response CCSD program.39 All SR and SU CC calculations reported in this work were obtained using the general-purpose programs developed in our group.jS In order to calculate properties, we equipped these programs with the FF option. Thanks to the OSA formulation of CC theory and a rather small size of studied systems, the computer cost of each individual calculation was relatively low. We were able to obtain highly accurate property values, since neither the convergence of CC solutions to a high accuracy nor the use of a relatively large number of field values constituted a serious problem. We also benefited from a newly designed algorithm for the solution of CC equations.94 Some recent improvements made this algorithm even more efficient.95 The required molecular integrals were calculated with the electronic structure package GAMESS96and with routines that form part of our linear-response CCSD codes.39 The latter routines were used to generate the one-electron integrals wt defining the perturbation operator W, eq 27. They were exploited in nonrelaxed CC calculations. All one- and twoelectron molecular integrals defining the unperturbed Hamiltonian H = H(0) in nonrelaxed CC calculations and all oneand two-electron integrals defining the perturbed Hamiltonian H(A) in relaxed CC calculations were computed using the routines that form part of the GAMESS CI system. GAMESS was also used for the initial SCF calculations and to obtain various limited or full CI results (limited CI for CH2 and FCI for the H4 models). All calculations were performed on our small Silicon Graphics Challenge L and Silicon Graphics 4D/ 380 GTX minicomputers.
IV. Zero-Order Restricted Hartree-Fock Properties
As mentioned earlier, the relaxed and nonrelaxed property values were calculated using the FF method. The required first and second energy derivatives at A = 0 (cf. eqs 31-33) were obtained by performing each time a series of energy calculations for up to 11 values of I,, such as A = 0.0, 10.001, f0.002, f0.003, 10.004, f0.005 au, and differentiating the resulting functions &(A) (Y = 1, 2) at A = 0 numerically using the Stirling interpolation formula.92 The relevant expressions for the highorder polynomial coefficients appearing in the Stirling interpolation formula were generated with the symbolic manipulation language MAPLE.93 The same procedure was employed for all property calculations and all methods considered in this study, including SCF, CI, SR CCSD, and SU CCSD. The SCF and nonrelaxed CI first-order properties [SCF and CI(N) multipole moments] for the ground IAl(C2J state were also calculated as expectation values, eq 32 (cf. ref 78). To facilitate CC
It is well-known that the degree of configurational quasidegeneracy in the H4 model can be continuously varied by changing the defining parameter a. In the a = 0 limit, the two model space configurations,
become exactly degenerate, while for a = 0.5 it is the reference RHF configuration I Q l ) , eq 39, that dominates the ground-state wave function IY 1),56,57,68 Although the potential energy curves for a E [0, 0.51 have been studied in great detail using various methods, little is known about the properties characterizing this model. To get an idea of the problems involved and, in particular, of the role played by the correlation effects, we first
J. Phys. Chem., Vol. 99,No. 42, 1995 15361
Coupled-Cluster Properties
O,O5I 0.04
Y
0.00
.o'ob.O
0.1
0.2
0.3
0.4
0.5
a
Figure 2. Dipole moment component p2 for the ground state of the DZP H4 model with a = 2.0 au as a function of the parameter a.The open squares and the triangles designate, respectively, the RHF and
FCI results. examine the studied properties at the uncorrelated, zero-order SCF (or RHF) level. The a-dependence of the pzdipole moment of the H4 ground state, obtained with the RHF wave function, is compared with the corresponding FCI values in Figure 2. Clearly, in both the a = 0 and a = 0.5 limits, the RHF and FCI dipole moments vanish thanks to the spatial symmetry of the corresponding nuclear configurations (D4h and D-h, respectively). Consequently, the general shapes of the RHF and FCI property curves are rather similar. Interestingly enough, the exact FCI p, values are invariably lower than the RHF ones (cf. Figure 2). However, the negative p y values in the vicinity of the square geometry are not faithfully approximated at the RHF level. For example, for a E [0.01, 0.021, p y ( 1 l A l ) -0.004 au, at least an order of magnitude larger absolute value than I p ~ ( 1 ' A I ) I< 0.0003 au [here and in the following text, nx(nR)designates the property n obtained with the method X for the nth electronic state of the symmetry R]. In fact, ptW(llAI)remains negative only for a < 0.02, whereas the exact FCI value p y ( 1'AI) is negative for a E (0, 0.05), so that in the region a 0.02-0.04 RHF gives a wrong sign of pz(1IAl). For a L 0.05 both dipole moments RHF 1 pz (1 A]) and pF(1IAl) are positive, reaching their respective maxima at a 0.24 and 0.26, while the difference between the RHF and FCI values remains substantial over the entire region of a. At a = 0.25 it reaches 0.007 763 au, which is nearly 17% of the corresponding FCI value p y ( 1 ' A I ) = 0.045869 au. This difference gradually decreases to zero when a approaches the nondegenerate a = 0.5 limit. Although the dipole moment of H4 is rather small (for both llA1 and 2'Al states), its value is extremely sensitive to the quality of the wave function. Other properties are equally sensitive. In the vicinity of the a = 0 limit, the correlation contribution to the polarizability component az,( l'A1) is nearly equal to the total RHF value a Y ( l 1 A I ) . In this region, the RHF method gives a wrong sign and a wrong order of magnitude, or at best extremely large (1100%) errors in the calculated values of Oyy(llAl)and Oz,(llAl). For the square geometry, we deal with an interesting case of symmetry breaking, where a simple relationship between the exact values of O y y ( l l A ~and ) Ozz(llAl),namely,
-
FCI 0, (a-0) = O:y(a=O)
(41)
which is a consequence of the D4h symmetry of the a = 0 H4 model, is violated at the SR SCF level, so that various correlated methods using RHF orbitals may also violate this relationship (see section VB for details). In fact, the correlation contribution
to each of the three components Oxix,,i = 1-3, remains large over the entire region of a , including the nondegenerate a = 0.5 limit, where the errors in the RHF values (relative to FCI) exceed 30%. It is thus interesting to find out to what extent, if at all, the SU CCSD formalism is capable of correcting these inadequacies of the RHF property values. In particular, the following two problems are of special importance: the ability of the tworeference SU CCSD theory to cope with the symmetry breaking problem occurring at the square geometry (cf. eq 41) and the quality of the SU CC property values for the first excited IAl(C2J state, which is relatively poorly described by our model space, spanned by configurations [ @ I )and 1@2), eqs 39 and 40, in the nondegenerate ( a = 0.5) limit (cf. refs 55-57 and 59). For similar reasons, it is interesting to examine the quality of the SU CC properties of CH2, whose ground and fist-excited 'A1 states require at least two-dimensional model space involving configurations
1
1)
= 1 (la1>2(2a1)2(lb2)2(3a1)2I
(42)
in order to obtain good energies and properties.55,63,64~69~70~83
V. Coupled-Cluster Results and Their Discussion Our results are presented in Tables 2-5. Tables 2-4 pertain to the H4 model, and Table 5 pertains to CH2. For the lowest 'A1 state we give the results for pz (nonrelaxed and relaxed values), Ox,,,,i = 1-3 (nonrelaxed values), and azz(nonrelaxed and relaxed values), and for the first excited 'AI state the results for pz and a,, (nonrelaxed and relaxed values). We did not calculate the relaxed values of Ox,, since in this case we require relaxed MOs and molecular integrals generated by GAMESS. However, GAMESS does not allow for other than dipolar fields. Obviously, the nonrelaxed SU CCSD values of quadrupole moments O,,,, for the first excited 'AI state were calculated along with the SU CCSD(N) values of O,,x,(llAl) (in this case we could rely on our own integral routines39), but we decided not to include them in Tables 2-5, since we could not compare them with the results of other (e.g. CI) calculations. Although the nonrelaxed CI calculations of O,,,, do not require energy differentiation, and the CI(N) results for Ox,, can be obtained as expectation values, eq 32 (cf. ref 78), we were unable to generate them for the first excited 'A1 state using GAMESS. For the similar reason, we do not include in Tables 2-5 the SU CCSD(N) values of the quadrupole polarizabilities Cua, Cyy,yy, and Czz,zz, which we obtained along with the quadrupole moments Qu, Qyy, and Q,, respectively. The SU CCSD(N) values of 0,,,(2'A1) and C,,,,,(nlA1) (i = 1-3; n = 1, 2) are available upon request. In all cases, the property results are accompanied by the corresponding energies. We begin our discussion with the H4 model (sections VA and VB). The SU CC properties of methylene are discussed subsequently. A. H4 Model. The reported SU CCSD results for H4 include the following six geometries: the limiting square (a = 0) and linear (a = 0.5) conformations,where the dipole moment vanishes (Table 2), two geometries corresponding to the strongly degenerate region (a= 0.01, 0.02), where RHF gives a wrong order of magnitude (a= 0.01,0.02) or a wrong sign (a = 0.02) of the dipole moment (Table 3), and two geometries from the intermediately degenerate region (a = 0.1, 0.25), where the differences between the RHF and FCI property values remain large, and where the two-configurational SU CC description of
15362 J. Phys. Chem., Vol. 99, No. 42, 1995
Piecuch and Paldus
TABLE 2: CC and FCI Energies and Properties (au) for the Two Lowest AI(CZJSinglet States of the DZP H4 Model with in the Square (D4h or a = 0) and Linear (D-h or a = 0.5) Conformations" 1' A I 2'Al
a = 2.0 au
a 0.0
0.5
methodb SCF FCI SR CCSD(N) SR CCSD(R) SU CCSD(N) SU CCSD(R) SCF FCI SR CCSD(N) SR CCSD(R) SU CCSD(N) SU CCSD(R)
e,,
0::
a::
0.365 95 0.418 98 0.427 70
-2.870 17 -0.209 49 -0.467 34
2.504 22 -0.209 49 0.039 64
C
C
c
0.417 02
-0.209 01
-0.208 01
c
C
C
-0.541 35 -0.407 93 -0.418 59
1.082 7 1 0.815 85 0.837 18
-0.541 35 -0.407 93 -0.41 8 59
C
C
C
-0.393 14
0.786 28
-0.393 14
C
C
C
14.990 7.894 8.385 8.366 7.860 7.905 1.566 1.566 1.568 1.567 1.567 1.567
a,,
E
-1.931 -2.063 -2.057 -2.057 -2.063 -2.063 -2.150 -2.232 -2.231 -2.231 -2.235 -2.235
750 112 604 604 798 798 368 700 948 948 079 079
E
a::
-1.981 434
8.904
-1.978 061 -1.978 061
8.857 8.892
-1.865 105
1.526
1.863379
1.538 1.534
-
- 1.863379
In this case, ,u,(llAl)= ~ ~ ( 2 ' A=l0. ) For a = 0, the CC methods employing ground-state RHF MOs, as well as the RHF method itself, break the relationship 0,!(l'Al) = O:-(llAI), which is satisfied by the exact (FCI) wave function. No such breaking occurs for a = 0.5, where @.rx(IIAl)= @::(l1A1) for all methods considered (see the text for details). Note that despite the higher symmetry for the square ( a = 0.0) and linear ( a = 0.5) structures, we use the C2Lnotation to label the electronic states. The C?
