The nature of Van der Waals bond

Z.B. Maksi6 and W.J. Orville-Thomas (Editors) Pauling's Legacy: Modern Modelling of the Chemical Bond

665

Theoretical and Computational Chemistry, Vol. 6 9 1999 Elsevier Science B.V. All rights reserved.

T h e nature o f V a n der W a a l s b o n d Grzegorz Chatasifiski a, Malgorzata. M. Szcz~niak b, and Stawomir M. Cybulski c aDepartment of Chemistry, University of Warsaw, Pasteura 1, 02-093 Warszawa, Poland bDepartment of Chemistry, Oakland University, Rochester, Michigan 48309, United States of America CDepartment of Chemistry and Biochemistry, Miami University, Oxford, Ohio 45056, United States of America

1. I N T R O D U C T I O N Linus Pauling has had a lasting impact on two areas of the theory of molecular interactions [1,2]. He was the first to recognize the importance of weak interactions, such as hydrogen bonding to biology, and he pioneered the theoretical investigations of these interactions. His second powerful idea was the role of molecular shape in such diverse areas as crystal packing and a complementarity of enzyme-substrate interactions. By his own admission Pauling never used a computer in his studies [3]. In a recent reminiscence he said: "I am sure that if I had been relying on a computer to make most of the calculations, some of these ideas, which have in fact turned out to be important, would not have occured to me." By contrast, today's applications of quantum chemistry to these systems are almost exclusively carried out on a computer, taking advantage of unprecedented advances in solving the Schr6dinger equation. The computational methods have become reliable enough so as to provide us with potential energy surfaces, which are suitable for simulations of the behavior of molecular clusters. The wealth of numerical data generated in this way enhances the need for simple and intuitive rationalizations, which were the hallmark of Pauling's work. In this contribution we intend to show that such a program can be most appropriately accomplished by dissecting the interaction between molecules into a few basic components. These components, electrostatic, induction, dispersion and exchange energies, are as important conceptually to the theory of weak interactions as the notions of ionic, covalent, and metallic bonding were to Pauling's approach to chemical bonding. We will demonstrate that any sensible modeling of intermolecular forces must rely on these four basic components. The primary concept of the theory of intermolecular interactions is the intermolecular potential energy surface (PES) [4]. The intermolecular force is defined as the negative gradient of this energy. In the quantum mechanical setting, the PES has its origin in the BornOppenheimer approximation as a potential energy for the motion of nuclei. In the case of molecular interactions, it is usually reasonable to treat the interacting species as rigid bodies in their equilibrium geometries. Such an approach provides the PES, which depends only on the intermolecular degrees of freedom. If the intermolecular forces are strong enough to cause a

666 sizable deformation of the monomers, the intramonomer degrees of freedom should also be included. The crudest approximation is to separate these two effects into the interaction of deformed monomers and the effect of monomer geometry relaxation. A rigorous approach should solve the equations for the nuclear motions on the PES, which parametrically depends on the inter- and (at least some) intramonomer degrees of freedom. At either level of complexity the ultimate task of theory is to interpret and predict the experimental measurements, such as the dissociation energies into free monomers, intra- and intermolecular vibrational frequencies in Van der Waals complexes, integral and differential cross sections for scattering experiments. In all such measurements the effect of the interaction is gauged with relation to relaxed monomers. 2. FUNDAMENTAL I N T E R A C T I O N ENERGY C O M P O N E N T S As the understanding of chemical bonding was advanced through such concepts as covalent and ionic bond, lone electron pairs etc., the theory of intermolecular forces also attempted to break down the interaction energy into a few simple and physically sensible concepts. To describe the nonrelativistic intermolecular interactions it is sufficient to express them in terms of the aforementioned four fundamental components: electrostatic, induction, dispersion and exchange energies. Typical closed-shell molecules feature an uneven distribution of charge density. This charge distribution may be described by permanent multipole moments, or in the simplest approximation, by point charges distributed over the molecules. When two molecules are far apart the moments interact via the Coulomb law, giving rise to a long-range part of the PES (the one, which decays as some inverse power of the intermolecular distance, R-n). This interaction is essentially exhibited in three different forms, primarily derived by London [5] (cf. also Refs. 6 and 7). First, it is the direct electrostatic interaction between multipole moments. Second, since the interaction perturbes the involved species, the multipoles are created and/or modified. These modifications give rise to what is usually referred to as the induction interaction energy. There is also a third type of interaction, called the dispersion interaction, that molds the long range shape of the potential. Whereas the first two interactions may be viewed as clearly related to classical electrodynamics, the dispersion energy emerges because of the quantum mechanical nature of the molecular world. Its semiclassical model, consisting of the interaction of "instantaneous multipoles" due to the fluctuating positions of electrons, may serve only as a simplified visualization. In fact, the dispersion energy is an electron correlation effect that takes place in the area between the interacting monomers. Bringing molecules closer to the Van der Waals minimum, roughly defined by the Van der Waals radii, generally leads to two additional effects. One is an alteration of all three interactions discussed earlier. Indeed, electrostatic, induction and dispersion effects can no longer be modeled by a multipole expansion. This is because the electron clouds begin to overlap giving rise to some exponential damping of the long-range interactions. This overlap brings about a secondary modification, which may easily be built into the classical models. A primary effect on the PES, however, is due to the appearance of a new, purely quantummechanical factor. It has been known under the name of the exchange repulsion or HeitlerLondon (HL) exchange energy. It is due to such quantum mechanical considerations as the delocalization of the electrons and their indistinguishability. It corresponds to "resonance"

667 integrals in Pauling's description of chemical bonds and is related to the Pauli exclusion principle. [8-11] Among closed shell molecules the fully occupied monomer orbitals repel rather than invite the electrons of the partner. This effect of blocking the electronic space has been related to such notions as atomic radii and molecular shapes - the concepts to the development of, which Pauling significantly contributed. It means that in the first approximation the atoms and molecules may be viewed as surrounded by rigid spheres or contours, which make impenetrable borders for partner monomers. In a better approximation, the shapes of monomers are soft, featuring an exponential rise and decay with R. The simplest picture of the intermolecular interaction in the complex is, thus, the following: the monomers stick together because they are glued by electrostatic forces in the form of direct electrostatic energy, the induction energy and the dispersion energy. They cannot collapse onto each other because of the exchange effect (the Pauli principle), which prevents the occupied orbitals from overlapping. This picture may be simulated in a variety of ways, from a very simple one (a direct electrostatic model restrained by hard spheres, which envelope molecules [12,13]), through combined ab initio, model and semiempirical approaches [14-16] to the most sophisticated based on the rigorous solutions of the Schr6dinger equation, in particular by means of the symmetry adapted perturbation theory (SAPT) [17-20]. SAPT, while preserving the backbone of the classical perturbation theory approach of London [5], supplements it with refinements and couplings of the electrostatic, induction, dispersion and exchange terms. This results in a beautiful theory, which blends mathematical rigor with the conceptual simplicity of common intuitions. Before proceeding any further, we will briefly outline the ab initio theory of intermolecular interactions.

3. AB INITIO A P P R O A C H TO I N T E R M O L E C U L A R F O R C E S All the information that is needed to describe a particular intermolecular interaction in mathematical and physical terms is included in the Schr6dinger equation for a system under consideration HW=EW

(1)

where H is the total Hamiltonian of the Van der Waals complex, and W and E are its wave function and energy, respectively. Unfortunately, except for a few systems, solving the Schr6dinger equation is a very difficult task demanding ingenious approximations and state-ofthe-art computational techniques. The first important approximation known as the BornOppenheimer approximation separates the motions of nuclei and electrons. Light electrons are assumed to instantly follow any infinitesimal change of the nuclear positions. This approximation brings about the crucial concept of the potential energy surface (PES) [4,15]. In general, one may visualize nuclei moving due to forces exerted on them by this potential field. In the case of intermolecular interactions, as long as monomers are kept rigid, we can view them as "moving" with respect to each other upon the intermolecular PES. In other words, the PES provides a playground for the interactions of molecules. Knowledge of a PES not only means the understanding of how two molecules interact, but also provides us with a means of simulating and predicting the properties of dimers, trimers and large aggregates of atoms and molecules.

668 The calculation of PESs is most easily accomplished by evaluating the interaction energy, Eint, which is defined as the difference between the energy of the dimer, EAB, and the energies of the monomers, E A and E B Eint = EAB - E A - E B

(2)

The procedure chosen to calculate Eint must ensure that electronic energies of the dimer and monomers are evaluated in a consistent manner [19,21-24]. It should be stressed that this requirement is absolutely crucial, as no method at present can in practice yield EAB, E A and E B energies with an absolute error smaller than Eint. Therefore, Eq. (2), which defines the interaction energy does not offer so simple a computational approach as might be expected at first glance. Two notorious inconsistencies to be alleviated in practice are: basis set inconsistency (same basis set expansion or numerical grid for A, B, and AB must be used, otherwise the basis set s u p e r p o s i t i o n error (BSSE) arises [21,23,24]) and the size inconsistency (a theory to describe AB must guarantee a correct dissociation into A and B, at the same level of theory [25]). Even if the algorithm is correct, it has the disadvantage of giving no direct insight into the nature of the interaction.Therefore, from the very beginning of quantum mechanics, scientists were inclined to solve the Schrrdinger equations through a perturbation procedure that would provide the interaction energy directly and would offer insights into its physical nature and functional form. Originally, the fundamental components related to the Coulomb interaction of permanent, induced, or instantaneous multipole moments were recovered with the aid of the Rayleigh-Schrrdinger (RS) perturbation theory [5-7]. The introduction of exchange effects occurred independently, and originated from the variational theory of the chemical bond formulated by Heitler and London. The development of a unified, rigorous treatment proved difficult because of the symmetry problem related to the indistinguishability of electrons. The classic RS formalism assigns electrons to either of two monomers, and this approximation ("polarization approximation") works as long as the monomers stay far apart [26]. In the intermediate and short-range, though, the antisymmetry of the total dimer wavefunction cannot be efficiently recovered by means of the RS perturbation theory, which becomes divergent. Attempts to deal with the increase-of-symmetry problem were eventually successful and led to a variety of SAPT formalisms [ 11,18,27]. Along with the rapid development of many-body techniques to cope with electron correlation, SAPT today provides us with a rigorous framework, as well as a detailed description of quantum theory of intermolecular forces. The details of SAPT are beyond the scope of the present work. For our purposes it is enough to say that the fundamental components of the interaction energy are ordinarily expanded in terms of two perturbations: the intermonomer interaction operator and the intramonomer electron correlation operator. Such a treatment provides us with fundamental components in the form of a double perturbation series, which should be judiciously limited to some low order, which produces a compromise between efficiency and accuracy. The most important corrections for two- and three-body terms in the interaction energy are described in Table 1. The SAPT corrections are directly related to the interaction energy evaluated by the supermolecular approach, Eq.(2), provided that many body perturbation theory (MBPT) is used [19,28]. Assignment of different perturbation and supermolecular energies is shown in Table 1. The power of this approach is its open-ended character. One can thoroughly analyse the role of individual corrections and evaluate them with carefully controlled effort and desired

