Structure and Reactivity in Aqueous Solution - American Chemical

Chapter 2

Application of Continuum Solvation Models Based on a Quantum Mechanical Hamiltonian

Downloaded by UNIV OF LEEDS on May 17, 2015 | http://pubs.acs.org Publication Date: September 29, 1994 | doi: 10.1021/bk-1994-0568.ch002

J. Tomasi Department of Chemistry and Industrial Chemistry, University of Pisa, Via Risorgimento 35, 56126 Pisa, Italy

This chapter discusses areas of study of chemistry in solution for which continuum models may be used profitably. A distinction among types of projects is introduced. The projects may be computational applications, in which existing computer codes are used to get numerical values of a desired property, or they may be methodological studies, addressed to the implementation of further elaboration of the basic procedure. Solvation energy, reaction mechanisms, energy derivatives, and local and large scale anisotropies are the main topics considered here. The exposition is mainly based on the past experience of our group, with the inclusion of some recent developments.

The main characteristic of continuum models of solvation is their simplicity in the description of the solvent structure. This quality is quite appealing for several reasons, in particular the possibility of extending such models to treat a large number of processes and systems, the ease of interpretation of experimental and computational results, and efficiency in routine calculations. The simplicity is tempered in quantum mechanical versions with a more accurate representation of the solute by means of ab initio or semi-empirical methods. We may draw from these general statements some indications about the most promising ways of using continuum models. For this purpose we shall consider some specific topics where modelling and elaboration of efficient computer codes have different importance: the evaluation of the solvation energy, the description and interpretation of chemical reactions, and the elaboration of more detailed continuum descriptions of the solvent. These examples, to which more could be added, also show when and how quantum continuum methods are useful. This choice has been suggested by our personal experience, and we

0097-6156/94/0568-0010$08.00/0 © 1994 American Chemical Society

In Structure and Reactivity in Aqueous Solution; Cramer, C., et al.; ACS Symposium Series; American Chemical Society: Washington, DC, 1994.

2. TOMASI

Use of Continuum Solvation Models

11

will mainly rely on the use of the polarized continuum model (PCM), a computational procedure we proposed years ago (1) and are continuously refining, as some previews in the following pages will show.


Models and routine computational codes Modelling of solvent effects via continuous models dates back to the beginning of the molecular description of solutions and continues into current years. Einstein (2) in 1906 used a continuous model to describe transport properties of solutes (continuous models are not limited to electrostatic effects); the fundamentals of the theory of ionic solutions are given by the 1920 continuum model of Born (3); the basic aspects of electron transfer reactions are described by the 1956 theory of Marcus (4) invoking solvent fluctuations in a continuum model; finer aspects of similar reactions may be interpreted using a quantum description of the fast component of the continuum electrostatic polarization, as Gehlen et al. suggested in 1992 (5). These few examples, to which many others could be added, show that the persistence in the use of continuum models is accompanied by a remarkable versatility of the basic concept that surely will continue being exploited in the future. The efficiency of some recent semi-empirical methods in routine computation of solvation energies of medium size molecules (6) is well documented and accompanied by a satisfactory quality of the results. The same methods may be used to describe reactions in solution: the thermodynamic balance of a reaction may be well reproduced (6), but description and interpretation of the reaction mechanism is limited by the approximations of the semi-empirical approach. Similar considerations hold for semi-classical computational codes. These codes are not based on quantum calculations in solution, but they may rely on quantum calculations in vacuo to get the necessary parameters. Extended Born methods yeld fairly accurate solvation energies (7). Further extensions and refinements of the semi-classical codes will profit from analyses of continuum quantum results. T h e Role of ab initio C o n t i n u u m Solvation Methods. We shall consider ab initio continuum methods as the primary object of this paper. There are two reasons for this choice. The first has a methodological character. A sound strategy to test qualities, defects, and potentialities of a model consists in examining the output of the model at its best, in the most complete and detailed form, and then reducing it to a more manageable level when the essential features to be preserved are well ascertained. Shortcuts are not advisable; they may lead to wrong conclusions. The theoretical chemistry literature is rich in wrong statements based on a hurried examination of a model; some of these statements refer to continuum solvation models. When the usual homogeneous continuum model is considered, the description of the solute charge distribution and of the mutual solute-solvent interaction effects should be done at the highest possible level of accuracy to detect limits and potentialities of the approach. Analogous considerations hold when the homogeneous medium is replaced by other continuum distributions.


