Multireference coupled cluster method based on the Brillouin-Wigner perturbation theory

q 2005 Elsevier B.V. All rights reserved. Theory and Applications of Computational Chemistry: The First Forty Years Edited by C. Dykstra et al.

465

CHAPTER 18

Multireference coupled cluster method based on the Brillouin– Wigner perturbation theory ˇ a´rsky1, Jirˇ´ı Pittner1 and Ivan Hubacˇ2 Petr C 1

J. Heyrovsky´ Institute of Physical Chemistry, Academy of Sciences of the Czech Republic, Dolejsˇkova 3, 18223 Prague 8, Czech Republic 2 Department of Chemical Physics, Faculty of Mathematics and Physics, Comenius University, 84215 Bratislava, Slovakia

Abstract This chapter deals with quantum chemical applications in which the electronic state of the system treated is degenerate or quasidegenerate. In chemistry this situation occurs in problems in which an accurate description of a large part of the potential energy surface is required to be treated on the same footing. Commonly used approaches such as the Møller – Plesset theory, configuration interaction, coupled cluster methods CCSD and CCSD(T), or DFT may fail in such cases. Instead, the multireference coupled cluster (MR CC) approach is recommended as the method of choice. A survey of MR CC methods of different kind is given, with a brief description of their main features, merits and drawbacks. Special problems of the MR CC methods, size extensivity and the intruder state problems, are discussed. A variant of the MR CC theory, based on the Brillouin – Wigner type resolvent, is presented. A few applications demonstrate its performance and utility in practical chemical applications. 18.1 INTRODUCTION Recent progress in the development of quantum chemical computational methods [1] and quantum chemical software provided a user-friendly tool to chemists, helping them to explain many problems met in a chemical laboratory. The most frequent task is to determine the optimum structure of molecules, relative energies of reaction components, heats of formation, vibrational frequencies, energies of activation and energies of ionization. In all these cases, a standard theoretical approach is to perform first References pp. 479 – 481

466

Chapter 18

a Hartree – Fock (self-consistent field, SCF) calculation for obtaining a set of molecular orbitals (MOs). In the next step, the main defect of the Hartree – Fock approach is eliminated, which is the neglect of dynamic electron correlation. This is done by assuming the wave function for the ground state (or a particular excited state) as a linear combination of singly, doubly and even higher excited configurations. The respective expansion coefficients are fixed in one way or another. The most popular options are the variational approach in a form of the configuration interaction (CI), the perturbational approach by means of the Møller– Plesset (MP) theory and coupled cluster (CC) expansion. Alternatively, if computational feasibility is favored over the rigor of strictly ab initio treatments, density functional theory (DFT) can be applied. In DFT the two-electron part of the Hamiltonian is replaced by an adjusted coulomb and exchange-correlation functional and the MOs so obtained give electron density which takes account of electron correlation. All these approaches can be referred to as ‘single-reference’ approaches because the starting (‘reference’) wave function is taken in the form of a single Slater determinant or a symmetry adapted configuration state function constructed from several Slater determinants. Besides the availability of the user-friendly software, also the progress in the chemical education of undergraduates should be noted. Basic quantum chemistry and its applications to chemical reactivity and spectroscopy is contained in a curriculum which gives an opportunity to students to acquire the necessary basic knowledge and skill for performing routine quantum chemical calculations on their own. This favorable state in the field of theoretical chemistry may lead to a false impression that all important problems in quantum chemistry have already been solved and that the reliability and feasibility of quantum chemical calculations will increase in parallel with the anticipated continuing progress in the construction of computers in the near future. Unfortunately, there is a type of applications that does not conform to this optimistic outlook. It concerns problems in which an accurate description of a large part of the potential energy surface is required to be treated on the same footing. If the respective part of the surface is associated with a chemical bond breaking, the singlereference methods may fail, unless an enormously large CI or CC expansion is used, which prevents applications to all but the smallest molecular systems. We show an example of such a failure in Section 18.2. A more accessible solution to this problem is provided by a class of ‘multireference’ methods. In Section 18.3, we give a survey of such most commonly used methods that are based on the coupled cluster expansion and comment briefly on their merits and drawbacks. A version of the multireference coupled cluster method based on the Brillouin–Wigner perturbation theory (MR BWCC) is a method developed by the present authors. In Sections 18.4 and 18.5, we describe briefly the essence of the method and its properties and present the working equations for its computer implementation. A few selected MR BWCC applications are presented in Section 18.6 and in Section 18.7 we summarize our experience with the method in its present form and note on the prospects for its next development. 18.2 SINGLE-REFERENCE VERSUS MULTIREFERENCE METHODS It is generally accepted that the family of coupled cluster methods, CCSD, CCSD(T) and CCSDT represents the most sophisticated class of single-reference methods that are

Multireference coupled cluster method based on the Brillouin –Wigner perturbation theory

467

commonly used for chemical applications. Unfortunately, even these sophisticated methods fail in some instances, as it has been demonstrated for potential curves of diatomic molecules on many occasions (see, for example, Ref. [2]). We present it in Fig. 18.1 for the F2 molecule because of our MR BWCC study of that molecule [3]. The potential energy curves shown in Fig. 18.1 were obtained with the cc-pVTZ basis set. The CCSD curve overestimates the dissociation energy by a factor of about 2. CCSD(T) gives excellent result for geometries near to the equilibrium bond distance but a typical failure of the method is seen at intermediate and large interatomic distances. Only CCSDT and MR BWCCSD give correct potential curves for the whole range of interatomic distances. At the equilibrium bond length the MR BWCCSD energy is closer to the CCSD energy than to the CCSDT energy. This is caused by the fact that in the wave function at the equilibrium F – F bond length a single configuration dominates and hence MR BWCCSD represents only a slight improvement over CCSD. Inclusion of higher excitations by means of T3 clusters in CCSDT is more effective than extension of the reference space to two configurations at the CCSD level. This is reflected in the obtained dissociation energies De of F2: 1.40 eV with MR BWCCSD and 1.60 eV with CCSDT, to be compared with the experimental value of 1.66 eV [4]. Unfortunately, the cost of CCSDT calculations prevents their application to larger molecules. By using the ACES program [5] and the Xeon 4 2400 MHz/2 GB RAM, computer time (no use of symmetry was made) for a single-point calculation was 3 min for two-reference BWCCSD and 14.6 h for CCSDT.

