1. Introduction

This review on continuum solvation models has been preceded in Chemical Reviews by others addressing the same subject. They are due to Tomasi and Persico¹ (published in 1994), Cramer and Truhlar² (published in 1999), and Luque and Orozco³ (published in 2000). These three reviews on the same topic in a journal covering all of the aspects of chemical research indicate the interest this topic has for a sizable portion of the chemical community.

Liquid solutions play in fact a fundamental role in chemistry, and this role has been amply acknowledged by Chemical Reviews since its very beginning. In the abundant number of reviews addressing different aspects of chemistry in the liquid phase, the number of those centered on the theoretical and computational aspects of the study of liquid systems has considerably increased in the past two decades. This reflects the increasing importance of computational approaches in chemistry, an aspect of the evolution of scientific research chemistry shares with physics, biology, engineering, geology, and all of the other branches of sciences, an evolution that is ultimately due to the widespread availability of efficient computers.

Computers have permitted the activation of many approaches in sciences that were dormant, or limited in their applications to the level of simple model, for the lack of appropriate computational tools; an example in chemistry is given by the activation of methods based on quantum mechanics for the description of isolated molecules, which have now reached the degree of accuracy in the description of molecular structures that chemists require.

Computers have also completely modified the ways of doing theoretical and computational studies of liquid systems, permitting the introduction of new approaches, new concepts, and new ideas. The most important innovations in this field are related to the use of computer simulations that directly, or indirectly, are the basis of our present understanding of condensed systems. Simulations, initiated about 50 years ago, have greatly evolved in the past 20 years and proceed now in covering all of the fields in which condensed matter occurs. One aspect of this evolution is of direct interest here: the merging of simulations with quantum mechanical (QM) descriptions of molecular structures. This merging, in progress for several years, has to overcome computational difficulties, due to the computationally quite intensive numerical procedures.

At this point of our rapid exposition we may introduce what we consider the second important innovation in the field of condensed systems made possible by the use of computers: the continuum models.

Continuum models were introduced more than a century ago, in very simplified versions, giving results of remarkable importance using computational instruments no more complex than a slide rule. These models have been used for more than 50 years, and in more detailed versions they are still in use, but the essential step, opening new perspectives for the study of solvent effects, has been its merging with the QM descriptions of molecules. Continuum models are in fact the ideal conceptual framework to describe solvent effects within the QM approach, as will be shown in section 2.

This merging was initiated more than 30 years ago by a small group of young researchers (Claverie, Rivail, Tapia, and Tomasi), working in French and Italian laboratories gravitating around Paris and Pisa. The initial stimulus was provided by the recognition that the QM description of the electrostatic potential generated by the charge distribution of a molecule could represent a valid analytic and interpretative tool to study intermolecular interactions.⁴

The perspectives opened by these findings were elaborated independently by the various laboratories and led to the alternative theoretical methods^5-8 that have been examined in detail in the preceding reviews. In the present one we abandon the historical perspective used in the first review¹ to present and analyze the methodological issues of the "modern" continuum solvation theory (modern in the sense that is essentially based on the QM description of the solute). The second² and third³ reviews give considerable space to methodological issues, but also pay attention to applications: detailed surveys on the results on chemical equilibria, spectra, and dynamics of reactions are reported in the former, whereas, in the latter, attention is centered on biomolecular systems, also including in the survey simulation methods with QM description of the solute.

The present review is again centered on methodological issues, with a perspective focused on the most recent developments. The body of the theory developed in the first 30 years of the "modern" solvation methods is summarized, giving emphasis to those aspects that are at the basis of all recent extensions and reformulations of the model. Some among them were already considered in the previous reviews, reflecting, however, the provisional state of the methodological elaboration, now accomplished; others correspond to new entries in the theory.

With respect to all of the previous reviews, here a much larger space is devoted to three aspects that we consider to be important in future applications of continuum models, namely, their use in studying phenomena involving time-dependent solvation, on the one hand, and molecular properties, on the other hand, and their coupling to discrete models. None of these three aspects is completely new, but only in the past few years have they acquired the necessary reliability and accuracy. In parallel, the extension of the continuum models to treat complex condensed phases (ionic solutions, anisotropic dielectrics, heterogeneous systems, liquid-gas and liquid-liquid interfaces, crystals, etc.) has allowed us to enlarge the range of possible applications and to consider phenomena and processes that have been until now the exclusive property of computer simulations. The table of contents should give a schematic but clear proof of this new trend in continuum solvation methods.

Throughout the text, we have also added remarks indicating other extensions of the continuum methods that are still in their infancy, or even in an earlier stage, but which seem to us to be possible and to promise a satisfactory reward. A further analysis of future prospects will be done in the last section dedicated to comments and conclusions.

We conclude this introductory section with a short exposition of the main features of the computational strategies based on continuum distributions. The considerations reported here do not claim originality, but we find it convenient to report them to put continuum approaches in the right perspective.

A continuum model in computational molecular sciences can be defined as a model in which a number of the degrees of freedom of the constituent particles (a large number, indeed) are described in a continuous way, usually by means of a distribution function.