12

STRUCTURE AND REACTIVITY IN AQUEOUS SOLUTION

The second reason is related to the intrinsic superiority of ab initio quantum methods in describing fine details of molecular charge distributions. For some problems a detailed ab initio description is absolutely necessary. The Solvation Free Energy One of the basic problems in the study of solutions is the development of models able to describe their thermodynamic properties, for example the free energy of solvation, A G i , of neutral solutes at infinite diluition. From the application of ab initio continuum methods we may state the following points: 1) Calculations of A G i (and of related quantities) with continuum models are able to yield results within the range of experimental errors. 2) These calculations require the use of a well shaped cavity. The use of cavities with simple shapes can be accepted only under limited conditions. 3) Electrostatic (Coulomb and polarization), dispersion, cavitation, and repulsion terms are all necessary: A G i = AG i+Gdis+Gcav+G epThe relative importance of these terms may be related to the bulk properties of the solvent and to the molecular properties of the solute. In making comparisons among solutes of the same class some terms may be neglected, but this choice must be based on the information derived from the decomposition of complete ab initio results. 4) The extra energy effects due to local variations induced in the solvent distribution by the solute (cybotactic changes) are of limited entity. 5) The contributions due to terms related to the vibrational, rotational and translational solute partition functions are not decisive for almost rigid solutes. Some corrections due to large-amplitude motions (hindered internal rotations, out-of-plane deformations) and to zero point contributions (stretching of M-H groups making hydrogen bonds with solvent molecules) may be easily introduced, when necessary. The semi-empirical procedures we have mentioned satisfy points 2, 3 and 4. The analysis of ab initio results justifies their formulation and suggests some minor improvements.


so

s o

8 0

e

r

N u m e r i c a l Results. We report here the results of a linear regression analysis of computed against experimental AGhyd values. The computed values refer to the P C M ab initio and to the AMSOL-SM2 semi-empirical procedures. The data reported in the Table are the coefficients of the regression line: AGhyd(exp) = a AGhyd(comp) + b (Kcal/mol), the regression coefficient R, the standard error σ and the number η of cases within each set. Set 1 refers to a sample of chemically similar compounds (esters) all described at the 6-31G* SCF level with geometry optimisation in water (8). Set 2, presented here for the first time, it is not chemically homogeneous (neutral organic solutes with heteroatoms) and is computed using the same 6-31G* basis set with geometries optimised in vacuo. Some of the


2. TOMASI

13

Use of Continuum Solvation Models 1

"experimental' values are drawn from the empirical formulas elaborated by Cabani et al.(9). Set 3 includes all the neutral polar solutes used by Cramer and Truhlar in their calibration procedure (β), but not the charged ones. The last set has been also computed by Cramer and Truhlar, but it was not used in the calibration. These results show that ab initio calculations may reach chemical accuracy and that semiempirical values are of a good level.


Table I. Model Predictions against Experimental Data Set

Method

a

b

R

σ

η

1 2 3 4

PCM PCM AMSOL AMSOL

1.03 0.98 0.90 0.98

0.33 0.35 -0.02 0.62

0.999 0.910 0.876 0.965

0.015 0.452 1.207 0.796

16 102 117 7