−199.2

−199.25

E , a.u.

−199.3

−199.35

2

3

4

5

6

R (F – F) , a.u. Fig. 18.1. Potential energy curves for the F2 molecule calculated by CCSD (dotted line), CCSD(T) (staggered line), CCSDT (dashed line) and MR BWCCSD (solid line) methods. The cc-pVTZ basis set was used.

References pp. 479 – 481

468

Chapter 18

The origin of such a failure is well understood. The ground-state configuration of F2 at the equilibrium bond length is (1pg)4(3sg)2. As the interatomic F –F distance increases, the energy gap between the (1pg)4(3sg)2 and (1pg)4(3su)2 configurations decreases to become zero at the dissociation limit. Evidently the (1pg)4(3sg)2 configuration becomes a poor representation of the electronic state of F2 at large interatomic distances. An obvious way how to treat this defect is to include highly excited configurations in the CI expansion or to include T3 and higher clusters in the CC expansion. Since an approximate inclusion of triples by means of CCSD(T) is not applicable in this case, and a rigorous treatment of T3 clusters by means of CCSDT is too costly, it is more profitable to apply a ‘multireference’ approach in which instead of the regular SCF procedure a multiconfiguration SCF (MC SCF) is used. In the case of F2 it means to include both the (1pg)4(3sg)2 and the (1pg)4(3su)2 configurations. This ensures that all parts of the potential curve are treated on the same footing. As Fig. 18.1 shows, extension of CCSD to two-configuration CCSD is successful in providing a realistic shape of the potential curve with the correct dissociation limit.

18.3 OVERVIEW OF MULTIREFERENCE CC METHODS Over the years, many different CC-based approaches for treatment of systems with static correlation have been proposed: 1. The treatments based on single-reference CC expansion as CCSDtq [6], CASCC [7,8], and the recent state-specific method [9], which are based on a single (possibly formal) Fermi vacuum and include higher amplitudes with part of their indices confined to an active space. This has the advantage of retaining the formal simplicity of the singlereference CC ansatz, while being able to describe processes like dissociation, where one needs to consider higher than double or triple amplitudes. Since the active space is typically very small, the number of the higher amplitudes can be kept manageable, resulting in a reasonable disc storage and computational cost. However, so far we are not aware of a really efficient implementation of these techniques, which would allow us to compute chemically interesting molecules. The disadvantage of these approaches is that they do not treat all reference configurations on an equal footing and the results may be dependent, e.g. on the reference chosen as the Fermi vacuum. 2. The externally corrected CCSD and reduced multireference CCSD methods developed by Paldus and coworkers [10 – 13] make use of external data, typically from a CI calculation, to correct the CC calculation for multireference effects. They can be subdivided into two basic classes: energy corrected and amplitude corrected [14]. The amplitude-corrected methods analyze the CI wavefunction and take a selected set of T3 and T4 amplitudes from it and include them in the CC calculation, while using the standard CC energy formula. On the other hand, the energy-corrected methods use socalled asymmetric energy formula [14] to correct a posteriori the energy obtained from a standard CC calculation. Recently, the authors developed a general scheme that covers the state-universal multireference as well as amplitude-corrected CCSD methods [15]. A general disadvantage of any externally corrected CC method is that a potential

Multireference coupled cluster method based on the Brillouin –Wigner perturbation theory

469

formulation of the analytic gradients for such a method becomes difficult, since the quantities which have been taken from the external source (like CI coefficients) depend on the molecular geometry and their derivatives would be needed. 3. The method of moments coupled clusters developed by Kowalski and Piecuch [16 – 18] belongs essentially also to the general scheme of energy-corrected methods described above, since it gives an a posteriori correction of a CC, EOMCC, or stateuniversal MR-CC energy based on an external wave function of CI or PT type. In principle, if the external wave function were full CI, this correction would recover the exact full CI energy, regardless on the level of approximation used as an input. In practice, the authors have shown a remarkable improvement of the CCSD(T) method, when applied to describe, e.g. bond breaking and dissociation of molecules like F2 or N2. As a disadvantage of this technique one can think of the lack of rigorous size-extensivity of the corrected energy. However, at least for some applications this should not be a serious obstacle, since on the basis of some test calculations the errors are expected to be small. 4. The spin-flip EOMCC method, suggested by Spilchenko and Krylov [19], describes a special class of systems exhibiting static correlation. It treats the state which has MS , S and has multireference character as an spin-changing excitation from the high spin MS ¼ S component, which is well described by a single reference. This technique is not limited to CC method; it has been implemented also in the DFT framework. It is particularly suited for treatment of diradicals or triradicals; probably it is not able to describe whole sections of potential energy surfaces, where the nature of the wavefunction changes from single reference to multireference one. 5. We should also mention the conceptually very interesting general two-body cluster expansion [20 – 22], which—if a practicable numerical treatment can be found—would include the static correlation implicitly in a FCI-like manner. The assumption on which this approach is based is still a subject of discussion. Recently, some arguments against its validity [23] have been given. Even if a rigorous solution could not be found, the idea still might be a source of some useful approximations. 6. The Fock space multireference CC methods and the intermediate Hamiltonian techniques (see e.g. Refs. [24 –29] and references therein), as well as closely related similarity transformed EOMCC [30 – 33] are methods particularly suited for calculation of excited/ionized states with a multireference character. Recently, a Brillouin – Wigner formulation of Fock space CC has also been derived [34]. As a disadvantage, we see particularly the very complex structure of the working equations resulting from these formalisms. This obstacle could be overcome by some automatized implementation tools like the ‘Tensor Contraction Engine (TCE)’ [35], however, at least some of the methods have also very big computational demands, since solutions for many sectors of the Fock space are required. 7. The Hilbert space multireference CC (see e.g. Refs. [36 –40]), based on the Jeziorski –Monkhorst ansatz for the wave operator [36]. This ansatz can be either combined with the standard (Rayleigh – Schro¨dinger) Bloch equation, or with the Brillouin –Wigner Bloch equation (cf. Section 18.4), or with a linear combination of both [41]. Recently, an important progress has been achieved by Li and Paldus [15,42,43], who generalized the Jeziorski – Monkhorst formulation to arbitrary incomplete model spaces. This brings two advantages—computational savings due to smaller number of References pp. 479 – 481