669 Table 1 Decomposition of two- and three-body supermolecular (S-MPPT) interaction energies. The contents of S-MPPT terms is described and the leading SAFF terms are indicated in square brackets. S-MPPT

SAFF

Physical interpretation Two-body

AESCF

AE(2)

Electrostatic energy between SCF monomers Exchange repulsion between SCF monomers

AESCF def

Mutual polarization restrained by exchange [13ind,rJ

0 t3( disp

Dispersion energy arising between SCF monomers (2nd order)

(20).

Electrostatic-correlation energy (2nd order). Intra-monomer correlation correction to e~1~ (2) AEdet

1. D e f o r m a t i o n - i n t r a - c o r r e l a t i o n .

e(22)

['~ind,r]

2. Deformation-dispersion. re (30) j1 t disp_ind rt(2)

Ar-,exch

1

(20) 1 9 Exchange-dispersion r[Eexch_dispJ

2. Exchange-intra-correlation AE(3), AE(4), etc.

Higher-order correlation corrections as for AE(2) Three-body

AE SCF

eeI~xch F

1. SE component: single exchanges between monomers 2. TE component: all monomers are involved in the exchange - (20)

e(30)1

SCF-deformation nonandditivity 113ind,r, ind,r I - (20)

1. Exchange-dispersion nonadditivity [Eexch_disp] 2. Exchange-intra-correlation nonadditivity

AE(2)

e(22) 3. Deformation-intra-correlation nonadditivity [Wind,r]

4. Deformation-dispersionnonadditivity rE (30) 1j t disp_ind AE(3 )

..(30) e-disp

1. Dispersion nonadditivity, accounts for triple-multipole terms 2. Higher-order correlation corrections as for AE(2)

AE(4), etc.

Higher-order correlation corrections as for AE(3)

670 precision. The price is involvement of many arbitrary choices and approximations, which at some point may render the method cumbersome and subjective. A proliferation of corrections with an increasing order of the perturbation treatment and additional approximations to circumvent a divergent character of theory may pose problems. Another shortcoming is the limitation to intermediate and long-range interactions. Fortunately, all these drawbacks may be overcome by using the relationship to the supermolecular energies [19,28]. The supermolecular interaction energies defined in Eq.(2) have the advantage of being free from such problems and represent the most reliable interaction energy values for adjusting the final potential. Only a concerted use of the direct supermolecular approach and SAPT may lead to very accurate interaction energies. It also leads to a deeper understanding of the origin and behavior of the components to the interaction energy. Below, we will describe some very recent advances.

3.1

Exchange repulsion versus molecular shape

"The Dutch physicist J.D. van der Waals found that in order to explain some of the properties of gases it was necessary to assume that molecules have a well def'med size, so that two molecules undergo strong repulsion when, as they approach, they reach certain distance from one another. [...] It has been found that the effective sizes of molecules packed together in liquids and crystals can be described by assigning Van der Waals radii to each atom in the molecule. The Van der Waals radius defines the region that includes the major part of the electron distribution function for unshared [electron] pairs." Cf. Fig. 1.A [2]. Strong repulsion described by Pauling in his "General Chemistry" textbook [2] originates from the electron exchange effect that arises between the electron clouds of interacting monomers. The quantum mechanical approach allows us to link the concept of the effective size of molecules or a molecular shape to the exchange effect. We may visualise this shape in terms of the exchange repulsion, which is sensed by a rare gas (RG) atom moving around a molecule. Several researchers have attempted to define alternative pictorial representations [29,30]. In Figs. 1.B and 2 we show the representation proposed by our group [31,32]. In the displayed diagrams, a contour is drawn in polar coordinates that were used to define the motion of a probing RG atom around the molecule. However, the distance R is replaced by the value of the exchange-repulsion energy. In that manner, the contour exhibits local decreases and increases of repulsion, which indicate local concentrations and depletions of electron density in the diffuse region. Several examples of such plots for different molecules are in Figs. 1.B and 2. Let us examine the contour for the C12 molecule in Fig 1.B. How does this compare with the drawing in Pauling's book? Both pictures reveal depletions in the middle of the C1-CI bond. However, the Van der Waals radii approximation fails to predict a characteristic flattening of the electron distribution at the chlorine ends along the internuclear axis. One can see in this and all other drawings that the assumption of a spherical distribution around an atom within a molecule is not justified. The deviation from a spherical distribution may be small, but it is sufficient to influence the shape of the PES, and to determine the equilibrium structures of Van der Waals complexes with non-polar species for, which induction effects are negligible. The PES of RG-C12 exhibits two kinds of minima, one for a T-shaped and one for a coUinear from. The PES for C1F has three minima, and two of them (at the C1 side and the T-shaped one) are due to a reduction of the repulsive effect, cf. Fig.2. The global minimum corresponds to the collinear form RG-C1F. The PESs of RG-HF and RG-HC1 systems feature two collinear minima, and it is clear from Fig. 2, that the one at the

671

Figure 1. A) A chlorine molecule, illustrating the difference between Van der Waals radius and covalent radius (from Pauling's book [2]). B) The exchange repulsion contour of C12, obtained for the He-C12 complex, and defined by two polar coordinates, measured from the center of mass: (Energy, O). The contour is the image of the chlorine molecule shape, detected by a rare-gas atom [32]. C) Relief map of the negative Laplacian of the charge density, -V2p, of the chlorine molecule. One can notice three depletions of the electron charge density: in the region perpendicular to the bond, and along the interatomic axis, for the colinear approach at the both ends.

672

F

F

H

Ct

C

Cl

0

H

Figure 2. The exchange repulsion contours for several molecules, obtained for interactions with rare-gas atoms, and defined by two polar coordinates measured from the center of mass: (Energy, O) [31,32]. The contours are the images of molecules' shapes, probed by structureless atoms. In contrast to plots that show isoenergetic regions, these contours reveal an enhanced anisotropy. Convex and concave regions indicate, respectively, the areas of increased and reduced exchange repulsion.

673 halogen atom end must be due to a reduction in repulsion. The RG-CO complexes are skew T-shaped, again in agreement with the depletion visible in Fig 2. In general, any non-polar species (such as methane, methyl group, N 2, etc.) will be preferably hosted in niches within the exchange repulsion, which are displayed by these contours. Finally, it is important to add that the regions of reduced (or enhanced) repulsion may be directly related to the regions of depletion (or concentration) observed in the diffuse part of the Laplacian of electron density plots proposed by Bader [33]. You may see this in Fig.2, where the negative of V2p for C12 has been drawn. In such a way, scrutiny of-V2p indicates the locations where we may expect reduced repulsion, and which are thus the most favorable sites for a nucleophilic partner to attach. It is important to stress that the analysis of the Laplacian must carefully distinguish between the relative charge concentration in short range and the relative charge concentration in the diffuse region - as they sometimes do not coincide, and only the latter is relevant. In quantitative modeling of PESs the description of the molecular shape as a superposition of atomic components remains an attractive approach, but it is clear from the earlier discussion that it must be extended to accommodate two important factors. The atomic shape is not a rigid, but rather a soft, exponentially decaying electronic charge cloud. In addition, it should be anisotropic with the anisotropy depending not only on the atom itself, but also on its partner in the chemical bond. 3.2

Dispersion as the intermonomer correlation effect The traditional view of dispersion effect introduces the notion of instantaneous multipoles. Such multipoles arise, when for an infinitesimal portion of time, electrons are frozen in their positions, and thus become oriented with respect to positively charged nuclei. A momentary multipole at, say, monomer A, is created, which induces a momentary multipole at monomer B. The resulting electrostatic interaction defines the dispersion effect. In this way fluctuations of negative electron charges of A become correlated with the electron charges of B. It is clear that this effect may be related to the intermonomer electron correlation phenomenon. This simple picture, however, is reasonable as long as we assume that the electrons of A and B do not mingle with one another, or more precisely, do not penetrate each other's occupied space. Mathematically, it means that the condition ra+rb
674

X2

b)

a)

4

"

i

"-----"

'

"-~','~

"

.'.'/1

"~ " ,

Xl

He1 _

~

He z -

' ~

- 2.8

;

2 "

OiO

'"""'

I ' , ' _ " > " ",;_/t~J// II \'. '_, ,' ,. .' ,, -_ ~' ." ~; '" /I \ ~ S / /

" 2.8

"if///// ~/// /

I o !- ', ",'" :-7,1~\~'"////// i

X2

li

I

9

9

!////

ill

," ~!111 \ ~

lilli

I

.

-

-4 x2

-2

0

2

0

2

Xl

4

x2

(''~ -~1

,

i

'

IIlll

d)[/~1/i;"

i

"''"""

4

Illl

~'~ 9 ' ,e.~.,t,1 m ,, :,,, ,,,.. Ii I ~ .

_

I', , ',,~ '-',,:I I ~ , ,'.', ,,.!41 _1 ', ', ,:....':~!~ hi-. t t ~ "dill'

#11 I I ! I

"\'..':: ~I<~ "-''

.,:;:":

2

l~

Ul

~

/ -4

, -4

, -2

~

, 0

x 1 -4~ 2

4

~ ".''dill

',," ,'7It [( -4

i -2

4

X1

Figure 3. Plots of the dispersion function for the equilibrium He dimer from Ref. 36. The electron coordinates, x 1 and x 2 are defined in fig. (a). The accurate dispersion function (calculated with Gaussian geminals) is shown in plot (b). If the dispersion term is reduced to the dipole-dipole component, one obtains plot (c). If every electron is allowed to use the complete dimer basis set as well as the bond functions one obtains plot (d).

675 correlation will now have to account for a sizeable contribution from the "cusp": a singularity predicted by the Coulomb law when two electrons occupy the same point of the configuration space. Correlated, but non-overlapping multipoles, are not capable of describing this effect properly. In Fig. 3 depicting two interacting He atoms, we show the dispersion ansatz to the electronic wavefunction: exact and reduced to the dipole-dipole interaction [36]. One can see that whereas the dipole-dipole interaction reproduces the general pattern of this ansatz (its "butterfly shape") it fails to reproduce a steep cliff, which arises in the region where the electron coordinates coincide. The related portion of the dispersion energy must be evaluated by some method, which does not refer to the multipole approximation. In addition, the rigorous approach should use explicitly correlated basis functions centered not only at the nuclei, but also off them, in particular along the Van der Waals "bond" line. In many more approximate calculations, however, it is enough to supplement a reasonable basis set (that is, the one, which properly describes the first few moments and polarizabilities of the monomers) with some small set of functions located in the mid-bond region [36-38,19]. They are termed bond functions. The bond functions efficiently account for the cusp in the region where it is the most important for inter-monomer correlation. The effect of bond functions is also shown in Fig. 3. One immediately notices a considerable improvement in the cusp region.

3.3

Induction, charge-transfer and SCF deformation

At first glance the induction energy is a very "classical" term, and conceptually a simple one. A multipole moment on A induces another multipole moment on B, and eventually they interact electrostatically. Different mechanisms behind induction interaction in a variety of molecular settings were elucidated and composed into an elegant theory [7]. It is this simple, however, so long as the condition of applicability of multipole approximation, ra+rb>R holds (cf. our discussion of the dispersion interaction). Not only does it keep the electrons of A away from the electrons of B (thereby obeying the Pauli exclusion principle), but also they are screened off the nuclei of the partner molecule (thus avoiding the singularity due to coinciding position of an electron and a nucleus). At close range, however, electrons of A mingle with B, and since the Pauli principle between monomers is not imposed, they may deeply penetrate each other's occupied space. There is nothing to prevent an energetically beneficial charge flow from one monomer into a manifold of symmetry-forbidden states of another monomer. Not restrained by the Pauli principle, such a flow may result in multiply occupied spinorbitals of this m o n o m e r - a process dubbed "unphysical charge transfer" [22,39]. A good example of the unphysical charge transfer was provided for the He-Li + interaction [39]. The induction interaction leads then to a transfer of the electrons from Li + to the occupied space of He. In effect, the induction interaction in the minimum region, although a well defined quantity, describes some huge, unphysical interaction. A proper remedy is enforcing the Pauli exclusion principle during the course of induction interaction, as in the routine supermolecular calculations for the dimer. The resulting induction effect is then accompanied and restrained by the exchange effects, and termed the SCF-deformation energy. However, dissecting the SCF-deformation energy into the induction and exchange-induction contributions would be largely dubious because, if the induction alone is not a physical quantity, neither would be the accompanying exchange effect. In addition, both are much larger than the SCF-deformation effect. From the preceding discussion it is clear that the induction energy may function as an approximation to the SCF-deformation, but not on its own. Therefore, one should carefully

676 monitor the legitimacy of this approximation in every case and for every geometry of a complex. Although a rigorous extracting of the classical induction effect cannot be performed, the induction multipole expansion may still be useful as a template for modeling the SCFdeformation effect and its calibration. Another consequence is the impossibility of a rigorous separation of the component termed "charge-transfer" energy. This type of energy appeared in early partitioning treatments of the interaction energy [40]. The charge-transfer effect referred to that part of the induction effect, which required a description of the charge cloud modification at A, not only with orbitals of A (that was referred to as the "true" induction), but also with the orbitals of B. In the language of the Valence Bond theory the latter would correspond to the ionic structures, e.g. A§ -. In fact, even for so called extended basis set calculations the induction component of the total electronic wave function is not efficiently reproduced unless "exchange" and "ionic" terms are allowed for. However, an obvious problem with this definition is that in the real dimer the borders of monomers are not defined and it is impossible to distinguish, which charge cloud deformation is already a charge-transfer and which is just a simple induction effect. One can do it for finite basis set treatments, but such a definition is strongly basis set dependent, and the "ionic" and "exchange" components of the wavefunction vanish in the limit of a complete basis set. These ambiguities led some researchers to propose that monomers in the dimer are defined through orthogonal orbitals. Then, however, the charge-transfer contribution grew very large and became a dominant contribution, absorbing portions of the electrostatic and exchange components [41]. In view of these problems the question is whether we really need a separate charge-transfer term. In a rigorous quantum mechanical treatment the induction, charge-transfer, and exchange effects are blended together in the SCF-deformation term and in the correlation corrections to it. Modeling is usually based on multipole expansion of the induction energy, restrained by exchange and penetration effects. This is a legitimate approach since asymptotically the deformation effect is correctly described by the induction multipole expansion. Some researchers, however, attempt to include a separate charge-transfer term, and argue that this results in a better numerical fit [42]. Indeed, ab initio calculations are often slowly convergent if ionic and exchange components are left out, and a similar problem may pertain to modeling. Despite all the aforementioned problems, the role of a separate CT term in the modeling of a potential should be further investigated.

3.4 Example 1. Ar-CO2: dispersion bound complex One reason for studying the RG-molecule complexes is that a RG atom serves as a structureless probe of the interaction proclivity of a molecule. In addition, such complexes are simpler than molecule-molecule interactions since they lack a strong long-range electrostatic factor, which may obscure other energy components. Not surprisingly, the RG-molecule systems provide convenient examples for studying weak intermolecular interactions for both the theory and the experiment. To describe the composition of the PES from well defined and physically sensible fundamental parts, we briefly describe our recent results for the Ar-CO 2 complex [43]. The coordinates for this complex are shown in Fig. 4.A The geometry of CO 2 was kept rigid. The interaction energy perturbation components are also shown in Fig. 4, for R = 3.7 A. The dominant short-range part, the HL-exchange energy, l~exch filL) has a strong angular dependence with a minimum at 90 ~ and maxima at 0 ~ and 180 ~ A wide niche in the HL-

677

A)

E, pE h 20000

'

'

I

'

'

_- ~e(HL) 15000 ~-.. ~, exch

i

'

'

i

'

'

i

'

'

I

c o [O"~'~.... R .... Ar

HL

1ooOO ooo

' /

~/ . ~

1

II

-~" <~-X---AE

/ -.54

~_AEScF

o~

'

~.-~..•

~- -x.._

---=5000 -10000

~'" F

1 "

,

E(: 0> I

0

,

E}20,>r ,

30

)

60

,

,

i

,

90 O, deg

,

I

,

,

I

120

,

,

180

150

B)

E, pE h 8000 ~ x ~ 2 ) 6000 4000 ~ . 2000

f

-/~i .,~....AE(Z)

117"~~.

-2000 F-

exch

. . o ~/ .-- .0~ *

-4000 ~ - ' ~ / o..
-6000 L , / ~ 1 7 6 0 30

,0

8(20)

/

_.---]

