Coupled Cluster Methods for Bond-Breaking - ACS Symposium Series

Chapter 5

Coupled Cluster Methods for Bond-Breaking Martin Head-Gordon, Troy Van Voorhis, Steven R. Gwaltney, and Edward F. C. Byrd

Downloaded by UNIV OF GUELPH LIBRARY on September 9, 2012 | http://pubs.acs.org Publication Date: August 14, 2002 | doi: 10.1021/bk-2002-0828.ch005

Department of Chemistry, University of California Berkeley, and Chemical Sciences Division, Lawrence Berkeley National Laboratory, Berkeley, CA 94726

A class of coupled cluster methods designed to describe bond-breaking processes are described. The approach taken involves first approximating the complete valence space Schrödinger equation to account for static electron correlations. For simplicity and tractability, the coupled cluster wave functions are restricted to double substitutions within this active space. Such reference functions are capable of correctly breaking only a single chemical bond, when developed within conventional coupled cluster theory. There are, however, at least two modifications to the theory that allow this limitation to be largely overcome. These modifications are described and discussed. A perturbative approach to correct active space coupled cluster methods for dynamical correlations is then described. The performance of these theories is tested on some model problems to illustrate their current capabilities and limitations.

© 2002 American Chemical Society

In Low-Lying Potential Energy Surfaces; Hoffmann, M., et al.; ACS Symposium Series; American Chemical Society: Washington, DC, 2002.

93

94


Introduction Single reference electronic structure methods are most widely used in quantum chemistry today. This is because of their ease of application; they require no molecule-specific customization unlike multireference methods where either reference configurations or an active space must be chosen for each problem at hand. However, computationally tractable single reference methods, ranging from simple mean-field Hartree-Fock theory to quite sophisticated methods such as CCSD(T) fail to correctly describe potential energy surfaces for bond-breaking. This is manifestly true for calculations based on restricted orbitals where incorrect dissociation products are obtained. The point is perhaps more debatable for unrestricted (spin symmetry-broken) orbitals where potential curves are qualitatively correct using single reference methods. However beginning at equilibrium with a pure spin state and finishing at dissociation with a statistical mixture is also unsatisfactory. B y contrast, multireference electronic structure methods are well-suited to the description of global potential energy surfaces, apart from their inherent difficulty of application. There is hence much incentive to develop electronic structure methods that retain the simplicity of use of existing single reference theories, but address their deficiencies for the description of bond-breaking processes. This chapter reports our recent progress towards this objective. Due to its short length, we make no pretense of reviewing the work of others, but instead seek to provide an overview of our recent work with some selected results. Our general strategy is as follows. We seek to partition electron correlation effects into nondynamical correlations associated with (relatively few) low-lying electron configurations, and dynamical correlations associated with the collective effect of (relatively many) high-energy configurations. We define the nondynamical correlation as being the difference between the mean field Hartree-Fock energy, and the energy obtained by solving the Schrôdinger equation in the space of valence orbitals only. This could be the valence minimum basis, or the perfect pairing active space that supplies one correlating virtual orbital for each occupied valence orbital. The dynamical correlation is then the difference between the valence space Schrôdinger energy, and the energy associated with solving the Schrôdinger equation in the full orbital basis. While the general strategy of dividing electron correlations into nondynamical and dynamical effects is similar to that of multireference methods (1) such as complete active space (CAS) methods (2,3), there is a crucial difference. Specifically, our insistence on the use of the full valence or perfect pairing active space means that there is no need to customize the number of active orbitals for each problem. In other words, techniques based on this procedure will meet the requirements of a theoretical model chemistry. Of course this means that much larger active spaces will be generated for a given


95


molecule than would typically be employed in a standard C A S calculation. As a result, exact solution of the valence space Schrôdinger equation (which scales approximately exponentially with the number of valence space electrons) is intractable with our approach. Instead, we must approximate the valence space Schrôdinger equation as the initial step, and then subsequently approximate the correction for dynamical correlations. Our approaches to each of these two steps are described in each of the following main sections, and are performed within the framework of coupled cluster theory. This ensures proper extensivity of the approximations with respect to the number of electrons in the system.

How good can a CCD wave function be? The trial function that we have chosen to focus on is of the coupled cluster doubles (CCD) type (4-6): φ

Ί*οοο) = ^| >

W

Here the excitation operator promotes pairs of electrons from filled to empty orbitals. We note in passing that a generalized T operator, which includes creation operators for filled orbitals, and destruction operators for empty orbitals is exact (7-9): however this form cannot yet be used tractably, so we do not consider it in this work. We anticipate that the orbitals in the reference single determinant | φ ) will 2