470

Chapter 18

amplitudes and decrease of the probability of occurrence of the intruder state problem. The Hilbert space MRCC methods are highly relevant to the subject of this review and will be mentioned in more detail in Section 18.4. 18.4 MULTIREFERENCE BRILLOUIN – WIGNER COUPLED CLUSTER METHOD Our interest in Brillouin – Wigner perturbation theory was stimulated by our finding [44] that this theoretical tool proved very useful in the scattering theory. The fundamental equation, known in the scattering theory as Lippmann – Schwinger equation, expresses the scattering operator as T ¼ V þ VB0 T

ð1Þ

which can be rewritten in the Brillouin – Wigner perturbation series as T ¼ V þ VB0 V þ VB0 VB0 V þ VB0 VB0 VB0 V þ · · ·

ð2Þ

This can be compared with the MP expansion, E ¼ kF0 lVlF0 l þ kF0 lVQ0 VlF0 l þ kF0 lVQ0 VQ0 VlF0 l þ kF0 lVQ0 VQ0 VQ0 VlF0 l 2 kF0 lVQ0 VlF0 lkF0 lVQ0 VQ0 VlF0 l þ · · ·

ð3Þ

well known in the electronic structure theory and commonly used for the evaluation of MP2, MP3, MP4, and higher order ground state MP energies. In Eqs. (1) – (3) V stands for the perturbation H ¼ H0 þ V

ð4Þ

and in the scattering theory it has the meaning of the interaction potential between the scattering particles. Eqs. (2) and (3) also differ by the form of the propagator B0 or Q0 : As already noted, in scattering theory use is made of the Brillouin – Wigner perturbation theory and therefore the propagator is defined as B0 ¼

1 2 lF0 lkF0 l 1 2 H0

ð5Þ

where 1 is the exact energy and H0 is the unperturbed Hamiltonian with ground state F0 : MP theory is based on the Rayleigh – Schro¨dinger perturbation theory and the propagator is defined therefore as Q0 ¼

1 2 lF0 lkF0 l E0 2 H 0

ð6Þ

where E0 is the unperturbed ground state energy. In contrast to the MP theory, interaction V in the Lippmann –Schwinger equation is strong and cannot be considered as a small perturbation. Hence, the scattering amplitude

Multireference coupled cluster method based on the Brillouin –Wigner perturbation theory

471

cannot be calculated by means of Eq. (2) order by order, as we are used to do so in the MP theory. Instead, Eq. (1) is used in the form T ¼ ð1 2 VB0 Þ21 V

ð7Þ

where contributions to all orders are summed up by means of the matrix inversion. We tried to exploit the merits of the Brillouin – Wigner for the evaluation of the correlation energy. On going from the Rayleigh –Schro¨dinger to the Brillouin – Wigner perturbation theory the ‘renormalization’ terms in Eq. (3) drop and we obtain E ¼kF0 lVlF0 l þ kF0 lVB0 VlF0 l þ kF0 lVB0 VB0 VlF0 l þ kF0 lVB0 VB0 VB0 VlF0 l þ · · ·

ð8Þ

with B0 defined by Eq. (5). Next we tried various approximations to Eq. (8) to obtain a practical method for calculations but all our attempts resulted [45] in one or another known version of the CI method [46 – 50]. We concentrated therefore on the development of the multireference coupled cluster method based on the Brillouin – Wigner perturbation theory. Derivation of working equations for the MR BWCC theory is beyond the scope of this chapter and therefore we refer the interested reader to our earlier papers [40,41,51 – 54]. Here we only present an outline of the method to explain its essence. As it is usual in the multireference coupled cluster theory we are using the concepts of the effective Hamiltonian and Bloch equation (see, for example, Ref. [55]). Let us consider a complete model space spanned by M reference configurations, so that the model function for the ground state is expressed as

C0P ¼

M X

m¼0

Cm Fm

ð9Þ

Such a wave function was used for the two-reference MR BWCCSD calculation of F2 in Section 18.2, where F0 was the (1pg)4(3sg)2 configuration and F1 the (1pg)4(3su)2 configuration. Higher number of reference configurations has been employed in the study of IBr [56] and oxygen molecule [57]. The superscript P is used to indicate that C0P is a ‘projected’ wave function. Its relation to the exact wave function C0 is provided by the (so far) unknown ‘wave operator’

C0 ¼ V0 C0P

ð10Þ

Our task is then to find an effective Hamiltonian which would give us exact energy from the following equation H eff C0P ¼ 10 C0P

ð11Þ

instead of solving rigorously the Schro¨dinger equation with the exact Hamiltonian and the exact wave function H C0 ¼ 10 C0 References pp. 479 – 481

ð12Þ

472

Chapter 18

As it is usual in the MR CC theory, we also used the wave operator in a form suggested by Jeziorski and Monkhorst [36]

V0 ¼

M X

m¼0

eTðmÞ lFm lkFm l

ð13Þ

which is then substituted on the right side of the Brillouin –Wigner analogue of the Bloch equation [58]

V0 ¼ 1 þ B0 V V0

ð14Þ