../-'/'~"

-J

~,~.. 0

q |

<>\ ~,,'" '-~ o

, , , , , , i , , ,-o 60 120 150 90 O, deg

,..~ 180

Figure 4. The system of coordinates and angular dependence of the interaction energy components in ArCO2 at R=7.0 ao; A) at the SCF level of theory, B) at the correlated level of theory.

678

B)

HL F_.,

E[~E] h

1250 _

'

'

'

'

I

'

-

'

'

'

I

'

'

'

'

exch

I

'

Ar

1200

'

'

'

I

'

'

'

'

I

'

'

'

'

z>

i

1150 1100 1 I'

O

1050 1000

,

-3

,

,

,

C

I

-2

,

,

~

,

I

-1

,

,

,

,

I

O

,

0 A[A]

,

,

,

I

1

,

,

,

,

I

2

,

-

,

,

,

3

Figure 5. A) Relief map of the negative Laplacian of the charge density, -V2p, of the carbon dioxide molecule. filL) as a function of A, the B) Repulsive effect of the HL-exchange term, •exch, coordinate of Ar atom moving along the molecular axis at R=3.75 A in a parallel fashion.

679 exchange plot at 60~ ~ is observed. As pointed out in Sec.3.1 [19] such regions of reduced repulsion are directly related to the electron density depletions in the outer, diffuse region of a molecule. The plot of the Laplacian of electron density [44] (cf. Fig. 5.A) indicates a depletion around the central C atom. At the same time, an Ar atom moving along the O=C=O molecule experiences a distinct reduction of exchange repulsion fight in the middle, Fig. 5.B. Towards the O atoms, beginning at 40 ~ the HL-exchange rapidly rises. Of particular interest is that no additional depletions or concentrations of electron density appear at the outskirts of the terminal O atoms, i.e., there is no evidence of the lone pairs at the O atoms. This is, again, in accord with the plots of Laplacian of electron density in Fig. 5.A. The anisotropies of the electrostatic (e~ls~ induction

- (20)) [eind,r, terms

9

SCF.

SCF-deformation (Z~def) and perturbation

are quite similar and qualitatively reciprocal to the anisotropy of the

HL-exchange energy, see Fig.4.A.

e~ls~ is reduced for this complex to the penetration part

only and has no long range component, ~ind,r' o(20) which includes the interaction of the quadrupole and higher CO 2 moments with induced moments of Ar (restrained by overlap effects and truncated at the second order), represents a poor quantitative approximation to 9

_SCF Al~de f .

This is

SCF

not unexpected since Al~def also encompasses exchange effects that prevent the Pauli principle violation in the course of interaction. The correlation components are shown in Fig 4.B. AE(2) is dominated by the dispersion ..(20) The electrostatic-correlation component, e~l,2r~ is practically negligible. The term e.disp. _(2) accompanying exchange corrections included in AEexch are small but not negligible in the vicinity of O atoms.

3.5 Example 2. Water dimer: introducing electrostatics The water dimer is the most important H-bonded system. The major attractive contribution to the interaction energy of the water dimer is the electrostatic effect. It dominates over other attractive terms, such as the induction and dispersion energies, and it is the most anisotropic. To discuss the properties of the fundamental components in the water dimer case we chose to demonstrate the angular dependence of various terms in the dimer geometry derived from the cyclic configuration of a trimer (see Fig 6). SCF HE e~10) and ~ d e f , are shown in The components of the SCF interaction energy, eexch, Fig.6. It is clear that the anisotropy is determined by the electrostatic contribution. This term generates both a minimum for the H-bond geometry and the barriers for the H-to-H and O-to-O configurations. The exchange energy behaves differently. It is greatly enhanced at the H-to-H structure, from which it gently rolls down as ct increases to reach a broad minimum in the region of the O-to-O form. The effect of the exchange term on the total anisotropy is consequently small, although quantitatively the attraction of the electrostatic energy is considerably reduced. The correlation contributions are of secondary importance. They primarily include the weakly anisotropic dispersion term (cf. Fig. 6) and a smaller but much more anisotropic electrostatic correlation term. In the PES modeling, the latter may be safely absorbed in the total electrostatic interaction. In the above description we need no reference to lone electron pairs. Nevertheless, one 9