Good results are also obtained with other semiempirical procedures. The results are somewhat scattered and we refer to a still unpublished review by Cramer and Truhlar (Cramer C.J.; Truhlar, D.G. Reviews in Computational Chemistry 1994, in press) for more information. Note, however, that this field is in rapid evolution and that the number of methods and the quality of the results are increasing rapidly. The reduction of the model may be carried further. Even semiclassical models give results of appreciable quality (see again the review by Cramer and Truhlar quoted above). A sequence of approximations starting from the quantum description and ending with atomic charges may be a guide to check this reduction of the model without shortcuts (10). The continuum methods compare well with other approaches. The first systematic comparison between ab initio continuum and M D based free energy perturbation calculations of AGhyd performed by Orozco, Jorgensen and Luque (11) gives an average error of 0.8 kcal/mol for the continuum P C M procedure and 1.5 kcal/mol for the M D - F E P technique with respect to the experimental values. Of course many refinements may be introduced in this picture, but details are not essential here: we conclude that computationally very convenient continuum methods may be used to obtain gas-liquid and liquid-liquid transfer thermodynamical properties. Chemical Reactions with the Continuum Model The evaluation of solvation energy is but a tiny part of the topics facing theoretical chemistry in solution. We shall consider now a more challenging subject, the description of chemical reactions. Here the ab initio formulation of the continuum model has an important role. Let us summarise the formal set-up of the approach. Continuum methods are able to give an evaluation of the free energy G of the solute in In Structure and Reactivity in Aqueous Solution; Cramer, C., et al.; ACS Symposium Series; American Chemical Society: Washington, DC, 1994.

14


solution (not to be confused with A G i ) as a function of the nuclear coordinates {R) of the solute ( as "solute" we consider the bare reacting molecules, supplemented when necessary with a restricted number of solvent molecules playing an active role in the reaction). The G(R) surface corresponds to the E(R) surface defined in vacuo. We apply to the G(R) surface (or, better, to the set of relevant G(R) surfaces) the same concepts originally derived for reactions in vacuo. First, we define geometry, energy and electronic structure of reactants, products, intermediates, and saddle points. Then, we define the reaction coordinate and the portions of the energy surface near the reaction path necessary for dynamical studies of the reaction. The chemical interpretation of the mechanism will be based on a scrutiny of the solute wave function, to be performed with suitable techniques. This simplified version of a rather formidable problem (a quantum study of a system composed of a large number of molecules) is corroborated by preliminary tests on the model. In particular the partition of the whole system into a "solute" and a medium is supported by the success in describing the energy profile of several significant reactions. This scheme must be supplemented when certain dynamical effects of the solvent are considered. It is known that in many important classes of reactions the dynamics cannot be properly described by relying on the G(R) surface alone: the typical case is the outer-sphere electron transfer model of Marcus (4), in which the dynamics is carried by a solvent coordinate alone, without intervention of the geometric coordinates of the solute. We are thus compelled to extend our definition of energy hypersurface G(R) by including some extra dynamical coordinates {S}, and to use, when and where necessary, a more general function G(R+S). This enlargement of the space is not equivalent to the addition of more solvent molecules in the "solute". We have thus recognised at least three important problems: the correct evaluation of G(R) and of its critical points, the description of the electronic structure at some significant points, and the definition of the additional {S} subspace, supplemented by the protocols for the use of G(R+S)


so

The Analysis of the G(R) Surface. Long experience in the evaluation of E(R) surfaces tells us that semiempirical methods give, in the most favourable cases, nothing more than a first order guess. Ab initio calculations, of good quality, are necessary. The analysis of the G(R) surface is made easier by the use of the gradient of G(R), grad G(R), supplemented by the diagonalization of the Hessian matrix H(R). The calculation of grad G and of H must be performed according to the conditions set in points 2 and 3 above, i.e. using a suitably shaped cavity and including in G(R) all the necessary contributions: G i(R) + Gdi (R) + G ( R ) + G ( R ) . Current ab initio continuum programs are not equipped for the analytical evaluation of gradG and H at this level of accuracy. We present here the computational scheme we have recently elaborated in the framework of the P C M (Cammi, R.; Tomasi, J., J. Chem. Phys., May 1994). e

s

c a v

rep


2. TOMASI

15

Use of Continuum Solvation Models

Iterative a n d Direct P C Method. The original iterative P C M version is not suited to get analytical derivatives of G i (the other terms of G are simpler to manage and may be treated separately). It is more convenient to resort to a matrix formulation of the electrostatic problem, exploiting the fact that in the PCM the solute is effectively replaced by a charge distribution σ on the cavity surface, and that this surface is divided into a finite number of tesserae. An accurate elaboration of the matrixP C M has been done by the Sakurai group (12). We have elaborated a similar procedure, more computationally effective and suitable to compute the first and second derivatives of G i with respect to parameters α and a, β, later indicated by G and G P respectively. This preliminary step highlights some points that deserve mention. In the iterative procedure the Schrodinger equation e