where V is the perturbation in the MP partitioning of the Hamiltonian (4) and B0 is defined by a multireference generalization of Eq. (5) which reads as B0 ¼

X lFq lkFq l q.M 10 2 Eq

ð15Þ

Once we know the wave operator through Eqs. (13) and (14), we are ready to evaluate the matrix elements of the effective Hamiltonian in the basis of reference configurations because the effective Hamiltonian is given by the following relationship (see, for example, Ref. [55]) H eff ¼ PH V0 P

ð16Þ

where P is the projection operator onto the model space P¼

M X

m¼1

lFm lkFm l

ð17Þ

We arrange the M2 matrix elements kFm lH eff lFn l in a square matrix, diagonalize it, and save the lowest root as the updated ground state energy to be used for a new resolvent (Eq. (15)). The updated resolvent is then used for M single-reference BWCCSD treatments for F1 ; …; FM for obtaining new Heff matrix elements. Actually, the theory yields slightly modified single-reference-like equations, in particular additional Brillouin –Wigner specific terms are included to their right-hand sides. Construction of a new Heff matrix completes a cycle of the iterative procedure for obtaining 10. Our approach may be referred to as MR BWCCSD because only T1 and T2 clusters are allowed to enter the exponential ansatz (13). The T1 and T2 amplitudes may be obtained by a small modification of the CCSD T1 and T2 equations [40]. 18.5 INTRUDER STATES AND SIZE EXTENSIVITY For avoiding the intruder state problem in BWCCSD we pay the price of losing size extensivity. Unfortunately, for the purpose of chemical application of MR BWCCSD, the size-extensivity error is not tolerable. For example, the MR BWCCSD value of the singlet –triplet gap in twisted ethylene exhibits an error of 3 kcal/mol with respect to full CI [40]. Even worse case is the dissociation of the F2 molecule, where MR BWCCSD

Multireference coupled cluster method based on the Brillouin –Wigner perturbation theory

473

Table 18.1 Size extensivity problem in different MR CC methods Method

Resolvent

Active space

Disconnected diagrams

Size extensivity

MR BWCCSD uncorrected MR BWCCSD uncorrected MR BWCCSD last iteration corrected MR BWCCSD last iteration corrected MR BWCCSD iteratively corrected MR BWCCSD iteratively corrected MR CCSD Kucharski – Bartlett [37] MR CCSD Kucharski – Bartlett [37] MR CCSD Jeziorski –Monkhorst [36] MR CCSD Jeziorski –Monkhorst [36] MR CCSD Mukherjee [38]

BW

Complete

Yes

No

BW

Incomplete

Yes

No

BW

Complete

Yes

(No)a

BW

Incomplete

Yes

(No)a,b

BW

Complete

No

Yes

BW

Incomplete

No

Yesb

RS

Complete

No

Yes

RS

Incomplete

No

Yesb

RS

Complete

No

Yes

RS

Incomplete

No

Yesb

RS

Complete

No

Yes

a

A small size-extensivity error can be expected. C-conditions for the incomplete model space [42,59] have to be employed.

b

overestimates the dissociation energy by 17 kcal/mol, but the error can be reduced to less than 0.5 kcal/mol by the size-extensivity correction (see below) [3]. Extent of the size-extensivity error depends on several factors: type of the resolvent, active space, derivation of working equations for amplitudes and the presence or absence of disconnected diagrams in the diagrammatic representation of working equations. In Table 18.1 we show some examples. As in the single-reference CC approaches, the working equations for amplitudes in the MRCC approach are obtained by projection of the Schro¨dinger equation on the manifold of excited configurations (for a review see, e.g. Ref. [1,60,61]). This can be done in two different ways. Assume the CCSD approach and denote the manifold of singles and doubles as SD. Then the Schro¨dinger equation for the ground state CC wave function HeT l0l ¼ EeT l0l

ð18Þ

can be projected on the manifold of singles and doubles as follows kSDlHeT l0l ¼ EkSDleT l0l References pp. 479 – 481

ð19Þ

474

Chapter 18

Alternatively, we can multiply by e2T from the left and do the projection in the Baker – Campbell – Hausdorff manner [1] kSDle2T HeT l0l ¼ kSDle2T EeT l0l

ð20Þ

For single-reference CC approaches the two equations, Eqs. (19) and (20), give the same set of diagrams and therefore the same set of equations for the amplitudes. It is profitable that only connected diagrams are obtained which in accordance with the linked cluster theorem guarantees size extensivity. With the MRCC approaches the situation is more complicated. First of all, for state-universal MRCC methods it is necessary to remove explicit occurrence of the energy E from the basic equation (analogy of Eq. (18)). One thus starts from the Bloch equation instead [55,60], which is, however, equivalent to the Schro¨dinger equation. The projection onto kSDle2T still gives only connected diagrams, whereas the plain kSDl projection gives both connected and disconnected diagrams and it therefore does not apparently ensure size extensivity. Hence the Baker –Campbell –Hausdorff approach may be viewed as the method of choice. There is a problem with this approach, however. The price one has to pay for removal of energy from the Bloch equation is the occurrence of a term quadratic in the wave operator, which leads to products of amplitudes as TðmÞTðnÞ; where the indices m and n denote different reference configurations. Explicit formulas for evaluation of these coupling terms have been recently derived by Paldus and his coworkers [42,62], however, they are rather complicated, their implementation tedious, and evaluation computationally expensive. Their code is able to treat any number of reference configurations and incomplete model space, but is limited to rather small molecules. It seems therefore preferable to develop methods based on an MR analog of Eq. (19), since they lead to coupling terms, which are substantially easier to implement. Kucharski, Balkova´ and Bartlett [37,63 –65] proceeded along that way. They identified and dropped disconnected diagrams and achieved size extensivity. An implementation of this method for open-shell singlet two-reference case has been coded in the ACES program by Szalay [66]. General implementation for more than two reference configurations has been developed recently by Pittner [41], also within the ACES2 program. Moreover, based on a certain formula derived by Kowalski and Piecuch [18] in the context of the method of moments coupled clusters [16], it can be shown [67] that a method based on Eq. (19) is equivalent to the Jeziorski –Monkhorst approach based on Eq. (20) in the case of complete model space, thanks to an exact cancellation of the disconnected diagrams arising from HeT and the coupling terms. In contrast, BWCCSD requires only one set of amplitudes for one reference configuration at a time and this feature makes the method easily amenable to tasks with more than two reference configurations and makes the extension of the method to connected triples (MR BWCCSDT) feasible. Elimination of the size-extensivity error in the MRBWCC theory was suggested by Hubacˇ and Wilson [68]. The Brillouin –Wigner perturbation theory has been known notoriously as a method not furnishing size extensivity. There was therefore every reason to believe that the size-extensivity error in MR BWCCSD originated from the use of the Brillouin – Wigner type resolvent (15) instead of the Rayleigh – Schro¨dinger