680

E, m E

O-to-O

H-bond

H-to-H

h

40 Z ~ scF

30

,s

--.

20

"

e"~ e$

,

"

. . . . . . . . .3. . .0. .

_

j 10

9

0~,'~-~

%

,,~

e (~2) ~ m B m m m ~

)

]~" 9

. .,~. ~

~ , . a,,~

A 1:;'ScF

a, .tY,- -

atsp

"-'~-" def

-20

,

0

,

m

% 1

30

,

,

1

.~, ,

,

60

I

90

,

,

I

120

1 150

(~, deg

Figure 6. Angular dependence of various two-body interaction energy components in the cyclic planar H20 trimer at R=3.0 ]k.

681

Figure 7. A). Relief map of the negative Laplacian of the charge density, -V2p in the plane perpendicular to the molecular plane [33]. -V2pexhibits two local maxima and two saddle points. The maxima may be identified with localized lone electron pairs. B) Electrostatic equipotential contours for the water molecule, in the plane perpendicular to the molecular plane [ 13]. The contour spacing is 2000 cm-1 (for a unit test charge), and the distance units for the horizontal and vertical axes are atomic units (au). C). The exchange repulsion contour of H20 derived from the He-H20 complex at R(He-O)=3.5 ,h,, defined by two polar coordinates: (Energy, o), and drawn in the same plane as the Laplacian. The contour is the image of the water molecule shape, detected by a rare-gas atom [31]. The regions of lone electron pairs are indicated with arrows, but no apparent sign of their presence is observed. The lone pairs electron concentrations are not diffuse enough to show up in the van der Waals minimum region.

682 may see local concentrations of the electron density in this region in the Laplacian of charge density plots, cf. Fig.7.A [33]. Similarly, one can notice local maxima of the electrostatic field along the directions of the lone pairs Fig.7.B [13]. On the other hand, no distinct local maxima are seen on the exchange-repulsion contour, cf. Fig.7.C [31 ]. The reason is that these features are fairly short ranged, and thus may affect the directionality of interaction only at shorter intermolecular separations. For example, in the water dimer the O--O distance is relatively small, and the familiar tilted structure of the dimer is rationalized by the OH group of one monomer pointing towards the lone pair of the other. On the other hand, when the moieties are farther apart, as in the Na+...water complex, the cation attaches itself along the C 2 axis of water and not along one of the lone pairs. A rare gas atom also seems to take little notice of a lone pair, cf. Fig.7.C. 3.6 General considerations The two preceding examples of the ArCO 2 complex and the water dimer well illustrate some general trends in shaping the PES following the behavior of its fundamental components. The short-range shape is always molded by the exchange-repulsion, primarily described by HL eexch. The long-range part of the PES is expected to follow the shape of the dispersion term if t;disp _(20) is the primary binding factor. If the complex binding is dominated by the electrostatic or induction interactions, then they assume the role of determining the long-range behavior. Interestingly, the anisotropies of all these fundamental components are only weakly dependent on the intermolecular separation - unlike the anisotropy of the total PES. The shapes of the long- and short-range parts of the PES are thus easily predictable. This is often not the case in the Van der Waals minimum region where we encounter a delicate balance of all components. The geometry of the global minimum of a Van der Waals complex may sometimes be difficult to a priori postulate unless quite accurate calculations are performed. In a majority of cases, however, the picture is simple. The dispersion term displays a relatively weak orientation dependence, whereas the exchange and electrostatic terms (and induction to some extent) are strongly anisotropic. Consequently, in the absence of long range electrostatic interactions the equilibrium geometry is determined by the exchange anisotropy. This is the case of the previously discussed ArCO 2. It is different, however, from interaction of ionic and polar species. If the electrostatic interactions dominate, their anisotropy determines the shape of PES in the minimum region. The latter fact was successfully exploited by the electrostatic model of Buckingham and Fowler [12] and other, similar models [13]. Between these two extremes one may encounter more complex combinations. For instance, the minima on the PES of rare gases with polar molecules often follow the anisotropy of the induction component. The interaction between a RG atom and a polar molecule such as HF can serve as an example [19].

4. MODELING OF PES AND ITS COMPONENTS There may be many different functional representations of PESs based on a variety of mathematical techniques (e.g. polynomials or splines) and different physical models (e.g. atoms-in-molecule approach). Theoretically, all of them should work adequately in simulations of scattering, spectroscopic or thermodynamic properties. However, we believe

683 that their usefulness can be greatly enhanced when a high numerical accuracy of the interaction energy is combined with a sensible representation of the fundamental components. Only then can we combine a full understanding of the nature of the interaction at the molecular level with the predictive power of the spectroscopic, scattering, and thermodynamic theories. The preparation of such potentials will be demonstrated using again ArCO 2 and (H20) 2 as examples. 4.1 Ar-CO 2 The total potential energy was expressed as a sum of the fundamental components, the short-range repulsive term (Vsr), the dispersion energy component (Vdisp) and the induction part (Vind) [43]

V = Vsr + Vdisp + Vind

(3)

The fourth term, the electrostatic interaction, was omitted since one monomer, the Ar atom, has no permanent multipole moments. The Vsr part accounts for the monomers' shape factor in the potential. It mainly includes the exchange effects. However, other effects, such as the overlap-dependent electrostatic interaction, may also be conveniently absorbed in Vsr. Following Buckingham, Fowler and Hutson [14], Vsr may be adequately modeled by a Born-Mayer exponential form: Vsr(R, O) = A exp[-~(O)(R- Rref(O)) ]

(4)

where Rref and 1~are expanded in the Legendre polynomials PL(cOsO)

Rref(O) = ~ Rref,L PL(COSO) L ~(0) = ]~ ~L PL(COSO)9 L

(5) (6)

This form of Vsr is particularly efficient for rod-shaped, elongated molecules. For such molecules a direct expansion of Vsr in terms of PL(COsO) would have a prohibitively slow convergence because of the strong radial-angular coupling of the PES. However, the Rref as a function of O effectively takes care of this problem through a relatively short expansion within the exponential function argument. The dispersion part of PES was represented using the following damped expansion:

r.,disp PL(COSO) even n-4 Dn(~,R)~n,L Vdisp =- ~ ~ Rn n=6 L=0

(7)

where ~disp "-'n,L are the Van der Waals dispersion coefficients, which are defined by the multipole expansion of the dispersion energy through n=10. Since the monomer electron distributions overlap, the resulting quantum mechanical effect is accomplished by the damping functions D n designed by Tang and Toennies [46]

684