be optimized to minimize the trial energy (10,11), so we shall neglect single substitutions. This also makes optimization within an active space possible (12). While the C C D wave function is well-known, and has desirable properties such as recovery of a substantial fraction of electron correlation at equilibrium geometries and proper extensivity with size, it has not traditionally been viewed as a suitable candidate for bond-breaking. This is because failures of conventional coupled cluster theory are well-known for bond-breaking (13). To investigate whether this is due to limitations of the C C D wave function, or the way in which the C C D energy and amplitudes are conventionally obtained, we have recently performed variational C C D calculations (14). Variational C C D calculations are a restricted form of full configuration interaction (with factorial cost), because the energy expectation value includes contributions from all orbital replacements as high as the number of electrons: ^var-CCD

=

^ C C D |#| ^ C C D ^ / ( ^ C C D | ^ C C D }

(2)

As an example of the results that are obtained in this way, Figures 1 and 2 illustrate the potential curves obtained for homolytically dissociating the two single bonds in the water molecule by full CI, conventional C C D , and variational C C D , in two basis sets.


96


25.0

R(O-H) [Angstrom]

Figure 1. Variational CCD (VCCD), conventional CCD, and full CI (FCI) calculations on the double dissociation ofH O in the STO-3G basis. The experimental bond angle (104.5°) and FCI Brueckner orbitals were used. 2

t

25.0

" 1.5

2.0

2.5 3.0 R(O-H) [Angstrom]

3.5

4.0

Figure 2: Same calculations as Figure l but in a double zeta (DZ) bas t

There are two principal conclusions that can be drawn from these results: While conventional C C D fails dramatically, the variational C C D calculations are well-behaved at long bond-lengths. The difficulties with conventional C C D for this 4 active electron problem are clearly a result of



97 solving for the energy in a nonvariational fashion. Note that conventional C C D is qualitatively correct for a 2 active electron process such as a single bond-breaking, while it fails even more dramatically in a 6 active electron process such as a triple bond-breaking. • The errors associated with variational C C D are significantly larger in the double zeta basis (Fig. 2) than in the minimal basis (Fig. 1). This simple trial wave function can do a near quantitatively accurate job of approximating the valence space (static) correlation (which is all that is present in the STO-3G basis), but is less successful in reproducing both static and dynamic correlation in the larger D Z basis. These results and conclusions suggest a possible path towards practical coupled cluster methods for bond-breaking. The simple C C D trial function is potentially useful, i f we are able to better approximate the variational C C D solution (conventional C C D is inadequate beyond single bond-breaking). Furthermore, the C C D trial function should be applied within a valence active space, because it performs far better in the minimal basis than the D Z basis. From the point of view of approximating the Schrôdinger equation, we are being led to first approximate the valence space (CASSCF) equation, and then subsequently the remaining dynamical eon-elation. We turn to the first of these tasks below.

How difficult must a CCD energy functional be? How might we approximate the variational C C D approach discussed above, without incurring its factorial computational cost? Let us briefly re-examine the conventional C C D energy ansatz as a prelude to discussing how we have chosen to go beyond it. In conventional C C theory the equations that determine the doubles amplitudes can be obtained by minimizing the following functional: (3) In this expression, the variables to be minimized are the usual doubles amplitudes (T ), the corresponding bra amplitudes (A ), and the reference orbitals. In Eq. (3), the T amplitudes are used to define a similarity transformed Hamiltonian: 2

2

2

(4) Since the bra amplitudes occur only as linear terms, the resulting C C D equations for the ket amplitudes (T ) are independent of A : 2

2

(5)


98 Indeed, the Λ amplitudes are not required to determine the C C D energy. They arise in analytical derivatives of the energy as the so-called Z-vector (75). Extended coupled cluster (ECC) theory (16,17) is a powerful existing framework for going beyond the standard coupled cluster energy ansatz, (3) In E C C theory (limited to double excitations, here), the energy is obtained by minimizing the following functional with respect to the bra (Λ ) and ket (T ) amplitudes: 2

2

=


-ECCD

2

Α2