Multireference coupled cluster method based on the Brillouin –Wigner perturbation theory

475

type resolvent. Hubacˇ and Wilson found the following identity for the two resolvents X lFq lkFq l X ðE0 2 10 ÞlFq lkFq l X lFq lkFq l ¼ þ 1 0 2 Eq E0 2 E q ð10 2 Eq ÞðE0 2 Eq Þ q.M q.M q.M

ð21Þ

where M is the number of reference configurations. Since the BW resolvent on the l.h.s. is not size extensive, and the RS-type resolvent in the first term of the r.h.s. is size extensive, the last term on the r.h.s. cannot be size extensive. The idea was to find all terms corresponding to that last term and eliminate them in the equations for amplitudes. In practice we did it in several ways. In the simplest one [52] the MRCC procedure is iterated as described in Section 18.4 without paying any attention to the presence of size-inextensive terms. These are only eliminated in the last iteration. The amplitudes so obtained are free from disconnected diagrams. From them the Heff matrix is constructed and its diagonalization gives the corrected energy. It should be realized that the size-extensivity error is not completely eliminated in this way. The amplitudes so obtained are not converged and they cannot be expected to give a rigorous result. A more sophisticated approach [41] is based on a continuous transition between Brillouin –Wigner and Rayleigh – Schro¨dinger perturbation theory. Computationally, a step-wise correction is performed in each CC iteration, finally leading to converged amplitude equations in the rigorously size-extensive Rayleigh – Schro¨dinger limit. Absence of disconnected diagrams in the working equations alone is not a sufficient condition for size extensivity. In addition to that a stronger condition of a complete active space (CAS) is required. By a CAS we mean all configurations that can be formed by a given number of N electrons in a given number of n orbitals. For two electrons in two orbitals the CAS is formed by four configurations, unless the contribution of two singly excited configurations is vanishing on symmetry grounds, in which case the CAS is represented by the ground and doubly excited configuration. For example, for carbenes the important reference configurations are those that represent two electrons in the active space consisting of HOMO and LUMO orbitals f1 and f2 located on the carbenic center. Four spin-unrestricted reference configurations can be formed: F1 ¼ ðf1 Þ2 ðf2 Þ0 ; F2 ¼ ðf1 Þ0 ðf2 Þ2 and F3;4 ¼ ðf1 Þ1 ðf2 Þ1 ; where the last two differ by the spin of the two electrons. However, except for systems of C1 symmetry only two of these, F1 and F2 ; are required for the description of the singlet ground state. There have been several attempts to avoid the CAS restriction and achieve size extensivity for an incomplete model space (see, for example, Refs. [69,70]). The methods based on abandoning the intermediate normalization condition resulted in an excessively complex formalism and as far as we know were actually never implemented. Only recently this problem has been attacked successfully by Li and Paldus, who introduced so called C-conditions for the amplitudes of internal excitations [15,42,43,71]. When these conditions are incorporated into the MRCC amplitude equations of the Jeziorski – Monkhorst method, the solutions with a general incomplete model space become also exactly size extensive. The C-conditions are, however, not limited to the original Jeziorski –Monkhorst formulation [36], but rather they apply to any Hilbert-space MRCC References pp. 479 – 481

476

Chapter 18

method, including BWCC, which can be corrected for size extensivity also in the incomplete model space [59]. As with the single-reference methods, size extensivity does not guarantee size consistency (a proper dissociation limit in the case of potential curves). The reference space must be carefully chosen to contain important configurations at any point of the potential curve as it was recognized long time ago [72]. But this is not a specific problem of MR BWCCSD. 18.6 PERFORMANCE OF THE MULTIREFERENCE BRILLOUIN – WIGNER CC METHOD AND APPLICATIONS Computationally, MR BWCCSD may be viewed as a set of weakly coupled singlereference CC calculations. Equations for the determination of amplitudes depend only on a set of amplitudes (from a previous iteration) of a single reference configuration and coupling between reference configurations is only ensured by the energy, which is eigenvalue of the effective Hamiltonian matrix. This simplicity of the computational scheme is beneficial for the feasibility of calculations. The computer time for a MR BWCCSD run with M reference configurations is about M times higher than it is for a standard CCSD run. This permitted us to perform MR BWCCSD calculations on a desktop PC even for molecules that have been of interest to organic chemists and biochemists. The largest molecules we treated were tetramethyleneethane (TME) [73] H2C

CH2 C

H2C

C CH2