Dn([3,R) = 1- exp(-~R) ~ [~(O) R] k k~ k=0

(8)

Finally, the induction interaction Vin d was cast into a similar form as the dispersion part with t-,ind "-'n,L describing the interactions between permanent and induced moments of particular ranks. e~n ~4 r Vind = -

PL(COSO)

~n,L n - 8 L--0

Rn

(9)

In fact, in complexes like ArCO 2 where multipole interactions are fairly weak, the induction interaction is significantly affected by the short-range exchange effects whose role is to prevent the Pauli-principle-violating redistribution of electrons. These effects are related to the molecular shape factor and may be absorbed in Vsr. The described overall procedure and the final form of the potential are applicable to any complex of an atom and a linear molecule. They also may be generalized to more complex systems, composed of polyatomic monomers, as exemplified in the next Section by the water dimer. How efficient is the described representation of the ArCO 2 potential? To answer this question the above PES along with a few empirical potentials have been used to derive a number of properties, such as the ground vibrational state and dissociation energy of the complex, ground state rotational constants, the mean square torque, the interaction second virial coefficients, diffusion coefficients, mixture viscosities, thermal conductivities, the NMR relaxation cross sections, and many others [47]. Overall, the ab initio surface provided very good simulations of the empirical estimates of all studied properties. The only parameters that were not accurately reproduced were the interaction second virial coefficients. It is important that its performance proved comparable to the best empirical surface 3A of Bohac, Marshall and Miller [48]. This fact must be greeted with satisfaction since n o empirical adjustments were performed for the ab initio surface. 4.2

Water

dimer

As a second model potential we shall briefly discuss the PES for the water dimer. Analytical potentials developed from ab initio calculations have been available since the mid seventies, when Clementi and collaborators proposed their MCY potential [49]. More recent calculations by Clementi's group led to the development of the NCC surface, which also included many-body induction effects (see below) [50]. Both potentials were fitted to the total energy and therefore their individual energy components are not faithfully represented. For the purposes of the present discussion we will focus on another ab initio potential, which was designed primarily with the interaction energy components in mind by Millot and Stone [51 ]. This PES was obtained by applying the same philosophy as in the case of ArCO 2, i.e., both the template and calibration originate from the quantum chemical calculations, and are rooted in the perturbation theory of intermolecular forces. Compared to the ArCO 2 complex there is an important new f a c t o r - the electrostatic

685 energy. A water molecule possesses sizeable dipole and quadrupole moments, and thus the electrostatic energy provides the dominant attractive contribution. Since the electrostatic interaction is strongly anisotropic, the orientational dependence of the interaction energy is governed by this term. The nonspherical symmetry of both monomers introduces an additional difficulty. The functional forms of the individual interactions should account for the mutual orientations of two water molecules. Customarily, two approaches are used to incorporate this angular dependence. The first distinguishes the single interaction centers at each molecule, which are connected by the intermolecular vector R=[R,c0], where co specifies its orientation. The orientations of the molecules in the space-fixed coordinate system are described by the two sets of Euler angles {tOa} and {COb}. If the multipole expansion within such a single-centered approach is to be convergent, the already mentioned convergence criterion, ra+rb
Electrostatic energy In the water-dimer potential of Ref.[51 ], each water molecule has three interaction sites centered on atoms. The long-range electrostatic energy assumes the functional form Ves=~., ~., QatTatubQub

(10)

a~A b~B

It includes the interactions of distributed multipole moments Q (up to a quadrupole) labeled t and u. The T matrix provides the Coulomb energy appropriate for particular multipoles and includes the distance between sites a and b and their relative orientations. The short range (penetration) component of the electrostatic energy, in a manner similar to the Ar-CO 2 case, can be absorbed into the exchange repulsion term.

Short-range e n e r g y The exchange-repulsion energy has been fitted to the following functional form [51,52] Vsh = E

E

exp{'tXab[Rab-Pab(~)]}

(11)

a~A b~B

with Rab the distance between sites a and b, ~ab a hardness parameter depending on the pair of sites, and an orientation-dependent parameter describing the effective size of atoms. The orientation dependence of Pab(f~) is given by a simple model that assigns to each site a shape described in terms of orientation-dependent radius p

686

pa(O,*) = Z P~k Clk(O,*) l,k Pab(~) = Pa(~) + Pb( ~ )

(12)

where Cm(O,.) are renormalized spherical harmonics and P~k the related expansion coefficients. Dispersion energy A reasonable multi-site representation for the dispersion energy was developed by Szcz~niak et al. [53]

Edisp = Z

Z

fab(R)C~b PCa6

aeA beB tl forR>ra+rb } fab(R) = Iex~ ra+rb R )12 forR
(13)

where C are Van der Waals site-site coefficients and f is a damping factor, ra and rb are the Van der Waals radii. The C coefficients were obtained by fitting to the second-order nonexpanded dispersion term. Since dispersion has a weaker anisotropy than other terms, it is also quite sufficient to use a single-centered expansion

10 Cn(ll,12,j ;kl,k2) ~klk2 Vdisp = - ~ Z Rn S1112j ((0a,0)b,fD)fn(a) n=6 11,12,j,k1,k2

(14)

where R is the distance between the centers of mass, Cn(ll,12,j;kl,k2) are the anisotropic ~klk2 dispersion coefficients, Sl112j (0~a,~b,C0) are the orientation functions first proposed by Stone [54], which include products of Wigner rotation matrices and 3j coefficients, fn(R) are damping functions of Tang and Toennies [46] (see Eq.(8)). The coefficients Cn(ll,12,j;kl,k2) have been derived by Rijks and Wormer [55]. The multi-center form of the dispersion term is more in the spirit of the other terms, and despite its simplicity proved to perform better in some applications, e.g. in calculations of tunneling splittings between nonsuperimposable forms of the water dimer [56,57]. Induction energy To describe the induction interaction the single-center approach with one-site polarizabilities was found to be the most efficient. Vind has been expressed as

Vind=21---~ ~ AQt~mb fn(Rab)l/2Qb A B#A

(15)

687 and AQtl = - Z

b b a'a Ttu ab fn(Ra,b)l/2 (Q~+AQu) (gt't

(16)

B;eA where T is the interaction function discussed in the context of Eq.(10), Q are multipole moments, and fn(R) are the damping functions. A distributed multipole moment Q at b induces multipole moment AQ at a, which corresponds the polarizabity tensor ot~:~. Due to its coupled form (i.e. each induced moment depends on the induced moments on all the other centers), the Eq. (16) should be solved iteratively. The missing element in this model is the accompanying exchange effect. As mentioned in the preceding section, this effect may be partly accommodated in the short-range repulsive component. Nevertheless, the modeling of the SCF deformation energy is still not quite satisfactory. The potential for the water dimer that was just described is not yet very accurate. Since, however, both the functional form and calibration were derived from ab initio calculations, there is room for well controlled improvements that would follow future more accurate ab initio data. This model provided the geometry of the global minimum in very good agreement with experiment, and a fair account of the second virial coefficient. It should be mentioned that the well-depth of this potential is smaller than the experimental value of 5.4 + 0.7 kcal/mol. However, the best ab initio calculations on the water dimer also consistently predict smaller values. This potential was also used in a number of bound-state calculations on the water dimer and trimer using either variational basis set approach [56] or the rigid-body diffusion Monte Carlo (RBDMC) method [57]. The predicted tunneling splittings in the dimer were in very good agreement with the experimental findings. Upon augmenting it by the three-body induction contribution this potential was also used by Gregory and Clary in the RBDMC simulations of the tunneling dynamics in the water trimer [58], tetramer, and pentamer [59].

5. TRIMERS AND NONADDITIVE EFFECTS So far, our discussion of intermolecular potentials has been limited to pair interactions. In clusters involving more than two monomers the three- and higher-body terms will appear. For example, the total energy of a trimer ABC may be expressed as follows

E(iAB)C =

Z E~) X=A,B,C

+ Z

AEO) X=A,B,C

+

Z