Continuum distributions are a very general concept. In the standard quantum mechanical description of a single molecule M based on the usual Born-Oppenheimer (BO) approximation (the cornerstone of the theory for molecular sciences), the electronic wave function, , is expressed in terms of one-electron wave functions, each depending on the coordinates of a single electron. From this single-particle description a one-particle distribution function is easily derived with an averaging operation, the one-electron density function (r), which contains a good deal of the information conveyed by the original wave function. According to the formal theory, one-electron (r) and two-electron (r,r') density functions collect all of the elements necessary for a full exploitation of the QM basic calculation. Electron density functions are endowed with many important formal properties, as they represent the kernels of integral equations from which properties can be derived. Actually, the formalism of the density matrices is completely equivalent to (and even more powerful than) the usual wave function formalism.

By tradition, in basic quantum mechanics, the emphasis is not placed on the continuity of the electronic distribution but on the discreteness of the molecular assembly. We are here interested in the application of the concept of continuum distribution functions to particles of different physical types, including electrons and nuclei as well. Continuum distributions of this kind, often supplemented by constraints acknowledging the existence of molecules, are of current use in statistical mechanics and find application, among others, in computer simulations of pure liquids and solutions.

The continuum models we shall consider are intermediate between the two extremes we have mentioned, a continuum model for the electrons of a single molecule, and a continuum model for a very large assembly of molecules. Our aim is to preserve the accuracy of the former in describing details of the molecule and the capability of the latter in strongly reducing the degrees of freedom of large molecular assemblies. To do it in a proper way it is convenient to introduce another concept, that of focused models.

In focused models the interest of the enquirer falls on a limited portion of the whole system. There is a large variety of systems and properties for which the focusing approach can be profitably exploited. A single solute molecule in a dilute solution is just an example, but many other examples can be cited: a defect inside a crystal, the superficial layer of a solid, the active part of an enzyme, the proton or energy transfer unit in a larger molecular assembly, and even a single component of a homogeneous system, as a single molecule in a pure liquid or in a gas.

The definition of a focused model presents specific problems for each different case; here, we shall highlight some general aspects referring to solvation models.

The concept of a focused model can be translated into a simple formal expression. The whole system is partitioned into two parts, which we define as the focused part F, and the remainder R. The Hamiltonian of the whole system may be written as

where {f} and {r} indicate the degrees of freedom of the F and R parts, respectively. To focus the model means to treat the F part at a more detailed level than the R part. An important parameter in focused models is the number of the degrees of freedom of R, which are not explicitly taken into account. In continuum solvation models the whole ^R(r) term is eliminated and the total Hamiltonian is reduced to an effective Hamiltonian (EH) for the solute in the form

In this approach in fact, there is no need to get a detailed description of the solvent, it being sufficient to have a good description of the interaction. This is surely a considerable simplification. In other solvation models ^R(r) is maintained and the focusing simplifications involve the number of freedom degrees within R. This is the case of the QM/MM methods, in which a limited number of the degrees of freedom for each molecule within R are explicitly considered, at least those involving position and orientation.

The elimination of the solvent Hamiltonian is not sufficient to eliminate the {r} degrees of freedom, because they appear in the interaction Hamiltonian. An almost complete elimination can be obtained by introducing an appropriate solvent response function that we indicate here with the symbol Q(,'), where () is no more the whole set of solvent coordinates, but just a position vector.

The solvent response function is similar in nature to the response functions that can be derived from the electron density function in the case of a molecule. The density here is that of the liquid system R, and Q(,') is, in the more complete formulations, expressed as a sum of separate terms each related to a different component of the solute-solvent interaction.

We shall be more specific in section 3, but for this general discussion we limit ourselves to the electrostatic contribution. The response to be considered in this term is that with respect to an external electric field. As in classical electrostatics, the polarization function of the medium is proportional to the external field

and the kernel response function we have to use is related to the function we have here introduced, namely, the permittivity, . The expression of may be quite simple, just a numerical constant, or more complex, according to the model one uses. This point will be fully developed in section 4, in which also the basic assumption of linearity introduced with eq 4 will be reconsidered, but here it is sufficient to remark that the whole set of solvent coordinates {r} is replaced by a function depending only on one parameter, the position vector , or by a couple of position vectors (,').

The quest for the accuracy requested in chemical applications of such focused models suggests the introduction of additional features, different according to the various cases. Within this very large variety of proposals a third concept emerges for its generality, the concept of layering.

Layering can be considered a generalization of focusing. The material components of the models are partitioned into several parts, or layers, because often these parts are defined in a concentric way, encircling the part of main interest. Each layer is defined at a given level of accuracy in the description of the material system and with the appropriate reduction of the degrees of freedom. There is a large variety of layering; for simplicity, they can be denoted with abbreviations, as, for example, QM/QM/Cont or QM/MM/Cont for a couple of three-layer models in which the inner layer is treated at a given QM level, the second at a lower QM or at a molecular mechanics level, respectively, and the third using a continuum approach. Some type of layering involving as chain end the continuum description will be examined in section 7.

2. Methodological Outlines of the Basic QM Continuum Model

This topic has been amply presented in the first review on the argument,¹ complemented in a second review,² and summarized in many other places, including textbooks.^9-12 For this reason we shall avoid an exhaustive presentation of the whole literature but shall only present the essential elements. By contrast, we shall pay attention to all topics that have been the subject of discussion in recent years and/or new methodological proposals.

In the next subsection we shall identify the characteristics that define the "basic QM continuum model": these in principle completely determine the model we are interested in but, in practice, they are seldom completely fulfilled. We thus suggest that the reader consider this definition more as a "literal" convenience and a research of formal clarity rather than as a description of a real methodology. As a matter of fact, in this section, there will be many occasions to present models that do not completely conform to this basic definition.

The cavity is a basic concept in all continuum models. The model in fact is composed of a molecule (or a few molecules), the solute, put into a void cavity within a continuous dielectric medium mimicking the solvent. The shape and size of the cavity are differently defined in the various versions of the continuum models. As a general rule, a cavity should have a physical meaning, such as that introduced by Onsager,¹³ and not be only a mathematical artifice as often happens in other descriptions of solvent effects. On the physical meaning of Onsager's cavity, see also the comments in ref 14. In particular, the cavity should exclude the solvent and contain within its boundaries the largest possible part of the solute charge distribution. Here, for convenience, we divide it into its electronic and nuclear components:

In continuum models, much attention has been paid to the portion of solute electronic charge outside the boundaries of the cavity; the terms "escaped charge" and "outlying charge" are often used to indicate this portion of charge. This subject will be treated in due detail in section 2.4.3. Here we will assume that all of the solute charge distribution lies inside the cavity, which in turn has a size not so large as to be in contrast with the solvent exclusion postulate.

The optimal size of the cavity has thus been a subject of debate, and several definitions have been proposed. The adopted definitions are the result of a tradeoff between conflicting physical requirements.

The shape of the cavity has also been the object of many proposals. It is universally accepted that the cavity shape should reproduce as well as possible the molecular shape. Cavities not respecting this condition may lead to deformations in the charge distribution after solvent polarization, with large unrealistic effects on the results, especially for properties. Here, once again, there is a tradeoff between computational exigencies and the desire for better accuracy.

Computations are far simpler and faster when simple shapes are used, such as spheres and ellipsoids, but molecules are often far from having a spherical or ellipsoidal shape.

Quantum mechanical calculations of the molecular surface can give a direct ab initio definition of the cavity.

An accurate description is based on the use of a surface of constant electronic density (isodensity surface).^15,16 Within this framework, one only needs to specify the isodensity level (typically in the range of 0.0004-0.001 au) and, thus, the cavity will be derived uniquely from the real electronic environment. Such a cavity has been inserted into the Gaussian computational package.¹⁷ Even if not largely used at the moment, in our opinion the isodensity surface represents an important definition of the cavity for continuum solvation models, and it will surely receive a renewed interest in the coming years.

A different technique pioneered by Amovilli and McWeeny¹⁸ and employed by Bentley and others^19-24 is based on the calculation of the interaction energy between a given molecule and an atomic probe (typically a rare gas atom, from He to Ar) placed at opportune positions in the outer molecular space. From these calculations, a set of three-dimensional isoenergy surfaces is determined. As described below, two kinds of surfaces are of interest, the solvent-accessible surface (SAS) and the solvent-excluded surface (SES). The first can be directly obtained from these calculations, whereas the second requires an additional assumption. According to Bentley the SES can be determined from the electronic density function of the system constituted by the molecule and the probe. Following the AIM topological analysis,²⁵ the points of the surfaces can be identified in the saddle points [also indicated as (3,-1) bond critical points] of such a function.

A connection of these surfaces with the thermal energy kT allows one to define T-dependent molecular surfaces and cavities. This technique is of potential interest as a benchmark for solvent calculations far from the ambient temperature; otherwise, the approach is too costly to be used in standard applications.

The generally adopted compromise between analytical but too simple and realistic but computationally expensive cavities is based on the definition of the cavity as an interlocked superposition of atomic spheres with radii near the van der Waals (vdW) values (the precise determination of such radii is related to the problem of the cavity size).

The most used set of vdW radii in the chemical literature is that defined by Bondi²⁶ (~5000 citations in the past 15 years). This set was confirmed as the recommended one some years ago,²⁷ after the examination of a quite large number of intermolecular contact distances drawn from the Cambridge Structural Database. For his tabulation, Bondi used data of other origin, mainly addressing the hard volume of the molecule and not the non-covalent contact distances. The data drawn from 28403 crystal structures confirm the Bondi values, with the only exception the hydrogen radius, set by Bondi at 1.2 Å, which is probably too high by 0.1 Å.

Another tabulation of vdW radii frequently used is that of Pauling,²⁸ available in many tabulations ofphysicochemical data, for example, in the CRC Handbook.²⁹

A third set of values has been inserted in the latest versions of the Gaussian package.¹⁷ It is drawn from the compilation of data for the universal force field (UFF),³⁰ and it covers the whole periodic table including groups not present in Bondi's or Pauling's tabulations.

The presence of small differences in the radius values currently used in continuum models and the consequent effect on the solvation energies deserve a few words of comment. Chemical applications of molecular calculations generally involve trends, in particular, comparisons of selected properties computed for different systems. In the particular case of solvation calculations, the comparisons may involve properties of different solutes in the same solvent or of the same solute in different solvents. Absolute values at a precision comparable with that of very accurate experiments are rarely requested. For this reason, the selection of the hard radii to use among the recommended tabulations seems to us to be not a critical issue. Attention has to be paid, however, to maintain a coherent choice for all of the calculations and to avoid comparisons among results obtained with different definitions of the radii.

Molecules often have an irregular shape, and the occurrence of small portions of space on their periphery where solvent molecules cannot penetrate is not a rare event. This intuitive consideration is at the basis of two definitions, those of solvent-excluding and solvent-accessible surface (SES and SAS, respectively).^31-35

Both introduce in the surface (and in the volume) changes to the vdW description related, in a different way, to the finite size of the solvent molecules. In both cases, the solvent molecule is reduced to a sphere, with a volume equal to the vdW volume (other definitions of this radius have been used, but this seems to be the most consistent definition). The positions assumed by the center of a solvent sphere rolling on the vdW surface of the solute define the SA surface, that is, the surface enclosing the volume in which the solvent center cannot enter. The same sphere used as a contact probe on the solute surface defines the SE surface, that is, the surface enclosing the volume in which the whole solvent molecule cannot penetrate (see Figure 1 for a schematic drawing of the different surfaces for the same molecule).

Figure 1 Solvent accessible surface (SAS) traced out by the center of the probe representing a solvent molecule. The solvent excluded surface (SES) is the topological boundary of the union of all possible probes that do not overlap with the molecule.

In the literature, the SES is also called "smooth molecular surface" or "Connolly surface", due to Connolly's fundamental work in this field. Indeed, the SE surface developed by Connolly^32,36 can be considered to be the prototype for the computational study of molecular surfaces. Visualization and handling of surfaces have given origin to a very large literature that cannot be reviewed here. The reader can be referred to Connolly's website³⁷ for a clear and concise review accompanied by a sizable selection of references, adjourned at 1996 (430 entries). Connolly's surfaces have been applied to a very large variety of problems, and they have been also used to computesolvation energies with continuum models (generally of classical type). The probe sphere divides the whole SE surface into pieces of three types: the convex patches in which the probe touches just one sphere of the hard vdW shape function, the toroidal patches in which the probe touches two spheres of the hard body, and the concave (reentrant) patches in which the probe touches three spheres. An analytical expression for this shape, easy to visualize on a computer screen with a probe provided with markers that put dots on the SES, has been given by Connolly within a short time from the first computer implementation of the procedure (1979-1981), and it is still in use, with some modifications. The analytical description presents some problems, among which we mention a few: the intersection between a torus and a sphere is described by a fourth-degree equation, for which the available solvers are not sufficiently robust; the SES may present singularities and cusps; these last problems are better treated with methods developed by Gogonea^38,39 and by Sanner et al.⁴⁰ The latter have developed a computational code called MSMS,⁴¹ standing for Michael Sanner's Molecular Surface, which has received much attention, especially among biochemistry-oriented computational researchers. MSMS computes, for a given set of spheres and a probe radius, the reduced surface and the analytical model of the SES. The MSMS algorithm can also compute a triangulation of the SES with a user-specified density of vertices.

Besides Connolly's SES, another SES-like surface will be reviewed here as of current use in QM continuum solvation methods. This alternative surface is defined following a different strategy, originally conceived in Pisa around 1984 and finalized in 1986 by Pascual-Ahuir in his Ph.D. research.⁴² This surface-building method, known as GEPOL,^43-45 is based on a sequential examination of the distance among the centers of each couple of hard vdW spheres and comparison of it with the solvent probe diameter. If the distance is such that the probe cannot pass between the two hard spheres, additional spheres are added. Only three cases are possible, each corresponding to a different positioning of the additional spheres, each with the opportune radius (position and radius are determined with very simple and unambiguous algorithms). The whole set of spheres, the original vdW spheres and those added, is subjected again to the same sequential examination to add new spheres (second-generation spheres) to smooth the surface. The program originally written by Pascual-Ahuir introduced thresholds and options to keep the number of additional spheres within reasonable limits.

In GEPOL, the final surface is thus always the result of the intersections of spheres and, in this sense, it can be seen as an alternative version of the SES made only by convex elements.

To complete the section on the definition of the cavity, we recall an alternative strategy to define van der Waals, solvent-accessible, and solvent-excluding molecular surfaces originally formulated by Pomelli in 1994-1995 for his master's thesis. This strategy, known as DEFPOL,^46,47 has never been implemented in publicly released computational packages, and thus its use is limited to a few examples; however, it still presents some aspects that are worth recalling here. The basic strategy consists of progressive deformations of a regular polyhedron with the desired number of faces (triangular faces are preferred) inscribed into a sphere centered on the mass center of the molecule. The sphere is deformed into the inertial ellipsoid and enlarged so as to have all shape functions of the molecule within it. Each vertex of the deformed polyhedron is then shifted along the line connecting the initial position with the origin of the coordinates until it lies on the surface of the shape function. The polyhedron is so transformed into a corrugated polyhedron topologically equivalent to the initial one, with faces still defined as triangles. The center of each triangle is then shifted along the axis orthogonal to the triangle until it touches the surface of the shape function. In the cases in which the volume of the tetrahedron defined by the three displaced vertices and the displaced center is larger than a given threshold, the procedure is repeated on a finer scale on the three triangles having as vertices the original ones and the triangle center. The final step consists of transforming the flat triangles into spherical triangles, each with the appropriate curvature.

The physics of the electrostatic solute-solvent interaction is simple. The charge distribution _M of the solute, inside the cavity, polarizes the dielectric continuum, which in turn polarizes the solute charge distribution. This definition of the interaction corresponds to a self-consistent process, which is numerically solved following an iterative procedure. It is important to remark that the corresponding interaction potential is the one we shall put in the Hamiltonian of the model. As this potential depends on the final value of _M reached at the end of this iterative procedure, the Hamiltonian (previously introduced as an "effective Hamiltonian", see eq 2) thus turns out to be nonlinear. This formal aspect has important consequences in the elaboration and use of the computational results (see section 2.4.4).

Reference is often made to the solvent reaction field for the interaction potential obtained with continuum models (and also with models using explicit solvent molecules). This label has a historical reason, being related to Onsager's seminal paper¹³ in which the solute was reduced to a polarizable point dipole and the electrostatic interaction between a polarizable medium and a dipole was expressed in terms of an electrostatic field, having its origin in the polarization of the dielectric. Actually, it is now convenient to speak in terms of the solvent reaction potential, because a potential is the term we have to introduce into the Hamiltonian.

The basic model requires the solution of a classical electrostatic problem (Poisson problem) nested within a QM framework. Let us consider the electrostatic problem first.

In our simplified model the general Poisson equation

can be noticeably simplified to



where C is the portion of space occupied by the cavity, is the dielectric function (actually a constant) within the medium, and V is the sum of the electrostatic potential V_M generated by the charge distribution _M and the reaction potential V_R generated by the polarization of the dielectric medium:

Observe that we have assumed that all of the real charges of the system (i.e., those described by _M) are inside the cavity. Therefore, this basic model is not formally valid for liquid systems having real charges in the bulk of the medium (as is the case, e.g., for ionic solutions) and also, strictly speaking, for the tiny portions of the electronic component of _M lying out of the cavity (see section 2.4.4 for a more detailed analysis). It is valid, however, for systems having a multiplicity of cavities C₁, C₂, ..., C_n, each containing a different _M1, _M2, ..., _Mn charge distribution.

Equations 7 and 8 are accompanied by two sets of boundary conditions, the first at infinity and the second on the cavity surface. At infinity we have

with finite values for and . These conditions ensure the harmonic behavior of the solution, but they have to be specifically invoked in just one of the methods we shall examine; in the other cases, the conditions are automatically satisfied in our basic QM model.

More important, from a practical point of view, are the conditions at the cavity surface . They may be concisely expressed as jump conditions:

The jump condition (eq 12) expresses the continuity of the potential across the surface, a condition valid also for other dielectric systems we shall consider later:

The second jump condition (eq 13) involves the discontinuity of the component of the field (expressed as a gradient of V) that is perpendicular to the cavity surface. In the model we are considering here, a cavity with a dielectric constant equal to 1 and an external medium with (a finite value >1), this condition leads to

where is the outward-pointing vector perpendicular to the cavity surface.

Equations 7-15 are the basic elements to use in the elaboration of solvation methods according to standard electrostatics. In the first review¹ the approaches in use were classified into six categories, namely, (1) the apparent surface charge (ASC) methods, (2) the multipole expansion (MPE) methods, (3) the generalized Born approximation (GBA), (4) the image charge (IMC) methods, (5) the finite element methods (FEM), and (6) the finite difference methods (FDM).

We maintain here this classification, but with the required modifications due to the important developments achieved in the past years in almost all of the categories. Only in the category of the IMC methods are no new important developments to be found, at least within the framework of molecular calculations, and thus this category of methods will be not considered in the present review: the interested reader is referred to our previous review.¹

2.3.1. ASC Methods

From the jump condition (eq 15) one may derive an auxiliary quantity that defines all of the ASC methods: an apparent surface charge (s) spread on the cavity surface. We are using the symbol s for the position variable, to emphasize that this charge distribution is limited to the surface .

The definition of the ASC is not unequivocal but changes in alternative versions of the model. In all cases, however, the ASC defines a potential over the whole space:

This potential is exactly the reaction potential V_R of eq 9. There are no approximations in this elaboration: the definition of V_R given with eq 16 is exact when (s) is defined according to the proper electrostatic equations.

The reduction of the source of the reaction potential to a charge distribution limited to a closed surface greatly simplifies the electrostatic problem with respect to other formulations in which the whole dielectric medium is considered as the source of the reaction potential. Despite this remarkable simplification, the integration of eq 16 over a surface of complex shape is computationally challenging. The solutions are generally based on a discretization of the integral into a finite number of elements. This technique may be profitably linked to the boundary element method (BEM), a numerical technique widely used in physics and engineering to solve complex differential equations via numerical integration of integral equations (see the Website http:// www.boundary-element-method.com/ for a global view on literature and applications of this method).

The cavity surface is approximated in terms of a set of finite elements (called tesserae) small enough to consider (s) almost constant within each tessera. With (s) completely defined point-by-point, it is possible to define a set of point charges, q_k, in terms of the local value of (s) on each of these tesserae times the corresponding area A_k. The integral of eq 16 is thus transformed in the following finite sum:

Actually the local value of the potential necessary to define q_k also depends on the whole set of the surface charges, and so the correct values of the surface charges, and the correct expression of the reaction potential, are to be obtained through an iterative procedure. These aspects will be examined in section 2.3.1.5; now we present the most important ASC models. In this presentation we shall not make use of the BEM version of the ASC equations involving the point charges q_k but instead the original ones in terms of a continuous surface charge (s). A description of the formal and practical aspects related to the use of the BEM approach for ASC methods will be given in section 2.3.1.5.

2.3.1.1. Polarizable Continuum Model (PCM): Original Formulation. PCM, the oldest ASC method, at present is no more a single code, but rather a set of codes, all based on the same philosophy and sharing many features, some specialized for some specific purposes, others of general use, but with differences deserving mention.

The original PCM version was published in 1981, after some years of elaboration,⁸ and subsequently implemented in local and official versions of various QM computational packages.⁴⁸ More recently, PCM was renamed D-PCM (D stands for dielectric)⁴⁹ to distinguish it from the two successive reformulations (CPCM and IEFPCM) that we shall present in the next sections. This acronym is not completely correct as also the other reformulations refer to dielectric media (directly or indirectly); however, we cannot forget that it has become of common use in these years, and, thus, in the present review, we will adopt DPCM to refer to the first version of the model, whereas PCM will be used to refer to the entire family of models.

DPCM, like all members of the PCM family, is able to describe an unlimited number of solutes, each equipped with its own cavity and ASC, interacting among them through the dielectric. In this way, DPCM permits an extension of the basic model to association-dissociation phenomena, molecular clustering, etc., and it can account for a continuous shift from a single cavity to two cavities during a dissociation and the merging of two or more cavities during association. In parallel, it permits an extension to models in which the medium is composed by a set of nonoverlapping dielectric regions at different permittivity, constant within each region.

The reason for this versatility is in the use of the ASC approach in an unsophisticated version. To better appreciate this point, let us look again at the basic electrostatics from a different viewpoint.

For systems composed by regions at constant isotropic permittivity (including systems composed of a single isotropic solvent with multiple cavities), the polarization vector is given by the gradient of the total potential V(r) (including also that deriving from apparent charges)

At the boundary of two regions i and j, there is an ASC distribution given by

The basic PCM definition may be derived from the general expression 18 by taking into account the facts that in the basic case _i = 1 and _j = 1 and that we have computed the gradient on the internal (in) part of the surface, namely

After its first formulation, the DPCM was revised many times both in its theoretical aspects and in its numerical implementation; among all of these revisions, a fundamental one was proposed in 1995 by Cammi and Tomasi⁵⁰ when a new, and more efficient, computational strategy was defined to solve the BEM equivalent of eq 20 (see section 2.3.1.5 for more details and comments on this strategy).

Since its first presentation in 1981, DPCM has been "adopted" by many groups, which have thus largely contributed to its diffusion and its development. Some of these extensions will be examined in following sections of this review, but here we cannot forget to cite the work of three groups that have given, and continue to give, important contributions to show potentialities of DPCM.

First we cite the extensive and systematic work done in Barcelona by Luque, Orozco, and co-workers. The work of this group has to be here mentioned for several reasons. Their reformulation of DPCM, known by the acronym MST from the names of the three authors of the first PCM paper, has been applied to almost all aspects of solvation problems, with special attention to organic and biological systems; the amount of results of remarkable quality is the largest given by a single group. The papers of this vast literature devoted to methodological innovations (the main theme of this review) are abundant, all deserving attention. For obvious reasons, only a few^51-54 are cited here, and the reader is referred to a recent review³ for many others. Barcelona's group has spent efforts during the past 10 years to analyze the performances of the method in many ways and to validate it by comparisons with other procedures: simulations of different type, other continuum models. Among them we cite a comparison with two methods we shall examine in the following pages, the MPE method of Nancy (section 2.3.2) and the GB method of Minneapolis (section 2.3.3), performed in collaboration with the authors of such methods.⁵⁵ In the following, we shall cite other papers by Luque and Orozco, but what stimulated us here is to manifest the appreciation for the considerable work in developing solvation methods, through the example of numerous impeccable applications, the work of analysis and validation, the suggestion of methodological innovations improving or extending the potentialities of the approach, and the activity of popularizing it by means of several short (or long) review papers.

Another important contribution to the development of the DPCM has been done by Basilevsky and co-workers; in their numerous works on solvation phenomena they started from the original definition of DPCM to develop a parallel model known by the acronym BKO,⁵⁶ standing for Born-Kirkwood-Onsager. During the years, this model has been applied to many different problems. Here it is not possible to report the vast literature, but in the following sections we shall analyze in detail the main developments introduced by Basilevsky and co-workers, namely, the reformulation of continuum models to nonlocal dielectrics, their application to studies of solvation dynamics, and their extension to mixed continuum/discrete approaches. In these sections, we shall also report all correspondent references.

A last contribution we quote here is that of the research team of the Gaussian computational package.¹⁷ This team soon realized the potentialities and computational interest of the ASC approach to solvation, in the form given by DPCM. A version of the original model was rewritten by Gaussian and inserted in some releases around 1990, but the official collaboration between the PCM group and Gaussian started only at the end of 1996, and it has continued since then: more details on the results of this collaboration will be given in the following sections.

The Gaussian team has also independently elaborated a different PCM-like model based upon the previously mentioned isodensity surfaces (see section 2.2). This version of PCM is generally known as the isodensity-PCM (IPCM) model.¹⁶ The same model has been further extended to allow the isodensity surface at each SCF iteration to be varied; that is, the cavity is not fixed, even once the geometry is fixed (as is the case in the standard PCM model), but it is relaxed to the isodensity of the solvated molecule at each SCF iteration. This further development is generally known as self-consistent isodensity (SCI-PCM).⁵⁷ Both IPCM and SCIPM models are available in the Gaussian computational package.¹⁷

2.3.1.2. Conductor-like Screening Model (COSMO). In this method, originally devised by Klamt and Schüürmann,⁵⁸ the dielectric constant of the medium is changed from the specific finite value , characteristic of each solvent, to = . This value corresponds to that of a conductor, and this change strongly modifies the boundary conditions of the electrostatic problem. The most important effect is that the total potential V(r) of eq 9 cancels out on the cavity surface. From this condition it follows that the ASC is determined by the local value of the electrostatic potential instead of the normal component of its gradient (as in eq 20). To recover the effects of the finite value of the dielectric constant of the medium, the ideal unscreened charge density, *, corresponding to = , is finally scaled by a proper function of , namely

In the original paper about COSMO, Klamt suggested k = 0.5,⁵⁸ with the remark that k seems to depend on the cavity shape and on the distribution of charges in the solute. Other versions of COSMO proposed in the following years adopted other values for k. Truong and Stefanovic, in their GCOSMO,^59,60 and Cossi and Barone, in their C-PCM⁶¹ (this program is a part of the PCM suite of programs), proposed k = 0 on the basis of an analogy with the Gauss law. Subsequently, Cossi et al.⁶² compared the solvation free energies for neutral and charged solutes in water and in CCl₄, and they found that the choice of k is irrelevant for water but important in CCl₄: for neutral solutes the best agreement is obtained for k = 0.5, whereas for charged molecules it is preferable to use k = 0. Pye and Ziegler⁶³ in their implementation of COSMO in the ADF computational package^64,65 use k = 0 as default, but leave the users free to select the value they prefer. Chipman, in a recent comparison of several continuum models,⁶⁶ gives results for both k = 0 and 0.5, remarking that k = 0 seems to be somewhat better.

We note that the corrections introduced by this factor are quite small for solvents with high dielectric constant, and the final output of the energy is almost insensitive to changes in the k parameter. Things are a bit more critical for solvents of low and, in particular, for the cases in which one has to separately compute the inertial and the electronic component of the polarization vector (see sections 5.1 and 5.2); in these cases attention has to be paid in the use of solvation codes of the COSMO family.

An extension of the COSMO model, called COSMO-RS (conductor-like screening model for real solvents), has been proposed by Klamt,⁶⁷ starting from the analysis that solvents do not behave linearly when we consider strong electric fields on the molecular surfaces of fairly polar solutes. COSMO-RS is a theory that describes the interactions in a fluid as local contact interactions of molecular surfaces, the interaction energies being quantified by the values of the two screening charge densities that form a molecular contact. Having reduced all interactions to local interactions of pairs of molecular surface pieces, one can consider the ensemble of interacting molecules as an ensemble of independently interacting surface segments. COSMO-RS has become a predictive method for the thermodynamic properties of pure and mixed fluids. In contrast to group contribution methods, which depend on an extremely large number of experimental data, COSMO-RS calculates the thermodynamic data from molecular surface polarity distributions, which result from quantum chemical calculations of the individual compounds in the mixture. The different interactions of molecules in a liquid, that is, electrostatic interactions, hydrogen bonding, and dispersion, are represented as functions of surface polarities of the partners. Using an efficient thermodynamic solution for such pairwise surface interactions, COSMO-RS finally converts the molecular polarity information into standard thermodynamic data of fluids, that is, vapor pressures, activity coefficients, excess properties, etc.;^68-72 the corresponding computational code is called COSMOtherm.⁷³

2.3.1.3. Integral Equation Formalism (IEF). With the IEF method^74-76 originally formulated by Cancès and Mennucci in 1997, we come back again to methods involving solutes in a liquid phase and using the characterization of the solvent given by macroscopic properties of the specific liquid (such as the permittivity, ). Actually, IEF is a method more general than DPCM, or COSMO, and the case of an isotropic dielectric is nothing more than one of its possible applications, as we shall show below. Together with the previously mentioned DPCM and CPCM, IEF is a member of the PCM suite of solvation codes implemented in Gaussian since its 1998 version and, in the latest G03 version,¹⁷ it has become the default PCM formulation.

As in the two other ASC methods we have described in the previous sections, also for IEF we start from the electrostatic system (eq 6) and we introduce the decomposition (eq 9) of the potential V in the sum of the electrostatic potential generated by the charge distribution _M in vacuo, V_M, and the "reaction potential", V_R.

In the IEF, however, the potentials are redefined in terms of the proper Green functions; we recall that the Green function of an electrostatic problem G(x,y) is the potential produced in x by a unit point charge located in y (here with x and y we indicate two general positions in the space, but for simplicity's sake we drop the vector symbol).

If we indicate by G(x,y) = 1/x - y the Green function corresponding to the operator -² (namely, its kernel), by G^s(x,y) the Green kernel of the operator -(), and G^R(x,y) = G^s(x,y) - G(x,y), the following relationships arise:

Exploiting some basic results of the theory of integral equations (see ref 76 for details), it can be proved that the reaction potential, V_R, which satisfies

We note that eq 27 is different with respect to that reported in the original papers^74,76 because here we have used a different form of the Poisson eq 6 including a 4 factor on the right-hand side. This is reflected in the Green functions, G_a, and thus in the corresponding operators S_a and D_a, which now do not include the 1/4 factor.

The operators defined in eq 28 are well-known in the theory of integral equations: they are three of the four components of the Calderon projector.⁷⁷ We recall some of their properties: the operator S_i is self-adjoint, and is the adjoint of D_i for the scalar product. Besides, S_i = D_iS_i.

Equation 26 may be further simplified using the equality (2 - D_i)V_M + S_i(V_M/n) = 0; in this way the expression for g in eq 27 can be rewritten as⁷⁸

Numerically, it is advantageous not only because the calculation of a single scalar function (the potential) is computationally less demanding than the parallel calculation of both the potential and the vectorial electric field, but also because the potential is less sensitive than the field to numerical instabilities that can appear when we introduce a BEM approach to solve the electrostatic equations. From a formal point of view, the reformulation of the IEF ASC in terms of only the potential is important because it represents an implicit correction of the error due to the fraction of solute electric charge diffusing outside the cavity (namely, the previously defined outlying charge). The details on this issue will be given in section 2.4.3. Here it is worth reporting the final equations for the IEF ASC in the case of isotropic solvents; these, in fact, will be useful in the following when we shall compare IEF to other ASC methods. Namely, using the relations, S_e = S_i/, D_i = D_e, we can transform the IEF operators A of eq 27 and g of eq 29 in

From eqs 24-29 it can be seen that the IEF approach is completely general, in the sense that it can be applied without the need of modifying either the basic aspects of the model or the basic equations (26-29) to all of those systems for which the Green functions inside and outside the cavity are known. As a matter of fact, as the interior form of G is always known, G_i(x,y) = 1/4x - y, the problem is shifted to the evaluation of the exterior, G_e(x,y). Analytical expressions of this function are available for standard isotropic solvents characterized by a constant and scalar permittivity [namely G_e(x,y) = 1/(x - y)], but also for anisotropic environments characterized by a tensorial but constant permittivity for ionic solutions described in terms of the linearized Poisson-Boltzmann equation (see section 2.3.1.6 for further details) and for sharp planar liquid-metal interfaces.⁷⁹ Obviously, it is not always possible to have analytical Green functions. However, in several cases the Green function can be effectively built numerically, and thus the IEF approach can be generalized to many other environments as, for example, a diffuse interface with an electric permittivity depending on the position.⁸⁰

As a final comment it is worth noting that IEF contains as subcases both the DPCM and the COSMO models. In the first case we have to consider an isotropic solvent with a scalar permittivity and use the electrostatic equation already introduced above to rewrite the IEF ASC eq 29 in terms of the potential, in the opposite way so as to keep (V^M/n) instead. To recover COSMO ASC, on the other hand, we have to consider that , and thus all of the terms in eqs 27 and 29 involving the S_e operator can be neglected (S_e is in fact proportional to 1/).

2.3.1.4. Surface and Volume Polarization for Electrostatic [SVPE and SS(V)PE]. From 1997 to date, Chipman and co-workers have developed a series of continuum models^81-85 that have, as a starting point, the consideration that unconstrained quantum mechanical calculation of solute charge density generally produces a tail that penetrates outside the cavity into the solvent region (the "outlying charge" mentioned in the previous sections). Exact solution of Poisson's equation in this situation would require invocation of an apparent volume polarization charge density lying outside the cavity in addition to the apparent surface polarization charge density lying on the cavity.⁸⁶ As a consequence, also the reaction potential heretofore written in terms of an apparent surface charge has to be supplemented with a term due to the apparent volume charge, namely, following the notation used by Chipman

The acronym used to indicate this formulation was SVPE, meaning that both surface and volume polarization for electrostatic interactions were included. Contrary to the previous ASC models, this method requires a discretization of both the cavity surface and the exterior volume; in particular, the volume polarization charge density was approximated by a collection of point charges located at various nodes on a series of layers covering the exterior region, as described in detail in ref 82.

The exact SVPE method is laborious to implement and time-consuming, because it utilizes a volume polarization potential arising from a discontinuous volume charge density. With a cavity surface that is adapted to the detailed nonspherical shape of a general molecular solute, this leads to difficult integrations over just part of the full three-dimensional space.

To avoid this large complexity a simpler approximate solution that involves only apparent surface charge distributions was subsequently introduced. Here, we limit ourselves to summarize the main aspects of this method denoted "surface and simulation of volume polarization for electrostatics", SS(V)PE.

The SS(V)PE method originated with the demonstration⁸¹ that all direct and indirect effects of the volume polarization charge density on the exact SVPE reaction potential can be exactly represented at all points inside the cavity (which is the most important region) by simulating the explicit volume charge density in terms of an additional surface charge density so that the total apparent surface charge satisfies the equation⁸⁴

The correspondence between SS(V)PE and IEF is not only of notation: the two methods in fact coincide when the isotropic IEF(V) eq 31 is considered; to show this, it is sufficient to multiply both sides of eq 33 by 2( + 1)/( - 1) and to use the relationship S_i = D_iS_i. The formal equivalence of the two methods was published⁸⁷ only one year after the publication of the SS(V)PE original paper,⁸⁴ and thus the two methods are often considered as two alternative ASC methods. As a matter of fact, when the whole solute charge is contained within the cavity, SS(V)PE is equivalent also to DPCM (besides IEF), but when a part of this charge lies outside this, equivalence with DPCM is no longer valid; on the contrary, the equivalence with IEF(V) (see the end of section 2.3.1.3) remains in all cases. The late recognition of the equivalence of the two approaches is also reflected in a further PCM approach proposed by Cossi et al.⁸⁸ In fact, this model, indicated as implicit volume charges PCM (or IVCPCM), was presented by the authors as the PCM reformulation of the SS(V)PE by Chipman but, in practice, it exactly coincides with the IEF version of PCM. Unfortunately, also in papers that appeared later, the authors never clarified this point, and thus the two methods continued to be described as different approaches.⁸⁹

2.3.1.5. Discretization of the Cavity Surface and the ASC-BEM Equations. In the presentation of the ASC models we have adopted the original formalism using a continuous surface charge . Here, we go a step farther and present a description of the formal and practical aspects related to the use of the BEM version of the ASC equations involving the point charges q_k. It is interesting to note that BEM is a relatively young technique (the acronym was first used in 1977⁹⁰), but now a very large literature, including books and dedicated international conferences (the principal conference on BEM conducts its 27th meeting in 2005⁹¹), can be found on this very powerful mathematical approach.

As introduced at the beginning of section 2.3.1, the definition of ASC elements uses the discretization of the surface into a finite number of elements (usually called tesserae), which is the cornerstone of BEM, and for this reason it may be considered an application of a BEM procedure (or, better, of the indirect BEM procedure; see also section 2.3.1.6). Actually it is not a trivial application, because the original BEM was limited to cases without internal sources, that is, based on the Laplace equation ²V(r) = 0 instead of the Poisson equation (7). Other features of BEM have been considered in more recent versions of ASC methods, but it must be said that the large wealth of suggestions available in the BEM literature has not yet been fully exploited in the continuum solvation methods.

The necessary preliminary step in BEM-ASC strategy is the generation of the surface finite elements. This step can be achieved in many different ways; here, in particular, we shall focus on the algorithm developed for the GEPOL cavity (see section 2.2).

In the original formulation, the GEPOL tesserae were generated by inscribing a polyhedron with 60 faces within each sphere, projecting the faces onto the spherical surface, and discarding those that were fully inside an intersecting sphere. If a polyhedron face happened to be partially inside a neighboring sphere, the tessera was cut into smaller triangles, and only those triangles that were completely exposed were kept; in this way the final tessera was represented by a complex polygon.

GEPOL has been considerably improved during the years; we mention here an improved scheme to fill the solvent-excluded space⁹² and the algorithm through which the exposed portion of a tessera partially inside other spheres is analytically cut by adding edges to the tessera being generated.⁹³ In this way, all of the tesserae are defined in terms of connected arc segments, and their surface areas can be analytically computed using the Gauss-Bonnet formula so that the analytical derivatives of all the tesserae can be computed. Another important improvement is the introduction of flexible tessellation⁹⁴ and symmetry-adapted tessellations.^95,96 This advanced version of GEPOL is available in GAMESS^97,98 and Gaussian¹⁷ computational codes. The last revision of Gaussian (Gaussian 03) includes an even more advanced version of GEPOL; this recent re-implementation has been realized by Scalmani et al.⁹⁹ and, among different important numerical improvements, it presents also a linear scaling computational cost, thus extending the range of applicability to very large molecular systems.

Recently, two independent novel schemes have been proposed^100,101 to sample the GEPOL molecular surface without an explicit discretization of the cavity surface, which may lead to unphysical discontinuities. In both of these schemes the sampling of the surface is done on each original sphere, and it is not changed after all of the spheres have been combined to construct the final cavity. Each sampling point is then associated with a "weight" (which in the simple scheme coincides with the area of the tessera) and to a switching function, the definitions of which are different in the two methods. In both cases, however, the main result is that no sampling points and weights will be discarded; their number remains constant also during a geometry optimization when changes in the solute geometry induce changes in the relative position of the atom-centered spheres.

Once the cavity surface has been partitioned in tesserae small enough to consider (s) almost constant within each tessera, it is possible to define a set of point charges q_k in terms of the local value of (s) on each of these tesserae times the corresponding area. In doing this the electrostatic equation presented in sections 2.3.1.1-2.3.1.4 for (s) within the various ASC methods can be rewritten as a set of T (with T equal to the number of tesserae) coupled equations, which can be recast in a matrix form of the type

The elaboration of the problem giving origin to an equation of type 34 has been done in a large number of ways, in practice one for each ASC method, and also in several different ways for the same method. There is no need to repeat here the various elaborations; instead, we provide a table with the expressions of the elements of the K and f matrices for the most important variants of ASC methods, namely, DPCM,CPCM, IEFPCM, and SS(V)PE. As far as concerns IEFPCM, two different sets of expressions are given, the first referring to standard isotropic solvents (characterized by a constant scalar permittivity ) and the second for anisotropic solvents (i.e., characterized by a constant but tensorial permittivity ) and for ionic solutions (in the limit of a linearized PB scheme, see eq 43). We also note that in all cases, IEFPCM equations have been rewritten in the IEF(V) form, that is, using only f = V and not the combination of V and E_n as in the first formulation of the method; this reformulation was originally proposed for the isotropic case only, but recently it has been generalized to anisotropic dielectrics;¹⁰² here, we present for the first time the parallel version for the linearized Poisson-Boltzmann scheme (see also section 2.3.1.6).

In Table 1 we have reported only the off-diagonal elements of the D and S matrices involved in the different versions, because different numerical solutions have been proposed for the diagonal elements. In particular, for S the following approximation is generally used

An alternative expression for D_ii has been adopted by Foresman et al.,¹⁶ by Chipman,^83,104,105 and by Cossi et al.⁸⁸ following the original proposal by Purisima;^106,107 this expression uses the sum of the off-diagonal elements, namely

The linear system (eq 34) can be solved either by matrix inversion or iteratively.⁵⁰ The cost of computing the ASC terms used to be usually negligible if compared to other steps in an ab initio calculation (e.g., Fock matrix building and diagonalization). Nowadays, this is no longer true as, in the field of ab initio methods, fast algorithms for which the cost grows linearly with the system size [e.g., the fast multipole method (FMM)] are beginning to be widely available and applied.^108-111 On the other hand, the significant research effort being devoted to hybrid quantum/classical approaches (see section 7) is driving toward the handling of large (partly) classical solutes, which could easily involve thousands of atoms.

Given these methodological and algorithmic advances, the cost of calculating the ASCs could easily become the computational bottleneck. Indeed, whenever fast and linear scaling algorithms, such as the FMM, are used to compute the electrostatic potential/field on the tesserae, the cost of evaluating the interaction of the polarization charges among themselves will grow more steeply and will rapidly become dominant as the number of tesserae increases.

Recently, an efficient linearly scaling procedure has been formulated¹¹² (and implemented in the latest revision of Gaussian) by exploiting the FMM for computing electrostatic interactions. Such a procedure has been formulated for DPCM, CPCM, and IEFPCM models, and it is based on an iterative solution of the linear system (eq 34). Here, we present the basic equations just for the DPCM case as more intuitive, but we refer the interested reader to the reference paper for the corresponding equations for CPCM and IEFPCM. By properly rewriting the system (eq 34), the iterative equation defining the apparent charges becomes

We conclude this section by mentioning an alternative approach to define the apparent charges formulated by York and Karplus¹¹³ within the COSMO formalism. The method uses smooth Gaussian basis functions of the form

For the appearance and disappearance of surface elements to occur smoothly as a function of geometrical changes, a switching layer around each atom is introduced. The switching layer serves to "turn off" or "turn on" the surface elements associated with other atoms as they pass into or out of the layer.

Very recently, this approach has been extended to the PCM suite of models by Scalmani and Frisch¹¹⁴ within the Gaussian code. In this new formalism the S and D matrices reported in Table 1 can be expressed in terms of two-center, two-electron integrals and their derivatives, for example

2.3.1.6. More on the BEM. When one wants to extend the basic model to solutions with nonzero ionic strength (i.e., to salt solutions), additional problems appear. One way to treat these problems consists of adopting other formulations of the BEM approach such as those reported here below. These methods can be applied also to zero ionic strength, and for this reason they are considered here, in a section of the review dedicated to isotropic liquids without dispersed ions. The basic electrostatic equations modify as

Equation 43, generally known as the linearized Poisson-Boltzmann (LPB) equation, is the result of the Debye-Huckel approximation applicable in the case of low potentials, a condition approached at low concentrations. The LPB equation is obtained by taking the linear term of the Taylor expansion of the hyperbolic sine function which, according to a treatment of statistical physics in the thermodynamic equilibrium approximation, represents the distribution of the mobile ions in the field of the electrostatic potential V.

The LPB equation can be solved analytically for simple cavity shapes such as a sphere.¹¹⁵ However, for a general cavity of more complicated shape (such as that adapted to the molecular solute), it must be solved numerically. Many different models have been developed so far. In this section, we are interested in the specific class of approaches, which determines the total electrostatic potential indirectly by using BEM techniques. To better analyze this class of approaches it is useful to introduce some further information about BEM that has not been presented yet.

To derive the BEM, one must replace the partial differential equation that governs the solution in a domain by an equation that governs the solution on the boundary alone. There are two fundamental approaches to the derivation of an integral equation formulation of a partial differential equation. The first is often termed the direct method, and the integral equations are derived through the application of Green's second theorem. The other method is called the indirect method and is based on the assumption that the solution can be expressed in terms of a source density function defined on the boundary.

Among the indirect methods, we can, for example, mention the IEF and the parallel SS(V)PE approaches described in the previous sections. Both of these two methods have been generalized to treat the LPB problem,^75,116,117 and, as for the basic isotropic solvent, only a single-layer (charge) distribution (s) is exploited. We note here that, for problems involving ionic effects, the use of indirect approaches is new and thus not largely used until now (to the best of our knowledge the original implementation of IEF in 1997 is the first example). A far more used BEM method is the direct one. To explain this method let us consider the electrostatic equations (42 and 43) governing the domain C bounded by the surface . By applying Green's second theorem, these equations can be replaced by integral equations; for example, eq 43 becomes

The surface integral on the right-hand side of eq 44 is called the single-layer potential, whereas the one on the left is called the double-layer potential. For this reason, the direct BEM approaches are also indicated as models employing both single-layer and double-layer distributions over the surface cavity.

The power of this formulation lies in the fact that it relates the potential V and its derivative on the boundary alone; no reference is made to V at points in the domain. To numerically solve the integral equations of the type of eq 44, the first step is to partition the surface into surface elements (i.e., to define an N-point grid on the surface).^118-124 Triangular and quadrilateral elements are commonly used for this purpose. As a second step, the unknown functions V(r) and V(r)/n are approximated by continuous trial functions over each boundary element. Generally, polynomials determined by their values at the N grid points (the nodes) are used. In this way, the integrals become sums of integrals over elements, which results in a set of linear equations for the polynomial coefficients. This set of equations can then be written as a single matrix equation of dimension 2N.^125,126 A fast multipole algorithm can be used to reduce computational costs,^107,127,128 or one can apply an iterative procedure on a pair of coupled N × N finite matrix linear equations in every iteration.¹²⁹

The seminal papers of Kirkwood^115,130 and Onsager¹³¹ have provided inspiration for various continuum solvation models based on a multipole expansion (MPE) of the solute charge distribution. Here, we present these models starting from the simplest one (a dipole inside a sphere) and passing to those using multipole expansions and/or cavities of more complex form.

Appreciable popularity was gained in past years by the Onsager-SCRF code, elaborated by Wiberg and co-workers^132,133 for the Gaussian computational code.¹⁰³ The Onsager model is the simplest version of the MPE approach. Solvation is described in terms of a dipole moment, drawn in an iterative way from QM calculations on the molecule. The appealing feature of Onsager-SCRF was that it permitted one to directly exploit almost all of the computational facilities of the Gaussian packages. For this reason, and for its very limited computational cost, it is still in use by people not requiring an accurate description of solvation effects but just a guess or a qualitative correction to the values obtained for the isolated molecule. Users must be aware of the limitations of the approach, of the unphysical deformation of the solute charge distribution it may induce, and of other shortcomings specific of the approach, such as the lack of solvation for solutes with zero permanent dipole.

A more general model going beyond the dipole approximation is that developed by Mikkelsen and co-workers.¹³⁴ This model starts from the idea proposed by Kirkwood to describe the interaction between a set of classical charges (described in terms of a multipole expansion) enclosed in a spherical cavity embedded in a structureless polarizable dielectric medium described by the macroscopic dielectric constant . Mikkelsen applies the model to QM charge distributions instead of classical point charges, as in the original Kirkwood paper, and the effects of the solvent polarization are described in terms of proper QM operators to be added to the Hamiltonian of the isolated system (i.e., adopting the same strategy introduced in section 2.3 to define the "effective Hamiltonians").

Even if the model is limited to spherical cavities, it presents very advanced features especially as far as concerns the QM aspects. The model has been implemented in the Dalton quantum chemistry program,¹³⁵ and it has been used to study solvent effects on a large number of molecular response properties; it has also been generalized to several QM approaches including multiconfigurational self-consistent field and coupled cluster methods^136-143 (see also section 2.4.6 and section 6).

The most complete, and also the most largely used MPE method, is that developed in Nancy by Rivail and co-workers.^6,144 Historically, this method represents the first example of QM continuum solvation methods, but over the years it has been continuously developed, and it continues its evolution.¹⁴⁵ In some cases the acronym SCRF (self-consistent reaction field) is also used to indicate this specific method, even if it also identifies the larger class of solvation methods implying the use of a reaction potential to be self-consistently solved together with the solute charge.

Here, we shall consider three steps in the evolution of this MPE method: (i) a one-center expansion for a spherical cavity (in this case it coincides with that developed by Mikkelsen); (ii) its extension to ellipsoid cavities; and, finally, (iii) its extension to multicenter expansion and to cavities of general shape.

(i) The expression of the multipole expansion for the potential of a charge distribution placed into a sphere of radius a, immersed into a continuum with dielectric constant , has been known for a long time and can be found in numerous textbooks. We report the expression given by Böttcher¹⁴⁶ for the potential inside the cavity, namely

The QM reformulation of the problem given by Rivail consists of using quantum mechanical approaches to compute the elements of the multipole expansion of _M and in introducing the corresponding V_r terms in the Hamiltonian to get a new solvent-modified _M, until convergence. This process, as in the ASC methods, is of iterative nature.

(ii) Rivail and co-workers rapidly abandoned the spherical cavity of the first formulation of the model to pass to ellipsoids. The shape of the cavity is a critical factor in continuum methods, and surely an ellipsoidal cavity, with its variable parameters, is more suited than a sphere to describe solvent effects on molecules. To do this extension, Rivail's group was compelled to develop a part of the mathematical theory on multipole expansions (in terms of ellipsoidal harmonics), which was almost completely neglected before. We shall use now the simplified notation introduced by Rivail and valid for both ellipsoids and spheres.

The solute-solvent interaction energy can be written in the form

(iii) The reaction field components, expanded over a set of centers I, J, ... have the following expression

In the latest version of the Nancy MPE method,¹⁴⁵ the ellipsoidal cavity has been abandoned, shifting to a molecular-shaped cavity such as that used in PCM. For this cavity, the number of points selected to fix the reaction factors is relatively large (usually of the order of 1500 or more), but this large amount of numerical values (quite simply individually computed, however) is efficiently compacted into linear equations, put into a matrix form, and solved with standard methods. The large number of points ensures rotational stability with respect to the selection of sampling points, a problem that disturbs other similar fitting procedures, such as the definition of the ESP atomic charges, that is, the charges fitted on the electrostatic potential, and some early versions of PCM.

Another point of general methodological interest in this latest version of the Nancy MPE code involves the selection of the multipole multicenter distribution. Because the partition is arbitrary, there is an infinite number of ways of defining the distribution, but they are not equivalent from a computational viewpoint. The literature on intermolecular interaction is rich with proposals and debates about the choice of the multipole distribution. Rivail and co-workers have selected and tested three different distributions. In the first, called GAD (Gaussian averaged distribution), each multipole is distributed over the N centers (the number of centers corresponds to the number of nuclei) with a Gaussian weight, with the translation expressed in terms of solid harmonics. In the second, called MSP (Mulliken-Sokalski-Poirier), use is made, for the electronic component, of the one-electron integrals of the pertinent _v elementary charge distribution for the multipole moment operator, namely

In the last expansion, called LSP, use is made of the Lowdin orthogonalization procedure applied to expression 51.

All three expansions are convergent, but GAD and MSP have a better convergence, MSP being the less time-consuming. The solvation code is accompanied by analytical derivatives with respect to the nuclear coordinates,^145,147 a feature essential for the performance of geometry optimizations in solution (see section 6.1.1).

We conclude this section by noting that this recent re-elaboration of the Nancy MPE code shows that previous criticisms about the unreliability of MPE codes for large molecules are not justified. Indeed, results are obtained that are of comparable quality to those of ASC codes and still with comparable computational times.

The name comes from an affiliation with the old Born model, which gives the interaction free energy of a spherical ion inside a dielectric¹⁴⁸

The passage from a single sphere to a multiplicity of interacting spheres, each with a charge at its center, required a considerable amount of ingenuity. Earlier versions of the GB approach were described in the 1994 review, and we start here from the definition given by Still and co-workers¹⁴⁹ of an empirical Coulomb operator 1/f_GB that has to be applied to all pairs q_i and q_j of atomic solute charges, namely, by generalizing the free energy (eq 52) to

We note that when eq 53 is applied to a two-dielectric system in which the set of charges occupies a molecular cavity with a dielectric constant _p surrounded by a uniform high-dielectric continuum environment with a dielectric constant _w, the prefactor becomes (1/_p-1/_w). An internal dielectric of _p = 1 is appropriate for simulations when dipole fluctuations occur explicitly as part of the model. A larger internal dielectric constant is more appropriate when the energies of minimized or averaged structures are evaluated.

A fundamental problem arising in the GB approach, and present in all of its successive formulations, is the definition of the appropriate Born radii _i.^151-155

For every charged atom i, the effective Born radius _i is related to the effective Born free energy of solvation, G_i, of a reference system through the Born formula (eq 52), where the reference system is, by definition, the same as the original solute/continuum solvent system, except that in the reference system atom i bears unit charge and all other atoms are neutral. Exact calculation of _i for every charged atom in a macromolecule through numerical solution of the Poisson-Boltzmann equation is time-consuming and of little practical use. The potential of GB models has come from the progress in the formulation of approximate methods to compute _i with sufficient accuracy and efficiency. An analytical formula proposed by Qiu et al.¹⁵⁶ forms the basis of several recent parametrizations and applications of the GB model for protein systems.^150,157-160 In Qiu's work, the parameters contained in the formula have been optimized to minimize the differences between the effective Born radii calculated by the finite difference Poisson-Boltzmann (FDPB) method (see section 2.3.5) and by formula 53. Besides the approach of Qiu et al., Hawkins et al. have proposed another approximate pairwise method to compute the effective Born radius analytically.^161,162

There are many different GB models, and still they are the subject of active research, especially in the simulation of macromolecules^150,163-165 such as nucleic acids and proteins. Here, we obviously cannot mention all of them and analyze their specificities. Moreover, most of these models are formulated within classical descriptions and thus are beyond the scope of the present section (see section 7.1 instead).

GB methods have been also successfully applied to QM descriptions of the solute.

The best known and more widely used QM GB method is that developed by Cramer and Truhlar in Minneapolis.¹⁶⁶ This method has been elaborated and distributed in a large variety of versions, with different names. At present, these models are collected under the collective name SMx, where x stands for an alphanumeric code indicating the version, with its specific features (the most used recent codes are SM5.42R and SM5.43R), but in the past other names were used such as AMSOL. In one of the most recent papers¹⁶⁷ 19 versions of the method are reported, but other versions are not present in this list.

The first versions known under the name AMSOL were all developed for various semiempirical methods (AMSOL still refers to a semiempirical code¹⁶⁸), whereas the SMx suite of programs permits ab initio QM calculations at several levels of the theory.^169-172 Some SMx programs are available in the GAMESSPLUS,¹⁷³ which is an add-on module to the GAMESS,^97,98 and in HONDOPLUS^174,175 program. The latest version of the method is called SMxGAUSS,¹⁷⁶ and it may use Gaussian output files or be connected with Gaussian 03¹⁷ to exploit the many features this last program contains. The objective of making the SMx method more accessible to users is not the only, nor the main, reason for so large a number of versions. The real reasons will be shown later in this subsection.

We have not yet examined how the solvation electrostatic problem is linked to the QM one in continuum methods. This will be done in detail in section 2.4, but it is convenient to anticipate here some aspects to better present the SMx models.

The Fock operator of the solute in vacuo (for simplicity we limit ourselves to the Hartree-Fock case, even if SMx codes can treat higher levels of the quantum molecular theory) is modified by a solute-solvent interaction potential describing the electrostatic interaction according to the GB model, that is, in terms of atomic charges q_i and of the modified Coulomb operator (see eq 53). The q_i charges are derived from the solute electronic wave function via a population analysis performed in different ways in the various versions of the method. First, Mulliken populations were used, then Löwdin charges, and finally charges that according to the Cramer and Truhlar definitions are of classes III and IV. The class IV charges¹⁷⁷ are derived from Mulliken NDDO charges with a mapping to better describe physical observables. Energy minimization cycles in the so-defined Fock equation lead to self-consistent modification of the set of q_i charges and to the introduction of solvent polarization in the solute (see also section 2.4.5 for further comments). The free energy of the electronic polarization, called G_P, can be so computed.

Another term of the free energy difference with respect to the isolated molecule is also computed. This term, called E_EN, corresponds to the change in the electronic and nuclear energy of the solute upon relaxation from the gas-phase geometry. These two terms are combined to give the electrostatic contribution to the solvation free energy:

In all stages of the evolution of the Cramer-Truhlar solvation code a complete parametrization, or reparametrization, has been performed on very large data sets, generally of some thousands of solute-solvent couples, with a coupled fitting for G_CDS and G_ENP. These parametrizations are quite specific in the sense that for each solvent they are specific for the type of QM calculation (AM1, PM3, HF, DFT, according to the functional, etc.), for the basis set in ab initio methods, for the definition of the q_i model (Mulliken, Löwdin, class III, class IV), for the parameters inside the coulomb operator 1/f_GB = g_kk' for the Born radius and for the definition of the solvent accessible surface (an alternative definition of SAS has been also given).¹⁷⁸ In addition, reparametrization has been also performed for specific classes of phenomena (solvent partition equilibria, solvatochromic effects, etc.). This remarkable attention to the fitting procedures explains the large number of different versions, each about a specific combination of the above-mentioned parameters of the computation. We note that the accurate documentation available for the users is a fundamental help for selecting the version of the programs best suited for a specific study.

2.3.4. Finite Element (FE) and Finite Difference (FD) Methods