(β Φ Η

(6)

Comparing Equations (3) and (6), it is evident that the usual C C D expression is the leading term in the E C C D energy, when exp(A ) is expanded as a Maclaurin series. Therefore E C C D and C C D will be similar when the amplitudes are small (when C C D works well). For strongly correlated problems such as bondbreaking, it is likely that the more even-handed treatment of bra and ket degrees of freedom would permit the E C C D energy to much better approximate the variational C C D energy than C C D . To our knowledge, molecular calculations at the E C C D level have not yet been reported. This is presumably because of the large increase in both algebraic and computational complexity that is inherent in the minimization of Eq. (6) Indeed, analysis of the E C C D functional (18) suggests that computational cost will scale as i V , far worse than the N scaling of C C D in a straightforward implementation. However, there is an interesting intermediate step between standard C C D and E C C D . This is a model that includes terms quadratic in Λ in the bra function. We term this energy functional the quadratic coupled cluster doubles (QCCD) method (19): 2

10

6

2

^QCCD

= ((1 + Α + 1 Α ! ) Φ Η Φ)

ω

2

The Q C C D energy is to be made stationary with respect to both the ket (T ) and bra (A ) degrees of freedom, as well as the orbitals. In this respect it is a generalization of the conventional optimized orbital C C D method (11). However, the extra term relative to C C D introduces new physics, such as the coupling of the A and T amplitudes. For example, the equations for the T amplitudes can be symbolically written as: 2

2

2

2

2

^(ΐ + Λ ) φ ^ ΗΦ^ 2

=

0

(8)

This coupling will be weak when the A amplitudes are small ( C C D will be recovered with only minor differences), but when the amplitudes are large as in bond-breaking, this may lead to energies that are not as prone to nonvariational collapse. It is beyond our present scope to discuss the Q C C D equations or their computational implementation in any further detail. This task is addressed in 2


99 both published (19) and pending (20) publications. However, let us note some general features. The Q C C D energy and associated amplitudes can be obtained in computational effort that scales with the 6 power of molecular size, and storage that scales as the 4 power of molecular size. These are the same scalings that conventional C C D exhibits, and therefore Q C C D calculations will generally be feasible on molecules for which C C D is feasible. With a valence active space, the cost ratio is a small factor that is independent of basis set size. In the full space, the extra computation increases as the basis set size increases because one additional step in Q C C D scales as V where V is the number of (active) unoccupied orbitals. B y contrast the rate-determining step in C C D scales as 0 V*. Still, Q C C D is dramatically less expensive than E C C D ! The question then is whether Q C C D is robust enough to be faithful to variational C C D in bond-breaking problems. For double bond-breaking in water (see Fig. 1), the calculated energy differences between Q C C D and V C C D do not exceed 5 microHartrees (19)\ On a plot like Fig. 1 such differences are not even visible, and Q C C D is thus an unqualified success. A second and more challenging test is triple bond-breaking in the nitrogen molecule, as shown in Fig. 3. These are also minimal basis calculations, because we know from the previous section that the performance of the C C D wave function for bondbreaking degrades when used out of the valence space. th

th


6

2

50.0

-200.0 1.0

1.5

2.0 R(N-N) [Angstroms]

2.5

3.0

Figure 3. Calculations of the triple bond dissociation ofN in the STO-3G basis, using full CI, conventional CCD, variational CCD, and quadratic CCD. 2


100 The results contained in Figures 3 are also quite encouraging. For N dissociation, the Q C C D results are qualitatively excellent, but there are noticeable differences between variational C C D and Q C C D , which amount to roughly 15 kcal/mol as dissociation is approached. The difference between this test case, and the water double dissociation discussed above reflects the role of A terms (which are neglected in QCCD) in the N problem (corresponding to

2


2

2

6-electron contributions), while they are insignificant in the double bondbreaking problem. Our conclusion at this stage is that Q C C D represents a very significant and yet tractable improvement upon C C D , when used with orbital optimization in the valence space. In this sense, our first generation valence optimized orbital C C D method (12) is clearly superceded by Q C C D in the valence space. We have implemented Q C C D (both in active spaces and in the full space) within the correlation module of the Q-Chem program (21). Both energies and gradients are available, as the latter can be immediately computed from the 1- and 2particle density matrices (and an energy-weighted density matrix) once the optimum orbitals and amplitudes have been determined.

How simple can a CCD wave function be? To this stage our cluster wave function has been left in the general form of Eq. (1), with the restriction that double substitutions are confined to an active space. Within the active space, the double substitution operator couples together all pairs of occupied orbitals with all pairs of virtual orbitals. Physically we expect this description to be more complicated than is essential for a qualitatively correct description of the main correlations in bond-breaking. As is evident in minimal basis dissociation of H , they are surely the alpha-beta bondantibond correlations necessary to permit homolytic bond separation. So, restricting ourselves to the 1:1 perfect pairing active space, the simplest possible version of a C C D wave function would be to retain only the linear number of excitations needed to provide alpha-beta bond-antibond correlations: 2

valence pairs i

This defines the C C D perfect pairing (PP) operator (22,23). Mathematically it converts the coupled cluster equations into the uncoupled cluster equations", because the amplitude equations for the valence pairs are no longer coupled to each other. The only coupling is indirect through the orbital optimization procedure. As a result of this decoupling, we anticipate (and it is true) that the


101 nonvariational catastrophes of conventional C C D will be eliminated. O f course C C D itself is exact for a single pair. Perhaps the next simplest form of trial C C D wave function would be to retain the 1:1 active space, but reintroduce some coupling between valence electron that are associated with different pairs. The simplest way to do this defines what we have termed the imperfect pairing (IP) operator (24): valence pairs t

p


n

=fr

+ Σ

(VK*r+

i*j We have used the unitary group generators, Ey* = a^a +ο£*αγ for conciseness. t

The IP form of the C C D operator is somewhat analogous to the G V B - R C I wave function, although there are important distinctions. For example, the cluster equations for this simplified C C D operator directly couple the different pairs (unlike PP), and accordingly the nonvariational failures associated with full C C D remain in the IP model (24). Given the much simpler form of Eq. (10), which involves only a quadratic number of amplitudes, relative to Eq. (1), it is also more straightforward to understand how and why conventional C C D exhibits difficulties for multiple bond-breaking. The simplest model problem to examine is a double bondbreaking, as in the separation of ethylene into 2 triplet methylene fragments, in a 2:2 active space comprising the 2 bond orbitals. For this problem, it is possible to show that the trial C C D operator, Eq. (10), does not yield a physically correct solution at dissociation. In particular, because there is a forced relation between double and quadruple excitations through the exp(F ) form of the wave function, we find a spurious ionic contribution at dissociation (25). This is a result of the restricted form of the cluster wave function, and is, in a practical sense, a fundamental problem with limited cluster wave functions for bond-breaking. While a detailed presentation of both the analysis of the model problem, and the design of a solution is beyond our scope here, we shall summarize the central result (25). Within the IP form of the C C D operator, we must modify the quadruples term to ensure that no ionic contribution is obtained at dissociation in the model 2:2 double bond dissociation problem. For this purpose, it is sufficient to eliminate terms in the quadruples which involve a repeated index 2

split between different amplitudes (for example the t\

2

t[

2

terms, as opposed

to t t terms). This procedure defines a modification of the IP model that we term the G V B - R C C model, by analogy with the G V B - R C I approach. Whether deleting these terms that prevent correct dissociation in the model double bond-breaking also helps to avoid the nonvariational catastrophes that occur in C C D and IP is unclear: it must be answered by numerical tests. However the answer appears to be affirmative. To illustrate the performance of n

22



102 the 3 models discussed in this section, Figure 4 shows potential curves for the PP method, the IP method and the G V B - R C C method for double bond-breaking in water relative to C A S S C F in the same perfect pairing active space. The results are clearly quite encouraging. PP is qualitatively correct, and G V B - R C C successfully overcomes the limitations of the IP model to obtain nearly quantitative agreement with C A S S C F for this problem. Other examples are presented elsewhere (25), and also support the likely value of G V B - R C C as a very inexpensive approximation to C A S S C F . The development of somewhat more complete pair approximations than Eq. (10) is underway, and promises to still further improve these results. With additional difficulty, it is also possible to apply these simplified T models to the Q C C D energy functional, which will be of considerable interest. 2

R(O-H) [Ang]

Figure 4. Double bond dissociation of the water molecule using the perfect pairing (PP), imperfect pairing (IP) and restricted pairing (GVB-RCC) local correlation models, compared to full configuration interaction (FCI) and Hartree-Fock theory in a minimal (STO-3G) basis.

How can a reference CCD energy be improved? To this stage we have established that the simple C C D wave function is capable of yielding good accuracy when applied in a valence space, with an



103

appropriately modified energy functional. However for quantitative accuracy it will be essential to also account for the neglected correlation effects associated with fluctuations of electrons into orbitals that are not in the valence space. This correlation is of the dynamical type. A s it involves primarily high energy fluctuations, we expect that it should be reasonably well accounted for by perturbation theory. Our problem therefore is to develop such a correction given a reference solution of the C C D type. Piecuch and co-workers (26-28) have recently considered an approach to this problem in which the denominators of perturbative corrections are stabilized by a renormalization term. Our approach, which we summarize below, has been described in a series of recent papers (2931). We begin by recalling (32) that the initially solved C C D problem can be written as a linear eigenvalue problem with the similarity transformed Hamiltonian, E q . (4), right eigenvector

/ ^ ° ^ = |θ),

and left eigenvector

0

(^ή ^ | = (θ| ^1 + A ) . This is true within the space defined by the reference plus 2

active double substitutions, which we shall term the primary space, |p). The remaining substitutions (singles, inactive doubles, triples, etc) comprise the secondary space, |q). Our objective is to perturbatively correct this initial solution towards exact solution of the similarity transformed Schrôdinger equation, H | R) = E | R), where, formally the exact right eigenfunction would be related

\^

to the full

configuration interaction (FCI) wave function by

) = cxp(f )\R).

FCl

2

The similarity-transformed Hamiltonian is partitioned for perturbation theory, with the zero order part recovering the C C D solution in the primary space, |p):

(0)

+

^ =|p)4(p| N)4(q| (1)

+

Coupled Cluster Methods for Bond-Breaking - ACS Symposium Series

Recommend Documents