AE(~+/.amAB CArt(i)

(17)

X>Y=A,B,C

where (i) denotes a particular level of theory, e.g. Hartree-Fock (HF) theory, an order of Moller-Plesset perturbation theory, or any other size-consistent treatment of correlation effects, such as Coupled Cluster (CC) Theory. The second, third and fourth terms describe, respectively, the one-, two-, and three-body contributions. The one-body term describes the effects of the geometry relaxation of the subsystem X in the trimer. A two-body term Art(i) term AE(~ describes the pairwise interaction between two monomers, and the z-U~ABC

represents the three-body contribution arising between the relaxed-geometry monomers

688 arranged in the same way as they occur in the complex. Many-body effects are known to affect directly and indirectly various bulk properties of liquids and solids. Binding energies of rare-gas crystals indicate 10% deviations from additivity [60] and are believed to be responsible for different lattice structures than expected from pairwise additivity [61,62]. Many-body effects are particularly important when hydrogen-bonds, ion-polar, or ion-ion interactions are present. For example, the simultaneous description of structural and thermodynamic properties of the three phases of water was found impossible without the explicit consideration of many-body contributions [63]. In ion-water interactions, the nonadditive effects were found to alter the coordination numbers of solvated ions [64], and structures of solvation shells [65]. They also proved to be a crucial binding factor in the beryllium clusters [66] and quartet state of alkali metal trimers [67]. Finally, in some alkaline earth-halogen crystals the induction nonadditivity has been linked to unexpected, layered crystal structures [68]. Although investigations of bulk properties offered a great deal of evidence of the importance of many-body effects, it is the study of clusters, which in recent years have led to a better understanding of these effects. Many-body interactions are much more difficult to retrieve from experimental data or to compute directly by ab initio methods than the pair interactions. Discerning different threebody terms may be accomplished within the framework of SAPT [69,70] and the relationship to supermolecular interaction energy terms has also been obtained [19,71] (cf. Table 1). Similar to the binary interactions, the nature of nonadditivity may be explained in terms of the fundamental nonadditive components: induction, dispersion, and exchange; the fourth one, the electrostatic energy, is always additive. While the absence of the electrostatic part reduces the number of elementary terms, the fact that in the trimer there are three different pair interactions leads to a greater variety of the mutual couplings of fundamental components than in the case of binary interactions. For example, in addition to a pure three-body induction term, there are also mixed terms, such as induction-dispersion, exchange-induction, or exchange-dispersion, which correspond to situations when one pair is involved in one type of fundamental interaction (e.g. induction) and another pair in a different type (e.g. dispersion) [ 19,71,72]. With the advances in cluster spectroscopy the field of Ar2-chromophore clusters has become a convenient laboratory for experimental and theoretical investigations of three-body interactions [73]. Due to the presence of a chromophore these complexes can be studied by IR or microwave spectroscopic techniques. Theoretical techniques are used to predict these spectra by assuming some form of the three-body interaction potential. This approach requires accurate two-body potentials to retrieve the three-body effects from the experimental data. Unfortunately, the need for accurate pair potentials greatly restricts the number of useful chromophores. At the time of this writing only a few pair potentials qualify for this label, the most important being those for: Ar-Ar, Ar-HF and Ar-HC1. For this reason, the related Ar 2HX clusters have been studied in great detail. Owing to the synergistic interaction between theory and experiment, the nature of nonadditivity is close to being well understood in these clusters.

5.1 Ar2-chromophore clusters: exchange and dispersion nonadditivity The three-body exchange nonadditivity, Eel~xch,is conveniently categorized as related to two general mechanisms depicted in Fig.8. The first, where permutation of electrons encompasses all three monomers and the second, where only two are involved. They are referred to as

689

1-e exchange

V

$

B

C

SE

TE

Figure 8. Diagrams representing the Heitler-London exchange terms, single exchanges (SE) and triple exchanges (TE). The electron exchange operators are symbolized with arrows and the interaction operator with a dashed line.

Ar

Ar

,

i

! i

,

! !

", : ,' ' l l !

I !

| ~ H ~ e.o.m X ~ 1~

~

H

' ! I

i ! !

I ! ! I !

! I !

! !

Ar

Ar X - F , C1

Figure 9. Definition of geometrical parameters in Ar2HX (X=F,C1). R is the distance between the center of mass of HX and the middle of Ar 2. O is the angle between the R vector and the HX axis. ~ is a dihedral angle between Ar-Ar and HX axes.

690 triple-exchange (TE) and single-exchange (SE) nonadditivities [72].

SE nonadditivity The exchange repulsion within dimer AB may be viewed as a distortion of charge densities of A and B; the electrons are pushed outward and effectively create an "exchange quadrupole." (This electrical effect will be the most pronounced if A and B are rare-gas atoms). Next, the "exchange quadrupole" interacts with the permanent moment of C, giving rise to the nonadditive effect [74,75]. In modeling, the SE nonadditivity is characterized by a different functional dependence with respect to different pair separation distances [72,75]. It has longrange dependence with respect to multipole interactions, but short-range exponential decay for exchange multipoles.

TE nonadditivity To better appreciate the behavior of this term let us consider two extreme configurations, an equilateral triangle cluster and a collinear arrangement. The exchange repulsion is often interpreted as a distortion brought about by the theoretical process of orthogonalizing occupied orbitals of interacting monomers. If two monomers are already orthogonalized, bringing around a third monomer in a perpendicular approach already requires a somewhat weaker orthogonalizing. The effective exchange repulsion in the trimer will thus be less than three times the exchange repulsion of a dimer. The nonadditive effect is attractive and may be modeled as Eexch= -A exp(-txrl2)exp(-13rl 3)exp(-Tr23)

(18)

An opposite situation arises for a collinear form. Orthogonalizing two monomers "gets in the way" of orthogonalizing the third monomer. Instead of cooperating, the monomers now compete. Eventually, the net repulsion is larger than the sum of all three pair repulsions. The nonadditive effect is repulsive.

Dispersion nonadditivity The mechanism of dispersion nonadditivity was proposed over 50 years ago by Axilrod and Teller [76] and independently by Muto [77]. It is referred to as the correlation of three instantaneous dipoles. To better appreciate the behavior of this term, let us consider the same two extreme configurations of a trimer as for the TE nonadditivity described earlier. In the equilateral triangle the three monomers cooperate in correlating with each other; i.e. when a third monomer gets close, it sees the other two conveniently "pre-correlated." In contrast, for the collinear approach of a third monomer this "pre-correlation" takes place in the wrong direction. Since pair dispersion interaction is attractive, the nonadditivity is repulsive for the equilateral trimer and attractive in the collinear form. The simplest model for the three-body dispersion energy is provided by the well-known triple-dipole expression VDDD = 3 ZDDD 1 + 3 cos01 cos02 cos03 R312R33 R31

(19)

where Rij describe the sides and 0 i the angles of a triangle formed by three atoms, while ZDDD

691 describes a triple-dipole dispersion coefficient. A more accurate representation should account for higher terms (like DDQ etc.) as well as for the overlap effects.

Ar2HC1 The structure of this complex is shown in Fig. 9. Ab initio calculations of the three-body potential determined the relative importance of each term. They showed [78,79] that the total three-body interaction is very anisotropic with respect to the in-plane and out-of-plane rotations of HC1 within the cluster (see Fig. 10). It is instructive to have a closer look at the composition of this three-body term. The dispersion component is only moderately anisotropic, and the induction component is slightly more anisotropic than dispersion. The former may be reliably approximated by the third-order DDD term. The latter is dominated by the interaction of multipoles induced on the Ar atoms by the HC1 dipole. The most anisotropic is the exchange nonadditivity. It may be split into two physically different components: the "pure" exchange effect, TE, and the electrostatic interaction of the Ar 2 "exchange quadrupole" with the HC1 dipole, SE. It is important to notice that both parts of the total exchange nonadditivity show the opposite behavior. As seen in Fig. 10 the SE term is strongly repulsive for the O=0 ~ geometry, slightly attractive for O=0 ~ angle, and again becomes slightly attractive for the O=180 ~ The TE contribution is attractive for O=0 ~ and cancels a large portion of the SE contribution. This different behavior is not surprising as SE follows the pattern of the dipole-quadrupole electrostatic interaction while TE may be related to the mutual orthogonalization of monomers discussed earlier. Overall, approximating the total exchange nonadditivity by only its SE contribution, leads to too anisotropic a model.

Ar2HF Ab initio calculations of the three-body potential of Ar2HF by Cybulski et al. [79] reveal considerable differences between this system and Ar2HC1 [78] (see Table 2). First of all, because of a much closer approach of HF to the center of mass of Ar 2, the HL-exchange nonadditivity is larger and much more anisotropic in Ar2HF. The difference between the maximal and minimal value of the interaction energy at a given distance (amplitude) amounts to about 42 ~tEh for Ar2HF vs. 22 ~tEh for Ar2HC1. Another difference is the induction nonadditivity; it is larger (since the dipole of HF is larger than that of HC1), and also much more anisotropic (the amplitude of 38 ].tEh for Ar2HF vs. 15 ~tEh for Ar2HC1). The dispersion nonadditivity is smaller in magnitude and slightly less anisotropic, mostly due to the smaller polarizability of HF. The decomposition of the HL-exchange effect into its SE and TE components is shown in Fig. 10. It is qualitatively similar to the Ar2HC1 case. A comparison of higher order effects also proves instructive: Ar2HF has a significant contribution from the induction-dispersion coupling, whereas Ar2HC1, from the exchangedispersion effect. In conclusion, the three-body effect in these systems represents a somewhat different blend, with the induction-type components being much more important for Ar2HF. The first analytical three-body potential for Ar2HC1 was proposed by Hutson and collaborators [75] on the basis of semiempirical considerations. It consisted of three terms: exchange, induction, and dispersion nonadditivities. The dispersion nonadditivity was represented by Eq.(19), and the induction nonadditivity as an interaction of multipoles induced

692

'

601

A)

E, ~tEh/

'

'

' ' ' A r2HCl

'

'

T'-- '1 ota

'

'

40

'

'

'

'

'

exchSE I

?.~,'",z~'t

20/,

,

b~

o~_~--+_._+

-~ ~ _

,-

~

e

disp

-'~.~ " , ~ ~ -~ _+ ~ - - - ~ ~ ~" ---.--" ~...

I

I --%-.._ .--~-. I -~ _~;-" - \

-20 ,[-

def

,

. "~

"--.- -- "~ ?~"

b--:.

-

exch

%..__..~

i-

.~""~I~"

.~:31

'>'...t, ~ ,i rT',Cl~ ,

,

- 180

70 BE h

,

,

,

- 120

'

'

B)

,

,

e x c , - , ~ ,

-60

I

'

'

0

I

'

'

,

'

~

50

'

,

,

,

120

I

'

'

I

180

'

'

exchSE

~/'~

+'~~..~ dee

30

-

,

Total

~

r

,

60

O, deg

1

A r2HF

,

exch

L'

n.~_ /

10

~- "+ * -~'~'.~

/7 f'" ~-tt

-10

--=-*-- +----,-- ---'U-. o ,

- 180

,

I

- 120

,

,

d isp

'.~'x,x~ .,.

_

t

.

" "*I

-60

,

exchTE ~

I

0

,

,r

-.-- *--' --- -,-----*--

-*" ,

I

60

,

+

1

120

,

,

1

I

180

O, deg Figure 10. The dependence of the three-body components upon the in-plane rotation of: A) HC1 in the Ar2HCI cluster, B) HF in the Ar2HF cluster. The following abreviations have been used: "exch" - eexch,rL."def" _ z-X~defA'~SCF;"disp" - "disp,~(30)""Total" - AEMP3; SE and TE denote single-exchange and triple exchange respectively.

693

by HC1 on two Ar atoms. The SE exchange nonadditivity was modeled by the interaction between the exchange-induced quadrupole moment on Ar 2 and the permanent moments of the HC1 molecule. Later, Ernesti and Hutson [80] introduced a number of modifications of the Ar2HF model by recognizing other mechanisms. For example, the interaction in the Ar dimer subsystem involves the dispersion effect, which induces a quadrupole moment of opposite sign to the exchange-quadrupole one. In addition, the induction effects in the trimer, due to the presence of HF rather than HC1, required a more accurate treatment. The nonadditive potential obtained in Ref. [80] was successfully employed to predict the change in the red shift of the HF stretch by Bali6 and collaborators [81].