and benzyne [74]. The purpose of the calculations on TME was to bring some theoretical evidence on the ordering of the lowest singlet and triplet states. From the EPR spectra of TME measured in a matrix it was concluded (for references to experimental work, see Ref. [73]) that TME has D2d structure at the triplet state. This was at variance with the gas-phase negative ion photoelectron spectra showing that the triplet state is about 2 kcal/ mol above the singlet state. The MR BWCCSD/cc-pVDZ calculations support the latter assignment. At the optimum triplet geometry (dihedral angle of 498) the calculated S – T energy gap is 1.2 kcal/mol, and the triplet is predicted to lie above the singlet at any dihedral angle. However, a firm interpretation of the observed triplet state of the matrixisolated TME would need more rigorous calculations, including triple excitations and using a larger basis set. In parallel, we also used the MR CCSD method by Kucharski and Bartlett [37]. The results obtained by the two methods were close in absolute value, showing that the deficiency of the MR BWCCSD method—its size inextensivity—may be eliminated by an a posteriori correction (correction done in the last CC iteration only). The calculations on benzyne [74] had a biochemical experimental background (see Ref. [74]). The natural products shown in Fig. 18.2 displayed high activity against a

Multireference coupled cluster method based on the Brillouin –Wigner perturbation theory O HO

sugar

NHCO2Me

H

MeSSS

O

HO MeSSS

sugar

477

NHCO2Me

H sugar

calicheamicin

esperamicin Me

OH

O

HN

COOH O OMe

OH

OH

O

dynemicin Fig. 18.2. Natural products with a high activity against a number of tumor cell lines.

number of tumor cell lines, both in vivo and in vitro, but unfortunately also showed unacceptable toxic effects in both animal and human trials. The active site is believed to be the hex-3-ene-1,5-diyne moiety, which yields 1,4-didehydrobenzene (benzyne) by the Bergman cyclization reaction (Fig. 18.3). Benzyne abstracts hydrogen from a saccharide phosphate backbone to form benzene, which denaturates the DNA and ultimately causes cell death. Benzyne is a diradical species and it is therefore difficult to calculate it by standard single-reference methods. The objective of MR BWCCSD calculations was to provide reliable data on the heat of reaction and enthalpy of activation for the Bergman reaction, which would be used as standards for less sophisticated (and less demanding) calculations, and to show the accuracy attainable by the present state-of-the-art techniques in designing new antitumor agents with an enediyne-like structure. The entries in Table 18.2 show that the MR BWCCSD method with the a posteriori correction for size extensivity improves greatly the CCSD results and it gives results in good agreement with the MR CI and experimental data. Good performance of CCSD(T) is probably fortuitous in this case.

H H

H

H

H

H H

H Reactant Fig. 18.3. Bergman reaction.

References pp. 479 – 481

TS

H

H

H

H

IR

478

Chapter 18 Table 18.2 Enthalpy of activation and heat of reaction of the Bergman reaction.a Method CCSD MR BWCCSD MR CI CCSD(T) Experiment a

DH†298 (kcal/mol)

DH298 (kcal/mol) 0

38.2 32.7 29.4 27.6 28.2

27.5 12.9 10.3 10.1 8.5

The cc-pVTZ basis set was used in all calculations; for details see Ref. [74].

We have also studied the automerization barrier in cyclobutadiene, where the transition structure has a diradical character [75] and the singlet – triplet gaps in alkylcarbenes [76]. Besides calculations of organic compounds, we have assessed the accuracy of the MR BWCCSD technique on benchmark systems, like the insertion of Be into hydrogen molecule [77], and compared its accuracy with other multireference methods. We have also employed the MR BWCCSD method for accurate treatment of diatomic molecules like F2 [3] or IBr [56], which required four reference configurations to span the complete model space. We have also extended calculations by our method to low lying excited states and up to eight reference configurations, used for calculation of the oxygen molecule [57]. The isoelectronic, but heteronuclear NF molecule has been investigated along the same lines [78]. A series of carbide diatomics, CaC, ZnC, BeC, and MgC, were another challenge for the MR BWCCSD to treat systems of multireference nature caused by near-degeneracy effects. The task was to examine theoretically [79,80] the competing 3S2 and 5S2 states. In Table 18.3 we present results for CaC. Single-reference methods CCSD and CCSD(T) give a wrong order of the two states, predicting the triplet to lie 10.0 and 0.4 kcal/mol above the quintet state. Experimentally, the ground state is not known but MR BWCCSD with the a posteriori size-extensivity correction and MR CI predicts consistently that the triplet is the ground state. The former gives Te ¼ 0:9 kcal=mol for the quintet and the later gives 3.0 kcal/mol.

Table 18.3 Spectroscopic constants of the CaC molecule Method ROHF CCSD CCSD(T) MR BWCCSD MR ACPF MR CI MR CI þ Q

De (kcal/mol)

