Theoretical chemistry in solution. Some results and perspectives of the continuum methods and in particular of the polarizable continuum model

Journal of Molecular Structure (Theochem), 234 (1991) 401-424 Elsevier Science Publishers B.V., Amsterdam

401

Theoretical chemistry in solution. Some results and perspectives of the continuum methods and in particular of the polarizable continuum model Jacopo Tomasi D~partwnento dl Chrmxa e Chrmza Industrtale, V&aRcsorgLmento 35, I-56126 Pua (Italy)

Rosanna Bonaccorsi Zstztutodr Chlmxa Quuntuttca ed Energetrca Molecolare CNR, Via Rlsorgrmento 35, I-56126 PLsa (Italy)

Roberto Cammi Zstztutodl Chtmzca Ftsrca, Vtale delk Scienze 1, Z-43100 Parma (Italy)

Francisco J. Olivares de1 Valle Departamento de Quimtca Fislca, Unrversrdad de Extremadura, 06071 Badaloz (Spa&n) (Received 5 October 1990)

Abstract A concise overview of the theoretical methods based on a continuous distribution of solvent molecules (the contmuum dielectric model is the most important example) is accompanied by the examination of a few sub]ects, selected as examples of the recent developments of the method. The subjects are: (1) electron correlation and dispersion; (2) optimization of the solute geometry; and (3 ) boundary surfaces in solutions.

1 INTRODUCTION

The expression “theoretical chemistry in solution” covers such a large variety of scientific interests, approaches and methods that it is hard to find a unifying perspective able to organize the basic trends of the present research in this field into a coherent scheme. Chemistry in solution has been, and still is, the main object of chemistry: more than a century of vigorously pursued investigations constitutes the background for theoretical investigations. Moreover, the energy range of chemical phenomena makes unavoidable (or, better, necessary) strong overlaps with research domains traditionally covered by physics, biology and engineering. Complex problems require the use of complex methodologies, and only an ap-

0166-1280/91/$03

50 0 1991 Elsevler Science Publishers B V All rights reserved

402

propriate combination of different approaches represents, in many cases, the strategy to follow in order to ensure substantial progress. This combination of different strategies is gaining importance in this area: until now the effort has been especially directed to improving and extending specific approaches. The present paper also considers only one specific method: the occasion is not appropriate to show in detail how this method may be combined with others to achieve a fuller understanding of the phenomena of a chemical or a physico-chemical nature occurring in solutions. We have considered it necessary to stress here the need of a multi-approach strategy in order to avoid misunderstandings about the context in which the present report must be placed. 2. AN ATTEMPT OF CLASSIFICATION