Ar2C02 In this cluster the three-body effect was detected via the observation of the asymmetric stretching frequency of CO 2 by Sperhac et al. [82]. Recent ab initio calculations confirmed the experimental predictions [83]. The exciting aspect of this cluster is that the nonadditive effect on the stretching frequency may be obtained directly with a very good accuracy. The reason is the well defined structure of the Ar2CO 2 cluster, shown in Fig. 11

O Ar ~

,-.~

()

., ,...

" "

II

""t~F

Figure 11. The T-shaped configuration of the Ar2CO 2 cluster. As long as the two Ar atoms are held in equivalent equatorial positions, the interaction with each of them should, in the pairwise additive approximation, result in the same incremental shifts of the asymmetric stretch of CO 2. In reality, a minute nonadditivity of shifts amounting to 0.042 cm-1 was observed by Sperhac et al. when the second Ar atom was added. Rak et al. [83] used a one dimensional model to calculate the ab initio estimate of three basic nonadditive components in the v=0 and v=l levels of this stretch (see Table 3). The three-body dispersion interaction affects the v=0 and v=l levels more than the other two terms. The effect of induction nonadditivity is about an order of magnitude smaller. Even smaller is the exchange nonadditivity. Interestingly, its SE and TE constituents are large, but quite accurately cancel one another. These results are sensible. Ar2CO 2 in the T-shaped configuration is predominantly dispersion bound and the induction effect should play a secondary role only. The signs of TE and SE terms are also easy to predict by means of the models outlined at the beginning of this Section: exchange-quadrupole interaction and distortion due to orthogonalizing, respectively. A different picture emerges when we analyse the effects of the three-body terms upon the

694 Table 2 Comparison of nonadditive terms in Ar2HC1 and Ar2HF. All values are in laEh. Three-body term

Ar2HC1

HL

8exch

Ar2HF

4.8

HL 8exch,SE

24.9

27.2 a)

44 .7 b)

-23.5 a)

- 19.8 b)

(30) ind,r

I 1.5 16.5

30.8 36 1

AEscF

16.3

55.7

AE(2)

8.3

-10.0

~(3o) disp

30.6

21.5

HL 8exch,TE AESCF def

AE(3) 26.3 AEMP3 50.9 a)Approximate calculation from Ref. [79]. b)Accurate calculation from Ref. [72].

20.2 65.9

Table 3 The effect of the three-body contributions upon the frequency shift of the CO 2 antisymmetric stretch (in ~tEh).

Three-body term

v=o

v=l

v=l-v=O

8exc hIlL HL SE Eexch,

0.774 10.849

0.629 10.864

-0.145 0.015

HL E Eexch,T

10.077

- 10.234

-0.158

(30) ind,r

4.376

4.815

0.438

AF~def SCF AFt(2)

1.542 2.321

1.970 2.254

0.428 -0.067

8(30) disp

27.170

27.069

-0.101

Sum:

HE +1~(30) + ..(30) ind,r e'disp Sum: EeI'~xch+AEdef SCF +AE(2)+ ~disp ..(30) s

Experiment [82]

0.192 0.125 0.192

695 frequency shift, i.e. the difference of v=l and v=0. Although the dispersion term is the largest, it differs very little for the the v=0 and v=l states states, and so do the SE and TE terms. The three-body induction term strongly differentiates between the two states, and thus has the most effect upon the shift. The dramatic change in induction effect can be rationalized in terms of the appearance of the dipole moment when CO 2 is deformed along the asymmetric stretching coordinate.

5.2 Water trimer: induction nonadditivity In the water trimer induction nonadditivity provides a dominant contribution, which effectively overshadows all the other terms. Its mechanism is simple. For instance, in a cyclic water trimer the multipoles of A inductively alter the multipoles at B, which, in turn, inductively alter the multipoles at C, which then alter those on A, and so on, until the selfconsistency is reached. Various formulations of this simple model were implemented in the simulations since the 1970s [84-87,63,64,50]. To include the many-body induction effects of point charges interacting with a set of polarizable atomic centers the following classical electrostatics equation is solved iteratively E

1 Fo poI'--2E ~l'i i i

(20)