e


a

A

|(Η0 + ν ) Ψ > = Ε | Ψ >

(1)

σ

may be solved with the traditional variational techniques because V, and Vq to which nuclear contributions must be added (1):

G = - 1 / 2 < Ψ | ν | Ψ > + ν Ν Ν + 1 / 2 υ . σ

σ

Ν σ

(2)

Here V N N is the intrasolute nuclear repulsion term also present in vacuo calculations and U N O is the interaction between the solute nuclei and the apparent surface charge σ. In the direct methods simultaneously optimizing the solute and solvent charge densities, the functional to be minimized is not the mean value of H° + ν , but rather the free energy fuctional G (13). When the problem is recast at the Hartree-Fock level, with expansion over a finite basis, the minimization of G is reduced to the solution of a pseudo-HF equation (12): σ

f

F C = ESC

(3)

with F» = h» +G'(P) = (h + 1/2(J + Y)) + (CKP) + X ( P » .

(4)

Here h and G(P) are the usual H F one-electron and two-electron integrals matrices, Ρ is the density matrix, J , Y and X(P) collect one-and twoelectron integrals involving interactions with the surface apparent charges. (To be more precise, J describes the interactions between the elementary solute electron charge distributions and the surface charges having as origin the solute nuclear charges, Y describes the interactions between the solute nuclei and the surface charges having as source the elementary In Structure and Reactivity in Aqueous Solution; Cramer, C., et al.; ACS Symposium Series; American Chemical Society: Washington, DC, 1994.

16


solute electronic distributions, X describes the interactions between elementary electron distributions and the surface charges they generate). The expression of F is different from that of the Fock matrix related to equation (1) and applied in the iterative procedure. The two approaches are formally equivalent (Cammi, R.; Tomasi, J . J. Comp. Chem., to be published) and must give the same values for G and for the coefficients C of the wave function.


1

The Renormalization of the Surface Charges. This formal equivalence is numerically verified only when a renormalization of the surface charge σ (or of the point charges q which describe σ) is introduced in the computational scheme. The integrated value of σ must in fact satisfy a simple relationship with the total charge Q M of the solute:

Jo(s)ds =(ε-ΐνε Q .

(5)

M

Similar relationships hold for the electronic and nuclear components of the surface charge (σ = o + σ ) having as sources the solute electron and nuclear charge distributions ( Q M = Q M + Q M )· These conditions are not satisfied by continuum quantum calculations (one reason is that a portion of the electronic charge distribution is spread out of the cavity, by definition). The effect of this lack of normalization is different in the iterative and direct methods and results without renormalization may differ. On this basis it has been stated that there are two alternative pictures for the description of the solvation energy. Actually, there is only one picture, and the choice is between two alternative but equivalent computational methods. The introduction of a suitable renormalization permits one to recover in actual computations the equivalence of the two methods. In addition, the renormalization permits to exploit the formal equivalence between J and Y (the first, as said, describes the interactions of the elementary electronic solute charge distributions with σ , the second describes the interactions of nuclei with 0 ) that is lost in un-normalized calculations, with a further reduction of the computational times. With this reformulation the ab initio direct method is almost as fast as the iterative method for solutes of small size. e

Ν

e

N

Ν

e

A n a l y t i c a l Derivatives w i t h the P C Method. Surface charge renormalization plays an even more important role in the analytical determination of derivatives. Without renormalization no meaningful analytical derivatives may be computed in the ab initio apparent surface charge methods. The expression of G , obtained using equation (3), is: 1

f

G = trPh + 1/2 trPG (P) + V N N

(6) N

where V ' N N = V N N + 1/2 U N N collects nuclear repulsions and nuclei-o contributions to G . From this definition of G we may derive formal expressions for its first and second derivatives with respect to the cartesian components α and β of the solute nuclear coordinates: In Structure and Reactivity in Aqueous Solution; Cramer, C., et al.; ACS Symposium Series; American Chemical Society: Washington, DC, 1994.

2. TOMASI

17

Use of Continuum Solvation Models f

G« = trPh*« + 1/2 trPG a(P) - trS

Structure and Reactivity in Aqueous Solution - American Chemical

Recommend Documents