A general classification of methods is very difficult, as we have already said. We put forward, however, a tentative classification based on three main domains (a) methods based on the elaboration of physical functions, (b) methods based on a molecular description of a portion of the whole liquid system and (c ) methods based on a molecular description of a restricted portion of the system accompanied by a less detailed description of the whole system. The classification is quite rough and neglects the mixing among the three domains. In domain (a) we include the perturbation theories which start from a simple reference system to obtain, via a suitable expansion, the physical (essentially thermodynamical) properties of interest. In the same domain we also include all the elaborations of the distribution function theories. This very large field is growing at a considerable rate, giving us an increasing amount of information about the physical aspects of the problem [l-3]. In domain (b) we collect supermolecule calculations and computer simulations. The philosophy of the two approaches is quite different, and the information gained with the two approaches is, to a large extent, complementary with a small overlap. In both cases, a detailed description of the disposition and structure of molecules composing the liquid system is searched. The supermolecule approach is more addressed to the chemical aspects of the problem: solvation sites, solvent effects on the solute properties, chemical reactions assisted by a limited number of solvent molecules, etc. [ 41. The full ab initio formulation of this approach gives us the most accurate (though static) description of phenomena occurring at the submolecular level; it is, however, limited to treating small solvation clusters, too small for a full appreciation of solvent effects (the same remarks apply to the related approach of a description of molecular interactions via the perturbation theory). Some positive aspects of the ab initio supermolecule approach are present in

403

more approximate methods which are able to consider a larger number of solvent molecules at the expense of the accuracy of the description. There is a hierarchy of methods with a progressive decrease in the sophistication of the description, ending with the molecular mechanics methods which now represent one of the most promising methods of treating the chemistry of biological molecules in solution. The computer simulation methods (Monte Carlo, Molecular Dynamics, Brownian Dynamics, etc.) are probably the most powerful tools available today [ 5-71. Limited in preceding years to the description of the structure of the liquid and to some information about the mean energy of the molecules in the condensed systems, they are able at present to treat real chemical problems, such as reactions, and to give information about a larger number of physicochemical quantities. Domain (c) collects methods in which there is a target subsystem, the “solute”, for which a full Hamiltonian, including an averaged solute-solvent interaction potential, may be written. There are several approaches to be classified into this domain. Some among them belong to the category of semiempirical methods and some play a role in the study of complex solutes (e.g. molecules of biological interest) [8-121. Others retain the formalism of the in vacua ab initio molecular calculations, including in the Hamiltonian, and in the corresponding Schriidinger equation, an explicit expression of the solute-solvent interaction potential VmtiThis is the approach we shall consider in more detail in the present paper. 3 THE CONTINUUM MODELS

The interaction potential, in the last type of approach, is related to a description of the solvent as a continuous distribution. In other words, the basic assumption is that the solvent molecules, S, or (better) the atoms, t, composing the solvent molecules, may be described by a continuous distribution function g,(FJ which may be averaged over time or may display a time dependence; these functions may reflect the effect of the solute on the solvent distribution (cybotactic regions) or may discard these effects. In general, g,(f;) corresponds to a thermal average of some of the motions of the solvent molecules. The origin and extent of the last average may be stated more precisely in the elaboration of the model. Once the solvent distribution function is given, one has to specify the interactions considered by the interaction potential V,,. The supermolecule approach (both in the variational and in perturbation theory formalism) give some hints for the definition of the interaction terms. In passing from simple two-body interactions to real liquid systems involving numerous molecules and a variety of boundary conditions (e.g. the occurrence of phase separations), the relative importance of the various interaction terms may change, and other

404

modifications occur when the chemical and physical definition of the condensed system is changed. Part of this shift in the emphasis is due to the occurrence of many-body interactions. When the interaction energy in the system is decomposed according to the number of bodies (i.e. molecules), each term presents its own problems. Let us consider the simple case of a very dilute solution with a solute M encircled by a large number of solvent molecules (S,) @NM,

S1, SP, ... S,) > = G&M,

S,) > + (=AE(S,,

+ (-=JE(M,

S,) >

S,, S,) > + <-=zdE(S,,

S,, S,) >

+ (CZZAE(M, S,, S,, Sk)) +...

(1)

where the sums run over the solvent molecules. The first term on the righthand side of eqn. (1) may be decomposed into a sum of two-body terms (zdE(M,

S,) > +C
Sz) >

(2)

for which the supermolecule approach gives us the tools for the analysis, e.g. via decompositions such as dE(M, S,; R) =&s(R)

+&L (R) +&x

(R) +&T (R) +&IS (R)

(3)

The properties of M in solution differ from those of M in vacua, and this difference influences the various components of dE(M, S,; R) according to different, but well understood, rules. Also, the second two-body term of expansion (1) may be analyzed in terms of instantaneous, not averaged, solvent-solvent interactions. The presence of M affects the average value valid for the pure solvent. The new versions of the continuum approach may take into account these changes when necessary. Three-body interactions introduce non-additive effects, mainly related to polarization and dispersion. Here again the simple clusters M-S,-S, and S-S,Sk represent the starting point for the elaboration of the model, but the final description should take into account both effects specific to the pure solvent and effects due to the presence of the solute. Non-additive three-body contributions are partially modified by four- and many-body interactions. These concise remarks should be sufficient to highlight the problem caused by the passage from models for two-body interactions to models for solutions. There are also interactions of another type, generally not considered in the study of two-body systems, which may be of some interest for solutions. Just to give an example, we quote the so-called external heavy atom (EHA) effect [ 131 of ubiquitous occurrence in solution which influences many spectral properties of the solute and which could be included in Vlnt via interaction operators not present in the usual electrostatic Hamiltonian (unpublished results).

405

The simplest realization of the continuum model considers the solvent as a polarizable isotropic dielectric, including in V,, only the mutual solvent-electron polarization effects. This very simple model (it should be called the primitive continuum model) is directly related to the old, and famous, continuum electrostatic theories connected with the names of Born, Kirkwood, Onsager. Actually, the new versions of the continuum model are the direct filiations of these models. Kirkwood [ 141 pointed out that the solute charge distribution used in these theories may be derived from quantum calculations, but at that time it was not possible to formulate practical methods to exploit that indication. We may use this point to make a distinction between two categories of continuum electrostatic models. The first category uses the charge distribution of the solute, even when obtained from quantum mechanical calculations, in a classical way. Either the charge distribution is not affected by the solvent, and there is only a classical electrostatic problem to solve, or it may be polarized by the solvent reaction field, and there is a more complex electrostatic problem to solve, with polarization effects described in terms of solute polarizabilities (often expanded over several centres) . The second category makes explicit use of the Schrijdinger equation (H& + V& Y’ = E’ !P’

(4)

and the problem is amenable to the contemporary solution of the electrostatic problem and of the quantum molecular problem [ 151. The distinction between the classical and non-classical approaches is not sufficient to characterize completely the continuum methods. Another point of some importance is how the interaction potential is represented and how the cavity is described. The historical tradition is to use a multipole expansion of the potential. In the general case, the expansion must be truncated and this fact may introduce a lowering in the quality of the results. Multipole expansions are currently accompanied by a description of the cavity in terms of simple shapes (a sphere or an ellipsoid): in such a way, the cavity surface corresponds to a fixed value of one of the coordinates of the reference frame, and the electrostatic problem becomes simpler but the model less realistic. The first exception has been the Claverie method (a classical method, with expansion of the interaction potential) which adopts realistic cavities, expressed in terms of intersecting spheres centred on the solute atoms [ 161. We present in Table 1 a short selection of methods related to the classical continuum electrostatic approach. The list is far from being complete, but is sufficient to show some trends in the evolution of the models. Recent methods adopt more realistic descriptions of the cavity; this is a consequence of the development of the computer facilities which permit the adoption of algorithms able to describe complex closed surfaces by numerical methods [27291. Advances in the description of molecular electrostatic potentials of quan-

406 TABLE 1 Classrcal versions of the contmuum electrostatic model for solutions Authors (year )

Shape of the cavity

Expression of the source potential

Ref.

Born (1920) Kirkwood (1934) Onsager (1936) Huron-Clavene (1972) Bevendge-Schnuelle (1976) Abraham-Lrszi (1979) Miertui-Tomasr (1982) Warwicker-Watson (1982) Shaw (1985) Rashin-Namboodiri (1987) Drummond (1988) Abraham et al. (1988) Zauhar-Morgan (1988)

Sphere Sphere Sphere General shape Concentric spheres Concentnc spheres General shape General shape General shape General shape General shape General shape General shape

Monopole Multipole-one-centre expansion Dipole Multipole-multicentre expansion Multipole-one-centre expansion Multipole-one-centre expansion Complete evaluation Distributed monopoles Distributed monopoles Distributed monopoles Dlstnbuted monopoles Distrrbuted monopoles Distributed monopoles

17 14 18 16 19 20 21 22 23 24 25 26 27

turn origin [30] via distributed monopoles or other many-centre expansions has given more confidence to the use of these descriptions of the charge distribution. The multipole expansion of the interaction potential is often replaced by other methods: the method of images [23,31], methods based on numerical integration over spatial grids [ 22,32,33], methods based on the use of apparent polarization charges defined on the cavity surface [l&21,24-27,34]. This last approach seems to us to be the most promising, at least for isotropic solvents. More information on the recent developments in classical continuum methods may be found in a review by Davis and McCammon [ 351. The numerical results obtained with these methods show that inadequacies and failures of the continuum model lamented in the past were due to poor implementations of the method rather than to intrinsic physical deficiencies. The change in the quality of the results is so sharp that we may consider the methodological progress as a change of generation in the models, the first generation corresponding to the venerable models elaborated in the 1930s with successive improvements, still useful in restricted areas, and the second generation to the models with more realistic descriptions of cavity, charge distribution, interaction potential, a generation initiated less than 10 years ago. It may also be considered as belonging to a second generation (in a different meaning) the models based on a quantum-mechanical approach, a short list of which is presented in Table 2. This list, as the preceding one, is far from being complete, especially for methods relying on semiempirical quantum descriptions. In addition, we have not considered here methods using a supermolecular approach, with minor

407 TABLE 2 Quantum me&an& Authors

versrons of the continuum

(year)

Rivail et al (1973-1983) Tapia et al (1975-1978) Hilton-McCreery et al. (1976) Tomasi et al. (1981-1983) Agren et al. (1987) Hoshi et al (1987) Karlstrom (1988)

electrostatic

model for solutions

Shape of the cavity

Expression potential

Sphere/ellipsoid Sphere Sphere General shape Sphere General shape Sphere

Multlcentre-multipole expansion Dipole-monopole expansion Multlcentre-multipole expansion Exact expression Multipole expansion Exact expression Mulkcentre-multipole expansion

Ref

of the source

36-39 40,41 42 15,43 44 45 46

TABLE 3 Aspects of the description contmuum model (PCM)

and of the phenomenology

of salvation

treated

with the polarizable

Physzcal znteractzons zncluded znto V,, Electrostatic Repulsion Dispersion Electron correlation Spin-orbit coupling Solvent dzstrzbutzon Homogeneous Locally modified by the solute Anisotropic at large scale With limiting boundaries Phenomena and propertzes Thermodynamical properties Vlbratrons Electron excitation Dynamical phenomena Motions in restricted volumes Transfer of solutes across boundaries Molecular conformations Non-covalent solute-solute mteractions Tautomeric equilibria Dlradical/zwitterion equilibria Ionic dissociation mtimate/solvent-separated Chemical reactions: electrophihc addition,

ion pairs nucleophllic substitutions

rearrangements,

etc

corrections from the continuum distribution, a family of methods which is giving quite interesting results. Here again the most efficient models combine cavities realistic shape, ac-

408

curate representations of the interaction potential and wavefunctions of good quality. The progresses in the field of quantum-mechanical calculations for isolated molecules may be easily transferred to the quantum-mechanical continuum model for solvation. It is a short-sighted view to consider quantummechanical ab initio models as limited to material systems of small size; on the contrary, advances in the field of application to large systems (e.g. biomolecules) will quite probably be prompted by ab initio investigations. In addition, the domain of physical phenomena occurring in solution which will be treated with this approach is increasing at a noticeable speed. We may talk here of a third generation of continuum models. We shall not examine the basic aspects of quantum-mechanical models of solution, including a continuum description of part of the solvent [ 471, nor the differences between the various methods, although an analysis of the evolution of the various approaches, of the mistakes, and of the corrections, would be quite instructive. We shall limit ourselves to presenting in Table 3 a selection of subjects treated with our version of the model (sometimes indicated by the acronym PCM). We have insufficient space to comment upon all the themes indicated in the table (a portion of which has been considered by other research teams or has also been treated with other versions of the solvation model) and we shall select a few subjects, preferring variety to a more systematic approach. 4. DISPERSION ENERGY AND ELECTRON CORRELATION

The Schrijdinger equation (4) is treated, in models of the second generation, at the Hartree-Fock level. This level of the quantum molecular theory does not permit electron correlation effects inside M and the physical phenomena due to the correlation between electrons of M and of the surrounding solute (dispersion effects) to be taken into consideration. Dispersion effects are often included, as in models of the first generation, as separate contributions evaluated via discrete models and making use of approximate formulas derived from perturbation theory for bimolecular systems. We have considered it convenient to explore two different approaches. The first approach starts from the elaboration of the problem given by Linder [ 481 who applied the reaction field concept to a fluctuating dipole. Generalization to the case of cavities of general shape, the insertion of this term in the Hamiltonian (4) and the reformulation of the iterative method [ 151 for the contemporary solution of the electrostatic and of the quantum molecular problem, including dispersion, has been performed at Badajoz [491. Elaboration of the method has been limited thus far to the dipolar component of the dispersion term, and to isotropic solvents, but its extension should not present prohibitive difficulties. (The dipolar dispersion term has also been treated by the Rivail group, again using the Linder formalism [ 501. )

An approach based on atom-atom dispersion contributions has been examined at Pisa [ 511. In this approach, the dispersion contributions are not inserted in the interaction operator of the Hamiltonian (4) and the formalism is thus less compact than the preceding one. It allows a greater flexibility and a lower computational cost. By a proper selection of the atom-atom potentials, one may gain a useful insight into the relative importance of the various terms of the dispersion energy expansion [ 521. A specimen of results obtained with different sets of potentials is reported in Table 4. This approach makes possible the use of more efficient computational strategies. There are, in fact, basis sets well balanced and able to provide a satisfactory solution for eqn. (4) in the electrostatic version, but unable to give a correspondingly good description of the oscillating dipole: electrostatic calculations with that basis set may be supplemented by an atom-atom dispersion contribution having the desired characteristics. This feature has been already exploited to obtain the appropriate definitions of the computational parameters to evaluate differences in the solvation energy in two solvents, and then partition coefficients, reaching the chemical accuracy [ 571. Even more important for the present discussion, the topic being related to the evolution of continuum models from the second to the third generation, is the application of the second approach for dispersion to solutions exhibiting anisotropies in the solvent distribution. The first kind of anisotropy we have examined is related to the self-organization of solvent molecules (water) around the solute. We have found excellent results for the evaluation of the dispersion energy using solute-solvent correlation functions derived from solvent simulations and from the application of the RISM version of the integral equation method [ 581. The second kind of anisotropy regards the occurrence of a separation surface in the system. We have explored, thus far, systems composed of two immiscible TABLE 4 Companson of GDIsvalues for selected solutes in water with different sets of dispersion coefficients (kcal mol-‘) Solute

W-L C,H,OH C,H,NH, “Ref 53 bRef. 54. “Ref 55. dRef 56

GDIS Set 1”

Set 2”

Set 3b

Set 4”

Set 5d

-872 -7 79 - 8.32

- 7.23 - 6.82 -738

-8.00 - 7.08 - 7.81

-9.20 -8.12 -8.71

-9.38 -8.61 -8.86

410

liquids separated by a flat surface. We shall discuss progress in this field later on, the topic not being limited to the dispersion contributions only. Our last remark on this subject regards the basic mathematical treatment of the information about the solvent distribution, which influences the interaction potential Vlnt.It is convenient to write the solvent distribution in terms of the spatial densities of the various types, Z,of atoms belonging to the solvent (5) The numerical density of the solvent, ps, is expressed as a function of the position, with respect to an arbitrary reference frame, to include the possibility that some portions of the space are filled by a different medium, or by the solvent, S, at a different density. The atom-atom correlation factors g,, ( Tml) refer to atoms, m, of the solute (but they may also belong to other portions of matter present in the system), placed in fixed position, which influence the spatial distribution of the solvent. The atom-atom correlation functions g,J TmI) , derived from computer simulations, are spherical averages of these functions. The mathematical algorithm we have adopted to get from eqn. (5) the expression of the interaction potential (or of the corresponding energy) is characterized by a simple integration over the cavity surface. This algorithm has also been used for the dispersion model developed in Badajoz and represents an important step in the development of continuum models, being of wide application, not limited to the dispersion effects. For more details see refs. 49, 51 and especially 52. Dispersion contributions to the energy are due to intramolecular correlation effects. Recently, we have developed a method to treat intermolecular correlation effects [59 1. The treatment of intramolecular electron correlation in solution is more complex than in vacua. In fact, the basic units of the independent electron model theory are modified in solution by the solvent reaction field. The reaction field is defined in terms of the solute charge distribution, which is influenced by electron correlation. To decouple these mutual effects, correlation and polarization, we have developed the theory at three separate levels. The energy for the system M in vacua may be expressed in terms of the expansion

(6) where ZD” is the correlated density matrix in vacua. The complete solution of the same problem for molecular solutes, denoted by the acronym (PTDE) , is

411

GsoL (PTDE) = C Gk(ZDn)

(7)

k=O

We have introduced free energy terms, instead of internal energy terms, since they are more appropriate to the case of a solution (see, for example, ref. 21). The correlated density matrix refers to the situation in solution, also taking into account the mutual effects between the solvent reaction potential and the fluctuation potential. In the formulation of the polarizable continuum model (PCM) at the Hartree-Fock level, the mutual interaction between the Fock potential and the solvent reaction potential has been taken into account by means of a double iterative procedure [ 15,431, giving as output an expression of Vlnt at the Hartree-Fock level. To get the PTDE expression we need a triple iterative procedure. An intermediate step is done by the consideration of the electron correlation effects with the solvent reaction potential expressed at the H-F level. The corresponding approximation is called the perturbation theory at the energy level (PTE) . G(PTE) = C G,(D”)

(3)

k

The functional dependence indicated in eqn. (8) demonstrates that the interaction potential is a function of the H-F first-order density matrix, modified by the solvent. The PTE approximation does not require the third iterative step in the calculations we have mentioned above. The third approximation is addressed to gain more information about the mutual interplay of solvent and electron correlation effects. This approximation, called the perturbation theory on the density distribution alone (PTD ), corresponds to a calculation of V,nt, taking into account correlation effects (and then three levels in the iterative procedure) but discarding the successive calculation of two- and more-body terms in the expansion of the independent electron model G(PTD) =G,(ZD”)

+G,(ZD”)

(9)

We have thus far examined the PTE approximation on several solutes, taking into account changes in the basis set and in the characteristics of the solvent [ 59,601 and exploring the potentialities of the approach for the description of the vibrational properties of a solute [ 611. We report an example of the results in Table 5. The numerical data are dissected into various components, according to the order (lower index) and to the number of bodies (upper index). The details may be found in ref. 59. A reduced specimen of the results obtained at the PTE and PTDE levels is reported in Table 6. A detailed comparison of the three approximations has not yet been performed.

412 TABLE 5 Components of the electrostatx free energy of solution at PTE theory level Numerical entry refers to some molecular solutes in water as a solvent and are relative to the 631G basis set. Free energy“

Hz0

H&O

Contrzbutzons wzth respect to the number of bodzes znvolved AGo - 19.373 - 24.289 AG’ 58.990 154.021 AG2 - 48.202 - 133.961 AG3 -0.181 -3.065 AC’ - 0.034 0.132

CHsOH

CH&&)

109.590 - 94.897 -1325 0 249

-21.205 44.390 - 30.352 -0.026 -0.012

- 22.296

Contrzbutzons wzth respect to thezr order AGO 39.618 AG, - 48.953 A& 0.503 & 0.034

129.732 - 140.040 3.766 - 0.620

87.295 -96.972 1.040 -0.041

23.186 -31.087 0.405 0.292

Total free energy wzth/wzthout correlation AG (PTE) - 8.798 AG (SCF)b - 9.335

- 7.162 - 10.308

- 8.679

- 7.204 -7.901

Total contrzbutzon of correlatzon to the free energy AAG,, 0.537 3.147 %” 5.8 30.5 APd 0.270 1.062

-9 678

0.999 10.3 0.518

“In kcal mol-’ 1 u=627.50959 kcal mol-‘. bAG (SCF) =AGo+AG1. ’ [ AAGJAG,, (SCF) ] x 100. dVariation of dipole moment: p (sol) - fl (vat ) , in Debyes.

TABLE 6 Salvation energy GELof Hz0 in water, computed at different levels of the theory Basis set 6-31G, values in kcal mol-‘. SCF PTE PTD PTDE PTDE + DISP

-

8.812 8.289 8.999 8.472 - 12.345

0.697

8.8 0.533

413 5 GEOMETRY OPTIMIZATION

To perform studies of chemical problems in solution we need an efficient algorithm to compute optimized geometries for isolated or interacting solutes. The complex iterative procedure necessary to obtain the energy, and then the potential energy surface (PES), even at the H-F level, prompted us first to examine computational strategies based on the numerical evaluation of the gradient. We have performed a systematic survey of different strategies, based on two methods for the search of a minimum on the PES and on three approximations for the definition of the cavity shape and of the interaction potential V,, in the intermediate steps of the optimization. The details may be found in the source paper [ 621. It will be sufficient to report here some conclusions. The selection of the most convenient strategy depends on the nature of the solute and also on the characteristics of the solvent. One of these strategies, which may also be used for intermediate steps in the cases in which more complex strategies are to be employed, can be set in a form for which analytical expressions of the gradient and of the Hessian may be elaborated. It is well known that the efficient analytical procedures for the search of TABLE 7 Linear H-bonded water dimer in solution 6-31G* equihbrium geometry (see Fig 1 for parameter definitions). Results relative to several solvents.

R

Vacua

Benzene

l,l-DCE

Water

2 971 32.520 117.500

2.944 35.200 121400

2.927 36.700 120.900

2.923 37.100 120.300

X

a, t’ * Y

Fig 1 Definition of the mtermolecular geometrical parameters of the water dimer considered m this work.

414

minima and saddle point have represented in the last years an important breakthrough for the study of chemical reactions in vacua. We feel that an analogous methodology will also be available shortly for chemical reactions in solution described via the polarizable continuum model. An example of the effect of solvent on the geometry of a “solute” (actually a hydrogen-bonded dimer) is given in Table 7. 6 MOLECULAR INTERACTIONS IN SOLUTION

The study of non-covalent molecular interactions is a prerequisite for the modelling of condensed systems and for the study of chemical reactions. Much attention has been paid in the past by our group to the interpretation of noncovalent interactions in vacua, and it is gratifying to note that, 20 years later, our hypotheses on the characteristics of non-covalent bonds of medium strength have been confirmed and accepted. Calculations (and interpretations) of non-covalent interaction energies are plagued by the so-called basis set superposition error (BSSE), present in all the calculations not performed with very large and especially elaborated basis sets. The present paper is not the appropriate place to comment upon the debate about the corrections to the BSSE. We have expressed our opinion on other occasions. Suffice it to say that the method proposed by some of us a few years ago to obtain an analysis of the interaction energy after the correction to the BSSE [63] has proved its validity in many applications, and that a new version of the total correction, based on the decomposition of refs. 63 and 64, has been proposed by one of us [ 651. This last proposal avoids several shortcomings of the previous definitions, and, according to the numerical checks thus far performed, gives results of noticeable interest. Also, calculations on interacting molecules in solution are affected by a BSSE. The experience gained in the study of the BSSE in vacua cannot be directly transferred to systems in solution. The occurrence of a solvent reaction field of considerable strength in the continuum models or the occurrence of other intermolecular interactions of similar nature in the discrete models introduces changes in the extent of the BSSE and on its distribution on the components of the interaction energy dEAB (R) which are not easily predictable. It is appropriate to recall here the accurate studies performed by Dykstra and coworkers on the effect of an external field on the BSSE and its correction [ 66,671. We have developed a method to introduce corrections for the BSSE to the energy of interaction acts A + B occurring in solution [ 681. As previously done for systems in vacua [63,64,69], the method starts from the decomposition of dEAB(R) elaborated by Kitaura and Morokuma [ 701 (or from a modified version elaborated by Bonaccorsi et al. [ 711) and applied to Hamiltonian (4)) including electrostatic interactions with the medium.

415

The energetic quantity of primary interest is again a free energy. We may thus write

AGAB

=GAB@)- [GA+GI

(10)

where GA and GBare the energies of the two separate systems, non-interacting and computed each with its appropriate cavity. The counterpoise procedure (CP) [ 721 for the correction of the BSSE requires calculations of the monomers performed under the same mathematical conditions set forth for the dimer. To do this we introduced an intermediate step in the description of the association process

(11)

= GD + Gm AGAB

where GPDis the partial desolvation (the calculation of the monomers A and B in the cavity suited for AB) and GINTis a quantity subjected to energy decomposition and CP corrections to its components. The quantity corrected with the counterpoise procedure are labelled with CP as superscript. As the final result, we have TABLE 8 Linear H-bonded water dimer in vacua. Partition of the interaction SCF energy (a) without and (b) with counterpoise corrections Numerical entries refers to the 6-31G* basrs set (a)

R(O-0)

2.40

E ES EPL

-

25.480 3.827 40.452 - 8.575 2.581 5.150

-

EEX

E CT E MIX AE

2.60

2.80

3.06

3.20

- 15.970 - 1.446 18.637 - 3.617 0.042 - 2.354

- 10.425 -0 731 8.458 - 2.138 -0.219 - 5.054

-7.186 - 0.439 3.795 - 1.569 -0 156 - 5.555

- 5.239 - 0.283 1.682 - 1.279 - 0.083 - 5.203

2.60

2.80

3.00

3.20

- 10.425 -0.731 8.642 - 1.864

- 7.186 - 0.439 3.934 - 1.337

- 5.239 - 0.283 1.779 - 1.007

- 3.987 0.391

- 0.392 4.636

-40 373

(b)

R(O-0)

2.40

E ES E

-

E%

40.744

Eg

-7.895

- 15.970 - 1.446 18.870 -3.191

6.824 3.283

- 0 1.062 676

A;S ECP

25.480

- 3.827

416 TABLE 9 Linear H-bonded water dimer in solutron. Partition of the interaction free energy (a) without and (b) with counterpoise correctrons Numerical entries refers to the 6.31G* basis set and to water as solvent. (a)

R(O-0)

2.40

2.60

2.80

3.00

3.20

GE,

-

- 16.860 - 1.021 17.924 -3 586 0.135 - 3.408 - 5.025 1.791

- 11.204 - 0.367 8.119 -2.194 -0.147 - 5.793 - 1.990 9.596

- 7.841 -0.127 3.630 - 1.647 -0.108 - 6.093 - 2 100 10.086

-5.776 -0 023 1.607 - 1.356 - 0.055 - 5.604 - 2 200 9.007

2.60

2.80

3.00

3.20

26 369 39.579 -3.262 -7 317

- 16.860 - 18.446 1.021 -2.918

-11.204 8.522 - 0.367 - 1.724

- 7.841 -0.127 - 3.931 1.258

- 5.776 - 0.023 - 0.959 1.820

2 976 5.606 1.479 12 692

0.545 -1808 1.791 -1825

0.298 - 4.476 1.990 -6.961

0.319 -4.977 2.100 - 7.853

0 313 -4 624 2 200 - 7.049

26.369 -3.262 38.925 - 8.282 2 560 3 571 8.622 1.479

GPL GEX GCT GMIX GINT d? G (b)

R(O-0)

2.40

GE,

-

G G% GE Gcr MIX GINT GPD

AG

=‘%I@) +Gxs(R) +&,(R) +GE:(R)+G:$(R) (12) The meaning of the subscripts may be found in the original paper on the interaction energy decomposition [ 701 or in the numerous papers on this subject. We have not yet applied the alternative scheme proposed in ref. 65 to interactions in solution; we recall that this correction scheme consists in deleting the correction to the GEx (R) and to two of the five components in which GMIx may be further partitioned. The results obtained thus far with eqn. (12) confirm, especially for hydrogen-bonded dimers, the conclusions drawn for in vacua systems [ 731. The relative importance of the classical (Gns + GPL) terms is increased, the results are less basis set dependent and, as a consequence, the validity of the assumptions made in passing from the many-body expression of the interactions in solution (eqns. (1) and (2 ) ) to the polarizable model is strengthened. We report, as an

417

example, a comparison between the water dimer interactions in vacua, Table 8, and in solution, Table 9, obtained with this method. An important description of the energetics and of the PES near the minimum for small clusters of molecules associated via hydrogen bonds may be of interest for many cases of chemical and physical relevance. In a recent paper on the vibrational properties of “free” H,O molecules in bulk water [ 611, we pointed out as problems accessible to the present level of elaboration of the polarizable continuum model the study of the components of the various structural models for water. Our feeling is also based on other evidence, which, for brevity, are not treated in this paper, but our experience on the CP corrections to the BSSE for dimers in solution gives support to our belief. 7 BOUNDARY SURFACES IN THE SOLUTION

In Sect. 4 we have briefly touched on a development of the continuum model which has given interesting results in preliminary attempts to set up the details of the method, and which we consider quite promising. This development regards the study of the physical and chemical behaviour of the molecular components when the solution presents a boundary surface. The variety of systems, and of problems, for which a knowledge of the physice-chemical properties of a solution near a boundary is required is immense. We leave to the reader to consider the number of effects in the field of chemistry, biology, physics and engineering based on the behaviour of a liquid near a separation surface. The basic formal aspects of the method are related to the transformation of volume integrals to surface integrals via the application of the Green theorem. This kind of transformation has also been exploited in the original derivation of the PCM [ 151 but applied on that occasion to a simple problem, the polarization of an isotropic continuum. Even for such a relatively simple problem, the use of this transformation has made it possible to give a satisfactory solution to the problem of using cavity boundaries of general shape. A modified version has been applied to study transformations in the DNA structure. The presence of numerous negative charges on the double helix, accompanied by the necessary amount of counterions, suggest the inclusion in the solvent model of a solvation layer at lower permittivity. In this case, the anisotropy of the solvent is reduced to the presence of two portions of liquids at different e. The model (in the classical version) has been applied to DNA specimens, reaching the maximum length of 1500 base pairs, in the normal geometry of the B form, twisted, bent and with partial opening of the double helix [ 74,751. The same model (in the quantum version) has been also used for the evaluation of the thermodynamical solvation functions (dG, LIH,AS) of charged solutes: some preliminary results are reported in ref. 76; a more complete report

418

regarding a unified model for a larger set of solvents is under elaboration. This last model seems to confirm the model elaborated by Abraham and Liszi a few years ago [ 771; an important difference is that the Abraham’s model is limited to spherical ions, which may be treated with solvent models described by concentric spheres, while our model may be applied to more general cases, for instance molecular cations having a large uncharged tail appended to the cationic head. This is one example more of the possibilities offered by the use of cavities of realistic shape. A different application of the models with two dielectrics has been attempted in ref. 78. In this case, the second dielectric represents a massive body, e.g. a biopolymer, and the paper is addressed to the study of the partial desolvation effect in the processes of molecular recognition. This application, related to processes of great importance in biological systems, will be further developed. A formal presentation of the PCM with the inclusion of a preselected partition of the space external to the cavity into many different non-overlapping dielectrics has been done by Hoshi et al. [ 451. We have directed part of our efforts to the study of the behaviour of solutes near the surface separating two immiscible liquids. The first example we have considered regards the water/benzene system, for which the modelling may be supported by the Monte Carlo studies performed by Linse [ 791. The free energy of solute depends on the distance and orientation with respect to the separation surface. We have thus introduced the concept of solvation free surface, dG (
419

without chemical groups with a strong local dipole pass into the bulk of the benzene component of the medium; molecules with important hydrophobic groups passes into the bulk water. The situation is intermediate for molecules such as C,H,OH and C2H5NH2, which prefer to stay in the vicinity of the surface. Analogous phenomena occur in other water/non-polar solvent systems: we have made some exploratory calculations for cyclohexane and carbon tetrachloride. The quantitative aspects of the phenomena are also governed from the cavitation and dispersion-repulsion contributions to the energy. A detailed examination of the AG changes during the crossing of the boundary shows some unexpected features. We give here an example drawn from an as yet unpublished wider study. The solvent system is water/benzene, the solute n-hexanol. As stated earlier, the minimum in the AG(?,co) energy surface corresponds to a buoyancy situation with the dipolar group in water. This result agrees with intuitive considerations based on electrostatic arguments. The surprises derive from the dissection of the solvation energy during the surface crossing performed at different orientations of the solute with respect to the surface. The role of the dispersion contributions is practically negligible; there are no large variations when the solute changes solvent, and there are no preferential effects, favouring the immersion of the alkyl tail in the hydrocarbon portion of the system (the role of GDIsis different, of course, for a polar solute). The essential role is played by the cavity contribution to AG (r’,o), which modifies the shape of the electrostatic contribution GEL(j?~) to AG($o) and determines the difference between different orientations. In Fig. 2 we report some numerical results. For simplicity, C&H,OH is kept at a fixed all trans geometry; o= 0 means an orientation with the long axis of the molecule perpendicular to the separation surface and the OH group pointing toward bulk water, w= 180 means the opposite orientation with OH pointing towards bulk benzene. The buoyancy state corresponds, of course, to w= 0 as the numerical values in the figure also indicate. The example shows that GEL+ GcAvpredicts at a good extent the energy trends: the other contributions mainly give a shift to the energy profile. For water/n-octanol systems, the situation could be somewhat different. The importance of the difference in the AGb values in water and n-octanol for QSAR studies [ 841 prompted us to study solvation processes in this two-liquid system. The description of a mathematical algorithm to obtain more precise AGL,sf (water/n-octanol) values with the PCM method has already been published [ 571 and other papers will follow. The water/n-octanol system is not the ideal system to study boundary effects. The solubility of water in n-octanol is far from being negligible, and it would be better to consider the n-octanol phase as a mixed solvent (this fact also has an influence on the true thermodynamic values of AGzransf(water/n-

420

.: j -2

0

2

4

ii

Fig 2. Free energy components for the passage of C6H50H through the surface of the water/ benzene system 0,O: total values for w = 0 and w = 180, respectively; 0, n *GEL+ Cc, for w = 0 and w= 180, respectively.

octanol) which cannot be directly correlated with the values obtained via partition coefficients [S&36] ). There are no Monte Carlo simulations on this system on which to base our model (remember that the water/benzene model has been based on the Linse’s calculations [ 791)) and we are obliged to rely solely on our previous experience, like that shown in the preceding example. We have considered a tentative model, according which the n-octanol molecules on the separation surface are, to a large extent, oriented. The orientation of the OH groups of the alcohol molecules seems to produce an energy barrier to the passage of a small dipolar molecule. It is advisable, however, to perform other checks before publishing a numerical appreciation of this barrier. Another important problem which remains unanswered is the description of G,,, for systems of this kind. We have indicated in the preceding example that G,,, may play an important role in the solvent transfer. The RSM model, in the Pierotti version [ 871, which we have used until now with some modifications to compute G,,, is hardly tenable for n-octanol-like molecules. We are now considering a modification of the RSM theory with the inclusion of ellipsoidal shapes for the molecules, but the internal rotational degrees of freedom and the presence of a polar head in the solvent molecules should be taken into account. With this discussion we have presented an actual example in which the combination of different methods, as stated in the Introduction, seems necessary.

421

The elaboration of integral functions, the computer simulation, and the continuum model must be combined into an unifying strategy to give a final answer to this problem. The continuum model may offer some hints, or partial solutions, to numerous other problems involving separation surfaces. We quote here some models for which we have only partial results at present. The passage of a molecule from the liquid to the gaseous phase: for hydrogen-bonded liquids (water, for example), the molecular details of the process seems to be quite complex. The shape of GFvlfor a liquid slab delimited by two other dielectrics at different E: the influence of external parameters (the thickness of the slab, the value of the two external dielectric constants, the addition of an external electric field at selected strengths, the degree of permeability of the two external portions, etc.) allows us to obtain an interesting insight on this material model, which simulates actual systems of interest in biology, physics and engineering. Solute distribution and chemical reactivity in small drops of solution (aerosols, foams): the chemical composition, the size of the drop, the temperature, the rate of increase or decrease of the drop size are among the parameters which we include in our models to study this kind of problem. A mere enumeration of models under study is not particularly interesting, and we cannot supplement this enumeration with worked-out examples, also for reasons of space. We hope that what has been reported above is sufficient to show the potential interest of the application of the PCM to problems involving liquid boundaries. 8. CONCLUSIONS

As was said before, we have chosen examples of the application of the PCM to problems regarding “chemistry”, in its wider acception, in solution. In this selection we have sacrificed both the most direct applications concerning the examination of solvent effects on the properties of solutes (a recent workedout example may be found in ref. 88) or chemical reactions (a recent paper [ 891 on a reaction extensively studied also by Monte Carlo techniques [ 901 may be mentioned here as an example ), and the dynamical problems on which we are at present engaged (nothing has been thus far published on this subject). The basic message we are trying to convey with this paper concerns the potentialities of the continuum approaches. Other versions, better than the PCM version we have developed, may surely be devised. It remains that the recent innovations have opened new perspectives to the glorious old continuum electrostatic model. The use of ab initio quantum molecular techniques seems to us a necessary requisite, at least in the first stages of the elaboration of models. Simpler formulations, based on semiclassical reformulations of the quantum models [ 911,

422

will surely follow. It seems to us inconvenient to persevere, at present, with the amelioration of old models. A comparison with the theory of isolated molecules may he instructive: the molecular models of 40 years ago are now almost completely forgotten. The advent of computational quantum chemistry has produced other models, more realistic and more effective (we quote as representative examples the models elaborated by Bader and coworkers [921, the models based on localized orbitals [91] and the models based on the new formulation of the VB theory [93] ). The process of renovation of models for molecules in solution has been slower, for the intrinsic difficulties of the problem. Only recently molecular simulations, based on quantum mechanical interaction potentials, and continuum methods, again based on a quantum mechanical molecular formalism, have been initiated to provide the necessary input for the elaboration of new simple models. The meaning of “simple model” has changed in recent decades. In the past, a simple model was equivalent to a simple analytical expression, sufficiently simple to give the required numerical output by the use of a slide-rule or an analogous device. At present, everybody has the capability of running programs on a small computer, and the “analyticity”, according to the meaning given above to this term, is no more an essential requisite. The “simplicity” of the model means a reduction of the mathematical machinery and of the physical aspects of the original model (see ref. 94 for a fuller analysis of these points) to the essential features, transparent and acceptable from the point of view of the physical interpretation of the phenomena, and not requiring excessive computational efforts. This criterion of simplicity is met, in our opinion, by most of the continuum solvation models employed today. ACKNOWLEDGMENT

The support of the CNR (Prog. Final&. acknowledged.

Chimica Fine) is gratefully

REFERENCES 1 2 3 4 5 6

J.A Barker and D. Henderson, Rev. Mod. Phys., 48 (1986) 587. J.P. Hansen and I.R. McDonald, Theory of Simple Liquids, Academic Press, New York, 1986. C.G. Gray and K.E. Gubbin, Theory of Molecular Fluids, Clarendon Press, Oxford, 1984 E. Clementi, Computational Aspects for Large Chemical Systems, Spnnger-Verlag, Berlin, 1980. W.F. van Gunsteren and P.K. Werner (Eds.), Computer Simulation of Biochemical Systems, ESCOM, Leiden, 1983. G. Crccotti and W.G. Hoover (Eds.), Molecular Dynamics Simulation of Statistical Mechanical Systems, North-Holland, Amsterdam, 1986.

423

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 51 52 53

M.P. Allen and D.J. Tildesley, Computer Simulation of Liquids, CIarendon Press, Oxford, 1987. G. Klopman, Chem. Phys. Lett., 1 (1967) 200. H.A. Germer, Theor. Chim. Acta, 34 (1974) 145. R Constanciel and 0. Tapia, Theor. Chim. Acta, 48 (1978) 75. R Constanciel and R. Contreras, Theor. Chim. Acta, 65 (1984) 1. R. Contreras and A. Amman, Int. J. Quantum Chem., 27 (1985) 293. M. Kasha, J. Chem. Phys., 20 (1952) 71 J.K Kirkwood, J. Chem. Phys ,2 (1934) 351. S Mlertul, E Scrocco and J. Tomasi, Chem. Phys., 55 (1981) 117. M J Huron and P. Claverie, J. Chem Phys,, 76 (1972) 2123; 78 (1974) 1853,1862. M. Born, Z. Phys., 1 (1920) 45. L. Onsager, J. Am. Chem. Sot., 58 (1936) 1486. D.L. Beveridge and G W. Schnuelle, J. Phys Chem., 79 (1975) 2562. M. Abraham, J. Liszi and L. Meszaros, J. Chem. Phys., 70 (1979) 2491 S Miertus and J. Tomasi, Chem. Phys., 65 (1982) 239. J Warwicker and H.C. Watson, J. Mol. Biol., 157 (1982) 671. P.B. Shaw, Phys. Rev. A, 32 (1985) 2476; 35 (1987) 2254. A.A. Rashm and K Namboodui, J. Phya. Chem., 91 (1987) 6003. A.A. Rashm, J. Phys. Chem., 94 (1990) 1725. M L.J Drummond, J. Chem. Phys., 88 (1988) 5014,502l. R J. Abraham, B.D. Hudson, M.W. Kermode and J.R Miles, J. Chem Sot. Faraday Trans. 2,84 (1988) 1911. R.J Zauhar and R S. Morgan, J. Comput Chem., 9 (1988) 171; 11 (1990) 603 M.L Connolly, J. Am. Chem. Sot., 107 (1985) 118. J.L. Pascual-Ahuir, E. Sda, J. Tomasi and R. Bonaccorsi, J Comput. Chem ,8 (1987) 778. E. Scrocco and J Toman, Adv. Quantum Chem., 11 (1978) 116. H L. Friedman, Mol. Phys., 29 (1975) 1533. W H. Orttung, Ann. N.Y. Acad Sci., 303 (1977) 22; J. Am. Chem. SOC,100(1978) 4369. M.K. G&on, K.A Sharp and B.H. Honig, J. Comput. Chem., 9 (1987) 327. J. Langlet, P Claverie, J. Caillet and A. Pullman, J. Phys. Chem ,92 (1988) 1617. M.E. Davis and J.A. McCammon, Chem. Rev ,90 (1990) 509. D. Rmaldi and J L. Rivail, Theor Chim Acta, 32 (1973) 57. J L. Rwail and D. Rmaldi, Chem. Phys., 18 (1976) 233 J L. Rival1 and B. Terryn, J. Chim. Phys., 79 (1982) 1 D. Rinaldi, M.F. Ruiz-Lopez and J.L. Rivail, J. Chem. Phys., 78 (1983) 834. 0. Tapia and 0. Goschmski, Mol. Phys., 29 (1975) 1653. 0. Tapia, F. Sussman and E Poulain, J. Theor Biol., 71 (1978) 49. J. Hilton-McCreery, R.E. Christoffersen and G.G. Hall, J. Am. Chem Sot., 98 (1976) 7198. R Bonaccorsi, R. Cimiraglia and J Tomasi, J. Comput Chem., 4 (1983) 567. H. &en, C.M. Llanos and K V. Mikkelsen, Chem. Phys., 115 (1987) 43. H. Hoshi, M. Sakurai, Y. Inoue and R. ChfijB, J. Chem. Phys., 87 (1987) 1107. G Karlstrom, J. Phys. Chem., 92 (1988) 1315. 0 Tapia, in H Ratalczak and W J. Orville-Thomas (Eds.), Molecular Interactions, Vol. III, Wiley, Chmhester, 1982, p. 47. B. Lmder, J. Chem. Phys., 33 (1960) 668; 35 (1961) 371; 37 (1962) 963. M.A Aguilar and F.J Olivares de1 Valle, Chem. Phys., 129 (1989) 439. D. Rmaldi, B.J. Costa Cabral and J.L. Rivail, Chem. Phys. Lett., 125 (1986) 495. F. Floris and J. Tomasi, J. Comput. Chem., 10 (1989) 616. F Floris, J.L. PascuaI-Ahmr and J. Tomasi, J. Comput. Chem., m press. C. Huiszoon and F. Mulder, Mol. Phys., 38 (1979) 1497; 40 (1980) 249.

424 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 82 83 84 85 86 87 88 89 90

91 92 93 94

Y.K. Kang and M.S. Jhon, Theor. Chim. Acta, 61 (1982) 41. A.I. Kltaygorodski, Tetrahedron, 14 (1961) 230. K. Mirsky, in H Schenk, R. Olthof Mazekamp, H. van Koningsveld and G.C. Baesi (Eds.), Computmg m Crystallography, Delft University Press, Delft, 1978. R. Bonaccorsi, F. Flons and J. Tomasi, J. Mol. Liquids, 47 (1990) 25. A. Tani, F. Florie and J. Tomasi, to be published. F.J Olivares de1 Valle and J. Tomaei, Chem Phys., 150 (1991) 139. F J. Olivares de1 Valle, R. Bonaccorsi, R. Cammi and J Tomaei, J. Mol. Struct. (Theochem), in press. M.A Aguilar, F J. Olivares de1 Valle and J. Tomasi, Chem Phys ,150 (1991) 151. R. Bonaccorzi, R. Camml and J. Tomasi, J. Comput. Chem., 12 (1991) 301. R Cammi, R. Bonaccorsi and J. Tomasi, Theor. Chim. Acta, 68 (1985) 271. R. Cammi and J. Tomasi, Theor. Chim. Acta, 69 (1986) 11. S. Tolosa, J.J. Esperilla and F.J. Olivares de1 Valle, J. Comput. Chem., 11 (1990) 576 S.K Luishin, S.-Y. Lm and C.E. Dykstra, J. Chem. Phys., 84 (1986) 2720. S.K. Lulshin and C E. Dykstra, J. Comput. Chem., 8 (1987) 81. R Camml, F.J. Ohvares de1 Valle and J. Tomasi, Chem. Phys., 122 (1988) 63. R. Bonaccorsl, R. Cammi and J. Tomasi, Int. J. Quantum Chem., 29 (1986) 373. K. Kitaura and K. Morokuma, Int. J. Quantum Chem., 10 (1976) 325. R. Bonaccorei, R. Cimlraglia, P. Palla and J Tomasl, Int. J. Quantum Chem., 26 (1983) 307. S.F. Boys and F. Bernardi, Mol. Phys., 19 (1970) 553. G. Alagona, C. Ghlo, R. Cammi and J. Tomael, in J. Maruani (Ed.), Molecules in Physics, Chemistry and Biology, Vol. II, Kluwer Press, 1988, p. 507. R. Bonacconi, E. Scrocco and J. Tomasl, Int. J. Quantum Chem., 29 (1986) 717. G. Alagona, R. Bonaccorsl, C. Ghlo, R. Montagnani and J. Tomasi, Pure Appl Chem., 60 (1988) 231. J. Tomasi, G. Alagona, R. Bonaccorsi and C Ghlo, in Z. MaksiC (Ed ), Modelling of Structures and Properties of Molecules, Elhs Horwood, Chichester, 1987, p. 330. M.H Abraham and J. Liszt, J. Chem. Sot. Faraday Trans. 1,74 (1978) 1604. R. Bonaccorsi, M. Hodogek and J. Tomasi, J. Mol. Struct. (Theochem), 164 (1988) 105. P. Lmse, J Chem. Phys., 86 (1987) 4177. A. Ben-Naim, J. Chem. Phys., 82 (1978) 792. R. Bonaccorsi, E. Ojalvo, P. Palla and J. Tomaei, Chem. Phya., 143 (1990) 243. R. Bonaccorsi, E. Ojalvo and J. Tomasi, Collect. Czech. Chem Commun., 53 (1988) 2320 N. Israelachvili, Intermolecular and Surface Forces, Academic Press, London, 1985 C Hanech and A. Leo, Substltuent Constants for Correlation Analysis m Chemistry and Biology, Wiley, New York, 1979. P. Berti, S Cabam, G. Conti and V. Molhca, J. Chem. Sot. Faraday Trans. 1,82 (1986) 2547. L Bemazzani, Thesis, University of Pisa, 1990. R.A. Pierotti, Chem. Rev., 76 (1976) 717. G. Alagona, C. Ghio, J. Igual and J. Tomasi, J. Mol. Struct. (Theochem ) ,204 ( 1990 ) 253. C. Alemgn, F. Maseras, A. Lledos, M. Duran and J. Bert&, J. Phys. Org. Chem., 2 (1989) 611. J. Chandrasekhar, S.F. Smith and W.L. Jorgensen, J. Am. Chem Sot., 107 (1985) 154 J. Tomasi, G. Alagona, R. Bonaccorsi, C. Ghio and R. Cammi, in Z.B Maknb (Ed.), Theoretical Models of Chemical Bonding, Vol. 3, Springer, Berlin, in press. R.F.W. Bader and T.T. Nguyen-Dang, Adv. Quantum Chem., 14 (1981) 63. D.L. Cooper, J. Gerratt and M. Raimondi, Adv. Chem. Phys., 69 (1986) 318. J. Tomasi, J. Mol. Struct. (Theochem), 179 (1988) 273.