where [.I,i is the induced dipole on center i and F ~ 1 is the electric field at center i arising from all other fixed charges in the cluster. The induced dipole moment ILl,i and the total electric field at the polarizable center j are evaluated self-consistently from the following expressions:

~l,i=o~iF i

(21)

and

o Fi=Fi + E TijlLtj j~-i

(22)

where ~i is the polarizability of the center (e.g. atom) i and Tij is the dipole-dipole interaction tensor closely related to the interaction matrix desribed earlier, cf. Eq.(10). Eqs.(20-22) SCF neglect the exchange effects arising in tA. . x127 X - , d e f . Yet, they fairly well approximate the latter over a wide range of configurations [45]. This approach has been successfully incorporated into the nonadditive molecular dynamics simulations [88]. The neglect of other nonadditive effects exchange and dispersion - appears to be justified for low energy configurations. They become important, however, for repulsive geometries, e.g. at barriers separating superimposable structures [45]. An ab initio model of nonadditive effects in water, which also includes SE and TE components of the exchange repulsion component, the dispersion term and perhaps some other, secondary terms should be investigated in the near future. So far the effect of nonadditivity in water has been studied in the context of various structural properties, vibrationally averaged structures, O-H frequency shifts [89], zero-point energies, rotational constants, cluster predissociation dynamics and tunneling splittings [58].

696 The importance of three-body effects in the determination of macroscopic properties has also been studied. Recently, nonadditive molecular dynamics simulation of water and organic liquids has been performed [86].

6. SUMMARY The nature of Van der Waals binding may be described in terms of four basic types of interactions: electrostatic, induction, dispersion and exchange. These interactions are useful to classify and understand the physical origin of intermolecular potentials, and the probable structures of Van der Waals complexes. In this context, they play a similar role to the concepts of covalent and ionic binding in strong chemical interactions. The fundamental interaction energy constituents are included in the rigorous quantum theory of weak molecular clusters. They can be calculated with a desired accuracy, and represented with analytical forms that are related to a variety of simple and physically sensible models. Today, both the ab initio calculations and the potential modeling can be performed for small- and medium-sized molecules, and provide reliable intermolecular PESs for a wide range of mutual intermolecular orientations. Future calculations should also include the intramonomer degrees of freedom, and incorporate them into the final potential form. Very little has been done in this area so far but the tools are already available. All these advances will be accompanied by simulations of the dynamics and vibrational averaging, with the ultimate goal of bridging the chasm between our understanding of the electronic structure of atoms and molecules, and macroscopic character of matter.

8. A C K N O W L E D M E N T S We thank Dr. Piotr Cieplak for reading and commenting on the manuscript, and Professor Richard Bader and Dr. Todd Keith for providing us codes, which draw the Laplacian of the electron density. Gaussian 92 codes [90] were used for electronic structure calculations. Support by KBN through the Department of Chemistry, University of Warsaw, within the Project BST/532/23/97 and by the National Science Foundation (Grant no. CHE-9527099) are gratefully acknowledged. The Interdisciplinary Center of Modeling, University of Warsaw is acknowledged for the computational grant.

REFERENCES 1. L. Pauling, The Nature of The Chemical Bond and the Structure of Molecules and Crystals; An Introduction to Modern Structural Chemistry, 3rd ed., Cornell University Press, Ithaca, N.Y., 1960. 2. L. Pauling, General Chemistry, 3rd ed., W.H. Freeman, San Francisco, 1970. 3. L. Pauling in: The Chemical Bond. Structure and Dynamics, A. Zawail (ed.) Academic Press, San Diego, 1992. 4. J.O. Hirschfelder, C.F. Curtiss, and R.B. Bird, Molecular Theory of Gases and Liquids, Wiley, New York, 1954.

697

5. F. London, Trans. Faraday Soc. 33 (1937) 8. 6. A.D. Buckingham in: Intermolecular Interactions: From Diatomics to Biopolymers, B. Pullman (ed.), Wiley, New York, 1978. 7. P. Piecuch, in: Molecules in Physics, Chemistry and Biology, P. Maruani (ed.), Kluwer, Dordrecht, 1988, vol. 2, p. 417. 8. H. Margenau and N.R. Kestener, Theory of Intermolecular Forces, Pergamon, Oxford, 1971. 9. J.O Hirschfelder and W.J. Meath, Adv. Chem. Phys. 12 (1976) 3. 10. J.N. Murrell, in: Rare-Gas Solids, M.L. Klein and J.A. Venables (eds.), Academic Press, London, 1976, p. 177. 11. B.Jeziorski and W. Kotos, in: Molecular Interactions, H. Ratajczak and W.J. OrvilleThomas (eds.), Wiley,New York, 1982, vol. 3, p. 1. 12. A.D. Buckingham and P.W. Fowler, J. Chem. Phys. 79 (1983) 6426; Can. J. Chem. 63 (1985) 1985. 13. C.E. Dykstra, Chem. Rev. 93 (1993) 2339. 14. A.D. Buckingham, P.W. Fowler, and J.M. Hutson, Chem. Rev. 88 (1988) 963. 15. A. van der Avoird, P.E.S. Wormer, and R. Moszyfiski, Chem. Rev., 94 (1994) 1931. 16. C. Bissonnette, K.G. Crowell, R.J. Le Roy, R.J. Wheatley, and W.J. Meath, J. Chem. Phys. 105 (1996) 2639. 17. D.M. Chipman and J.O Hirschfelder, J. Chem. Phys. 59 (1973) 2838. 18. B. Jeziorski, R. Moszyriski, and K. Szalewicz, Chem. Rev., 94 (1994) 1887. 19. G.Chalasifiski and M.M. Szcz~niak, Chem. Rev., 94 (1994) 1723. 20. I.C. Hayes and A.J. Stone Mol. Phys. 53 (1984) 83; ibid. 53 (1984) 69. 21. S.F. Boys and F. Bernardi, Mol. Phys. 19 (1970) 553. 22. G. Chatasifiski and M. Gutowski, Chem. Rev. 88 (1988) 943. 23. F.B. van Duijneveldt, J.G.C.M. van Duijneveldt-van deRijdt, and J.H. van Lenthe Chem. Rev. 94 (1994) 1873. 24. M. Gutowski, G.Chatasifiski and M.M. Szcz~niak, Chem. Phys. Lett. 241 (1995) 140. 25. R.J. Bartlett and J.F. Stanton in: Reviews of Computational Chemistry, K.B. Lipkowitz and D.B. Boyd (eds), VCH Publishers, New York, 1994, Vol. 5 p. 65. 26. J.O. Hirschfelder, Chem. Phys. Lett. 1 (1967) 325. 27. K. Szalewicz and B. Jeziorski, Mol. Phys. 38 (1979) 191. 28. G.Chatasiriski and M.M. Szcz~niak, Mol. Phys. 63 (1988) 205. 29. A.J. Stone and C.-S. Tong, J. Comput. Chem. 15 (1994) 1377. 30. C. Amoviolli and R. McWeeny, J. Mol. Structure (Theochem), 227 (1991) 1. 31. B. Kukawska-Tarnawska, G. Chatasifiski, and M.M. Szcz~niak, J. Mol. Structure (Theochem), 297 (1993) 313. 32. B. Kukawska-Tarnawska, G. Chalasifiski, and K. Olszewski, J. Chem. Phys., 10t (1994) 4964. 33. R.F. Bader, Atoms in Molecules, Clarendon Press, Oxford, 1994. 34. J.C. Slater and J.G. Kirkwood, Phys. Rev. 37 (1931) 682. 35. L. Pauling and J.Y. Beach, Phys. Rev. 47 (1935) 686. 36. R. Burcl, G. Chatasifiski, R. Bukowski, and M.M. Szcz~niak, J. Chem. Phys. 103, (1995) 1498. 37. M. Gutowski, J. Verbeek, J.H. van Lenthe, and G. Chatasifiski, Chem. Phys. 111 (1987) 271.

698 38 (a) F.-M. Tao, and Y.-K. Pan, J. Chem. Phys. 97 (1992) 4989. (b) F.-M. Tao, and Y.K. Pan, Chem. Phys. Lett. 194 (1992) 162. 39. M. Gutowski and L. Piela, Mol. Phys. 64 (1988) 337. 40. K. Morokuma and K. Kitaura, in: Molecular Interactions, H. Ratajczak and WJ. OrvilleThomas (eds.), Wiley, New York, 1982, vol. 1, p. 21. 41. L.A. Curtiss, A.J. Pochatko, A.E. Reed, and F. Weinhold, J. Chem. Phys. 82 (1985) 6833; E.D. Gladening and A. Streitwieser, J. Chem. Phys. 100 (1994) 2900. 42. A.J Stone, Chem. Phys. Lett., 211 (1993) 101. 43. P.J. Marshall, M.M. Szcz~niak, J. Sadlej, G. Chatasifiski, M.A. ter Horst, C.J. Jameson, J. Chem. Phys. 104 (1996) 6569. 44. R.F. Bader and T.A. Keith, J. Chem. Phys. 99 (1993) 3685. 45. G. Chatasifiski, M.M. Szcz~niak, P. Cieplak, and S. Scheiner, J. Chem. Phys., 94 (1991) 2873. 46. K.T. Tang, J.P. Toennies, J. Chem. Phys., 95 (1992) 5918. 47. M.A. ter Horst and CA. Jameson, J. Chem. Phys. 105 (1996) 6787. 48. E.J. Bohac, M.D. Marshall, and R.E. Miller, J. Chem. Phys. 97 (1992) 4890. 49. O. Matsuoka, E. Clementi, and M. Yoshimine, J. Chem. Phys. 64 (1976) 1351. 50. U. Niesar, G. Corongiu, E. Clementi, G.R. Keller, and D.K. Bhattacharya, J. Phys. Chem. 94 (1990) 7949. 51. C. Millot, AJ. Stone, Mol. Phys. 77 (1992) 439. 52. A.J. Stone, The theory of intermolecular interactions, Clarendon Press, Oxford 1996. 53. M.M. Szcz~gniak, Rd. Brenstein, S.M. Cybulski, and S. Scheiner, J. Phys. Chem. 94 (1990) 1781. 54. A.J. Stone in: Theoretical Models of Chemical Bonding, Z.B. Maksi6 (ed.), Springer Verlag, Berlin, 1991, Vol. 4, p. 103. 55. W. Rijks and P.E.S. Wormer, J. Chem. Phys. 90 (1989) 6507; ibid 92 (1990) 5754. 56. S.C. Althorpe, D.C. Clary, J.Chem.Phys., 101 (1994) 3603. 57. J.K. Gregory, D.C. Clary, J.Chem.Phys. 101 (1995) 7817. 58. J.K. Gregory, D.C. Clary, J.Chem.Phys. 103 (1995) 8924. 59. J.K. Gregory, D.C. Clary, J.Chem.Phys. 105 (1996) 6626. 60. W.J. Meath and M. Koulis, J. Mol. Structure (Theochem) 226 1 (1991). 61. G. Dotelli and L. Jansen, Physica, A234 (1996) 151. 62. V. Lotrich and K. Szalewicz, Chem. Phys. 106 (1997) 9688. 63. P. Cieplak, P.A. KoUman, and T.P. Lybrand, T.P., J. Chem. Phys., 92 (1990) 6755. 64. P. Cieplak, and P. Kollman, J. Chem. Phys. 92 (1990) 6761; P. Cieplak, T.P. Lybrand, and P. Kollman, J. Chem. Phys. 86 (1987) 6393. 65. L. Perera and M.L. Berkowitz, M.L.J. Chem. Phys. 100 (1994) 3085. 66. W. Kotos, F. Nieves, and O. Novaro, Chem. Phys. Lett. 41 (1976) 431. 67. J. Higgins, C. Callegari, J. Reho, F. Stienkemeier, W.E. Ernst, K.K. Lehmann, M. Gutowski, and G. Scoles, Science 273 (1996) 629. 68. M Wilson, and P.A. Madden, J. Phys. Condensed Matter, 6 (1994) 159. 69. R. Moszyriski, P.E.S. Wormer, B. Jeziorski, A. van der Avoird, A. J. Chem. Phys., 103 (1995) 8058. 70. V. Lotrich and K. Szalewicz, J. Chem. Phys. 106 (1997) 9668. 71. M.M. Szczg~niak and G. Chatasifiski, in: Molecular Interactions, S. Scheiner (ed.), Wiley, Chichester, 1997, p.45.

699 72. G. Chatasiriski, J. Rak, M.M. Szcz~niak, and S.M. Cybulski, J. Chem. Phys., 106 (1997) 3301. 73. M.J. Elrod and R.J. Saykally, Chem.Rev., 94 (1994) 1975. 74. L. Jansen, Adv. Quantum Chem. 2 (1965) 119. 75. A.R. Cooper and J.M. Hutson J. Chem. Phys., 98 (1993) 5337. 76. B.M. Axilrod, E. Teller, J. Chem. Phys. 11 (1943) 299. 77. Y. Muto, Proc. Phys. Math. Soc. Jpn. 17 (1943) 629. 78. M.M. Szcz~niak, G. Chatasiliski, and P. Piecuch, J. Chem. Phys., 99 (1993) 6732. 79. S.M. Cybulski, M.M. Szcz~niak, and G. Chatasiriski, J. Chem. Phys., 101 (1994) 10708. 80. A. Emesti and J.M. Hutson Phys. Rev. A, 51 (1995) 239. 81. P. Niyaz, Z. Ba~,i6, J.W. Moskowitz, and K.E. Schmidt, Chem. Phys. Lett. 252 (1996) 23; Z. Ba~,i~, American Conference on Theoretical Chemistry, 1996. 82. J.M. Sperhac, M.J. Weida, D.J. Nesbitt, J. Chem. Phys. 104 (1996) 2202. 83. J. Rak, M.M. Szcz~niak, G, Chatasiriski, and S.M. Cybulski, J. Chem. Phys. 106 (1997) 3301. 84. F.H. Sillinger, C.W. David, J. Chem. Phys. 69 (1978) 1473. 85. P. Barnes, J.L. Finney, J.D. Nicholas, and, J.E Quinn, Nature 282 (1979) 459. 86. J.A. Rullman and P.T. van Duijnen, Mol. Phys. 63 (1988) 451. 87. M. Sprik and M.L. Klein, J. Chem. Phys. 89 (1988) 7556. 88. J.W. Caldwell, L.X. Dang, and P.A. Kollman, J. Am. Chem. Soc., 112 (1991) 9144; J.W. Caldwell and P.A. Kollman, J. Phys. Chem. 99 (1995) 6208. 89. J.G.C.M. van Duijneveldt-van de Rijdt and F.B. van Duijneveldt, Chem. Phys. 175 (1993) 271. 90. M.J. Frisch, G.W. Trucks, M. Head-Gordon, P.M.W. Gill, M.W. Wong, J.B. Foresman, B.G. Johnson, H.B. Schlegel, M.A. Robb, E.S. Replogle, R. Gomperts, J.L. Andres, K. Raghavachari, J.S. Binkley, C. Gonzalez, R.L. Martin, D.J. Fox, D.J. Defrees, J. Baker, J.J.P. Stewart, and J.A. Pople, Gaussian 92, Gaussian, Inc., Pittsburgh PA, 1992.