˚) re (A

ve (cm21)

ve xe (cm21)

1.2 35.5 46.1 45.8 48.0 48.6 49.3

2.838 2.368 2.383 2.357 2.365 2.364 2.364

464.0 464.2 468.3 465.1 462.4 461.3

3.67 2.39 2.06 2.65 2.37 2.51

Multireference coupled cluster method based on the Brillouin –Wigner perturbation theory

479

18.7 SUMMARY We do not advocate MR BWCCSD as the method of choice for general use. Instead we wanted to show that in applications where a system with a quasidegenerate electronic state is to be calculated, it may be profitable to use a multireference approach as a more advantageous alternative to the single-reference treatment. This applies particularly to cases where the number of reference configurations needed is low. The computer time for MR BWCCSD increases about linearly with the number of reference configurations. Hence, if a regular single-reference CC calculation is feasible for the particular system, also the MR BWCCSD would be likely affordable. Of course, as with any other multireference method, some experimentation with the number of reference configurations would be necessary. In its present form and at the level of MR CCSD, the MR BWCC method seems to be slightly less accurate than large scale MR CI calculations with some sort of a posteriori correction [52,68] (such as the averaged coupled pair functional) for retaining size extensivity. However, MR BWCCSD has good prospects for its further development. Its inherent feature of being not size extensive is now well understood [41] which is a good starting point for the development of more accurate sizeextensivity corrections. Also, development of the MR BWCCSDT computer program is in progress. First results obtained using an approximation of MR BWCCSDT for calculation of vibrational frequencies and anharmonicities of O2 molecule [81] show that inclusion of T3 clusters into the multireference CC expansion greatly improves the accuracy of the method. A profitable feature of MR BWCCSD is also the similarity of the MR BWCCSD working equations to those in CCSD, which enables to formulate analytically the MR BWCCSD energy gradient along the same lines as in CCSD. However, even in its present form the applications performed so far show that MR BWCCSD can be taken as a more economic alternative to MR CI for treatments of systems with a quasidegenerate electronic state. 18.8 ACKNOWLEDGEMENTS This work has been supported by the COST D23 action (grant OC D23.001 of the Czech Ministry of Education). We acknowledge also the support by the Grant Agency of the Czech Republic (grant No. 203/04/0425) and the Grant Agency of the Academy of Sciences of the Czech Republic (Grants No. A4040401 and 1ET400400413). We also acknowledge the support by the Academy of Sciences of the Czech Republic (Project No. K4040110). 18.9 REFERENCES 1 2 3 4

P.v.R. Schleyer, N.L. Allinger, T. Clark, J. Gasteiger, P.A. Kollman, H.F. Schaefer, III and P.R. Scheiner (Eds.), The encyclopedia of computational chemistry, Wiley, Chichester, 1998. X. Li and J. Paldus, J. Chem. Phys., 108 (1998) 637. J. Pittner, J. Sˇmydke, P. Cˇa´rsky and I. Hubacˇ, J. Mol. Struct. (THEOCHEM), 547 (2001) 239. K.P. Huber and G. Herzberg, Molecular spectra and molecular structure, Vol. 4, Van Nostrand Reinhold, New York, 1979.

References pp. 479 – 481

480 5

6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50

Chapter 18 ACES II is a program product of Quantum Theory Project, University of Florida, J.F. Stanton, J. Gauss, J.D. Watts, M. Nooijen, N. Oliphant, S.A. Perera, P.G. Szalay, W.J. Lauderdale, S.R. Gwaltney, S. Beck, A. Balkova´, D.E. Bernholdt, K.-K. Baeck, P. Rozyczko, H. Sekino, C. Huber and R.J. Bartlett. Integral packages included are VMOL (J. Almloef and P.R. Taylor); VPROPS (P. R. Taylor); and ABACUS (T. Helgaker, H.J.Aa. Jensen, P. Joergensen, J. Olsen and P.R. Taylor). P. Piecuch, S.A. Kucharski and R.J. Bartlett, J. Chem. Phys., 110 (1999) 6103. L. Adamowicz, J.-P. Malrieu and V.V. Ivanov, J. Chem. Phys., 112 (2000) 10075. V.V. Ivanov and L. Adamowicz, J. Chem. Phys., 113 (2000) 8503. M. Ka´llay, P.G. Szalay and R. Surja´n, J. Chem. Phys., 117 (2002) 980. X. Li and J. Paldus, J. Chem. Phys., 107 (1997) 6257. X. Li and J. Paldus, J. Mol. Struct. (THEOCHEM), 547 (2001) 69. X. Li and J. Paldus, J. Chem. Phys., 115 (2001) 5759. X. Li and J. Paldus, J. Chem. Phys., 115 (2001) 5774. J. Paldus and X. Li, Collect. Czech. Chem. Commun., 68 (2003) 554. X. Li and J. Paldus, J. Chem. Phys., 119 (2003) 5334. K. Kowalski and P. Piecuch, J. Chem. Phys., 113 (2000) 18. K. Kowalski and P. Piecuch, J. Chem. Phys., 115 (2001) 2966. K. Kowalski and P. Piecuch, J. Mol. Struct. (THEOCHEM), 547 (2001) 191. L.V. Spilchenko and A.I. Krylov, J. Chem. Phys., 117 (2002) 4694. M. Nooijen, Phys. Rev. Lett., 84 (2000) 2108. T.V. Voorhis and M. Head-Gordon, J. Chem. Phys., 115 (2001) 5033. P. Piecuch, K. Kowalski, P.-D. Fan and K. Jedziniak, Phys. Rev. Lett., 90 (2003) 113001. E.R. Davidson, Phys. Rev. Lett., 91 (2003) 123001. J.-P. Malrieu, P. Durand and J.-P. Daudey, J. Phys. A, 18 (1985) 809. L. Meissner, J. Chem. Phys., 108 (1998) 9227. L. Meissner, P. Malinowski and A. Nowaczyk, Chem. Phys. Lett., 381 (2003) 441. L. Meissner, P. Malinowski and J. Gryniako´v, J. Phys. B, 37 (2004) 1. A. Landau, E. Eliav, Y. Ishikawa and U. Kaldor, J. Chem. Phys., 113 (2000) 9905. N. Vaval, S. Pal and D. Mukherjee, Theor. Chem. Acc., 99 (1998) 100. M. Nooijen, J. Chem. Phys., 104 (1996) 2638. M. Nooijen and R.J. Bartlett, J. Chem. Phys., 106 (1997) 6441. M. Nooijen and R.J. Bartlett, J. Chem. Phys., 107 (1997) 6812. M. Nooijen and V. Lotrich, J. Chem. Phys., 113 (2000) 494. N.D.K. Petraco, L. Horny´, H.F. Schaefer, III and I. Hubacˇ, J. Chem. Phys., 117 (2002) 9580. P. Sadayappan, G. Baumgartner, D.E. Bernholdt, R.J. Harrison, S. Hirata, M. Nooijen, R.M. Pitzer and J. Ramanujam: The Tensor Contraction Engine, http://www.cse.ohio-state.edu:/~gb/TCE B. Jeziorski and H.J. Monkhorst, Phys. Rev. A, 24 (1981) 1668. S.A. Kucharski and R.J. Bartlett, J. Chem. Phys., 95 (1991) 8227. U.S. Mahapatra, B. Datta and D. Mukherjee, J. Chem. Phys., 110 (1999) 6171. I. Hubacˇ, in: A. Tsipis, V.S. Popov, D.R. Herschbach and J.S. Avery (Eds.), New methods in quantum theory. NATO ASI series, Kluwer, Dordrecht, 1996, pp. 183–202. ˇ a´rsky, J. Ma´sˇik and I. Hubacˇ, J. Chem. Phys., 110 (1999) 10275. J. Pittner, P. Nachtigall, P. C J. Pittner, J. Chem. Phys., 118 (2003) 10876. X. Li and J. Paldus, J. Chem. Phys., 119 (2003) 5320. X. Li and J. Paldus, J. Chem. Phys., 119 (2003) 5346. ˇ a´rsky, V. Hrouda, V. Sychrovsky´, I. Hubacˇ, P. Babinec, P. Mach, J. Urban and J. Ma´sˇik, Collect. P. C Czech. Chem. Commun., 60 (1995) 1419. ˇ a´rsky, Mol. Phys., 88 (1996) 1137. V. Sychrovsky´ and P. C Z. Gershgorn and I. Shavitt, Int. J. Quantum Chem., 2 (1968) 751. G.A. Segal and R.W. Wetmore, Chem. Phys. Lett., 32 (1975) 556. L.E. Nitzsche and E.R. Davidson, J. Chem. Phys., 68 (1978) 3103. T.H. Dunning, Jr., Chem. Phys., 42 (1979) 249. J.J. Diamond, G.A. Segal and R.W. Wetmore, J. Phys. Chem., 88 (1984) 3532.

Multireference coupled cluster method based on the Brillouin –Wigner perturbation theory 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81

481

J. Ma´sˇik and I. Hubacˇ, in: R. McWeeny, J. Maruani, Y.G. Smeyers, S. Wilson (Eds.), Quantum systems in chemistry and physics: Trends in methods and applications, Kluwer, Dordrecht, 1997, pp. 283– 308. ˇ a´rsky, J. Chem. Phys., 112 (2000) 8779. I. Hubacˇ, J. Pittner and P. C I. Hubacˇ, J. Ma´sˇik, P. Mach, J. Urban and P. Babinec, in: J. Leszczynski (Ed.), Computational chemistry. Reviews of current trends, Vol. 3, World Scientific, Singapore, 1999, pp. 1–48. ˇ a´rsky, I. Hubacˇ, P. Mach, J. Pittner and S. Wilson, in: J. Maruani, R. Lefebvre and E. Brandas (Eds.), P. C Advanced topics in theoretical chemical physics, Progress in theoretical chemistry and physics, Vol. 12, Kluwer, Dordrecht, 2003, pp. 71–118. I. Lindgren and J. Morrison, Atomic many-body theory, Springer, Berlin, 1982. J. Pittner, O. Demel, P. Cˇa´rsky and I. Hubacˇ, Int. J. Mol. Sci., 2 (2002) 281. ˇ a´rsky and I. Hubacˇ, Int. J. Quantum Chem., 90 (2002) 1031. J. Pittner, P. C J. Ma´sˇik, P. Mach and I. Hubacˇ, J. Chem. Phys., 108 (1998) 5671. J. Pittner, X. Li and J. Paldus, Mol. Phys., in press. J. Paldus, in: S. Wilson, G.H.F. Diercksen (Eds.), Methods in computational molecular physics. NATO ASI series, Plenum, New York, 1992, pp. 99 –194. J. Paldus, in: G.L. Malli (Ed.), Relativistic and correlation effects in molecules and solids. NATO ASI series, Plenum, New York, 1994, pp. 207 –282. J. Paldus, X. Li and N.D. Petraco, J. Math. Chem. 35 (2004) 215. A. Balkova´, S.A. Kucharski and R.J. Bartlett, Chem. Phys. Lett., 182 (1991) 511. A. Balkova´, S.A. Kucharski, L. Meissner and R.J. Bartlett, J. Chem. Phys., 95 (1991) 4311. A. Balkova´ and R.J. Bartlett, J. Chem. Phys., 101 (1994) 8972. P.G. Szalay and R.J. Bartlett, J. Chem. Phys., 101 (1994) 4936. J. Pittner, K. Kowalski and P. Piecuch, to be published. I. Hubacˇ and S. Wilson, J. Phys. B, 33 (2000) 365. L. Meissner, S.A. Kucharski and R.J. Bartlett, J. Chem. Phys., 91 (1989) 6187. L. Meissner and R.J. Bartlett, J. Chem. Phys., 92 (1990) 561. J. Paldus and X. Li, J. Chem. Phys., 118 (2003) 6769. G.C. Lie and E. Clementi, J. Chem. Phys., 60 (1974) 1288. ˇ a´rsky and I. Hubacˇ, J. Phys. Chem. A, 105 (2001) 1354. J. Pittner, P. Nachtigall, P. C ˇ a´rsky, P. Stampfuß and W. Wenzel, Collect. Czech. Chem. Commun., 68 O. Rey-Puiggros, J. Pittner, P. C (2003) 2309. ˇ a´rsky and I. Hubacˇ, J. Chem. Phys., 112 (2000) 8785. J.C. Sancho-Garcı´a, J. Pittner, P. C O. Demel, J. Pittner, P. Cˇa´rsky and I. Hubacˇ, J. Phys. Chem. A, 108 (2004) 3125. ˇ a´rsky, Chem. Phys. Lett., 386 (2004) 211. J. Pittner, H. Valdes-Gonza´lez, R.J. Gdanitz and P. C ˇ a´rsky and A. Mavridis, Int. J. Quantum Chem., Int. J. Quantum Chem., in S. Kardahakis, J. Pittner, P. C press. ˇ a´rsky, A. Mavridis and I. Hubacˇ, J. Chem. Phys., 117 (2002) 9733. I.S.K. Kerkines, J. Pittner, P. C V.I. Teberekidis, J. Pittner, P. Cˇa´rsky, A. Mavridis and I. Hubacˇ, Int. J. Quantum Chem., in press. J. Pittner and O. Demel, J. Chem. Phys., in press.

References pp. 479 – 481