Quantum Mechanical Continuum Solvation Models
Jacopo Tomasi,*,† Benedetta Mennucci,† and Roberto Cammi‡
Dipartimento di Chimica e Chimica Industriale, Universita` di Pisa, Via Risorgimento 35, 56126 Pisa, Italy, and Dipartimento di Chimica, Universita` di
Parma, Viale delle Scienze 17/A, 43100 Parma, Italy
Received January 6, 2005
Contents
1. Introduction 3000
1.1. Generalities about This Review 3000
1.2. Generalities about Continuum, Focused, and
Layered Models
3001
2. Methodological Outlines of the Basic QM
Continuum Model
3003
2.1. Definition of the Basic Model 3003
2.2. Cavity 3003
2.3. Solution of the Electrostatic Problem 3005
2.3.1. ASC Methods 3006
2.3.2. MPE Methods 3015
2.3.3. Generalized Born (GB) Approaches 3016
2.3.4. Finite Element (FE) and Finite Difference
(FD) Methods
3018
2.4. Solution of the Quantum Mechanical Problem 3019
2.4.1. Intuitive Formulation of the Problem 3019
2.4.2. Electrostatic Operators 3020
2.4.3. Outlying Charge 3021
2.4.4. Definition of the Basic Energetic Quantity 3023
2.4.5. QM Descriptions beyond the HF
Approximation
3024
3. Some Steps beyond the Basic Model 3025
3.1. Different Contributions to the Solvation
Potential
3025
3.2. Use of Interactions in Continuum Solvation
Approaches
3027
3.2.1. Cavity Formation Energy 3028
3.2.2. Repulsion Energy 3033
3.2.3. Dispersion Energy 3035
3.2.4. Charge Transfer Term 3036
3.2.5. Definition of the Cavities in the
Calculation of Solvation Energy
3037
3.2.6. Contributions to the Solvation Free
Energy Due to Thermal Motions of the
Solute
3038
4. Nonuniformities in the Continuum Medium 3039
4.1. Dielectric Theory Including Nonlinear Effects 3039
4.2. Nonlocal Electrostatic Theories 3040
4.3. Nonuniformities around Small Ions 3041
4.3.1. (r) Models 3041
4.3.2. Layered Models 3041
4.3.3. Molecular Cluster Models 3042
4.4. Nonuniformities around Neutral Molecules 3043
4.5. Nonuniformities around Systems of Larger
Size
3043
4.6. Systems with Phase Separation 3044
5. Nonequilibrium in Time-Dependent Solvation 3046
5.1. Dynamic Polarization Response 3047
5.2. Vertical Electronic Transitions 3047
5.3. Solvation Dynamics 3049
5.4. Spectral Line Broadening and Solvent
Fluctuations
3053
5.5. Excitation Energy Transfers 3054
5.6. Time-Dependent QM Problem for Continuum
Solvation Models
3056
6. Molecular Properties of Solvated Systems 3058
6.1. Energy Properties 3059
6.1.1. Geometrical Derivatives 3059
6.1.2. IR and Raman Intensities 3061
6.1.3. Surface-Enhanced IR and Raman 3062
6.2. Response Properties to Electric Fields 3063
6.2.1. QM Calculation of Polarizabilities of
Solvated Molecules
3064
6.2.2. Definition of Effective Properties 3064
6.3. Response Properties to Magnetic Fields 3066
6.3.1. Nuclear Shielding 3066
6.3.2. Indirect Spin−Spin Coupling 3067
6.3.3. EPR Parameters 3068
6.4. Properties of Chiral Systems 3069
6.4.1. Electronic Circular Dichroism (ECD) 3069
6.4.2. Optical Rotation (OR) 3069
6.4.3. VCD and VROA 3070
7. Continuum and Discrete Models 3071
7.1. Continuum Methods within MD and MC
Simulations
3072
7.2. Continuum Methods within ab Initio Molecular
Dynamics
3074
7.3. Mixed Continuum/Discrete Descriptions 3075
7.3.1. Solvated Supermolecule 3076
7.3.2. QM/MM/Continuum: ASC Version 3076
7.3.3. ONIOM/Continuum 3077
7.3.4. (Direct) Reaction Field Model 3078
7.3.5. Langevin Dipole 3078
7.4. Other Methods 3079
7.4.1. ASEP-MD 3079
7.4.2. RISM-SCF 3080
8. Concluding Remarks 3081
* Author to whom correspondence should be addressed (e-mail
tomasi@dcci.unipi.it).
† Universita` di Pisa.
‡ Universita` di Parma.
2999Chem. Rev. 2005, 105, 2999−3093
10.1021/cr9904009 CCC: $53.50 © 2005 American Chemical Society
Published on Web 07/26/2005

9. Acknowledgment 3084
10. References 3084
1. Introduction
1.1. Generalities about This Review
This review on continuum solvation models has
been preceded in Chemical Reviews by others ad-
dressing the same subject. They are due to Tomasi
and Persico1 (published in 1994), Cramer and Tru-
hlar2 (published in 1999), and Luque and Orozco3
(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 acknowl-
edged by Chemical Reviews since its very beginning.
In the abundant number of reviews addressing dif-
ferent 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 computa-
tional 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 de-
scription 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 indi-
rectly, 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
Jacopo Tomasi received his “Laurea” degree in chemistry in 1958
(University of Pisa) discussing a thesis on the theoretical determination
of the intensities of vibrational overtones. Since 1980 he has been a full
professor of physical chemistry at the University of Pisa. His research
interests cover several aspects of theoretical chemistry with propensity
to the formulation and elaboration of models based on ab initio quantum
chemistry with a special emphasis on the exploitation of the interpretation
of the phenomenon to obtain computational codes of easy use. This
approach has been applied to molecular interactions, chemical reaction
mechanisms, photochemical processes, solvent effects on molecular
response properties, and other related subjects. He has authored, or
coauthored, more than 300 scientific papers.
Benedetta Mennucci was born in Lucca, Italy, in 1969. She received her
“Laurea” degree in chemistry in 1994 and her Ph.D. degree in chemistry
from the University of Pisa in 1998 discussing a thesis on theoretical
models and computational applications of molecular phenomena involving
the environment effect. In the same year she became assistant professor
in the Department of Chemistry of the University of Pisa. Since 2002 she
has been associate professor of physical chemistry at the same institution.
Her research interests focus on the elaboration of theoretical models and
computational algorithms to describe molecular systems in condensed
phase with particular attention to molecular properties and time-dependent
phenomena. She has authored, or coauthored, more than 80 publications.
Roberto Cammi was born in Busseto, Italy, in 1954. In 1979, he was
awarded the degree of “Dottore in Chimica” at the University of Parma,
discussing a thesis in theoretical chemistry. Since 1983 he has been a
researcher with the Institute of Physical Chemistry of the University f
Parma. In 2000 he became an associate professor in the Department of
General and Inorganic Chemistry, Analytical Chemistry, Physical Chemistry
of the University of Parma, and since 2002 he has been a full professor
of theoretical chemistry of the University of Parma. He teaches physical
chemistry and theoretical chemistry. His research field is theoretical and
computational chemistry, mainly the developments and applications of
quantum mechanical continuum methods to the study of solvent effects
on molecular processes and properties. He has published more than 80
research papers and 8 chapters of collective books.
3000 Chemical Reviews, 2005, Vol. 105, No. 8 Tomasi et al.

is of direct interest here: the merging of simulations
with quantum mechanical (QM) descriptions of mo-
lecular structures. This merging, in progress for
several years, has to overcome computational dif-
ficulties, 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 computa-
tional 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 electro-
static potential generated by the charge distribution
of a molecule could represent a valid analytic and
interpretative tool to study intermolecular interac-
tions.4
The perspectives opened by these findings were
elaborated independently by the various laboratories
and led to the alternative theoretical methods5-8 that
have been examined in detail in the preceding
reviews. In the present one we abandon the historical
perspective used in the first review1 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 second2 and third3 reviews give consider-
able 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 method-
ological issues, with a perspective focused on the most
recent developments. The body of the theory devel-
oped 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, het-
erogeneous 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.
1.2. Generalities about Continuum, Focused, and
Layered Models
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 continu-
ous way, usually by means of a distribution function.
Continuum distributions are a very general con-
cept. In the standard quantum mechanical descrip-
tion 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 func-
tion is easily derived with an averaging operation,
the one-electron density function FM
e (r), which con-
tains a good deal of the information conveyed by the
original wave function. According to the formal
theory, one-electron FM
e (r) and two-electron FM
e (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
Quantum Mechanical Continuum Solvation Models Chemical Reviews, 2005, Vol. 105, No. 8 3001

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 in-
termediate between the two extremes we have men-
tioned, 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 Hamilto-
nian 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 Hö 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 solva-
tion models Hö R(r) is maintained and the focusing
simplifications involve the number of freedom de-
grees 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(rb,rb′), where
(rb) 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 mol-
ecule. The density here is that of the liquid system
R, and Q(rb,rb′) is, in the more complete formulations,
expressed as a sum of separate terms each related
to a different component of the solute-solvent inter-
action.
We shall be more specific in section 3, but for this
general discussion we limit ourselves to the electro-
static 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 PB 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 rb, or by a couple of
position vectors (rb,rb′).
The quest for the accuracy requested in chemical
applications of such focused models suggests the
introduction of additional features, different accord-
ing 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.
Hö FR(f,r) ) Hö F(f) + Hö R(r) + Hö int(f,r) (1)
Hö eff
FR(f,r) ) Hö F(f) + Hö int(f,r) (2)
Hö eff
FR(f,r) ) Hö F(f) + Vö int[f, Q(rb,rb′)] (3)
PB )  - 1
4ð
EB (4)
3002 Chemical Reviews, 2005, Vol. 105, No. 8 Tomasi et al.

2. Methodological Outlines of the Basic QM
Continuum Model
This topic has been amply presented in the first
review on the argument,1 complemented in a second
review,2 and summarized in many other places,
including textbooks.9-12 For this reason we shall
avoid an exhaustive presentation of the whole litera-
ture 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 char-
acteristics 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.
2.1. Definition of the Basic Model
For the basic QM model we intend models with the
following characteristics:
(1) The solute is described at a homogeneous QM
level. Computational procedures based on semiclas-
sical or classical descriptions of the solute will be
considered as derivation of the basic QM model and
briefly examined. Other models based on layered QM
descriptions of the solute (including layers treated
all at a QM level or mixed QM and semiclassical
descriptions) will be considered in section 7.
(2) The solute-solvent interactions are limited to
those of electrostatic origin. Other interaction terms
exist, and they must be taken into account to have a
well-balanced description of solvent effects (see sec-
tion 3). This point has to be emphasized; our choice
of paying attention first to the electrostatic model,
dictated by convenience of exposition, should not lead
the reader to the false conclusion that only electro-
statics is important in solvation. Often the opposite
happens, and in addition some problems arising in
treating the electrostatic term are greatly alleviated
by the consideration of other solute-solvent interac-
tions.
(3) The model system is a very dilute solution. In
other words, it is composed of a single solute molecule
(including, when convenient, some solvent molecules,
the whole being treated as a supermolecule at a
homogeneous QM level) immersed in an infinite
solvent reservoir.
(4) The solvent is isotropic, at equilibrium at a
given temperature (and pressure). Possible exten-
sions beyond the isotropic approximation will be
considered mainly to show new potentialities of the
most recent solvation models.
(5) Only the electronic ground state of the solute
will be considered. Extensions to other electronic
states will be considered in section 5.
(6) No dynamic effects will be considered in the
basic model. Under the heading of “dynamic effects”
there is so large a variety of important phenomena
that it would require a separate review. The main
aspects of these phenomena will be considered in
sections 5 and 6.
Once we have better defined what we intend with
the expression “basic QM continuum model”, we can
consider some of its essential elements.
2.2. Cavity
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 differ-
ently 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 On-
sager,13 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:
Obviously these requirements are in contrast with
the description of the whole system given by any QM
level. The electronic charge distribution of an isolated
molecule, in fact, persists to infinity. In a condensed
medium the conditions on FM
e at large distances are
less well-defined, but in any case there will be an
overlap with the charge distribution of the medium,
not explicitly described in continuum models but
existing in real systems.
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 condi-
tion may lead to deformations in the charge distribu-
tion 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 el-
lipsoids, but molecules are often far from having a
spherical or ellipsoidal shape.
FM ) FM
e + FM
n (5)
Quantum Mechanical Continuum Solvation Models Chemical Reviews, 2005, Vol. 105, No. 8 3003

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 environ-
ment. Such a cavity has been inserted into the
Gaussian computational package.17 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
McWeeny18 and employed by Bentley and others19-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 func-
tion of the system constituted by the molecule and
the probe. Following the AIM topological analysis,25
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 molec-
ular surfaces and cavities. This technique is of
potential interest as a benchmark for solvent calcula-
tions far from the ambient temperature; otherwise,
the approach is too costly to be used in standard
applications.
The generally adopted compromise between ana-
lytical but too simple and realistic but computation-
ally 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 Bondi26 (5000 citations
in the past 15 years). This set was confirmed as the
recommended one some years ago,27 after the exami-
nation of a quite large number of intermolecular
contact distances drawn from the Cambridge Struc-
tural 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 struc-
tures confirm the Bondi values, with the only excep-
tion 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,28 available in many tabulations of
physicochemical data, for example, in the CRC
Handbook.29
A third set of values has been inserted in the latest
versions of the Gaussian package.17 It is drawn from
the compilation of data for the universal force field
(UFF),30 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 com-
puted 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 periph-
ery 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, respect-
ively).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 mol-
ecule).
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 Connolly32,36 can be consid-
ered 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 website37 for a clear and concise
review accompanied by a sizable selection of refer-
ences, 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 compute
3004 Chemical Reviews, 2005, Vol. 105, No. 8 Tomasi et al.

solvation 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 ex-
pression 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 imple-
mentation 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 ro-
bust; the SES may present singularities and cusps;
these last problems are better treated with methods
developed by Gogonea38,39 and by Sanner et al.40 The
latter have developed a computational code called
MSMS,41 standing for Michael Sanner’s Molecular
Surface, which has received much attention, espe-
cially 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 algo-
rithm 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 sur-
face is defined following a different strategy, origi-
nally conceived in Pisa around 1984 and finalized in
1986 by Pascual-Ahuir in his Ph.D. research.42 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 examina-
tion 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 cur-
vature.
2.3. Solution of the Electrostatic Problem
The physics of the electrostatic solute-solvent
interaction is simple. The charge distribution FM of
the solute, inside the cavity, polarizes the dielectric
continuum, which in turn polarizes the solute charge
distribution. This definition of the interaction corre-
sponds to a self-consistent process, which is numeri-
cally solved following an iterative procedure. It is
important to remark that the corresponding interac-
tion potential is the one we shall put in the Hamil-
tonian of the model. As this potential depends on the
final value of FM 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).
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.
Quantum Mechanical Continuum Solvation Models Chemical Reviews, 2005, Vol. 105, No. 8 3005

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 paper13 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 equa-
tion
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 VM generated by the charge distribution FM
and the reaction potential VR 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 FM) 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 FM 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 C1, C2, ..., Cn, each containing
a different FM1, FM2, ..., FMn 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 R 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 nb 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 review1 the ap-
proaches in use were classified into six categories,
namely, (1) the apparent surface charge (ASC) meth-
ods, (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 develop-
ments 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 con-
sidered in the present review: the interested reader
is referred to our previous review.1
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 VR of
eq 9. There are no approximations in this elabora-
tion: the definition of VR given with eq 16 is exact
- rB[(rb)rVB(rb)] ) 4ðFM(rb) (6)
-r2V(rb) ) 4ðFM(rb) within C (7)
-r2V(rb) ) 0 outside C (8)
V(rb) ) VM(rb) + VR(rb) (9)
lim
rf∞
rV(r) ) R (10)
lim
rf∞
r2V(r) ) â (11)
[V] ) 0 on ¡ (12)
[@V] ) 0 on ¡ (13)
[V] ) Vin - Vout ) 0 (14)
[@V] ) (@V@n)in -  (@V@n)out ) 0 (15)
Vó(rb) ) s¡
ó(sb)
jrb - sbj
d2s (16)
3006 Chemical Reviews, 2005, Vol. 105, No. 8 Tomasi et al.

when ó(s) is defined according to the proper electro-
static 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 simpli-
fication, 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, qk, in terms
of the local value of ó(s) on each of these tesserae
times the corresponding area Ak. The integral of eq
16 is thus transformed in the following finite sum:
Actually the local value of the potential necessary to
define qk 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 qk but instead the original ones in
terms of a continuous surface charge ó(s). A descrip-
tion 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 philoso-
phy 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,8 and subsequently
implemented in local and official versions of various
QM computational packages.48 More recently, PCM
was renamed D-PCM (D stands for dielectric)49 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 clus-
tering, etc., and it can account for a continuous shift
from a single cavity to two cavities during a dissocia-
tion 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 permit-
tivity, 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)
where i is the dielectric constant of the region i.
At the boundary of two regions i and j, there is an
ASC distribution given by
where nbij is the unit vector at the boundary surface
pointing from medium i to medium j. The basic ASC
model is so transformed into a similar system, with
several ASCs that must be treated all on the same
footing.
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
where nb indicates the unit vector perpendicular to
the cavity surface and pointing outward.
After its first formulation, the DPCM was revised
many times both in its theoretical aspects and in its
numerical implementation; among all of these revi-
sions, a fundamental one was proposed in 1995 by
Cammi and Tomasi50 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 develop-
ment. Some of these extensions will be examined in
following sections of this review, but here we cannot
Vó(rb) = ∑
k
ó(sbk)Ak
jrb - sbkj
) ∑
k
qk
jrb - sbkj
(17)
PBi(rb) ) -
i - 1
4ð
rVB(rb) (18)
óij ) - (PBj - PBi)ânbij (19)
ó(s) )  - 1
4ð
@
@nb
(VM + Vó)in (20)
Quantum Mechanical Continuum Solvation Models Chemical Reviews, 2005, Vol. 105, No. 8 3007

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, Ga, and thus in the
corresponding operators Sa and Da, 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.77 We
recall some of their properties: the operator Si is self-
adjoint, and Di
/ is the adjoint of Di for the scalar
product. Besides, SiDi
/ ) DiSi.
Equation 26 may be further simplified using the
equality (2ð - Di)VM + Si(@VM/@n) ) 0; in this way
the expression for g in eq 27 can be rewritten as78
and thus the surface charge ó depends only on the
potential VM (and no longer on a normal component
of the field) exactly as in the COSMO approach. This
simplification is important from both numerical and
formal points of view.
Numerically, it is advantageous not only because
the calculation of a single scalar function (the poten-
tial) 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 instabili-
ties 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, Se ) Si/, Di
) De, we can transform the IEF operators A of eq 27
and g of eq 29 in
and thus eq 26, defining the ASC, reduces to78
where, once again, we have used the relationship Si
Di
/ ) DiSi. In the following, this version of IEF will
be denoted IEF(V) to indicate that only the solute
electrostatic potential (VM) is required to determine
the ASCs.
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, Gi(x,y) ) 1/4ðjx - yj, the problem is shifted
to the evaluation of the exterior, Ge(x,y). Analytical
expressions of this function are available for standard
isotropic solvents characterized by a constant and
scalar permittivity  [namely Ge(x,y) ) 1/(jx - yj)],
but also for anisotropic environments characterized
by a tensorial but constant permittivity E for ionic
solutions described in terms of the linearized Pois-
son-Boltzmann equation (see section 2.3.1.6 for
further details) and for sharp planar liquid-metal
interfaces.79 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 permit-
tivity depending on the position.80
As a final comment it is worth noting that IEF
contains as subcases both the DPCM and the COS-
MO 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 (@VM/@n)
instead. To recover COSMO ASC, on the other hand,
we have to consider that  f ∞, and thus all of the
terms in eqs 27 and 29 involving the Se operator can
be neglected (Se 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 models81-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 “outly-
ing 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.86 As a conse-
quence, 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
where the integration is over the whole volume
excluding the molecular cavity.
The acronym used to indicate this formulation was
SVPE, meaning that both surface and volume polar-
ization 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
VR(x) ) Vó(x) + Vâ(x)
Vâ(x) ) sext
â(y)
jx - yj
dy ) - ( - 1 )sext
FM(y)
jx - yj
dy
(32)
g ) [(2ð - De) - SeSi
-1(2ð - Di)]VM (29)
A ) (1 - 1/)[2ð( + 1)/( - 1) - Di]Si
g ) (1 - 1/)(2ð - Di)VM
(30)
[2ð( + 1 - 1) - Di]Sió ) - (2ð - Di)VM (31)
3010 Chemical Reviews, 2005, Vol. 105, No. 8 Tomasi et al.

Recently, two independent novel schemes have
been proposed100,101 to sample the GEPOL molecular
surface without an explicit discretization of the cavity
surface, which may lead to unphysical discontinui-
ties. 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 con-
stant within each tessera, it is possible to define a
set of point charges qk in terms of the local value of
ó(s) on each of these tesserae times the corresponding
area. In doing this the electrostatic equation pre-
sented 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
where K is a square matrix T  T collecting the
cavity geometrical factors (the tesserae representa-
tive points sk and the corresponding areas) and the
dielectric constant of the medium and q and f are
column matrices, the first containing the unknown
charges and the second the values of the proper
electrostatic quantity, namely, the normal component
of the electric field En or the electrostatic potential
V, calculated at the tesserae. Equation 34 represents
the electrostatic BEM problem we have to solve.
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 elabora-
tions; instead, we provide a table with the expres-
sions 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 permit-
tivity ) 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 En 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;102
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 solu-
tions have been proposed for the diagonal elements.
In particular, for S the following approximation is
Table 1. Matrices To Be Used in Equation 34 To Get the Apparent Charges in Various ASC Methods
K f
DPCM (2ð + 1 - 1 A-1 - D*)-1 with Dij ) (sbi - sbj)ânbjjsbi - sbjj3 98ifj Dij/ En
COSMO [CPCM] S-1 with Sij )
1
jsbi - sbjj
V
IEFPCM(iso) [IVPCM] {[2ð( + 1 - 1)A-1 - D]S}-1[2ðA-1 - D] V
SS(V)PE {[2ð( + 1 - 1)A-1 S - A-1 (DAS + SAD*)2 ]}-1[2ðA-1 - D]
IEFPCM {[2ðA-1 - De]S + Se[2ðA-1 + D*]}-1{[2ðA-1 - De] -
SeS
-1[2ðA-1 - D]}
V
aniso
[De]ij )
(sbi - sbj)ânbj
xdet E[(-1â(sbi - sbj))â(sbi - sbj)]
3/2
[Se]ij )
1
xdet E[(-1â(sbi - sbj))â(sbi - sbj)]
1/2
ionic
[De]ij )
exp[-kjsbi - sbjj][1 + kjsbi - sbjj](sbi - sbj)ânbj
jsbi - sbjj
3
[Se]ij )
exp[-kjsbi - sbjj]
jsbi - sbjj
q ) -Kf (34)
3012 Chemical Reviews, 2005, Vol. 105, No. 8 Tomasi et al.

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 inte-
grals and their derivatives, for example
where
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 dis-
persed ions. The basic electrostatic equations modify
as
where k accounts for the ion screening: its inverse
is known as the Debye screening length (1/k), namely,
k2 ) 8ðe2I/kT (I is the ionic strength ∑i zi2ci/2 for a
dissolved salt containing two or more types of ion i,
each having charge zie and an average bulk concen-
tration ci).
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 treat-
ment of statistical physics in the thermodynamic
equilibrium approximation, represents the distribu-
tion 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.115 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 applica-
tion 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 ap-
proaches 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 involv-
ing 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
where the function G is the Green function exp(-kjx
- yj)/jx - yj.
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 Trian-
gular 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 ele-
ment. 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.129
Sij ) 〈ijj〉 )
erf(œijjsbi - sbjj)
jsbi - sbjj
Dij )
@Sij
@sbj
nˆj )
(sbi - sbj)nˆj
jsbi - sbjj
3 [erf(œijjsbi - sbjj) -
2
xð
œij
j
sbi - sbj
j
e-(œijjsbi-sbjj)
2]
(41)
œij ) œiœj/xœi
2 + œj
2
-r2V(rb) ) 4ðFM(rb) within C (42)
- (r2 - k2)V(rb) ) 0 outside C (43)
s¡
@G(x,y)
@ny
V(y) dy + 1
2
V(y) )
s¡
G(x,y)
@V(x,y)
@ny
dy
(44)
3014 Chemical Reviews, 2005, Vol. 105, No. 8 Tomasi et al.

seminal paper, but a simple effective expression was
given
where Rij ) (RjRj)0.5 is the geometrical mean of the
pertinent pair of the so-called generalized Born radii
Ri (which may be interpreted roughly as the distance
from each atom to the dielectric boundary). The
exponent rij
2/2Rij)2 is a damping factor between the
charges placed on atoms i and j, leading smoothly to
the Born formula, when the spheres are separated
and placed at very large distances among them, and
to the normal Coulomb operator, when two spheres
merge.150
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 pre-
factor 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 struc-
tures are evaluated.
A fundamental problem arising in the GB ap-
proach, and present in all of its successive formula-
tions, is the definition of the appropriate Born radii
Ri.151-155
For every charged atom i, the effective Born radius
Ri is related to the effective Born free energy of
solvation, ¢Gi, 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 refer-
ence system atom i bears unit charge and all other
atoms are neutral. Exact calculation of Ri 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 Ri
with sufficient accuracy and efficiency. An analytical
formula proposed by Qiu et al.156 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 ap-
proach 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 macromolecules150,163-165 such as nucleic
acids and proteins. Here, we obviously cannot men-
tion 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.166 This method has been elaborated and
distributed in a large variety of versions, with dif-
ferent 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
papers167 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 code168),
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 GAMESS-
PLUS,173 which is an add-on module to the
GAMESS,97,98 and in HONDOPLUS174,175 program.
The latest version of the method is called SMx-
GAUSS,176 and it may use Gaussian output files or
be connected with Gaussian 0317 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 electro-
static interaction according to the GB model, that is,
in terms of atomic charges qi and of the modified
Coulomb operator (see eq 53). The qi 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 Lo¨wdin charges, and
finally charges that according to the Cramer and
Truhlar definitions are of classes III and IV. The class
IV charges177 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 modi-
fication of the set of qi 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 GP, can be so computed.
Another term of the free energy difference with
respect to the isolated molecule is also computed.
This term, called ¢EEN, 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
fGB ) xrij
2 + Rij
2 e-Dij (54)
Quantum Mechanical Continuum Solvation Models Chemical Reviews, 2005, Vol. 105, No. 8 3017

solutions. We give here some reasons supporting this
statement.
BEM is capable of giving values of the unknown
(in our case the sources of the reaction electrostatic
potential) in the interior of the domain in a pointlike
form, thus avoiding the interelement continuity
problems that disturb FEM applications to solvation.
The reduction of the spatial problem from three to
two dimensions greatly facilitates mesh refinement
studies, and it leads to a system of algebraic equa-
tions much smaller than in an equivalent FEM
formulation. This advantage is partly compensated
by the sparcity of FEM linear equations, which is
exploited in the PBF method described above.
Comparison between the two approaches could lead
to different conclusions in some specific solvation
cases. BEM is less efficient than FEM for systems
with one or two spatial dimensions extremely small
with respect to the others but dimensionally effective,
as is the case of bicontinuous liquid systems. BEM
has difficulties in treating problems exhibiting rapid
changes in the physical properties of a domain. This
last problem has been excellently solved in solvation
chemistry for the inhomogeneities and anisotropies
of the molecules inside the cavity, but, for example,
it remains a problem for the solvent anisotropies due
to nonlinear dielectric saturation.
Although BEM and FEM aim at giving a numerical
solution of the integral equations, the FD methods
aim at solving the differential equations. The dif-
ferential equation to solve is, once again, the basic
electrostatic PB problem (eq 43); we note that, with
FD methods, also the nonlinear version of such an
equation can be solved.
The FD methods have a long history in continuum
solvation, and they are very popular, especially in
classical version making no explicit use of QM
methods.184-186
In these approaches, a discrete approximation to
the governing partial PB differential equations is
obtained through a volume-filling grid. Ideally, a
boundary-conforming grid (i.e., one that does not
intersect the molecular surface) is preferred, but such
a grid is difficult to generate for a complex molecular
shape. Therefore, in most implementations, such as
the widely used UHBD187 (University of Houston
Brownian Dynamic)188-190 and DelPhi191 code (now
included in Insight II package),192,193 a regular lattice
is laid over the molecule and cells are allowed to go
over the molecular surface. A regular lattice arrange-
ment is also useful for an efficient multigrid solution
technique for solving the algebraic equations result-
ing from the discretization of the partial differential
equations.194
A further FDM approach has been developed by
Bashford, and it has been implemented in a compu-
tational code known as MEAD (macroscopic electro-
statics with atomic detail).195 The MEAD electrostatic
model has been coupled with QM techniques: an
inner region containing the active-site atoms is
treated quantum mechanically, whereas the sur-
rounding region is treated as a classical electrostatic
system that generates a reaction field and, possibly,
a field due to permanent atomic charges of the
protein. This technique has been used to calculate
redox properties196 and pKa values.197
There are many similarities between FE and FD
methods. In both cases, a grid of points covering the
space is defined. In FE procedures, these points have
to be considered as representative points of a finite
three-dimensional domain, and in the FD procedures,
nodes in which the differential equation has to be
solved. In addition, the finite element method differs
from the finite difference approach in that variational
principles rather than finite difference approxima-
tions are used to derive the discrete equations.198,199
A major reason for adopting a finite element ap-
proach is that it accommodates unstructured grids,
which offer improved geometric flexibility and vari-
able mesh spacing, albeit at considerably higher per-
node storage and CPU costs than a regular lattice
arrangement.
2.4. Solution of the Quantum Mechanical Problem
As described in section 1, the Hamiltonian appear-
ing in the Schro¨dinger equation of the basic model
is the effective Hamiltonian of eq 3, namely
This effective Hamiltonian is composed by two terms,
the Hamiltonian of the solute HM
0 (i.e., the focused
part M of the model) and the solute-solvent interac-
tion term Vö int (i.e., the solvent reaction potential);
The definition of Vö int depends on the method em-
ployed to set the electrostatic problem, which has to
be solved within the framework of the QM eq 58.
To describe this procedure in more detail let us first
consider an intuitive formulation of the problem.
2.4.1. Intuitive Formulation of the Problem
Being interested in a QM problem, let us introduce
the currently used Born-Oppenheimer (BO) ap-
proximation, which relies on a partition of the solute
variables into electronic and nuclear coordinates.
The solute charge distribution is conveniently
divided into electronic and nuclear components FM-
(r) ) FM
e (r) + FM
n (r). We recall that the QM procedure
does not affect the nuclear component of FM, although
it will modify the electronic component with respect
to a starting provisional definition. This modification
is done in an iterative way under the action of Vö int,
which is in turn modified in the iterative cycle.
Let us look now at the structure of Vö int in eq 59.
This operator can be divided into four terms having
a similarity with the two-, one-, and zero-electron
terms present in the solute Hamiltonian.
To show it in an intuitive way, we consider the
solute-solvent interaction energy Uint given as the
integral of the reaction potential times the whole
charge distribution FM.
The interaction potential has, as sources, the two
components of FM, and thus it is composed of two
terms, one stemming from the electronic distribution
of the solute M and one from its nuclear distribution.
Hö effj¾〉 ) Ej¾〉 (58)
Hö eff ) Hö M
0 + Vö int (59)
Quantum Mechanical Continuum Solvation Models Chemical Reviews, 2005, Vol. 105, No. 8 3019

Uint can thus be partitioned into four terms
where Uxy corresponds to the interaction energy
between the component of the interaction potential
having as source FM
x (r), namely, Vö int
x , and the charge
distribution FM
y (r).
As in all QM applications to molecular systems, the
solution of the Schro¨dinger eq 58 is based on an
expression of the unknown wave function ¾ in terms
of molecular orbitals expressed over a finite basis set
{ł}. Within this framework, it is convenient to
consider the FM
e charge distribution expressed as
sums of contributions due to the elementary charge
distributions ł/íłv, whereas FMn (r) is
Following this formalism, three different QM op-
erators appear, namely, Vö nn, Vö ne (it may be shown
that Une and Uen are formally identical), and Vö ee.
These correspond to zero-, one-, and two-electron
operators appearing in HM
0 , respectively. We note
that the zero-order term gives rise to an energetic
contribution Unn, which is analogous to the nuclei-
nuclei repulsion energy Vnn and, thus, it is generally
added as a constant energy shift term in HM
0 . The
conclusion of this analysis is that we have intuitively
defined four operators (reduced in practice to two,
plus a constant term) which constitute the operator
Vö int of eq 59.
Before proceeding further in this presentation, we
note that this analysis has a direct link with the ASC
approaches, but it is valid also for the MPE methods,
with minor changes, and for the FE and FD methods,
as well. Note that the QM versions of the last two
approaches describe the interaction potential in
terms of cavity apparent surface charges, and thus
they exploit the same features of ASC methods. By
contrast, for GB methods the formal setting is dif-
ferent. These methods reduce, in fact, the solute
density function to atomic charges that are treated
in terms of the generalized Coulomb operators we
have already introduced (see section 2.2.3).
2.4.2. Electrostatic Operators
The next step consists of solving the QM problem
with the determination of the electrostatic interaction
operators nested within it. It is convenient now to
make explicit use of the apparent surface charges,
because they simplify the exposition; as stated above,
this does not represent a real restriction as almost
all of the QM solvation methods follow the same basic
pattern.
The most naive (and transparent) formulation of
the process of mutual interaction between real and
apparent charges is that used in the first version48
of DPCM. We recall it here, even though it was
already presented in the 1994 review, as helpful in
the understanding of the basic aspects of the mutual
polarization process.
One starts from a given approximation of FM
e (let
us call it FM
0 ) that could be a guess, or the correct
description of FM
e without the solvent, and obtains a
provisional description of the apparent surface charge
density, or better, of a set of apparent point charges
that we denote here {qko,o}. These charges are not
correct, even for a fixed unpolarized description of
the solute charge density, because their mutual
interaction has not been considered in this zero-order
description. To get this contribution, called mutual
polarization of the apparent charges, an iterative
cycle of the PCM eq 20 (including the self-polarization
of each qk) must be performed at fixed FM
0 . The result
is a new set of charges {qko,f}, where f stands for
final. The {qko,f} charges are used to define the first
approximation to Vint, and a first QM cycle is per-
formed to solve eq 58. With the new FM
1 the inner
loop of mutual ASC polarization is performed, again
giving origin to a {qk1,f} set of charges. The procedure
is continued until self-consistency.
We remark that, in this formulation, we have
collected into a single set of one-electron operators
all of the interaction operators we have defined in
the preceding section, and, in parallel, we have put
in the {qk} set both the apparent charges related to
the electrons and nuclei of M. This is an apparent
simplification as all of the operators are indeed
present.
It is interesting here to note that this nesting of
the electrostatic problem in the QM framework is
performed in a similar way in all continuum QM
solvation codes, including GB methods, in which
there is an iterative updating of the atomic charges
during the QM cycle.
We report here the Schro¨dinger equation of the
ASC version of the basic model with the introduction
of a new formalism to make the exposition more
general:
With the superscript R we indicate that the corre-
sponding operator is related to the solvent reaction
potential and with the subscripts r and rr′ the one-
or two-body nature of the operator. The convention
of summation over repeated indexes followed by
integration has been adopted. The emphasis given
in this formula to the density operators is related to
a formal problem arising in focused models: the main
component of the system, M, interacting with the
remainder does not have a wave function satisfying
all of the rules of canonical QM. Only the density
function has a more definite status (see ref 200 for
comments on this point). For this reason we use only
density operators and expectation values in expres-
sion 62.
The Fˆr
e Vö r
R operator describes the two components
of the interaction energy we have previously called
Uen and Une. In more advanced formulations of
continuum models going beyond the electrostatic
description, other components are collected in this
Uint ) Uee + Uen + Une + Unn (60)
FM
n ) ∑
R
nucl
ZRä(r - RR) (61)
Hö effj¾〉 ) [Hö M0 + Fˆre Vö rR + Fˆre Vö rr′R 〈¾jFˆr′e j¾〉]j¾〉 )
Ej¾〉 (62)
3020 Chemical Reviews, 2005, Vol. 105, No. 8 Tomasi et al.

contraction of j¾HF
S 〉〈¾HFS j, and thus it is nonlinear.
It must be added that this nonlinearity is of the first
order, in the sense that the interaction operator
depends only on the first power of FM
e . The model is
in fact based on an approximation motivated by
physical considerations. The solvent reaction field
has, in general, a relatively low strength, and thus a
linear response regime is sufficient to describe its
effects. For this reason, the expansion of the interac-
tion potential in powers of the reaction field is limited
to the first order in almost all of the continuum
methods. In other words, it may be said that these
methods treat the linear aspect of a nonlinear
problem. True nonlinear solvation effects are, how-
ever, present. Generally, they are treated in other
ways, for example, making reference to phenomeno-
logical concepts, such as dielectric saturation or
electrostriction, and directly using higher terms in
expansion of the interaction potentials (more details
on these aspects will be treated in section 4).
Some comments about nonlinearities in the Hamil-
tonian are here added. This subject is rarely treated
in an explicit way in quantum chemistry, despite the
large use of the Fock operator, which is nonlinear,
and, in the case of ionic solutions, of the use of the
Poisson-Boltzmann approach, which also leads to a
nonlinear problem. The case we are considering here
is called scalar nonlinearity (in the mathematical
literature it is also called “nonlocal nonlinearity”):
204 this means that the operators are of the form P(u)
) 〈Au,u〉Bu, where A and B are linear operators and
〈.,.〉 is the inner product in a Hilbert space. The
literature on scalar nonlinearities applied to chemical
problems is quite scarce (we cite here a few
papers7,205-210), but the results justified by this ap-
proach are of universal use in solvation methods.
Let us now come back to eq 74. Attentive readers
will have noticed that we have changed the symbol
for the energy. The symbol G emphasizes the fact
that this energy has the status of a free energy. The
explicit identification of the functional (eq 74) with
the free energy was first done by Yomosa,205 to the
best of our knowledge. In the 1994 review alternative
justifications for the factor 1/2 in the expression of the
energy were given, deduced from the perturbation
theory, by statistical thermodynamics and by clas-
sical electrostatics, all valid for a linear response of
the dielectric. We report here only a consideration
based on classical electrostatics. Half of the work
spent to insert a charge distribution (i.e., a molecule)
into a cavity within a dielectric corresponds to the
polarization of the dielectric itself, and it cannot be
recovered by taking the molecule away and restoring
it in its initial position. This half of the spent work
is irreversible, and it has to be subtracted from the
energy of the insertion process to get the free energy
(or the chemical potential).
2.4.5. QM Descriptions beyond the HF Approximation
In the past few years, a large effort has been
devoted to extending solvation models to QM tech-
niques of increasing accuracy. In this way, models
originally limited to the Hartree-Fock level can be
now used in many post-SCF calculations. This com-
putational extension has been accompanied by a
reformulation of the various QM theories so as to
include the specifics of the solvation model. Such a
requirement is particularly pressing for those models
in which an effective Hamiltonian is introduced.
Electron correlation is more commonly introduced
into these techniques using CI or MCSCF and DFT
methods.
DFT does not require any specific modification with
respect to the formalism presented in the previous
section, as far as concerns the solvation terms, and
thus there are numerous examples of DFT applica-
tions of continuum solvation models. As a matter of
fact, continuum solvation methods coupled to DFT
are becoming the routine approach for studies of
solvated systems.
Applications of continuum solvation approaches to
CI or MCSCF wave functions have quite a long
history. In 1988 Mikkelsen et al.134 presented a
MCSCF version of their multipole approach (see
section 2.3.2) first limited to equilibrium solvation
and then extended to nonequilibrium descriptions137
(see section 5) and to response methods136 (see section
6). A similar history is presented by the family of
PCM approaches, which have been first applied to
equilibrium211 (including gradients212) and nonequi-
librium78,213,214 and subsequently to the response
method.215
Also, configuration interaction (CI) approaches
soon attracted the interests of researchers working
in the field of solvation methods.205,216-218 In this case,
however, some delicate issues not present for isolated
systems appear.
As remarked by Houjou et al.,219 an important point
to analyze is the physical meaning given to the
excitation energy of solvated systems when obtained
directly from the CI calculation. Karelson and Zern-
er220 have indicated that the solvent effect is incor-
porated ‘‘automatically” into the CI only by adding
solute-solvent interaction term(s) to the Fock matrix
for the in vacuo state, whereas the resulting excita-
tion energy must be corrected for the electronic
relaxation and for the costs required for the polariza-
tion of solvent. Similarly, Klamt derived the correc-
tion term corresponding to the electronic relaxation,
denoted ‘‘electronic screening correction” in ref 221,
and the diagonal correction in the CI matrix.
In addition, the nonlinear dependence of the sol-
vent reaction field operators on the solute wave
function introduces a new element of complexity. In
a CI scheme, this nonlinearity is solved through an
iteration procedure; at each step of the iteration, the
solvent-induced component of the effective Hamilto-
nian is computed by exploiting the first-order density
matrix of the preceding step. The main consequence
of a correct application of this scheme is that each
electronic state, ground and excited, requires a
separate calculation involving an iteration optimized
on the specific state of interest. This procedure has
been adopted by Mennucci et al.78 in their IEFPCM-
CI approach as well in the successive extension to
multireference CI wave functions (including pertur-
bative corrections)222 within the framework of the
CIPSI method.223-225
3024 Chemical Reviews, 2005, Vol. 105, No. 8 Tomasi et al.

models. Differences between the two approaches will
be commented on at the appropriate place.
The basic information on the nature of molecular
interactions giving origin to the solvent effect derives
from the theoretical analysis and interpretation of
the interactions occurring in molecular clusters.
There is an extremely large literature on this
subject, covering a span of more that 50 years, that
cannot be summarized here. Studies have progressed
from simple systems, dimers composed of two very
simple molecules, increasing step by step the number
of components of the aggregate and the internal
complexity of the molecules. Several aspects of these
studies are pertinent to our subject, but we shall limit
ourselves to the basic information derived from
bimolecular associations, complemented by some
aspects derived from the analyses of trimers and
larger aggregates.
The interaction energy between two molecules A
and B can be expressed as the difference between the
energy of the supermolecule and that of the two
separate partners:
This energy has been analyzed in a significant
number of cases using a variety of techniques based
on perturbation theory approaches as well as varia-
tional techniques. The details and the variants
among the different approaches are not of interest
here. What is of interest is to examine a partition of
¢E(AB) into terms each having a physical interpre-
tation and the possibility of being individually com-
puted. This partition reads
The Coulomb term, Ecoul, is defined as the electro-
static energy among the two unperturbed (i.e., not
polarized) partners. This term is very sensitive to the
mutual orientation of the two partners and strongly
anisotropic at short distances. The interaction may
be positive or negative and, for a given A-B system,
it may change sign according to the orientation when
at least one of the two molecules has a (local) charge
distribution with a dipole. The Coulomb energy
continues to be numerically important also at large
monomer separations; actually, it is the term having
the highest long-range character. The typical distance
dependency is R-1.
The polarization term, Epol, has its origin in the
mutual electrostatic polarization of the two charge
distributions, and it is everywhere negative. At short
monomer separations, it is moderately anisotropic,
being sensitive to the mutual orientation of the two
partners. Epol has a relatively long-range character:
the typical distance dependency is R-4. Epol(AB) can
be easily decomposed into two terms, the polarization
of A under the influence of B, Epol(A r B), and the
polarization of B under the influence of A, Epol(B r
A).
The exchange term, Eexc, is related to the presence
in the QM formulation of intermonomer antisymme-
trizer operators, which ensure the complete antisym-
metry of the electrons in the AB dimer. Complete
antisymmetry means to introduce repulsion forces
between the electron distribution of the monomers
related to the Pauli exclusion principle. Eexc is posi-
tive for every orientation of the partners, roughly
proportional to the mutual overlap of the two charge
distributions, and rapidly decaying with the distance.
The typical distance dependency is dictated by a R-12
term.
The dispersion term, Edis, can be related to dy-
namical effects, as the formation of an instantaneous
dipole moment in a monomer due to the motion of
the electrons. This dipole polarizes the charge dis-
tribution of the other partner and gives origin to an
interaction energy similar in some sense to Epol. It is
a distinct phenomenon, however, and must be sepa-
rately treated. Dispersion energies are negative,
relatively weak, and decaying as R-6.
The charge transfer term, ECT, is not present in the
standard analyses of ¢E based on perturbation
theory, whereas it appears in the decompositions
based on the variational approach [i.e., starting from
the variational calculation of the E(AB) energy] in a
form that seems to be a bit artificial, being present
also in symmetric AA dimers under the form of
monomer contributions which cancel each other. In
chemistry, charge transfer effects are quite important
in large classes of chemical reactions; these, when
considered under the viewpoint of a theory address-
ing the general characteristics of solutions, belong to
the category of “rare events”, and they will not be
treated here. However, we also note that small
amounts of charge transfer, not related to chemical
reactions, can play a role in assessing solvent effects
on some molecular properties, and in the elaboration
of solvation models we have surely to pay attention
to these interactions even if, at present, they have
attracted very limited attention (see section 3.2.4).
This very schematic presentation of the decomposi-
tion of dimeric interactions is far from being com-
plete, and our presentation of decomposition results
should be extended to the more complex analyses
performed on trimers and larger aggregates. As a
matter of fact, here we shall be even more concise,
reporting only the conclusions that seem more inter-
esting for liquids.
(1) There are no new types of interactions to
include in the model. Specific interactions, often
invoked or described with empirical ad hoc formulas,
as, for example, hydrogen bonds, are in this approach
described in terms of the contributions we have
already presented.
(2) All of the contributions, with the exception of
Ecoul, are nonadditive, and the weights of such non-
additivity ares very different for the different terms.
The larger nonadditivity, both in weight and in
range, is that due to the polarization energy. Mol-
ecules in an assembly are strongly influenced by
mutual many-body polarization effects. The nonad-
ditivity of exchange is in se very short-ranged. The
nonadditivity of dispersion is moderate and plays an
important role in some special cases only, as, for
example, in clusters composed of rare gas atoms in
which classical electrostatic terms are zero. The
nonadditivity of charge transfer has been less stud-
¢E(AB) ) E(AB) - E(A) - E(B) (76)
¢E(AB) ) Ecoul + Epol + Eexc + Edis + ECT (77)
3026 Chemical Reviews, 2005, Vol. 105, No. 8 Tomasi et al.

ied. It should be a relatively short-range effect,
especially for the energy.
3.2. Use of Interactions in Continuum Solvation
Approaches
We have now the elements to give a more detailed
appraisal of the use of these interactions in the
continuum description of solute-solvent interactions.
Among the various solvation methods, reference will
be made mostly to the PCM, which is the most
complete in this field, and to GB methods (see section
2.3.3) that have paid much attention to the consid-
eration of nonelectrostatic terms.
Continuum methods start from the definition of a
liquid, the pure solvent, equilibrated at the given
macroscopic conditions (temperature and pressure,
for instance). The solute molecule must be inserted
within this liquid. This specific aspect of continuum
models gives rise to a problem that can be solved by
introducing in the definition of the overall intermo-
lecular interaction energy (see eq 77) the concept of
the energy of formation of a suitable cavity within
the liquid. We shall call it Gcav.
The other contributions to the interaction energy
can be collected into three terms, named electrostatic,
repulsion, and dispersion. There are some changes
with respect to the partition given in eq 77. Coulomb
and polarization terms are grouped together; the
separate calculations of the two contributions can be
done (and it has been done in the first DPCM
version), but this separation is of no practical help
in the QM calculations. In any case, it may be used
for more detailed analyses about the origin of the
interactions. The repulsion term takes into account
the presence of the cavity; the physical origin of the
term is unchanged, but all solvent molecules are kept
out of the cavity. The dispersion term, too, is left
unchanged, and it involves interactions with solvent
molecules out of the cavity.
We have now the elements to write the interaction
energy decomposed into four terms, the three we
have just mentioned and the cavity formation term.
Before doing this, we have, however, to consider two
additional points.
The first is related to the use of ab initio QM
methods. These methods give the molecular energy
measured with respect to a quite peculiar reference,
that is, noninteracting electrons and nuclei. The
second point to consider is related to the existence
of free energy contributions due to other degrees of
freedom which are present in the focused model but
are neglected in the BO approximation we have used
so far. These degrees of freedom are related to the
thermal motions of nuclei. We shall call Gtm the
corresponding free energy contribution. More com-
ments will be made below, but we may now give an
explicit definition of the energy in continuum ab initio
models and of the solvation free energy, two related
but not identical concepts.
Definition 1: Free Energy of the Solution in ab
Initio Continuum Models. The energy of a system
composed by a liquid and a solute is defined making
reference to an ideal system composed of the pure
liquid at equilibrium under given values of the
macroscopic parameters (pressure, temperature) and
by the appropriate number of noninteracting elec-
trons and nuclei necessary to describe the solute and
at zero kinetic energy. Putting the energy of the
reference system equal to zero, the free energy of the
solution at the equilibrium is given by the following
expression:
The definition of the solvation free energy given
below can be considered as a supplementary defini-
tion drawn from eq 78, but it is the primary definition
for computational methods, which do not allow an
evaluation of the solute molecular energy. This is the
case for computational procedures based on semiem-
pirical QM methods or the use of classical models.
Definition 2. Solvation Free Energy. The solvation
free energy, measured with respect to a system
composed of the pure unperturbed liquid at equilib-
rium and by the solute molecule(s) in a separate
phase considered as an ideal gas, is given by the
following expression:
If we define ¢Gel ) Gel - G0, where G0  E0 is the ab
initio energy of the isolated molecule in the BO
approximation and Gel is the analogous quantity
computed in solution, we may obtain the first term
of eq 79 from ab initio calculations. ¢Gel can be
computed directly in all non ab initio calculations.
The term ¢Gtm of eq 79, being a difference between
two analogous terms related to the molecule in
solution and in the gas phase, may be simpler to
compute than the Gtm term in eq 78 because any
partial contributions cancel each other at a good
extent. This fact is often exploited when calculations
address solvation energy problems.
In definition 2 we have paid a bit more attention
to the thermodynamic definition of the reference
state. The reason the P¢V term was not inserted in
definition 1 will be mentioned later (section 3.2.6).
The second comment involves the internal geom-
etry of the solute. In passing from the gas phase to
the solution there is always some change in the
internal geometry, which often cannot be neglected
when the definition ¢Gel ) Gel - G0 is used. This
problem does not exist for definition 1 because no
reference is made to molecules in the gas phase in
eq 78. Semiempirical and classical models should
consider the energy term due to changes in the
equilibrium geometry, but this sometimes is not done.
The first three terms in the right-hand side of eq
78 can be computed all at the same time with a
unique QM calculation based on the appropriate
Schro¨dinger equation, or they can be computed
separately with semiclassical models. Gcav must be
computed separately; the same holds for Gtm. It is
worth remarking that contributions to the free energy
due to thermal motions of nuclei are separately
treated also in the QM descriptions of isolated
molecules; the same standard techniques (largely
simplified with respect to a full QM treatment, rarely
G ) Gel + Grep + Gdis + Gcav + Gtm (78)
¢Gsol + ¢Gel + Grep + Gdis + Gcav + ¢Gtm + P¢V
(79)
Quantum Mechanical Continuum Solvation Models Chemical Reviews, 2005, Vol. 105, No. 8 3027

The cavity formation energy, in the meaning used
by Uhlig in 1937,243 plays an important role in the
theory of hydrophobicity. We shall not review hydro-
phobicity theories here, because they encompass
continuum models, and we shall consider hydropho-
bic solvation only as a specific case of the more
general field of solvation phenomena. This remark
is not out of line here, because the majority of
available calculations of Gcav via computer simula-
tions246-270 have been motivated by hydrophobicity
studies, and this introduces a bias to the definition
of the simulation, inasmuch as the interests of
hydrophobicity are not coincident with those of
general solvation.
In fact, from simulations we would have informa-
tion about a large number of topics, necessary to
model methods for the calculation of Gcav for the
whole domain of solvation problems. The topics of
interest range from the chemical nature of the solvent
and of its physicochemical properties to the physical
macroscopic characteristics of the liquid (tempera-
ture, pressure, and also externally applied factors,
as electric fields and hydrodynamic flows), to the form
of the cavity that must faithfully reproduce the
extreme variety of sizes and shapes solutes have.
It would be extremely convenient to have a satis-
factory analysis of these factors before we begin to
develop Gcav computational methods. Parenthetically,
we remark that we need analyses and not tables of
specific values, because in the chemical applications
addressed by continuum models there is often the
need of examining the whole solvation free energy,
and hence Gcav, for a sizable number of solute
geometries. It would be impractical and in contradic-
tion with the basic philosophy of continuum models
to use computer intensive methods, as simulations
are, to compute a single element of the solvation
energy at each geometry of the solute. Simulations
are to be mainly used as validation tests of simpler
expressions describing Gcav and are better if they are
of analytical type.
Of course, this complete analysis has not been
performed, but the available simulations have been
of great help in assessing the relative merits of the
various methods we shall examine in the remainder
of this section.
There are at present a sizable number of Gcav
values from simulations for small cavities, in general
spherical but with a few examples of nonspherical
cavities,258,262 in various solvents. These latter simu-
lations are mostly due to Hofinger and Zerbetto,269-271
who have recently added to the half score of solvents
considered in the past 12 additional solvents, all of
chemical relevance. Some of these simulations also
give information about temperature and density
dependence. These results validate the guesses made
years ago (and discussed in the following pages)
about the general validity of some simple models,
dispelling doubts and criticisms and ending incon-
clusive discussions.
One may wonder why after more than 30 years of
computer simulations only now is there a satisfactory
situation (at least for small cavities). One reason is
that only in recent years has progress in computa-
tional machineries permitted a more efficient evalu-
ation of Gcav values. In 1998 Levy and Gallicchio272
remarked that the intensive computational demand
of simulations raised obstacles to the systematic
study of solvation in the preceding years. Now the
situation is greatly improved, and we are expecting
from simulations additional help, especially for larger
cavities and for liquids far from the standard ther-
modynamic state.
We cannot wait, however, and to proceed we have
to base our decisions on other criteria such as indirect
analyses of the large wealth of Gcav values computed
with simpler models.
The most extensive analyses performed thus far
are, to the best of our knowledge, those performed
by the Luque and Orozco group in Barcelona. We
refer interested readers to one273 of their papers as
it reflects the evolution of the cavity formation
methods over the past 10 years.
In our 1994 review1 we analyzed a variety of
methods; many have been abandoned and, in fact,
only two approaches have been considered in succes-
sive papers: those based on statistical mechanics and
those based on the area of the solute exposed to the
solvent. Here we shall consider both approaches.
3.2.1.2. Statistical Mechanics and Gcav Values.
The elaboration of PCM has strongly relied on
statistical methods to supplement the basic electro-
static code.274 For the cavity formation two ap-
proaches have been presented, one based on the
cavity surface and the other, which is the subject of
this section, based on a statistical mechanics ap-
proach, namely, the scaled particle theory (SPT),
introduced and elaborated over the years by Re-
iss.275,276 We may anticipate that the analytical
expressions of Gcav obtained from the SPT version
give values in good agreement with numerical simu-
lations.
This is not the appropriate place to fully explore
the great merits of the statistical mechanics methods
for our understanding of the properties of liquids. A
short comment is, however, necessary to put SPT,
and the use we are making of this theory, in the
appropriate context.
Statistical mechanics methods encountered dif-
ficulties in passing from gases to liquids, and to
overcome these difficulties very simplified models
were adopted. What is of interest here involves the
reduction of the fluid particles (atoms or molecules)
to hard spheres, neglecting the portion of the interac-
tion called the soft portion (i.e., the electrostatic and
dispersion terms which, in this review, we are more
interested in). Hard sphere models were extensively
used to study the various aspects of the theory for
fluid systems. Actually, the detailed and satisfactory
description of the properties of the hard sphere fluids
was found to be unable to provide numerically
acceptable estimates of the properties of real fluids.
The connection between the worlds of hard spheres
and of real liquids was given by the SP theory, a
model based on hard spheres again, but with radii
suitably modified with a scaling procedure to satisfy
some macroscopic experimental quantity. SPT is a
rigorous theory (in the sense that there are no
Quantum Mechanical Continuum Solvation Models Chemical Reviews, 2005, Vol. 105, No. 8 3029

drawn from experiments as well as from calculations:
The assumption that dispersion energy in more
complex systems can be described as a sum of atom-
atom terms has been expressed, and accepted as a
reasonable approximation, for many years (see ref
335 and references cited therein). Accepting this
approximation and indicating by m the atoms of the
solute and by s the atoms of the solvent molecules,
the dispersion energy can be written as a sum of
volume integrals, which in turn may be reduced to
surface integrals with the aid of appropriate vector
functions ABdis,ms
(n) (rbms) and reduced to a finite summa-
tion of elements defined for each tessera of the cavity
surface. The procedure exactly parallels that in use
for the repulsion term, and there is no need to repeat
what has been summarized in the preceding section
and in the 1994 review.1 Here, it is sufficient to say
that the form of the dispersion contribution to the
solvation free energy is exactly equivalent to that
reported in eq 96 for the repulsion interaction, but
this time the vector functions ABdis,ms
(n) (rbms) become
The serious limits of the use of atom-atom terms for
the description of interactions within small clusters
composed of polyatomic molecules are well-known.
In passing to solutions the constraints imposed by
the model (a fixed cavity, a continuous distribution
of the solvent molecules) a part of the artifacts
specific to the atom-atom model are eliminated.
However, the final values of dispersion and repulsion
energies depend on the values of the coefficients used
in the potential. The availability of QM versions for
Gdis (as well as for Grep) makes it possible to define a
computational strategy to get independent values for
these coefficients.336
3.2.4. Charge Transfer Term
Among the various components of the intermolecu-
lar interaction energy in small clusters we have
mentioned in section 3.1, there was a charge transfer
term that was not included in the decomposition of
the solvation energy we presented in section 3.2.
Indeed, charge transfer (CT) is an important term
in the decomposition of the variational definition of
the interaction energy ¢E, even if it is not included
in the standard definition of this quantity performed
with the aid of perturbation theory.
We have not put a CT term in the decomposition
of ¢Gsol because CT effects are quite specific and
directional and, thus, not appropriate for description
by continuum models that are essentially based on
an averaged distribution of solvent molecules. How-
ever, ionic and some strongly polar solutes tend to
accumulate electronic charge from the media or
distribute it to the media, irrespective of the local
orientation of the nearby solvent molecules. Envisag-
ing processes of just this type, Gogonea and Merz337,338
have proposed an addition to QM ASC methods
including CT effects. The procedure starts from a
standard continuum BEM calculation giving the
apparent charge distribution, called by the authors
óRF to signify that it depends on the reaction field
(RF). The CT process is mimicked by introducing an
additional CT charge, óCT, spread over the same
cavity surface. This charge is related to a new concept
introduced in the QM formulation of the theory, the
fractional electron pair (FEP). The FEP is added to
the LUMO, or deleted from the HOMO of the solute,
with the continuum medium acting as a generator
or a sink of electronic charge, respectively.
The whole calculation consists of an iterative
scheme which can be schematized in three steps.
First, the dielectric is charged by distributing óCT
as a charge density onto the BEM grid. Then RF
charge densities are obtained from the BEM equa-
tions
where óRF and óCT are vector quantities containing
the set of RF and CT charge densities, respectively,
and 0 and i are the solvent and solute dielectric
constants. In eq 106 T is a vector having elements
that are normal components of the electric field
generated by the solute charge density on the bound-
ary elements, and the Kij terms depend only on the
shape of the dielectric interface and have the form
where i and j are boundary elements and ai is the
corresponding area.
Finally, the RF and CT charge densities are fed
into the Hamiltonian, and the solute is (dis)charged
by adjusting the density matrix during the SCF
procedure (FEP creation) such that the total charge
on the solute is less qCT
where k is either the HOMO or the LUMO orbital,
which is determined by the direction of charge flow,
and qCT is the pointlike equivalent of the óCT at each
boundary element, with
The whole procedure has to be iteratively repeated
to get convergence on the SCF energy. The final
expression of ¢Gsol has an additional term, called GCT,
whereas the other nonelectrostatic terms are de-
scribed with the standard procedures we have docu-
mented in this section.
Vdis(rms) ) ∑
n)6,8,10
-
dms
(n)
rms
n
(104)
ABdis,ms ) -
dms
(6)
3rms
6
rbms (105)
(2ðI - i - 0i + 0 K)óRF ) i - 0i + 0[T - 10(2ðI - K)óCT]
(106)
Kij ) sai
(rbi - rbj)ânˆ
jrbi - rbjj
3
dai (107)
Píî ) 2 ∑
i
occ
cíicîi + qCTcíkcîk (108)
qCT ){<0 charge flow: solute f dielectric (k)HOMO)>0 charge flow: dielectric f solute (k)LUMO)
3036 Chemical Reviews, 2005, Vol. 105, No. 8 Tomasi et al.

lent ionic compounds, the same authors347 suggested
that the solute cavity be reduced with respect to
neutral species: the optimum cavity for the hydration
of ions was defined by scaling the atomic radii by a
factor of 1.10-1.15. This reduction was justified by
the fact that due to strong electrostatic interactions
solvent molecules are in general closest to the atoms
of charged molecules rather than those of the neutral
ones. Very recently, a further refinement has been
introduced;348 instead of reducing the whole cavity
of the solute, which might be appropriate for small
size ions, the scaling down is now limited to the heavy
atom (and the hydrogen atoms attached to it) that
bears the formal charge.
In 1997 Barone et al.349 presented a set of rules
for determining the atomic radii of spheres used by
GEPOL to build the molecular cavities within the
DPCM solvation model. The procedure was applied
to compute the hydration free energy for molecules
containing H, C, N, O, F, P, S, Cl, Br, and I at the
Hartree-Fock computational level with a medium
size basis set [6-31G(d) with eventual diffuse func-
tions for anionic species]. This new definition of the
cavity, which was referred to as the united atom
model for Hartree-Fock (UAHF), was based on
chemical considerations, basically hybridization, for-
mal charge, and first-neighbor inductive effect. The
UAHF parameters were obtained by imposing two
initial conditions; the hydrogens do not have indi-
vidual spheres, but they are included in the spheres
of the heavy atoms to which they are bonded (united
atom model), and the elements of each periodic table
row have the same ‘‘basic” radius, modified by the
molecular environment. As previously suggested by
Orozco and Luque, ions present atomic radii different
from (and smaller than) those of the corresponding
neutral molecules.
Due to its implementation in the Gaussian code,17
the UAHF method has been used in many studies;
it is, however, important to note here that this
method works well when all of the conditions used
in the original fitting procedure are maintained,
namely, the use of DPCM solvation model and the
HF/6-31G(d) level of calculation. Extensions to other
more recent versions of PCM (CPCM or IEFPCM) or
to other QM levels (for example, DFT), or only to
other, larger, basis sets, in fact, cannot reproduce the
good results of the original version. It is also impor-
tant to note here that, due to this limitation, the new
version of Gaussian code (G03) no longer considers
UAHF cavities as the recommended (or default)
cavities for PCM (which, in turn, is no longer DPCM
as default but IEFPCM) and, thus, careful attention
has to be paid by all users.
We terminate this brief excursion on different
proposals to define cavity radii by quoting a paper
by Dupuis and co-workers350 in which they describe
a new protocol to define the radii of the spheres to
be used for oxoanions, XOmn-. These radii are deter-
mined by simple empirically based expressions in-
volving effective atomic charges of the solute atoms
that fit the solute molecular electrostatic potential
and a bond-length-dependent factor to account for
atomic size and hybridization.
3.2.6. Contributions to the Solvation Free Energy Due to
Thermal Motions of the Solute
We have so far considered the contributions to the
solvation energy that may be computed at fixed
geometry of the solute (BO approximation). These
calculations can be extended to the whole space of
the nuclear conformations, a 3N-6 space, that will
be indicated by {R}. The energy G(R) (definition 1,
see eq 78) defines the PES surface of the focused
system in solution. We remark that to act as PES,
the ¢Gsol(R) function (definition 2, eq 79) must be
supplemented by an E0(R) function, which has to be
independently computed.
The G(R) function replaces the corresponding E(R)
function defined for the isolated system. The differ-
ence is in the physical nature of the two surfaces (a
component of the free energy versus an internal
energy) and not in their use. We remark, however,
that for the study of some phenomena, it is conve-
nient to expand the {R} space to an {RxS} space
including dynamical solvent coordinates.351-357 The
definition of “solvent coordinates” is indeed a very
important subject to which continuum solvation
methods have been successfully applied. In this
review, we have decided not to consider this topic,
not because of less interest but just because the
subjects to treat are so many that we had necessarily
to make a selection. Therefore, here, we shall not
consider dynamical couplings of solute and solvent
coordinates and we shall limit ourselves to analyze
G(R).
In analogy with what is done for the E(R), each
value of G(R) has to be supplemented with additional
contributions due to the internal and external mo-
tions of the nuclei. We are here only interested in
the contributions permitting the calculation of the
free energy. For E(R) these contributions depend on
the temperature; for G(R) the temperature has been
already fixed in the definition of the PES (the
dielectric constant  as well other solvent properties
used in the calculations depend on T). The calculation
of the free energy is often limited to the critical points
of the PES corresponding to local minima and to
transition states (saddle points of the first type: SP1)
for which the calculation of the internal motions
contributions (i.e., vibrations) is computationally
simpler and similar in vacuo and in solution.
For the contributions due to the external motions
(translation and rotation of the whole molecule), the
procedure in solution is not so simple as for molecules
in the gas phase. It is convenient to make explicit
reference to Ben-Naim’s definition of the solvation
process.358 According to Ben-Naim, the basic quantity
is the free energy change ¢Gsol
/ for transferring the
solute from a fixed position in the ideal gas phase to
a fixed position in the liquid phase. The constraint
of fixed position in the liquid is released later, with
the addition of a further contribution to the energy,
called “liberation free energy”, related to the momen-
tum partition function ⁄M
3 and to the number den-
sity FM. The explicit formula is FM⁄M
3 . An analogous
release of the fixed position is introduced for M in
3038 Chemical Reviews, 2005, Vol. 105, No. 8 Tomasi et al.

this model, the value of the dielectric constant in the
layer, and the thickness of the layer, which have been
numerically optimized in the context of the cited
studies.399-401 In the dual solvent region model, the
dielectric constant of the layer corresponds, in a good
approximation, to the optical dielectric constant of
the solvent and its thickness to the solvent molecular
diameter. The numerical optimization leads to a
description of the phenomenon acceptable on intui-
tive grounds; the first solvent layer is immobilized
and dipole-oriented, contributing to the polarization
with its electronic term only. Of course, there is
compensation of effects, but these findings reinforced
the opinion, already expressed in the literature, that
pure continuum models could exploit the difference
between the electronic and total permittivity to
introduce effects due to electric saturation.
Considerations of this type have led Basilevsky and
co-workers to formulate a layered model without
constraints on the symmetry of the two surfaces (it
is based on an ASC approach).411-413 The approach,
called the frequency-resolved cavity model (FRCM),
is based on the use of the electronic permittivity in
the first zone and on the static value in the outer
space; more details on this method will be given in
section 5.3.
4.3.3. Molecular Cluster Models
A third way to describe dielectric saturation effects
around ions consists of replacing the continuous
dielectric in the vicinity of the ion with solvent
molecules. We are here at the border between two
approaches, the continuum approach, which adds
some discrete molecules in critical positions, and the
discrete model, which replaces distant molecules with
a continuous distribution. We shall return to this
subject, with a discussion of where to put the border
between the two approaches (see section 7). Here we
are in the continuum domain with the addition of
discrete elements.
Attempts at inserting some discrete solvent mol-
ecules in the Born-like description initiated very early
to try to eliminate the “failures” of this model; we
report here a selection of some old references.390,414-419
By comparing these references with those reported
in section 4.3.1 (more recent), it appears that the
main attention has shifted from solvation energies
to the structure of the liquid around the ion. These
references are also testimony to the beginning of the
use of QM calculations in the study of condensed
systems. Actually, the discrete solvent molecules used
in these, and in other successive papers, were mod-
eled at a level that was practically unable to describe
the electronic solute polarization considered to be
essential in the models of the preceding subsections.
We have already described several shortcomings
of the use of semiclassical discrete solvent models for
simulations. In the present case, the need for using
simplified models is less justified, because the num-
ber of molecules is small.
The models to be used here should combine the
features necessary to describe molecules in a liquid,
which are different from those of molecules interact-
ing in the gas phase, and the features specific for
molecules feeling the strong field of the ion. An
attempt to formulate a procedure able to introduce
both features in semiclassical models was elaborated
by Floris et al.420-424 According to this procedure, ab
initio effective pair potentials can be obtained using
a continuum solvation model (the PCM in the cited
papers). In this way, many-body effects are included
in the ion-solvent potential at the same time keeping
the computational convenience and simplicity of two-
body functions in the development process. These
potentials are then used in fully discrete computer
simulations. Molecular dynamics in particular is
appropriate for the study of the dynamical features
of the system, such as the persistence time of a
solvent molecule in the first or second solvation shell,
the determination of the diffusion coefficients, the
spectra of the solvation shell, and the dynamics of
exchange of solvent molecules of different nature in
mixing solvents. All of the aspects we have here
mentioned have been studied for a variety of cations
in the above-cited papers, to which we add here a
more recent one on Cd(II) in aqueous solution.425
Other semiclassical potentials including polariza-
tion have been employed. We cite among them the
effective fragment potential (EFP) developed by Gor-
don and collaborators,426 which has given good results
in modeling solvent molecules, especially for neutral
solute (see also section 7.3.2). The EFP description
presents shortcomings when there is remarkable
charge transfer such as in anions and in some
cations.427,428
The use of a full QM description of the cluster,
however, is advisable when the property under
examination is suspected to be sensitive to layering
effects. It is not yet fully ascertained when an
accurate and well-balanced description of the solva-
tion cluster is compulsory, but surely it is convenient
to have available descriptions of this type. An ap-
preciable number of studies have been performed
on isolated metal-solvent clusters of the type
[MSn]k+,429-433 but often, when the number n of
solvent molecules reaches or exceeds the number
necessary for the first solvation shell, the geometry
does not correspond to that found in solution. The
concept of a “hydrated ion” as a new species formed
in solution by the ion and a limited number of water
molecules is now widely accepted.434 A simplified
strategy is often employed to study these systems.
It consists of performing a geometry optimization of
the hydrated ion in the gas phase followed by a
refinement of the geometry with the hydrated ion
within the continuum solvent.435,436 It is noteworthy
that these geometry optimizations of the clusters are
sensitive to details in the continuum solvation model.
This fact was pointed out by Sanchez Marcos and co-
workers,436 who noticed that there was ambiguity in
the results according to the procedure used. In
particular, it was shown that the incorrect lengthen-
ing of the metal-oxygen distance found in some
optimizations with continuum models was due to the
spherical shape of the cavity used in those models;
with a molecular shaped cavity the correct shortening
is in fact recovered.
3042 Chemical Reviews, 2005, Vol. 105, No. 8 Tomasi et al.

4.4. Nonuniformities around Neutral Molecules
Some of the early attempts to model dielectric
saturation for ions in solution were also extended to
small neutral molecules, reduced to a point dipole.
However, the electric field produced by a dipole is
decidedly smaller than that produced by a localized
charge, and thus the saturation effect on the solva-
tion energy is far smaller; this induced the consid-
eration of ions only. More recently, the progress in
computational methods and the availability of data
on time-dependent phenomena have renewed the
interest in larger and more complex molecular sys-
tems. For example, the method of Basilevsky and co-
workers411-413 considered in section 4.3.2 is not
limited to atomic ions, but it can be extended to
molecular systems, charged or uncharged.
A method addressing equilibrium solvation ener-
gies of polyatomic molecules has been recently pro-
posed by Edholm et al.437 The description of the
dielectric nonlinearity is based on a modification of
the Langevin-Debye model,438 with a finite differ-
ence method applied to calculate the spatial depen-
dence of the electrostatic potential using a quite fine
grid borrowed from DFT methods and limited to a
layer of 1.2 nm (the solvent outside this layer is
described with a linear continuum). The method also
contains the calculation of the nonelectrostatic terms,
van der Waals and cavity formation, according to the
schemes discussed in section 3. The solute is de-
scribed at fixed geometry, not influenced by the
solvent, and with atomic charges drawn from the
OPLS-AA force field.439,440 In addition, solvent polar-
ization effects are neglected. The hydration free
energy values for around 200 neutral organic com-
pounds with their decomposition into electrostatic,
van der Waals, cavity formation, and nonlinear
electrostatic terms are reported. The nonlinear con-
tribution is decidedly smaller that the cavity forma-
tion and van der Waals terms, but comparable to or
even larger than the sum of these two terms, which,
in water, compensate each other to a good extent. The
authors remark that the extensive parametrization
used in the Cramer and Truhlar SMx models (see
section 2.3.3) should implicitly include nonlinear
effects and that PCM (see section 2.3.1) uses adjusted
cavity radii, which will quench the saturation effects.
In our 1994 review, we explained the origin of the
scaling factor, applied both to polar and nonpolar
molecules, and we confirm here that neither dielectric
saturation nor electrostriction (the opposite term
related to nonlinearity) have been considered in
fixing this scale factor, but we accept the remark
about the quenching.
Sandberg and Edholm have expressed two other
formulations of the nonlinear response solvation. One
already cited paper398 involves spherical ions and
gives an analytical formulation exploiting the sim-
plifications made possible by the symmetry of the
problem. The other441 uses for the charge distribution
a dipole, derived from the GROMACS force field442
and specialized for amino acids. These methods seem
to be still in evolution, and it is convenient to await
future developments to assess their validity to com-
pute solvation energies.
4.5. Nonuniformities around Systems of Larger
Size
In extending our survey to larger molecules we
enter into a field dominated by computations on
biomolecules. The solvation methods we have de-
scribed also apply to biomolecules, and we do not
intend to open here a chapter dedicated to this very
important domain of computational chemistry, for
which there is a wealth of specialized reports. Our
attention will be limited to the very specific aspect
of the description of nonuniformities arising from the
large size of biomolecules.
There is a long tradition of using molecular me-
chanics (MM) methods to describe large molecules.
Molecular mechanics, we remind, is based on empiri-
cal force fields in which classical terms related to
intermolecular energy deformations are accompanied
by site-site non-covalent interactions of van der
Waals and electrostatic type. These last are ex-
pressed in terms of point charge interactions, which
are not screened because electronic polarization is
missing in standard MM models. This defect of the
model has been empirically corrected by introducing
in the calculation of Coulombic interactions a dielec-
tric constant, with values ranging from 2 to 20 (the
definition of a proper value of  to be used with
biomolecules has given rise to considerable contro-
versy).443-446 The fact that a “constant” value of  was
not sufficient, at least for groups bearing a net
charge, was soon noticed. In 1923 Bjerrum447 pro-
posed the use of a distance-dependent dielectric
function to describe the ionization constants of bi-
functional organic acids. Bjerrum’s proposal was later
elaborated by Kirkwood and Westheimer448 and
thereafter amply used, especially for the calculation
of pKa of titrable groups in proteins. Numerous
checks have been done, and several modifications
have been proposed. There is in the literature a
sizable number of variants of a linear form of the
distance-dependent dielectric function: (r) ) kr, with
k determined in different ways.
However, a linear relationship of this kind is not
appropriate to describe the saturation effects deter-
mining the nonlinearity (see section 4.1). Many
proposals have been suggested449-462 to take into
account the physical origin of the problem using one
among the various models we have reviewed, or
defining, through other ways, (r) functions with a
sigmoidal shape.
Here we cannot discuss all of these proposals, but
we can point out some authors whose work has had
a larger impact on the literature. Warshel et al.463
used experimental pKa values and redox potentials
to derive a distance-dependent dielectric function,
and Mehler and Eichele464,465 introduced a sigmoidal
distance-dependent dielectric permittivity function to
calculate the electrostatic effects in water accessible
regions of proteins. For metal ion-nucleic acid in-
teractions, Hingerty et al.466 proposed a modification
of Debye’s distance-dependent dielectric function,
which is sigmoidal and rises rapidly with distance.
An alternative mathematical expression for this
function was proposed by Ramstein and Lavery.467
Guarnieri et al.468 proposed an extension of Mehler
Quantum Mechanical Continuum Solvation Models Chemical Reviews, 2005, Vol. 105, No. 8 3043

and Eichele’s work, in which the position dependence
of the screening is partly accounted for by using
different values of a parameter for buried and for
exposed atoms. More recently, in an attempt at a
more generally applicable function, Sandberg and
Edholm469 proposed a modification of the function of
Warshel et al., with parameters depending on the
mean distance of the atoms from the protein surface.
Hassan et al.,470 following a formulation originally
developed by Bucher and Porter471 and subsequently
by Ehrenson,472 obtained a reaction-field-corrected
form of the radial-dependent permittivity profile.
In our rapid exposition we have simplified several
aspects of the methodological exposition of the use
of “screened” Coulomb potentials (for example, no
mention has been made of the use of models without
and with solvent cavity). Interested readers are
referred to the above-cited papers to which we add
here an older but rich review by Mehler,473 but
further comments as well as additional references
will be reported in section 7.1.
4.6. Systems with Phase Separation
We consider another class of systems to which the
concepts related to nonuniformities in the local
properties of the liquid can be applied, namely,
systems with phase separation boundaries. This is a
very large class of systems, collecting elements of
quite disparate nature. We anticipate that in some
cases the explicit consideration of local deviations of
the liquid properties with respect to the bulk is
compelling, whereas in others the description of these
changes is optional, in the sense that they need to
be considered only when one aims at a more accurate
description of fine details not given by simpler
descriptions of the medium.
Local changes involve several aspects of the liquid,
the density, the local concentration of dissolved
species (an aspect of particular relevance in ionic
solutions), the viscosity, and other parameters related
to the motion of solvent molecules, and, of course, the
dielectric function to which our attention will be here
limited. Before entering into details, we need to
introduce additional characterizing elements in our
set of systems, thus far defined only in terms of the
presence of phase boundaries. In this way, in fact,
we can obtain subclasses that are more homogeneous
with respect to the methodological tools to use. We
shall limit ourselves to a rough level of classification,
finer classification being necessary only for more
specific reviews.
An important feature, useful for our attempt of
classification, is the length scale for the boundaries.
This feature has been proposed to define ‘‘complex
liquids” in which the appropriate scale length is in
the mesoscopic range.474 This definition of complex
liquids collects, for example, colloidal suspensions,
micellar solutions, microemulsions, closed bilayers,
foams, gels, and other variants of the basic types we
have here mentioned as bicontinuous liquids. We
shall consider the complex liquids no further in this
review, but we note that a description of complex
liquids coherent with continuum models can be given,
even if not in a straightforward way.
The local effects due to the surface play a far more
important role in systems with multiple mesoscopic
phase boundaries than in systems with a single phase
boundary, and in general there will be the need of
extending the model to include mutual interactions
among the mesoscopic domains, an extension that
was not yet considered in the QM continuum models.
We have here mentioned multiple mesoscopic sys-
tems because in our opinion they represent one of
the most important perspectives of development of
continuum models
Other subclasses with macro- or mesoscopic length
scale for the boundary may be characterized by the
number of phases, by their geometry, and by the
physical and chemical nature of the components.
Some among them (liquid surfaces, pores, thin films,
etc.) are the subject of intensive studies, deserving
separate reviews.
Let us here consider some examples of two-phase
systems, namely, systems with a liquid/gas, liquid/
liquid, or liquid/solid surface. The continuum model
we have examined in the preceding sections can be
applied to these systems with minor changes. We
thus return to models in which the focused compo-
nent continues to be a single molecule or a small
cluster, whereas the remainder has a more complex
definition, possibly including local modifications near
the boundary (for the definition of focused models see
section 1).
The scientific and practical interest of two-phase
systems has led to a detailed experimental charac-
terization of at least some among them, and it has
spurred a considerable number of attempts of theo-
retical modeling also including continuum models.
These models exhibit specific features as the char-
acteristics of two-layer models may be quite different.
For example, the characteristics to model in processes
occurring at the electrode surface are remarkably
different from those necessary in the study of the
effects of a metal on the spectroscopic properties of a
chromophore, and both have little to do with those
required to model processes and equilibria in aerosols
of environmental interest. For this reason we limit
ourselves to a few selected examples, covering the
three types of two-phase systems mentioned above
and bearing a strict connection with the methods
presented in sections 2 and 3.
As remarked above, in some cases the explicit use
of local corrections to the bulk properties is not
compulsory, and we start to examine some methods
of this simpler type.
An ASC method, based on PCM and aiming at
studying systems of this type, has been elaborated
years ago by Sakurai and co-workers in Japan.475-477
The description of the medium is given in terms of
nonoverlapping portions of a homogeneous con-
tinuum with different values of the dielectric constant
within each portion. The code was later improved219,478
and used for a number of studies, mostly of organic-
biologic interest, from which the versatility of the
approach may be appreciated.479-483
Among the successive variants of Sakurai’s proce-
dures, one deserves explicit mention in this re-
view.484,485 It involves a QM/MM method in which the
3044 Chemical Reviews, 2005, Vol. 105, No. 8 Tomasi et al.

feature of separate dielectric domains is exploited to
allow instantaneous differential polarization in the
protein fragments. In this model, called the polariz-
able mosaic model (PMM), every covalent bond in the
MM region is replaced by a cylindrical stick made of
dielectric, the permittivity of which depends on the
nature of the original bond. As a consequence, a
protein molecule is represented by a mosaic consist-
ing of cylinders with various polarizabilities. This
model offers a perspective of solution to two kinds of
problems. First, other bodies can be described in a
first approximation as a whole region, introducing
then local anisotropies with a similar proceduresas,
for example, a relatively homogeneous membrane
with local elements involved in some local processes.
Second, the introduction of mutual polarization of the
classical elements in QM/MM simulations is a com-
putational challenge, and this approach is a possible
way of solving this problem.
The first version of PCM, namely DPCM, was also
used to study problems characterized by the presence
of a boundary. The first paper, published in 1986,486
presented a model addressing energy changes in long
DNA fragments subjected to large-scale deformations
and to local opening of the double helix. The problems
posed by this model are representative of those
occurring in “complex liquids” and in related systems.
Around the double helix (a polyion) there is a portion
of solvent showing large nonlinear dielectric effects,
combined with partial counterion condensation, and
both effects change when the double helix is bent or
open. In this model, the local nonlinear dielectric
effects were described in terms of the layered model
described in section 4.3.2. Another application was
on the effects of partial solvation in modifying the
molecular electrostatic potential of a substrate placed
in the cleft of the active site of a model enzyme.487
Other studies involved polar solutes (with hydropho-
bic tail) placed near the flat separation surface of two
immiscible liquids (or at a liquid/vacuum separa-
tion).488-491
In all of these DPCM studies, use was made of a
tessellation of the boundary (of planar infinite,
spherical, or irregular shape, according to the case).
More recently, the same systems were reformulated
within the IEF version of PCM.80 Within this formal-
ism, by introducing novel Green-type operators of
analytical form, tessellations of the infinite plane (or
the spherical boundary) separating the two phases
are no longer required, and thus the dimension of the
problem is identical to that of a homogeneous and
isotropic solution. We note that full analytical Green
operators are relatively easy to obtain for surfaces
with a regular shape but not for more complex
shapes. In these cases it is thus useful to introduce
Green operators of semianalytical form (requiring a
partial numerical integration) or to combine the IEF
formalism with the original BEM formulation based
on a tessellation of the boundary surface.
The use of the IEFPCM formalism also allows us
to introduce an (r) function for the liquid portion.
The present version of the computational code can
deal with (r) functions limited to a single spatial
variable (in practice, all of the models considered in
section 4.3.1); these kinds of functions, for example,
are well suited for planar surfaces and spheres. The
extension to some other cases, such as cylindrical
surfaces or liquid with a simple surface combined
with a tensorial permittivity, is immediate, but not
tested yet. The extension to more anisotropic per-
mittivities exhibiting a dependence on all three
spatial coordinates would require a more intensive
computational elaboration, but it is surely feasible.
This model, which eliminates all of the divergences
occurring in simpler image charge models with a
sharp surface, has thus far been applied to the
crossing of liquid/vapor and liquid/liquid surfaces by
neutral and charged molecules, namely, a two-surface
model mimicking a membrane has been studied,80 as
well as the changes of the polarizability of the halide
anions passing from water to air.492 In both cases the
function, (r), has been modeled on the corresponding
density profile, F(r), derived from simulations, as-
suming proportionality between the two functions.
Another continuum electrostatic model to treat
solute at interfaces has been elaborated by Benjamin,
especially to compute the spectral line shape at liquid
interfaces and surfaces.493 In this model based on the
image charge approach, the interface is described as
a sharp boundary separating two bulk media, each
characterized by different optical and static dielectric
constants, and the chromophore is modeled as a pair
of spherical cavities in which two equal and opposite
point charges are embedded.
As a matter of fact, Benjamin’s contribution to the
knowledge of interfaces and their modelization goes
well beyond this model.493-499 His main activity has
in fact involved molecular dynamics (MD) simula-
tions and has been focused on the study of electronic
properties of chromophores as a function of their
location at the interface.
Let us now pass to another kind of system with
boundaries, namely, that represented by lipid bilayer
membranes. These systems have been treated within
the framework of continuum electrostatics by con-
sidering a slab of low dielectric material intended to
represent the hydrocarbon tails of the lipid molecules,
surrounded by high dielectric material intended to
represent the headgroups and aqueous solution. As
a matter of fact, a solvated lipid bilayer considered
at the nanometer length scale is an example of a
system nonuniform in one dimension, that is, along
the axis normal to the bilayer. In general, the
dielectric permittivity for such a system will be a
function of the position along one axis. Furthermore,
it will not be a scalar, as has been assumed in
virtually all applications of continuum electrostatics
to molecular systems, but will be a second-rank
tensor.
Several authors have presented methods for cal-
culating the local dielectric constant for a nonuniform
system. King, Lee, and Warshel500 and Simonson and
co-workers450,501 have presented methods for deter-
mining the dielectric constant of a spherically sym-
metric but nonuniform system and applied them to
proteins solvated in water. These methods are based
on a partition of the system into two or more
concentric spherical shells. Once an estimate of the
Quantum Mechanical Continuum Solvation Models Chemical Reviews, 2005, Vol. 105, No. 8 3045

dielectric constant in the outer shells is made, the
dielectric constant in the inner shell can be calculated
either from fluctuations of the total dipole moment
for that region or by direct observation of the polar-
ization response to a uniform field. A dielectric profile
for a lipid bilayer in solution was presented earlier
by Zhou and Schulten.502 Their approach was to
divide the membrane into slices and compute the
permittivity for each slice with the usual expression
for a uniform, isotropic system. The implicit assump-
tion is that the slices do not interact, which seems
hard to justify. More recently, Stern and Feller503
presented a method for calculating the static dielec-
tric permittivity profile for a system nonlinear in one
dimension under the periodic boundary conditions
used in computer simulations.
As a last example of systems with phase separa-
tion, we consider liquid/metal systems with particular
attention to the surface enhancing (SE) effects of a
metal on the photophysical properties of a molecular
chromophore placed in the vicinity of the surface.
Some of these studies will be examined in more detail
in sections 5 and 6, and here we limit ourselves to
presenting some methodological aspects related to
the theme of this section.
The two models we consider describe the metal on
the same footing as the solvent in homogeneous
solutions, in the sense that the metal contributions
to the proprieties of the focused part (the chro-
mophore) are treated in terms of a continuous
response function without explicit consideration of
the atomistic structure of the metal itself. The
definition of the dielectric function for a metal speci-
men is more complex than for a liquid, and thus
specific modifications have to be introduced.
In the first model79,504-509 we cite (which belongs
to the family of PCM methods), the metal is consid-
ered to be a perfect conductor for time-independent
external electric fields. To study photophysical prop-
erties, however, one needs to introduce a frequency-
dependent electrical perturbation and thus, in this
case, the metal is considered as a dielectric with a
frequency-dependent dielectric function. The starting
point to model this basic quantity may be given by
the Drude theory presented in all textbooks of solid
state physics,510 which gives an expression of the
dielectric function suitable for the bulk of the metal.
The metal specimens used in SE measurements,
however, often have a finite size, and are all subjected
to other sources of electric potential, static as well
as frequency dependent, related to the presence of a
chromophore of molecular size. The finite size of
metal specimens first affects the relaxation time,
which surely has not the value of the bulk used in
the Drude theory. To account for these effects, a
correction based on the Fermi velocity of the electrons
has been introduced in the model. In addition, the
dielectric response of the electrons in the metal has
a nonlocal character, in the sense that the polariza-
tion vector at a given point depends on the value of
the electric field over all other points. This nonlocal-
ity, completely neglected in the original Drude for-
mula, has been partially described in the model506
making use of a hydrodynamic model511 to introduce
a correction in the Drude expression, or, more ac-
curately,508 using the Lindhard theory512 with Mer-
min corrections for the finite relaxation time.513
Another correction introduced in the model involves
the effect of surface roughness, found to be important
to realistically describe the numerical output of some
experiments. To this end a modified version of the
Green function proposed by Rahman and Maradu-
din514 has been employed.508
The second model515-517 we cite assumes that the
molecular system is enclosed in a half-spherical
cavity embedded in a linear, homogeneous, isotropic
dielectric medium and adsorbed on the surface of a
perfect conductor. The effects of the media on the
properties of the molecular system are described
using the method of image charges. According to this
method, a charge located in the cavity induces three
charges in the surrounding dielectric environment:
one charge in the metal surface, due to the perfect
conductor/vacuum interface, a second charge due to
the half-spherical environment between the vacuum
and the dielectric medium, and a third charge that
is the image in the metal due to the interactions
between the perfect conductor and the induced
charges in the solvent.
5. Nonequilibrium in Time-Dependent Solvation
In the sections 3 and 4 we have presented exten-
sions of the basic model but with the common feature
of retaining the condition of equilibrium between
solute and solvent. In this section we instead review
the physics we have to change and/or add to this basic
model to properly describe time-dependent (TD), or
nonequilibrium, processes.
TD processes embrace a large variety of phenom-
ena, which span an impressive range of time scales
(from 10-15 s for linear and nonlinear spectroscopy
to 10-2 s for diffusion-controlled reactions in dilute
aqueous solution). Obviously this broad range in-
volves many physical aspects of solvation. Here, we
limit ourselves to consider only those related to the
most recent progress in the description of the dy-
namic effects in QM solvation models.
Rigorously, the dynamic evolution of the molecular
solute should require that we postulate the proper
TD QM equations, which generalize the time-inde-
pendent effective Schro¨dinger equation introduced in
section 2.4 for the “static” basic model.
As a matter of fact, not all dynamic processes
require an explicit TD description. For example, if
the time dependence of the solute-solvent interac-
tions is sufficiently slow (as in the case of diffusive
relaxation processes of the slow degrees of freedom
of the solvent), it is possible to apply the adiabatic
approximation and the state of the solute at time t
can be obtained as the eigenstate of H(t) using time-
independent formalism. In contrast, the dynamic QM
equation for the solute must be explicitly introduced
in describing its response properties to TD electro-
magnetic perturbations. This aspect will be consid-
ered in section 6 devoted to solvent effects on solute
response properties.
In the presence of a time evolution of the solute
we have also to extend the description of the solvent
3046 Chemical Reviews, 2005, Vol. 105, No. 8 Tomasi et al.

specific analysis when in the presence of a polarizable
solvent. In this case, in fact, it becomes fundamental
to properly account for the dynamic response of the
solvent.
This class of phenomenon is of very general occur-
rence; here, in particular, we shall focus on the
vertical electronic transitions, but electron and en-
ergy transfers, or nuclear vibrations, are also ex-
amples of the same class.
Historically, in continuum models large use has
been made of the approximation according to which
it is sufficient to decompose the response function
into two terms. Within this approximation, also
known as the Pekar or Marcus partition, the polar-
ization vector PB described in the previous section
becomes
where fast indicates the part of the solvent response
that always follows the dynamics of the process and
slow refers to the remaining inertial term. Such
splitting in the medium response gives rise to the so-
called “nonequilibrium” regime. Obviously, what is
fast and what is slow depends on the specific dynamic
phenomenon under study.
In a very fast process such as the vertical transition
leading to a change of the solute electronic state via
photon absorption or emission, Pfast can be reduced
to the term related to the response of the solvent
electrons, whereas Pslow collects all of the other terms
related to the various nuclear degrees of freedom of
the solvent.
The operative partition of the total polarization can
be performed using two alternative, but equivalent,
schemes, which we prefer to indicate without using
Marcus and Pekar names because in the literature
there is an annoying confusion about these names.
One scheme, which we shall call partition I, was
proposed by Marcus in 1956 in his landmark paper521
on the nonequilibrium polarization in a classical
continuum solvation framework. The Marcus parti-
tion has been extended to QM continuum models by
several authors78,137,213,221,238,239,522,523 to properly de-
scribe excited/de-excited electronic states after a
vertical (Franck-Condon) transition. In this scheme,
eq 123 is rewritten as
where superscripts E and A refer to the “early” (or
initial) and the “actual” time in relation to the
chronological order of the transition process. In this
partition, the previous slow and fast indices are
replaced by the subscripts or and el referring to
“orientational” and “electronic” polarization response
of the solvent, respectively.
In the second scheme, which we shall indicate as
partition II, eq 123 has to be rewritten as
where the subscripts in and dyn refer now to an
“inertial” and a “dynamic” polarization response of
the solvent, respectively. The use of eq 125 within
the continuum solvation methods has been applied
both to classical and to QM approaches.56,214,216,219,524
The differences between the two schemes are
related to the fact that, in partition I, the division
into slow and fast contributions is done in terms of
physical degrees of freedom (namely, those of the
solvent nuclei and those of the solvent electrons),
whereas in partition II, the concept of dynamic and
inertial response is exploited. This formal difference
is reflected in the operative equations determining
the two contributions to PB as, in II, the slow term
(Pin) includes not only the contributions due to the
slow degrees of freedom but also the part of the fast
component that is in equilibrium with the slow
polarization, whereas, in I, the latter component is
contained in the fast term Pel (see Table 2 in which
we compare the electrostatic equations identifying
the two schemes).
Table 2. Basic Equations within Partitions I and IIa,b
partition I partition II
PBA ) PBor
E + PBel
A PBA ) PBin
E + PBdyn
A
PBor
E )
łor
ł
PBE(0) PBin
E ) PBE(ł) - PBdyn
E (łel)
PBel
A ) łelEB[FA,PBor
E , PBel
A] PBdyn
X (łel) ) łelEB[FX,PBdyn
X ] with X ) A, E
-r2V ) 4ðFA -r2V ) 4ðFA
-∞r2V ) 0 -∞r2V ) 0
Vin - Vout ) 0 Vin - Vout ) 0
(@V@n)in - ∞ (@V@n)out ) 4ðóorE (@V@n)in - ∞ (@V@n)out ) 0
Gneq ) G
∞
or-el +
s
FAVor
E dr3 - 1
2 s
FEVor
E dr3 -
1
2 s
Vor
E (óel
E - óel
A) dr2
Gneq ) G
∞
in-dyn +
s
FAVin
E dr3 - 1
2 s
FEVin
E dr3
G
∞
or-el ) E0
A +
1
2 s
FAVel
A dr3 G
∞
in-dyn ) E0
A +
1
2 s
FAVdyn
A dr3
a ł ) (0 - 1)/4ð, łel ) (∞ - 1)/4ð and łor ) ł - łel. b Even if not indicated, the potential V has the same sources of the
corresponding electric field EB.
PB  PBfast + PBslow (123)
PBA ) PBor
E + PBel
A (124)
PBA ) PBin
E + PBdyn
A (125)
3048 Chemical Reviews, 2005, Vol. 105, No. 8 Tomasi et al.

To better illustrate this aspect, let us consider two
polarizable dipoles having centers fixed at r and r′,
both placed in a homogeneous external field. This
electric field determines the mean orientation of the
dipole at r′. The induced dipole in r has two contribu-
tions, the first resulting from the external field and
the second originating from the interaction with the
mean orientation of the dipole in r′. If the external
field suddenly changes so that the mean orientation
of the dipole in r′ has not enough time to relax, the
induced dipole in r will change only as far as concerns
the component related to the interaction with the
external field, whereas the second component will
maintain the value corresponding to the initial field.
As a consequence only the first component presents
a dynamic character in response to the changes of
the external field.
It is important to stress that the two schemes I and
II differ in the form of the slow and fast components
of the polarization, but they are practically and
physically equivalent (in the limit of a linear response
regime) in the sense that they give the same value
of the total reaction potential and, thus, the same
effects on the solute as well as the same interaction
energies. This equivalence has not been always
recognized in the literature,221,525 and thus some
confusion exists on the use of these two partitions.
However, as shown by Cammi and Tomasi522 and,
more explicitly, by Aguilar,526 such confusion disap-
pears when the proper expression for the nonequi-
librium free energy functional is used in each ap-
proach.
The use of one partition instead of the other is
dictated only by the requirements of the specific
computational code in which the model is imple-
mented; for example, in Gaussian,17 the nonequilib-
rium version of PCM has been implemented in the
framework of partition I, exactly as the nonequilib-
rium version of the MPE implemented in Dalton.135
In addition to the proper description of the solvent
response, another aspect must be considered when
we simulate electronic excitations of solvated mol-
ecules. This additional aspect, however, appears only
when a QM description is used.
Within the continuum solvation framework, as in
the case of isolated molecules, it is practice to
compute the excitation energies (or better excitation
free energies) with two different QM approaches: the
state-specific (SS) method and the linear-response
(LR) method. The former has a long tradition and is
related to the classical theory of the solvatochromic
effects (see references in our 1994 review1). The latter
was introduced later in connection with the develop-
ment of the linear response theory for continuum
solvation models,136,138,140,215,527-529 which will be de-
scribed in section 5.6. Here, however, we anticipate
a relevant aspect of the theory when compared to the
SS method for the calculation of transition energies
of solvated molecules.
The state-specific method solves the nonlinear
Schro¨dinger equation for the state of interest (ground
and excited state) usually within CI or MCSCF
descriptions, and it postulates that the transition
energies are differences between the corresponding
values of the free energy functional, the basic ener-
getic quantity of the QM continuum models. As
discussed in section 2.4.6, the nonlinear character of
the reaction potential requires the introduction in the
SS approaches of an iteration procedure not present
in parallel calculations on isolated systems. At each
step of this iteration, the additional solvent-induced
component of the effective Hamiltonian is computed
by exploiting the first-order density matrix of the
preceding step, and the resulting energy, E )
〈¾jHeffj¾〉, is corrected for the work required to
polarize the dielectric G ) E - (1/2)〈¾jVRj¾〉 (we note
that the inclusion of the nonequilibrium needs some
further refinements, as indicated by all of the most
detailed papers on this subject mentioned above). An
important consequence (and disadvantage) of this
iterative approach is that one has to solve a different
problem for each different state instead of obtaining
the energies of all states in a single run as for isolated
systems.
A different analysis applies to the LR approach (in
either Tamm-Dancoff, RPA, or TDDFT version) in
which the excitation energies, defined as singularities
of the frequency-dependent linear response functions
of the solvated molecule in the ground state, are
directly determined by avoiding explicit calculation
of the excited state wave function (see section 5.6).
In this case, the iterative scheme of the SS ap-
proaches is no longer necessary, and the whole
spectrum of excitation energies can be obtained in a
single run as for isolated systems.
Although it has been demonstrated that for an
isolated molecule the SS and LR methods are equiva-
lent (in the limit of the exact solution of the corre-
sponding equations), a formal comparison for mol-
ecules described by QM continuum models shows
that this equivalence is no longer valid. By using a
nonequilibrium approach in both frameworks, it can
be established that these intrinsic differences are
related only to the interaction of the solute with the
fast response of the solvent medium, and they are
connected with the basic assumption of a Hartree
partition between solute and solvent.530 The analysis
of this aspect is quite new, and thus further inves-
tigations will surely appear in the future.
5.3. Solvation Dynamics
Solvation dynamics is the time-dependent, non-
equilibrium, relaxation of the solvent due to an
instantaneous change in the interactions between the
solvent molecules and a different, newly created
electronic state on the solute. Upon excitation of the
solute, the surrounding solvent molecules, which
were initially equilibrated with the ground electronic
state of the solute, are in a nonequilibrium state and
will adjust their positions and momenta in order to
relax to equilibrium with the new interactions as-
sociated from the new electronic state. This relax-
ation is driven by a decrease in the solvent system
interaction energy, and as a result, any measurement
that probes this energy will mirror the underlying
solvation dynamics of the system.
Among the techniques used to investigate solvation
dynamics, the time-resolved fluorescence Stokes shift
Quantum Mechanical Continuum Solvation Models Chemical Reviews, 2005, Vol. 105, No. 8 3049

dependent solvation energy of a generic distribution
of point charges (representing the chromophore) due
to the dielectric is reduced to an integral over the
chromophore molecular surface by introducing a
time-dependent generalization of the BEM. As an
application, a comparison with experimental mea-
surements of solvation relaxation functions for cou-
marin 343- (C343-) in water and C153 in methanol
and acetonitrile was reported. The required data for
the calculations are the atomic charge distribution
change of the chromophore from the electronic ground
state to its first electronically excited state, the
chromophore’s molecular surface, and the frequency-
dependent dielectric constant of the solvent.
More recently, the approach developed by Marcus
and co-workers has been reformulated by Ingrosso
et al.544 within the IEFPCM continuum model and,
for the first time, coupled to a QM description of the
solute electronic states. In this case the time-depend-
ent solvation energy, Esolv(t), is the electronic com-
ponent of the solute-solvent electrostatic interaction
(the nuclear part is a constant independent of time
if we assume a negligible role of solute nuclear
relaxation)
where Vel and qel are the vectors collecting the solute
potential at the surface and the solvent apparent
charges. By comparing eq 131 with the corresponding
one (eq128) in the Onsager-like model, it is evident
that the solute property of interest is no longer the
dipole moment but the electrostatic potential, and the
solvent response is here represented in terms of the
time-dependent apparent charges. Thus, assuming
for V the same evolution used in the Marcus model
to define í(t) in terms of the step function ı(t) and
applying the same mathematical framework of Fou-
rier and anti-Fourier transformations, one obtains
the IEFPCM equivalent of the right-hand side of eq
130. The model has been applied to the calculation
of the STCF for coumarin 153 in various polar
solvents (water, acetonitrile, dimethyl sulfoxide, and
methanol) described in terms of a generalized disper-
sion relation (in the form of a Cole-Davidson or
Cole-Cole equation) for (ö) in the low-frequency
limit and a fit of experimental data at high frequen-
cies. Using a configuration interaction approach
including all singly substituted determinants (CIS)
for the calculation of the coumarin excited states, all
of the major features observed in experimental
results of S(t), the initial inertial decay and the
multiple exponential long time relaxation, are repro-
duced almost quantitatively by the calculations.
Subsequently, the same model has been extended
to describe a general time-dependent relaxation of
the solvent. This extension has been initially pro-
posed within the DPCM scheme by Caricato et al.545
and then applied to the IEFPCM.546 Once again the
starting assumption is that the change in the solute
molecule occurs instantaneously (at time t ) 0) and
the corresponding variation of the electrostatic po-
tential V between the initial (0) and the final (fin)
state is a step change. The time-dependent solvent
polarization charges at a generic time t can thus be
written as
where the vector q0 collects the polarization charges
in the initial solute-solvent equilibrium.
Under the assumption of linear response of the
solvent the variation of the polarization charges äq
at time t due to a change in the electrostatic potential
at time t ) 0, one obtains
where R now indicates the PCM solvent response
function. As done before in eq 130, this expression is
transformed in a numerical procedure by passing
from the time domain to the frequency domain by
using Laplace transformations. This change is re-
quired as the dielectric response of the solvent is
described in terms of its complex dielectric permit-
tivity as a function of a frequency (ö). The use of
Laplace-transformed equations to pass from the time
to the frequency domain has the effect of simplifying
the formalism and of allowing the straightforward
use of the function (ö). Starting from eq 133, in fact,
the variation of the charges becomes
where the definition of the matrix K for the different
PCM schemes can be found in Table 1.
To find the proper form for the frequency-depend-
ent response function R(ö), it is useful to resort to
the iterative form of the PCM equations (see section
2.3.1.5), obtaining
where
The integral in eq 136 can be solved analytically if
we use the Debye expression for (ö); however, the
approach can be applied to any other functional form
for the complex dielectric permittivity (like, for
example, a multiple Debye, a Davidson-Cole, or a
Cole-Cole form) as well as to a combined scheme
including a fit of experimental data for the high-
frequency portion of (ö). The only practical differ-
ence is that in the latter case, the integral in eq 136
is solved.
Another continuum model used to describe solva-
tion dynamics is the frequency-resolved cavity model
(FRCM) formulated by Basilevsky and co-workers.547
Esolv(t) ) -Vel
† qel(t) (131)
q(t) ) q0 + äq(¢V,t) (132)
äq(¢V,t) )
s-∞
t
dt′R(t - t′)ı(t′)¢V (133)
äq(¢V,t) = äq′(¢V,t) + K¢V
äq′(¢V,t) ) 2
ðs0
∞ dö
ö
Im[R(ö)]cos(öt)¢V
(134)
äq′i
(n) )
2ð
Sii
gi(t)[bi + 12ðzi[aj(Sjjäq′j(n-1) +
yj[äq′k
(n-1)])] - 1
2ð
zi[aj](Siiäq′i
(n-1) + yi[äq′i
(n-1)])] -
Sii
-1yi[äq′j
(n-1)] (135)
gi(t) =
2
ðs0
∞ dö
ö
Im[( 4ðˆ(ö)ˆ(ö) - 1)-1] cos(öt) (136)
Quantum Mechanical Continuum Solvation Models Chemical Reviews, 2005, Vol. 105, No. 8 3051

The FRCM model coincides with the DPCM scheme
for the electrostatic equations of the apparent charges
(it is known by the acronym BKO as described in
section 2.3.1.1), but differs from it in the use of two
surfaces instead of a single one. The two surfaces are
here introduced to separate the inertialess (high-
frequency) response of the medium from the inertial
(low-frequency) one. Between the two surfaces the
medium is modeled by the inertialess high-frequency
dielectric constant ∞; outside the outer cavity the
medium is modeled by the static dielectric constant
0. The layer between the two surfaces corresponds
to the first solvation shell.
Calculations in this scheme amount to simulta-
neously solving two equations describing the surface
charge densities on the two cavities, namely, ó∞ (on
the inner surface) and ó (on the outer). In general, a
simplification is applied to the FRCM model, namely,
the charge density on the outer surface is assumed
to be sufficiently small for its influence on the inner
charge density to be neglected. Using this approxi-
mation, the equation describing ó∞ may be separated
from that determining ó, and thus it is identical to
the ordinary BKO equation but with a dielectric
constant now given by ∞. ó∞ then enters as a
parameter into the equation describing ó, which
effectively becomes a BKO equation with scaled
parameters.
Dynamics calculations are performed by replacing
0 with the complex number (ö), thereby splitting
the equation in ó into two equations describing its
real and imaginary parts. Solving these equations
proceeds in a similar way as for the standard BKO
case. Once again, the STCF C¢E(t) is expressed in
terms of ¢E(t), the inverse Fourier transform of the
reorganization energy E(ö), which can be considered
as a solvation energy for the charge distribution
difference between ground and excited state.
The one-cavity BKO and the two-cavity FRCM
approaches have been also extended to dynamical
nonlocal models for the dielectric function by using
the equation
where the spatial dispersion is determined in terms
of a Lorentzian model [¢(k) ) 1/(1 + ì2k2)]. This
nonlocal model when applied to the calculation of the
solvation time correlation function C¢E(t) gives a long-
time asymptote which differs markedly from that of
the local theories and of the experiments. In this
region, corresponding to the low-frequency limit ö f
0 in the frequency domain, the STCF behavior is
completely determined by the size difference between
the solvent and solute particles. A detailed analysis
of these findings led the authors to observe that the
dynamics for the nonlocal spherical model reduces
to the single-exponential dynamics of the simplest
dynamical continuum theory, and thus more realistic
models of the nonlocal dielectric functions should be
introduced.
The examples we have described above demon-
strate that reasonable modifications of the conven-
tional continuum solvent model enable one to over-
come the common belief that continuum medium
models are incapable of adequately treating solvation
dynamics. The improved versions of the continuum
theories are in fact able to qualitatively monitor the
same picture of polar solvent dynamics as do molec-
ular theories. Then, by properly adjusting the model
parameters (a feature that is common for both
molecular and continuum models), experimental de-
pendencies can be well reproduced. The key step in
all continuum treatments is to use the complex-
valued solvent frequency-dependent dielectric func-
tion (ö) as a phenomenological characteristic of
solvent dynamics.
Indeed, continuum treatments present an impor-
tant advantage with respect to standard molecular
simulations, and this could become decisive soon:
when combined with modern quantum chemical
calculations, they allow a precise description of
electronic structure and charge distributions in chemi-
cally nontrivial solute species, whereas current mo-
lecular treatments still invoke crude solute models
and concentrate on a description of the details of
solvent structure.
We finish this section by noting that almost all of
these detailed studies on time scales and mechanisms
of solvation dynamics have been limited to polar
solvents. For nonpolar solvation, with the notable
exceptions of the works from Berg and co-
workers548-550 and from the Maroncelli group,551,552
few experimental studies exist. In addition, the
feedback between theory and experiment that has
been so productive for polar solvation has not yet
been achieved for nonpolar or, more generally, non-
dipolar systems, where nonpolar systems refer to
systems that interact via weak nonelectrostatic forces,
such as dispersive and repulsive interactions, whereas
non-dipolar systems include higher electrostatic mul-
tipole (e.g., quadrupole or octupole) moments in
addition to the nonpolar interactions. There are a
number of reasons why non-(di)polar solvation has
been less deeply analyzed. First, there is no obvious
method to extract the dynamical time scales in non-
dipolar liquids from existing frequency domain data
such as dielectric relaxation measurements, because
non-dipolar systems do not exhibit significant dielec-
tric signals. Second, the assumptions that greatly
simplify the understanding of dipolar solvation,
linear response and the relative insensitivity of the
time scales to the specific nature of the probe
molecule, may not hold for nonpolar or non-dipolar
solvation.539 In other words, the short-range nature
of the interactions implies that the molecular proper-
ties of the solvent-solute combination play a more
important role in the solvation dynamics of non-
dipolar systems than in the dipolar case.
Despite this increased difficulty, there are some
models that can successfully predict the long time
dynamics observed in nonpolar systems.
For instance, Berg553,554 formulated a novel theory
to describe nonpolar solvation in which the response
of the nonpolar solvent is not based on dipolar
interactions, but instead treats the solvation dynam-
ics as a re-equilibration of the solvent to a changing
volume of a solute. This model is an analogue of
continuum models of polar solvation, which predict
(k,ö) ) 
∞
+ ((ö) - 
∞
)¢k(k) (137)
3052 Chemical Reviews, 2005, Vol. 105, No. 8 Tomasi et al.

solvation dynamics from the frequency-dependent
dielectric constant, but here the solvation reponse
function is calculated from the time- or frequency-
dependent mechanical moduli of the solvent. The
hypothesis is the assumption that a change in the
solute’s effective size or shape upon electronic excita-
tion is the primary interaction driving nonpolar
solvation. Following excitation, the solvent must
move to allow the solute cavity to expand. This
movement is modeled by treating the solvent as a
viscoelastic (VE) continuum. This VE model predicts
two distinct time scales that can be associated with
the inertial and diffusive mechanisms and time
scales. Comparisons with experimental data show
this model fits long time solvation dynamics well.
5.4. Spectral Line Broadening and Solvent
Fluctuations
The progress in solvation dynamics we have sum-
marized in the previous section has also stimulated
the application of continuum models to the study of
the line shape of bands corresponding to absorption/
emission processes and of the solvent fluctuations.
Such fluctuations in the structure of the solvation
shell surrounding the chromophore555 and the con-
sequent variation in the local electric field lead to a
statistical distribution of the energies of the electronic
transitions. This phenomenon is called inhomoge-
neous broadening and, in most cases, its extent is
much greater than that of homogeneous broadening
due to the existence of a continuous set of vibrational
sublevels in each electronic state. For obvious reasons
we focus here only on the inhomogeneous broadening
by giving a short review on the main issues related
to its modelization.
A widely used approach to the analysis of solvent
spectral broadening is based on the model formulated
by Mukamel and co-workers556,557 in which the sol-
vent is represented as a bath of Brownian oscillators.
This model provides convenient analytical expres-
sions for the line broadening function g(t). Imple-
mentation of the Brownian oscillator model requires
only two adjustable parameters: the frequency and
the variance of solvent-induced frequency fluctua-
tions, which in turn is related to the solvent reorga-
nization energy.558,559 The application of the Brown-
ian oscillator model to solvent-induced line broadening
typically assumes that the solvent modes coupled to
the electronic transition are overdamped, and thus
the time correlation function for solvent-induced
frequency fluctuations follows an exponential relax-
ation. However, as we have described in the previous
section, theory and experiments on solvent relaxation
suggest that at sub-picosecond times this relaxation
is a Gaussian function of time, resulting from inertial
solvent motion, followed by exponential decay char-
acteristic of diffusional motion on the picosecond time
scale.
To account for this inertial part of solvent relax-
ation, one can start from a phenomenological cor-
relation function, which is a Gaussian function of
time. In the limit of linear solvent response, in which
the solvent dynamics are assumed to be independent
of the electronic state of the chromophore, the
frequency-frequency correlation function (FFCF) is
defined by
where äöeg(t) is the solvent-induced fluctuation in the
frequency öeg of the electronic transition of the
chromophore and the angle brackets indicate an
equilibrium average over the solvent coordinates. We
note that in this linear solvent response limit the
FFCF of eq 138 is related to the function S(t)
determined from time-resolved Stokes shift measure-
ments as shown in eq 127 (namely, the previously
defined STCF is the FFCF normalized to unity at t
) 0).
Several very general relationships between optical
observables can be obtained for such linearly re-
sponding (Gaussian) solvents interacting with a fixed
charge distribution of the solute. More precisely,
these spectroscopic observables can be related to the
solvent reorganization energy, ìs.560 This is in fact
the only parameter needed to characterize the ener-
getic intensity of the solvent nuclear fluctuations
coupled to the solute.
The relationship between ìs and the solvent-
induced absorption, and emission, spectral width
(óabs/em), can be written as561,562
where â ) 1/kT. This relationship is exact for any
fixed charge distribution of the solute in a linearly
responding solvent.
Here, we do not want to enter in the large field of
theoretical models to determine the solvent reorga-
nization energy, which followed the original Marcus
formulation of electron transfer521 (some reviews can
be found in refs 563 and 564). Here, it is important
to recall that this concept is strictly connected to the
difference in solvation free energy between the initial
nonequilibrium and the final equilibrium relaxed
solvent configurations originated by any fast process
of solute charge redistribution such as the vertical
electronic transition described in section 5.2.240,412,545,565
The effect of relaxing the assumption of a fixed
charge distribution of the solute allowing the elec-
tronic density to flow between the initial and final
states (electronic delocalization) is to alter the solute
electric field, thus modifying the solute-solvent
coupling and the line shift, and, even more, the
inhomogeneous broadening and the observed optical
width. This is because, in an electronically delocalized
solute, the charge distribution changes instanta-
neously with solvent configuration, so that each
configuration “sees” a different charge distribu-
tion.566-568 This self-consistent action results in es-
sentially nonlinear features of the solvent effect.569-571
We conclude this section by mentioning a different
class of approaches to the study of the line shape in
solvated systems, which explicitly account for the
positional fluctuations of solvent molecules around
the solute. For example, this is the approach devel-
oped by Kato and co-workers572 within the RISM-SCF
model (see section 7.4.2). In this framework the
spectral bandwidths are estimated by calculating the
C(t) ) 〈äöeg(t)öeg(0)〉 (138)
ìabs
w ) âóabs
2 /2 and ìem
w ) âóem
2 /2 (139)
Quantum Mechanical Continuum Solvation Models Chemical Reviews, 2005, Vol. 105, No. 8 3053

[( + 2)/3]2, where  is the optical dielectric constant
of the medium. Agranovich and Galanin obtained the
same prefactor in 1982 from basically the same
considerations in a classical theory.579 More recently,
two QM models for EET including solvent effects
have appeared in the literature. Both models are
based on a time-dependent variational approach even
though they follow two completely different paths.
In the first model by Tretiak et al.580,581 a Hartree-
Fock approximation is combined with semiempirical
Hamiltonians, whereas in the second, by Hsu et al.,582
a DFT description is used. To account for solvent
effects, both methods introduce a simplified con-
tinuum description based on two important ap-
proximations, one on the form of the molecular cavity,
which is always limited to a sphere, and the other
on the solute electronic densities, which are repre-
sented in terms of a dipolar (or a multipolar) expan-
sion.
More recently, a new model has been presented by
Iozzi et al.583 This is still based on a continuum
description of the solvent, exactly as in the models
described above, but it differs from them in two
fundamental aspects. On the one hand, it introduces
the solvent effects in all of the steps of the QM
calculation of the EET process without requiring any
simplification in the representation of the molecular
charge distributions or the introduction of fixed
prefactors, but defining proper operators to be added
to the Hamiltonian and to the resulting response
equations. On the other hand, it describes the mo-
lecular cavity on the real three-dimensional structure
of the interacting solutes without imposing an arti-
ficial spherical boundary. In this framework, the EET
process is described through the time-dependent
density response theory (or TDDFT, see section 5.6)
originally proposed in the paper by Hsu et al.,582
whereas the introduction of the solvent effects is
realized in terms of the IEFPCM described in section
2.3.1.3.
The model considers two solvated chromophores,
A and D, with a common resonance frequency, ö0,
when not interacting (this condition is trivially
fulfilled when A  D, but it is still realistic in the
case of either large or condensed systems, i.e.,
whenever high density of nuclear states can be
found). When their interaction is turned on, their
respective transitions are no longer degenerate. In
contrast, two distinct transition frequencies ö+ and
ö-, appear. The splitting between these defines the
energy transfer coupling, JDA ) [ö+ - ö-]/2, which
can be evaluated by computing the excitation ener-
gies of the D x A system, as, for example, through a
proper TDDFT scheme (see eq 152 in the following
section).584-589 The alternative strategy used in the
model is based on a TDDFT-perturbative approach,582
which considers the D/A interaction as a perturbation
and includes solvent effects defining an effective
coupling matrix, which couples transitions of D with
those of A in the presence of a third body represented
by the IEF dielectric medium. Within this frame-
work, the first-order coupling J becomes a sum of two
terms, one of which is always present (i.e., also in
isolated systems) and another that has an explicit
dependence on the medium, namely
and
where FT(r) indicates transition densities of the
solvated systems D and A in the absence of their
interaction. In particular, J0 describes a solute-
solute Coulomb and exchange correlation interaction
(through the exchange correlation functional of the
proper DFT scheme), corrected by an overlap contri-
bution. The effects of the solvent on J0 are implicitly
included in the values of the transition properties of
the two chromophores before the interaction between
the two is switched on. These properties can in fact
be significantly modified by the “reaction field”
produced by the polarized solvent. In addition, the
solvent explicitly enters into the definition of the
coupling through the term JIEF of eq 142, which
describes a chromophore-solvent-chromophore three-
body interaction (we note that this term is originated
from the explicit solvent term defined in the TDDFT/
IEFPCM eq 155 of the following section).
The JIEF term can account for effects due to
dynamic solvation by introducing the dynamic (ö)
instead of the static permittivity in the definition of
the transition apparent surface charges q(sk; ö,FT).
These latter generalize the concept of ASC described
in section 2.3.1 to transition densities as the source
of polarization, but their definition is exactly equiva-
lent to that given when the source is represented by
a standard state density.
There are two the fundamental advantages of the
perturbative TDDFT/IEFPCM description of the
coupling: one of a physical nature and the other
numerical. The conceptual aspect, unique to this
approach, is that the coupling is computed within a
scheme coherently accounting for different kinds of
interactions (namely, direct and indirectsi.e., solvent-
mediated Coulomb, exchange, possibly including cor-
relation effects and overlap interactions). In this way,
both the reaction and the screening effects due to the
solvent are introduced in a unified model, which
merges Fo¨rster- and Dexter-type approaches. The
numerical advantage of the perturbative approach
with respect to the standard one in which the TD
scheme for the supermolecule (D x A) is solved is that
only the properties of the solvated chromophores (D
and A) are required to get the coupling.
This feature of the TDDFT-PCM perturbative
approach becomes particularly important when one
applies the model to the calculation of excitonic
splittings590 in conjugated organic materials formed
by equivalent polymeric chains in a three-dimen-
sional array.591 In this case, in fact, one can introduce
both short-range and long-range interchain effects
simply by computing the properties of a single chain
embedded in a continuum medium characterized by
J0 )
s
dr
s
dr′ FD
T*(r′)( 1
jr - r′j
+
gxc(r′,r,ö0))FAT(r) - ö0 s dr FDT*(r)FAT(r) (141)
JIEF ) ∑
k (s dr FDT*(r) 1jr - skj)q(sk; ö,FAT) (142)
Quantum Mechanical Continuum Solvation Models Chemical Reviews, 2005, Vol. 105, No. 8 3055

an anisotropic dielectric tensor (the two different
components of the tensor correspond to the dielectric
permittivity along and perpendicular to the main axis
of the polymer chain). This extension to anisotropic
dielectrics is made possible by the use of the IEFPCM
formalism presented in section 2.3.1.3.
Processes of energy transfer are also responsible
for another interesting phenomenon strongly de-
pendent on the environment. It is well-known that
the presence of a metal body can strongly affect the
response properties of a molecule placed in its close
proximity. In section 6.1.3, we shall examine the well-
known surface-enhanced Raman scattering (SERS)
for which enhancement factors of >10 orders in the
magnitude of the Raman scattering of molecules close
to metal particles have been reported. As far as
concerns the excited state properties, the presence
of the metal can have also the opposite effect and is
responsible, in many cases, for a decrease of the
molecular responses. This result is due to the fact
that the excitation energy of the molecule can be
efficiently transferred to the metal body (for example,
via the resonant energy transfer we have just de-
scribed) and can undergo, inside this medium, several
dissipation processes.
The great technological interest in these phenom-
enona has led to the formulation of many theoretical
models. Most of them consider the molecule as a
polarizable pointlike dipole.592-595 For metal-mol-
ecule distances of less than a few nanometers (a
situation that occurs in many experimental situa-
tions) this approximation can be too rough and a QM
description of the whole molecule would instead be
necessary. On the other hand, for many of the
physical systems of interest, the metal shape is one
of the most important factors, and a proper descrip-
tion of this aspect is compulsory.
To account for both of these demands, Corni and
Tomasi79,504,506 have formulated a QM method for the
calculation of dynamic response properties of a
molecule in close proximity to metal bodies and
possibly in the presence of a solvent. Such a model
(see also sections 4.6 and 6.1.3) treats the molecule
at a QM level and can explicitly consider metal
particles of complex shape. The metal (and the
solvent, if present) is described as a continuous body
characterized by its (dielectric) response properties
to electric fields, both those imposed on the system
from outside and those arising from the molecular
charge distribution. The metal-molecule and solvent-
molecule interactions are treated within the PCM
approach in its D and IEF versions. Recently, this
procedure has been generalized for the calculation
of the contribution to the molecular deexcitation rate
due to the metal;508,509 this generalization is based
on the response theory (applied at a TDHF or TDDFT
method, see the following section), in which excitation
energies and lifetimes can be obtained from the real
and imaginary parts of the poles of the linear
response function. Strictly speaking, the lifetimes
obtained in this way are related to the width of the
absorption peaks of the molecule.
5.6. Time-Dependent QM Problem for Continuum
Solvation Models
We conclude section 5 devoted to time-dependent
(TD) solvation by considering in a more detailed way
all of the QM issues implicitly or explicitly cited in
the previous subsections more focused on the physical
and modelistic aspects of the phenomenon.
Once again, the key point is the nonlinear charac-
ter of the effective Hamiltonian defining the molec-
ular system in the presence of the solvent (see section
2.4.5).
Here, in particular, such an effective Hamiltonian
becomes time dependent. There are a large variety
of processes in solution that can give rise to this time
dependency. Here we limit ourselves to two cases: (a)
the solvent dynamics relaxation and (b) the presence
of an external time-dependent perturbing field. In
general, we can write
where Hö M
0 is the Hamiltonian of the unperturbed
solute, Vint the solute-solvent interaction term, and
W(t) a general time-dependent perturbation term
that drives the system.
In case a discussed in section 5.3 we have that W(t)
) 0 and the time dependence of the solute-solvent
interaction Vint(t) originates from dynamical pro-
cesses involving inertial degrees of freedom of the
solvent. The time scale of these processes is orders
of magnitude higher than the time scale of the
electron dynamics of the solute, and an adiabatic
approximation can be used to follow the electronic
state of the solute, which can be obtained as eigen-
state of the time-dependent effective Hamiltonian (eq
143).
This approach is a direct extension of the QM
problem for the basic model (see eq 58) to the effective
Hamiltonian Heff(t) and, as shown in section 5.3, has
been used by various authors377,544,545,596 to study the
time correlation function for the time-dependent
Stoke’s shift following a vertical excitation and to
study electron-transfer processes of molecular sol-
utes.
For case b, which will be considered in more detail
in section 6, the external perturbing fields of interest
are usually optical field, and the study of response
properties of a molecular solute to the field requires
an extension of the basic model to describe the time-
dependent nonlinear Schro¨dinger equation for the
solute, namely
This extension has been formulated by Mikkelsen et
al.136 for the MPE solvation model (see section 2.3.2)
and by Cammi and Tomasi597 within the framework
of the PCM model (see section 2.3.1).
Mikkelsen and co-workers have considered the
response of a reference variational (HF/MCSFC)
136-138,598 state by using the Frenkel variation prin-
ciple for the TD-QM problem in the form of the
Ehrenfest equation, which describes the time evolu-
Heff(t) ) HM
0 + Vint(t) + W(t) (143)
i
@¾(t)
@t
) [HM
0 + Vint(t) + W(t)]¾(t) (144)
3056 Chemical Reviews, 2005, Vol. 105, No. 8 Tomasi et al.

tems;624,625,628 in these cases, as no analytical deriva-
tives are available, studies of geometry effects have
been limited to a few important motions such as
twisting and wagging modes in push-pull aromatic
molecular systems.629 In parallel, linear and nonlin-
ear response methods (see section 5.6) have been used
to calculate solvent effects on excited state properties
such as dipole moments and (hyper)polariza-
bilities.215,630-632
6.1. Energy Properties
Modern quantum chemistry has been revolution-
ized by the ability to calculate analytical derivatives
of the wave function and of the energy with respect
to various parameters. A reference book by Yama-
gouchi et al., published in 1994,633 is rightly titled A
New Dimension to Quantum Chemistry. Analytical
Derivative Methods in ab Initio Molecular Electronic
Structure Theory. The appreciation of the role of
analytical derivatives expressed in this title was
speculative when formulated in 1993, but now may
be considered a reality. Nowadays, practically all ab
initio methods are provided with analytical deriva-
tives, and this availability has had a remarkable
impact on computational strategies.
In the past decade there has been a similar
evolution in QM solvation models, and in this section
we describe how the methodologies originally formu-
lated for isolated systems have been generalized to
evaluate energy properties such as IR and Raman
spectra.
6.1.1. Geometrical Derivatives
The characterization of the free energy hypersur-
face of solvated molecules, G(R), is based on free
energy gradients, refined by second-order derivatives
with respect to the nuclear coordinates. The defini-
tion of the equilibrium geometry is based on the
search of critical points corresponding to local minima
on G(R), the definition of molecular motions leading,
for example, to chemical reactions are again based
on the examination of the gradient, supplemented by
the examination of the sign of the local Hessian
matrix to characterize transition states and minima.
The diagonalization of the Hessian is also necessary
to have the elements to characterize the molecular
vibrations, and thus to study all of the properties that
depend in some way on nuclear motions.
These geometrical derivatives can be computed as
finite differences as currently done in classical and
semiempirical continuum models, in which repeated
calculations of the energy are not excessively costly.
For ab initio methods the use of numerical deriva-
tives strongly reduces the size of the systems sub-
jected to scrutiny and the quality of the QM formu-
lation.
The elaboration of analytic QM methods to com-
pute derivatives of the free energy with respect to
nuclear coordinates in continuum solvation models
started at the beginning of the 1990s. The evolution
of this elaboration, which started from methodologies
developed in the 1970s for ab initio methods in vacuo,
has encountered additional specific problems. Most
of these problems have been solved, but the develop-
ment in the algorithms is still in progress.
Some among the first codes for gradients in con-
tinuum methods adopted the approximation of keep-
ing the cavity rigid (see, for example, ref 634). This
is a minor approximation for crude cavities, but it
represents a strong approximation when the cavity
is accurately modeled in terms of spheres centered
on nuclei as done in the most recent versions of the
continuum methods. The infinitesimal displacement
of a nucleus, considered in a given partial derivative,
entails an analogous infinitesimal motion in the
corresponding atomic sphere, which in turn produces
an infinitesimal change in the intersection of this
sphere with the others. The neglect of this infinitesi-
mal motion of the cavity gives origin to apparent
forces acting on the surface, which can have an
important effect, as it results from the comparison
of the derivatives computed with or without this
approximation.
These considerations hold for all methods using a
cavity modeled on the molecular volume. Methods
belonging to the ASC family, or making otherwise
use of a discretization of the cavity (see section
2.3.1.5), have other specific problems. For example,
using a cavity formed by interlocking spheres, surface
elements at the intersection of two spheres in general
have an irregular form, depending of the position of
the intersection circle with respect to the mesh. A
correct procedure for analytical derivatives must be
able to give a satisfactory solution to this problem.
In addition, in some of the algorithms used to build
the cavity (for example, the GEPOL approach43-45,92
used in the PCM class of models) an additional
problem appears. To avoid crevices in which solvent
cannot enter, additional spheres not centered on the
nuclei can be automatically introduced by the algo-
rithm; the positions of these added spheres depend
on the infinitesimal displacements of nuclei, and thus
their contribution to the analytical derivatives must
be correctly taken into account. Analogous consider-
ations hold for cavities making use of re-entrant
portions of the surface (the Connolly surface, for
example).
Another point deserving mention is the use of
derivatives for solutes of large size. Methods for
analytical derivatives in molecules now scale linearly
with the size of the system. Linear scaling in con-
tinuum methods is being actively pursued, and the
solutions to this problem require notable changes in
the formulation of analytical derivatives. Some ex-
amples in the direction of more efficient computa-
tional strategies to compute geometrical derivatives
of large solvated systems have recently been pre-
sented within the framework of PCM.62,635
Following the formalism used in section 2.4.2,
analytical derivatives for solvated systems can be
obtained by directly differentiating the expression of
the free energy G given in eq 75, taking into account
the SCF stationarity constraints. In this way we
Quantum Mechanical Continuum Solvation Models Chemical Reviews, 2005, Vol. 105, No. 8 3059

obtain for the first and second energy deriva-
tives636,637
where R and â are nuclear coordinates.
In eqs 156 and 157 we have introduced the short-
ened notation GR ) @G/@R and GRâ @2G/@R@â, omitted
the superscript 0 for the matrices h0 and G0, and
added superscripts on the matrices to denote deriva-
tives of the corresponding integrals over the basis set.
The solvent-induced terms, hR
RXR(P) and XRâ(P),
are obtained by differentiating once or twice the
matrices of eqs 64-67, for example
where the second term in brackets involves deriva-
tives of the apparent charges.
In eqs 156 and 157 we have introduced the ad-
ditional matrix, W ) PFP, and its derivatives, for
instance
where F is the Fock matrix including solvent terms
(see eq 63) and F(â) its derivative, namely
We note that the derivative of the density matrix
P is not necessary for the first derivative GR because
the SCF stationary condition leads to the equivalence
tr PRF ) tr SRPFP as for the isolated molecule. In
contrast, density matrix derivatives are necessary for
the second- and higher order derivatives. In particu-
lar, the second derivatives require the first deriva-
tives of P (one for each parameter), which can be
obtained by solving the appropriate coupled per-
turbed equations. These equations can be put in the
general form
with the usual orthonormality condition:
The elements of the matrices hR
R,XR(P) and
hR
Râ,XRâ(P) contain the derivatives of the surface
apparent charges with respect to nuclear displace-
ments. These derivatives depend on the specific
equation used to compute q (see eq 34), but in any
case they require that we differentiate all of the
geometrical factors which define the ASC matrix K.93
Cance`s and Mennucci638,639 proposed an alternative
way to proceed within the IEFPCM framework; such
a procedure avoids computing any geometrical de-
rivative of the surface charges. Their approach
consists of differentiating first the basic electrostatic
equation
where ì indicates a solute nuclear coordinate and -
(ì) ) 1 inside the cavity and (ì) )  outside, and
obtaining @V/@ì as a solution of the differentiated
equation. This quantity is then inserted in the
formula of the derivatives of the solute-solvent
interaction energy. Within this approach, the only
term coming from the dependence of the cavity on
the motions of the nuclear coordinates is reduced to
a simple sum of the apparent charges q(sk) lying on
the moving part ¡ of the cavity itself, namely
The function U¡ assumes a very simple form for
standard cavities given as unions of spheres, each of
them centered on a solute nucleus. In this case we
can write
when ì is the Rth coordinate of the lth nucleus. In
the aforementioned expression, (e1, e2, e3) is the
orthonormal basis of the real space and ¡l is the part
of the cavity belonging to the sphere centered on the
lth nucleus.
Going back to the QM expressions, the free energy
derivative (eq 156) reduces now to
where the fifth and sixth terms are those collecting
the explicit solvent terms. The cavity-dependent term
ô is shown in eq 164, whereas D, which denotes the
interaction energy between the variation of the solute
charge distribution @F/@ì and the apparent charges,
is expressed as
From eqs 164-167 it appears evident that no
derivatives of the apparent charges q are required.
The same method has been extended to second
derivatives;640 also in this case no derivatives of the
charges are required.
GR ) tr PhR + 1
2
tr PGR(P) + tr PhR
R +
1
2
tr PXR
R(P) - tr SRW + Vnn
R +
1
2
Unn
R (156)
GRâ ) tr PhRâ + 1
2
tr PGRâ(P) + tr PhR
Râ +
1
2
tr PXRâ(P) - tr SRâW + 1
2
Unn
Râ + Vnn
Râ +
tr PâhR + tr PâGR(P) + tr PâhR
R + tr PâXR(P) -
tr SRWâ (157)
Xíî
R (P) ) ∑
k
[@Víî(sk)@R qe(sk) + Víî(sk) @qe(sk)@R ] (158)
Wâ ) PâFP + PF(â)P + P[G(Pâ) + XR(P
â)]P +
PFPâ (159)
Fâ ) hâ + hR
âP + Gâ(P) + XR
â (P) (160)
FPR + FRP - PRF - PFR + FPSR - SRPF ) 0
(161)
PPR + PRP + PSRP ) PR (162)
-div[(ì)rV(ì)] ) 4ðF(ì) (163)
ô )
4ð
 - 1
∑
k
q2(sk)
ak
[U¡(sk)ân(sk)] (164)
U¡(sk) ) {0 if sk ∉ ¡leR if sk 2 ¡l (165)
GR ) tr PhR + 1
2
PGR(P) - tr SRW + Vnn
R +
D(@F/@ì,q) + 1
2
ô (166)
D(@F/@R,q) ) ∑
k
@V
@Rjsk
q(sk) (167)
3060 Chemical Reviews, 2005, Vol. 105, No. 8 Tomasi et al.

of the incident electric field due to the whole metal
aggregate, that due to the interaction of the molecule
with the local hot spot combined with the effect of
interaction of the molecule with the liquid phase (if
present), that due to the interaction of the light
emitted by the molecule plus hot spot system with
the rest of the metal cluster, and that due to chemical
interaction of the molecule with the metal.
The effects on electromagnetic fields due to the
presence of metal bodies represent a unifying element
in the origin of the SEIRA and SERS phenomena
even if each phenomenon has its own features.
Because SERS has been the most studied, here we
shall focus on that only.
A considerable amount of theoretical work on this
subject has been done since the discovery of SERS,
but most of these, although quite refined in the
description of the metal,668 represent the molecule
just as a polarizable point dipole. Another approach
used on some occasions669,670 consists of a whole ab
initio description of the molecule interacting with a
metal cluster of a few atoms. This model accounts
also for chemisorption effects, but the small metal
cluster clearly cannot reproduce the behavior of, for
instance, a nanosized metal particle, which is com-
posed of thousands of atoms.
If we limit ourselves to the study of molecules
physisorbed on a metal, that is, molecules which are
not chemically bound to the metal surface, the
problem of a proper description of such systems
strongly resembles that for solvated molecules. In
fact, a molecule in solution weakly interacts with a
great number of other electrons and nuclei, feeling
all of the fields of external as well as internal origin,
as a molecule physisorbed on a metal surface does.
Following the theoretical lines of the PCM approach,
an extension of such a model to treat molecules close
to metal particles (which, in turn, can be immersed
in a dielectric medium) has been presented79,504,507
(see also section 5.5). In such an extension, the metal
is described as a continuous body characterized by
electric response properties only, which behaves as
a perfect conductor for static fields and as a dielectric
for time-dependent fields. Infinite planar metal sur-
faces, complex shaped nanoparticles, and aggregates
of nanoparticles have been considered, but in any
case the only interactions considered among the
various portions of the complex system (the molecule,
the metal, and, possibly, the solvent) are electrostatic
in origin.
The basic idea underlying this model is that it is
not necessary to treat the whole system (metal
particle plus molecule) to the same degree of ac-
curacy, but it is better to focus the attention on the
portion of the system of which we want to study the
properties.
Subsequently, the model has been improved507 by
partitioning the system into three “layers,” described
at different levels of accuracy. The whole metal
cluster (the first layer) is treated as a collection of
spherical metallic particles described as polarizable
dipoles to be able to recognize where the hot spots
are at a given excitation wavelength and how strong
the electric field locally acting on them is. The metal
particles in a given hot spot (the second layer) can
be considered as a unique complex shaped particle
and treated by using the original model, which
neglects the detailed atomic and electronic structure
of the metal but which takes into account effects on
the electric field due to the complex shape (union of
interlocking spheres). Finally, the molecule (the third
layer) is described at the ab initio level, taking into
account electrostatic interaction with the metal par-
ticle and properly treating the electric field acting on
it. By exploiting this model, huge enhancement
factors have been found for Raman spectra of mol-
ecules close to different metal particle aggregates,
which compare well with experiments.
6.2. Response Properties to Electric Fields
The investigation of linear and nonlinear optical
(NLO) properties of solvated molecules and liquids
has been one of the most important developments of
the QM continuum solvation models in recent years.
Such progress has been prompted by the success of
modern quantum chemical tools in predicting re-
sponse properties for molecules in the gas phase.
The experience gained for isolated molecules has
in fact represented a fundamental background as
concerns theoretical issues as well as more physical
aspects. Following these studies, for example, con-
tinuum solvation models have been extended to treat
not only the leading electronic term of the optical
properties but also (even if in a less widespread way)
the vibrational contributions. Due to the far larger
number of studies present in the literature, in the
following sections we shall focus on the electronic
contribution only. However, it is interesting to add
some notes on the vibrational part also because
studies appeared so far seem to show that vibrational
effects can become very important (if not dominant)
for NLO properties of solvated molecules.671,672
Vibrational contributions are generally divided into
two components,673-676 the “curvature” related to the
field dependency of the vibrational frequencies (i.e.,
the changes in the potential energy surface in the
presence of the external field) and including the zero-
point vibrational correction and the “nuclear relax-
ation” arising from the field-induced nuclear relax-
ation (i.e., the modification of the equilibrium geometry
in the presence of the external field).
The nuclear relaxation is generally the dominant
contribution, especially when in the presence of
applied static fields. The most largely used method
to account for these effects is based on series expan-
sion of the energy and of the electric dipole with
respect to both the nuclear coordinates and the
external field components. As shown in refs 671, 672,
and 677, the same procedure can be applied to
solvated systems when described through continuum
solvation models; this equivalence is made possible
if expansions of the free energy functional are used
and the related dipole and polarizability derivatives
with respect to normal coordinates are computed in
the presence of the solvent.
The equivalence between condensed phase and gas
phase systems concerning the calculation of optical
properties, however, is not complete as the presence
Quantum Mechanical Continuum Solvation Models Chemical Reviews, 2005, Vol. 105, No. 8 3063

of the environment introduces more critical issues
than for isolated systems, but it also offers the
possibility to compare the computed microscopic
properties with experimental macroscopic quantities.
Formal correspondences and physical differences will
be analyzed in the following two subsections, respec-
tively.
6.2.1. QM Calculation of Polarizabilities of Solvated
Molecules
In section 5.6 we have presented the general time-
dependent QM theory for solvated systems; here, we
come back to such a theory but specialize it to the
case of a perturbation represented by the combina-
tion of a static (E0) and an oscillating electric field
(Eö). In this case, the operator W(t) in eq 143 becomes
where í is the electronic dipole moment operator of
the solute.
By using the Frenkel variational principle (eq 146)
within the restriction of a one-determinant wave
function with orbital expansion over a finite atomic
basis set, we obtain the following time-dependent
Hartree-Fock (or Kohn-Sham) equation
with the proper orthonormality condition; S, C, and
 represent the overlap, the MO coefficients, and the
orbital energy matrices, respectively.
In eq 173 the prime on the Fock matrix indicates
that terms accounting for the solvent effects are
included, that is
where the first four terms on the right-hand side of
the equation represent the Fock operator for the
solvated molecule (see section 2.4.2 for more details)
and m represents the matrix containing the dipole
integrals; we note that in eq 174 we have discarded
the superscript “0” used in the previous sections to
indicate the Hamiltonian terms of the isolated mol-
ecule (here the same superscript is used to indicate
static fields) and the superscript R for the solvent-
induced terms, j and X(P).
The solution of the time-dependent HF or KS eq
173 can be obtained within a time-dependent coupled
HF (or KS approach). We first expand eq 174 to its
Fourier components, and subsequently, each compo-
nent is expanded in terms of the components of the
external field. The separation by orders leads to a
set of CPHF equations for which the Fock matrices
can be written in the form
where a,b, ... indicate the Cartesian components of
the field and öx are the frequencies related to the
external fields (eventually static and thus öx ) 0).
We note that the elements of the solvent-reaction
matrices XöT (see eqs 67 and 68) depend twice on the
frequency-dependent nature of the field, in the
perturbed density matrices Pab...(ö1,ö2,...), and in the
value of the solvent dielectric permittivity (öT) at
the resulting frequency öT ) “öx (see section 5.2).
Once solved, the proper CPHF equations, the
dynamic polarizabilities of interest, can be expressed
in the following forms:
We note that an alternative (but equivalent) method
can be used to get Rab(ö); namely, we can resort to
the linear response theory presented at the end of
section 5.6 and use eq 151 with B and A being the
dipole moment operators, that is, 〈〈íˆa;íˆb〉〉ö (to get the
analogue equation for â and ç a quadratic and a cubic
response theory have to be invoked).
Applications of the QM theory for the calcula-
tion of linear and nonlinear optical properties of
solvated molecules have been carried out by
Willets and Rice,678 Yu and Zerner,679 Mikkelsen
and co-workers,139,516,680-684 and the PCM
group,599,655,656,671,685-689 all applying different solvent
models as well as different levels of the quantum
theory.
6.2.2. Definition of Effective Properties
A key factor in the study of response properties is
the relationship between the molecular and the
material properties, that is, the distinction between
the microscopic field at a molecule and the macro-
scopic field applied to the material.690 The macro-
scopic field is determined directly from the experi-
mental conditions, but the microscopic field has to
be calculated using a suitable theory.
The classical approach is to use the Onsager-
Lorentz model and write the measured susceptibili-
ties (the macroscopic equivalent of the linear and
nonlinear optical, NLO, molecular properties) in
terms of the gas phase molecular values multiplied
by so-called local field factors (see also section 6.1.2).
In recent years different models have been proposed
to improve such a description. Most of them are based
on the concept of a solute in a cavity in a dielectric
continuum provided by the rest of the material.
Within the framework of semiclassical descriptions,
important generalizations of the Onsager-Lorentz
model for NLO properties have been made by
Wortmann and Bishop691 and by A¡ gren and co-
workers,677,692-695 whereas a QM approach has been
formulated within the PCM approach described in
section 2.3.655,656,689
The PCM extension to get macroscopic susceptibili-
ties that can be directly compared with experiments
has been achieved by introducing so-called effective
Wö (t) ) íˆâ[EBö(eiöt + e-iöt) + EB0] (172)
F′C - i @
@t
SC ) SC (173)
F′ ) h + G(P) + j + X(P) +
mâ[Eö(eiöt + e-iöt) + E0] (174)
F′a(ö) ) ma + m˜
ö + G(Pa(ö)) + Xö(P
a(ö))
F′ab...(-öT;ö1,ö2,...) ) G(P
ab(-öT;ö1,ö2,...)) +
XöT(P
ab(-öT;ö1,ö2,...))
(175)
Rab(-ö;ö) ) -tr[maP
b(ö)]
âabc(-öT;ö1,ö2) ) -tr[maP
bc(-öT;ö1,ö2)] (176)
çabcd(-öT;ö1,ö2,ö3) ) -tr[maP
bcd(-öT;ö1,ö2,ö3)]
3064 Chemical Reviews, 2005, Vol. 105, No. 8 Tomasi et al.

only difference is that now perturbed electron density
Pa(ö) [or Pab(ö)] is different due to the presence of
the new perturbing term (eq 184) in the related Fock
matrices.
We conclude by noting that to extend the QM
calculations of effective polarizabilities to pure liq-
uids, a further element has to be introduced. In that
case, the optical radiation at the fundamental fre-
quency ö produces in the liquid a macroscopic
polarization density at the output frequency, which
acts as source of an additional perturbing term.700
6.3. Response Properties to Magnetic Fields
In this section, we present a research field rela-
tively new for QM continuum solvation models,
namely, that of the study of solvent effects on
magnetic molecular properties, which determine
nuclear magnetic resonance (NMR) and electronic
paramagnetic resonance (EPR) phenomena. Due to
this “young existence”, many aspects of both formal
and physical nature have yet to be deeply understood,
but the research is very active in many laboratories.
To date, almost all of the studies in this field have
involved two classes of continuum models, that based
on a multipolar expansion (MPE) of the reaction field
and that introducing an apparent surface charge
(ASC). In particular, we note that greatest progress
has been achieved within the ASC approach [in its
different formulations, CPCM, IEFPCM, and SS(V)-
PE]. On the one hand, generalizations to include
more detailed aspects of solvation, such as specific
solute-solvent interactions, through couplings with
discrete pictures have been presented using solute-
solvent clusters obtained through QM minimization
methods or taken from MD simulations. On the other
hand, extensions to more complex environments have
been also realized; we recall here, as an example, the
recent extension of the IEFPCM for anisotropic
dielectrics to the study of NMR properties of solutes
immersed in liquid crystalline solvents.102,701 Both of
these extensions indicate new perspectives of ap-
plications of continuum models toward systems of
increasing complexity, heretofore considered as com-
pletely prohibitive for a continuum description.
To better quantify the importance that QM studies
of solvent effects on magnetic properties are rapidly
achieving, in Table 3 we try to summarize recent
papers in which continuum solvation models have
been used.
6.3.1. Nuclear Shielding
The effects of solvent on nuclear magnetic shielding
parameters derived from NMR spectroscopy have
been of great interest for a long time. In 1960
Buckingham et al.743 suggested a possible classifica-
tion in terms of various additive corrections to the
shielding arising from (i) the bulk magnetic suscep-
tibility of the solvent, (ii) the magnetic anisotropy of
the solvent molecules, (iii) van der Waals interac-
tions, and (iv) long-range electrostatic interactions.
In the original scheme, strong specific interactions,
such as those acting in intermolecular hydrogen
bonds, were not specifically dealt with but just
mentioned as a possible extreme form of the electro-
static, or, more generally “polar”, effect; in the
numerous applications that followed Buckingham’s
classification, however, this further effect has been
always included as a separate contribution.
On the basis of such a classification an empirical
approach based on the so-called solvent empirical
parameters was formulated to evaluate solvent ef-
fects on nuclear shieldings. In brief, this approach,
originally proposed by Kamlet, Taft, and co-work-
ers,744 does not involve QM or other types of calcula-
tions but introduces a numerical treatment of ex-
perimental data obtained for a given reference system
to obtain an estimate of solvent effects on various
properties. Its generalization to the study of the
solvent effect on nitrogen nuclear shielding in various
compounds has been proposed by Witanowski et
al.745-749
In parallel to these semiclassical analyses, recently
we have observed growing attention to an explicit
description of the electronic aspects of the solvent
effects on NMR properties and in particular on the
nuclear shielding. This change of perspective has
been made possible by the large development of QM
solvation models we have described in section 2,
which have been coupled to QM methodologies ini-
tially formulated for isolated systems.
The theory of chemical shielding was originally
developed many years ago,750,751 but only later have
ab initio methods752-756 and density functional theo-
ries (DFT)757,758 been reliably used for the prediction
of chemical shielding for molecular systems.
Because these recent developments greatly im-
proved the accuracy of the calculation of the nuclear
shieldings, the following step has been the inclusion
of solvent-induced shifts. This has been achieved by
extending the cited theoretical methods to continuum
solvation models. Such extension began by redefining
the components of the shielding tensor ó, as second
derivatives of the free energy functional defined in
section 2.4712
where Ba and Mb (a,b ) x, y, z) are the Cartesian
components of the external magnetic field B and of
Table 3. NMR and EPR Properties (Nuclear
Shielding, Magnetizability, Spin-Spin Coupling,
EPR) Calculated with Continuum Solvation Models
(or Their Hybrid Discrete/Continuum
Generalizations)
solvation model authora refs
MPE Mikkelsen 702-708
Pecul 709, 710
Pennanen 711
ASC Cammi/Mennucci 102, 701, 712-722
(CPCM, Barone/Cossi 723-734
IEFPCM, SS(V)PE) Manalo 735, 736
Chipman 104
Cremer 737, 738
Zaccari 739, 740
Ciofini 741
Rinkevicius 742
a Reference name, not necessarily the first author.
óab
X )
@2G
@Ba@Mb
X
(185)
3066 Chemical Reviews, 2005, Vol. 105, No. 8 Tomasi et al.

organic compounds in solutions were first observed
by Biot, at the beginning of the 19th century. Biot
also introduced the modern definition of “specific
rotatory power” or “specific rotation,” [R], of a liquid
as R/Fl, where R is the measured optical rotation in
degrees, l is the optical path length (or sample cell
thickness) in decimeters, and F is the density of
liquid. For an isotropic dilute solution of a chiral
molecule, the specific rotation [R] (usually measured
at the sodium D line of 589.3 nm) is generally written
as785-787
where k is a numerical constant, M is the molar mass
in g mol-1, and î is in cm-1. In eq 195, f(î) indicates
the local field factor of eq 193 calculated at the proper
frequency î, and â(î) is equal to one-third the trace
of the frequency-dependent electric dipole-magnetic
dipole polarizability tensor having components de-
fined as
In eq 196, 0 and k label ground and excited electronic
states and íel and ímag are the electronic electric and
magnetic dipole moment operators. This sum-over-
states has been cast in the framework of time-
dependent linear response theory for SCF (either HF
or DFT), multiconfigurational SCF, and coupled
cluster wave functions. For SCF wave functions, the
explicit evaluation of the sum-over-states can be
avoided by rewriting âRâ in terms of electric and
magnetic field derivatives of the ground state elec-
tronic wave function, namely788-792
where @¾0/@E and @¾0/@H are frequency-dependent
derivatives of the ground state electronic wave func-
tion with respect to the electric and magnetic fields,
E and H, respectively.
In eq 195, the Lorentz local field factor f(î) of eq
193 has been introduced to account for the solvent-
induced changes in the microscopic (or local) electric
field acting on the chiral molecule with respect to the
macroscopic electric field of the light wave; â(î) also
depends on the solvent, exactly in the same way we
have outlined in the previous section for the CD
rotatory strength. Once more, to properly account for
solvent effects on â(î), a QM description is required.
The first QM study of the solvent effects on â(î)
has been presented by Mennucci et al.600,793 within
the framework of the IEF version of the PCM
solvation model (see section 2.3.1.3). In this study, a
frequency-dependent DFT/GIAO response approach
has been used to compute the electric dipole-
magnetic dipole polarizability: the mixed nature of
â(î), however, requires two response procedures
containing an electric and a magnetic perturbation,
respectively. For the electric perturbation, the inclu-
sion of the additional solvent terms in the coupled
perturbed equations is exactly equivalent to what
was seen in the case of polarizabilities (see section
6.2.2). Due to the imaginary nature of the magnetic
perturbation, solvent-induced terms do not appear
explicitly but only contribute to the first-order expan-
sion term of the Fock operator exactly as for the
nuclear shielding (see section 6.3.1). In this study,
solvent contributions to the electric perturbation were
obtained in terms of solvent charges calculated using
the value of the dielectric constant at the frequency
of the external field (nonequilibrium solvation). In
the general case, this is the sodium D line frequency,
and thus the value for (î) coincides with the so-called
optical dielectric constant, opt, defined as the square
of the refractive index.
The model was first applied to some conformation-
ally rigid chiral organic molecules for a set of polar
and apolar solvents. The predicted variation in
specific rotation for polar solvents was found to be
in excellent agreement with experiment for all of the
molecules. For apolar solvents the agreement was
much poorer and thus, because only electrostatic
solute-solvent interactions were included in the
model, it was concluded that nonelectrostatic effects
may have some importance in determining specific
rotation, at least for aromatic and chlorinated sol-
vents such as carbon tetrachloride, benzene, and
chloroform. As already noted for the ECD, the varia-
tions predicted by the classical Lorentz factor (eq 193)
were found to be inconsistent, both qualitatively and
quantitatively, with those calculated in terms of the
QM IEF model. In subsequent papers, the same
model was successfully applied to the study of optical
rotation of more complex and flexible molecules such
as glucose in aqueous solution794 and paraconic acid
in methanol;795 in both cases it becomes fundamental
to account for the solvent effects both in the property
and in the relative population of the various conform-
ers that contribute to the final OR value.
An alternative methodology to account for solvent
effects on the OR has been presented by Mikkelsen
and co-workers796 using a coupled-cluster approach
combined with a nonequilibrium version of the single-
center multipole based on the spherical cavity dielec-
tric continuum model described in section 2.3.2.
6.4.3. VCD and VROA
Chiral molecules exhibit optical activities also in
their vibrational spectra. Specifically, vibrations as-
sociated with the chiral molecular structure produce
spectral intensities that are dependent on the circular
polarization direction of the excitation light and the
detection optics. This optical activity can be observed
in infrared spectra, as vibrational circular dichroism
(VCD), and in Raman spectra, where it is referred
to as Raman optical activity (ROA).
Both VCD and ROA were discovered in the early
to mid-1970s (see refs 797, 798 and 799, 800, respec-
tively), and since then, they have evolved into mature
fields of research and application. Both can be used
to determine the absolute configuration and solution
conformation of chiral molecules. Each has advan-
tages relative to the other that parallel the relative
advantages of ordinary infrared absorption (IR) and
Raman scattering.
[R] ) kf(î)â(î)î2M-1 (195)
âab(î) )
c
3ðh
Im[∑k*0〈0jíaeljk〉〈kjíbmagj0〉îk02 - î2 ] (196)
âab )
hc
3ð
Im[〈@¾0@Eaj@¾0@Hb〉] (197)
3070 Chemical Reviews, 2005, Vol. 105, No. 8 Tomasi et al.

Instrumentation for measuring VCD spectra and
software for calculating VCD spectra from QM first
principles have been commercially available since
1997. In contrast, commercial instrumentation for
measuring ROA spectra became available only some
years later, and user-friendly software for calculating
ROA spectra is not yet available, even if it will
become so very soon. As a result, there are many
more publications describing the use of VCD for the
determination of absolute configuration than there
are for ROA and, therefore, this section will mainly
focus on VCD.
Application of VCD to the assignment of absolute
configuration and identification of dominant solution
conformations entail comparison of observed spectra
with spectra calculated for a specific configuration
and conformation of the molecule or a suitable
fragment of the molecule. Thanks to the work of
Stephens and co-workers,801-803 these calculations,
which are typically obtained at the DFT level, are
nowadays routinely performed for isolated molecules
and exhibit an accuracy similar to that of other
spectroscopic quantities (see, for example, ref 804).
Less studied is the prediction of VCD spectra for
solvated systems: only a few papers are present in
the literature. These resort to the IEFPCM714,715,805
of section 2.3.1.3 or to the simple Onsager continuum
model (eventually with explicit solvent molecules)806,807
to treat the solute-solvent electrostatic interaction.
To better understand how solvation models have
been applied to the study of solvent effects on VCD
we briefly present the general aspects of their QM
formulation.
VCD intensity is proportional to the rotational or
rotatory strength, the scalar product of the electric
dipole and magnetic dipole transition moments. In
the harmonic approximation, the rotational strength
R is proportional to the scalar product between the
derivative of the electric dipole moment of the
molecule with respect to normal mode displacement
(Qa) and the derivative of the magnetic dipole mo-
ment of the molecule with respect to the nuclear
velocities of the normal mode, expressed in terms of
the conjugate momentum (Pa), namely
Quantum chemistry programs calculate these in-
tensities by taking derivatives of the energy of the
molecule with respect to electric field or magnetic
field and vibrational normal mode displacement or
momentum. The calculations first require optimiza-
tion of the geometry and determination of the normal
modes of vibration. Coupled perturbed Hartree-Fock
or Kohn-Sham methods are used to obtain the
energy derivatives. In the case of the magnetic dipole
derivatives, the method that is incorporated into
commercial software such as Gaussian17 utilizes the
magnetic field perturbation method of Stephens et
al.,801 with GIAO orbitals to eliminate the origin
dependence in the magnetic field derivatives.808
VCD spectra are affected by the solvent in different
ways.
First, solvent can affect geometries and thus vi-
brational frequencies: for this effect the same com-
ments reported in section 6.1.2 for IR spectra are
valid here. Then, solvent can modify VCD intensities.
The calculation of VCD intensities in the presence
of a solvent still relies on the same equations
formulated for isolated molecules, but with some
refinements. As already remarked with regard to
infrared intensities for molecules in solution (see
section 6.1.2), the electric dipole in eq 198 has to be
replaced by the sum of the dipole moment of the
molecule and the dipole moment arising from the
polarization induced by the molecule on the solvent.
The latter takes into account effects due to the field
generated from the solvent response to the probing
field once the cavity has been created (the so-called
“cavity field”). In principle, the magnetic dipole
should be similarly reformulated. However, by as-
suming the response of the solvent to magnetic
perturbations to be described only in terms of its
magnetic permittivity (which is usually close to
unity), it is reasonable to consider that the magnetic
analogue of the electric “cavity field” gives minor
contributions to the rotational strength. A formula-
tion of this solvent-specific term has been given in
ref 805 within the IEF version of the PCM model:
in this case an additional set of charges spread on
the solute cavity surface is introduced in addition to
the PCM apparent charges representing the reaction
field. These additional charges are those produced by
the (Maxwell) electric field of the radiation in the
medium, exactly as in the PCM treatment of the IR
and Raman intensities shown in section 6.1.2.
A further important solvent effect on VD spectra
is implicit, that is, through the possible changes in
the relative energies of different conformers of the
solute molecule when passing from gas phase to
solution. In fact, if the solute can exist in different
conformations that are differently stabilized by the
solvent, the VCD spectra must be obtained by com-
bining population-weighted spectra of all conformers
where the weights are evaluated in solution.805
7. Continuum and Discrete Models
In section 1 we have anticipated one of the most
important extensions that continuum solvation mod-
els have presented in recent years. This extension,
which we denoted layering, has been propelled by two
opposing reasons: namely, the necessity to reduce
the complexity of calculations using discrete ap-
proaches on systems composed by a large number of
molecules and the increasing accuracy required in the
study of solvated systems (especially their molecular
properties).
The layering can be considered a generalization of
focusing. The material component of the models
partitioned into several parts, or layers, and each
layer is defined at a given level of accuracy in the
description of the material system and with an
appropriate reduction of the degrees of freedom.
There are many possibilities for layering, which can
be simply abbreviated 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
Ra ) p
2( @íb@Qa)â(@mb@Pa) (198)
Quantum Mechanical Continuum Solvation Models Chemical Reviews, 2005, Vol. 105, No. 8 3071

combined approaches mixing different QM or QM and
MM descriptions will be reviewed in the next sec-
tions.
7.3.1. Solvated Supermolecule
The coupling of the supermolecule approach with
a continuum is an attractive technique that has
proved to be extremely useful in studying solvent
effects on electronic and vibrational spectra and
response molecular properties as well in studying
reaction mechanisms and energetics in solution; in
section 6 we reported several references to which we
add here a few more.842-848
One of the main reasons to explicitly treat some
solvent molecules is to get an accurate description
of strong specific interactions (e.g., hydrogen bond-
ing); however, such an approach is not straightfor-
ward, and its correct realization requires a detailed
analysis on the nature and the strength of such
interactions.
In the case of very strong solute-solvent interac-
tions such as solutes with strong H-bond acceptor
centers in water or in alcohols, a sufficient description
is often obtained by including all of the solvent
molecules necessary to saturate such strong acceptor
centers (in general, a few solvent molecules are
enough) in the QM optimization procedure. The
resulting structure is then used to compute the
corresponding vibrational spectrum or the properties
of interest or just to compute the thermal free energy
necessary to get the energetics of a reaction. In all of
the steps of this strategy an external continuum is
used to account for long-range (or bulk) effects.
When we are in the presence of weaker interac-
tions, however, a representation of the supermolecule
in terms of a single rigid structure obtained as the
minimum of the potential energy surface of the
cluster cannot be adopted. The real situation is in
fact dynamic and a variety of different representative
structures can and do occur. A possible way to get
such a picture is to consider structures derived from
either classical or ab initio MD shots taken at
different simulation times: for each of these selected
structures the number of solvent molecules to be used
in the supermolecule calculation is determined by a
threshold imposed in the distance between the solute
H-bond acceptor/donor center and the corresponding
donor/acceptor center in the solvent molecules. Each
of the resulting clusters (generally involving more
than the solvent molecules really needed to saturate
H-bonding centers) is then embedded in the external
continuum and such a solvated supermolecule is
studied at the proper QM level. In the most recent
examples, the PCM method (either in the IEF or
CPCM version, see section 2.3.1) is exploited to
include bulk effects,716,720-722,729,731,733 whereas older
studies849-852 adopted an image reaction field plus
exclusion model to represent the polarization of the
dielectric material embedding the clusters.853
In this approach, no geometry optimizations of the
clusters are required; this of course makes the
procedure easier from a computational point of view
but not necessarily faster as an average among all
of the selected structures has to be done to get a
reliable description. A disadvantage of this procedure
is that it cannot be used to directly compute vibra-
tional spectra; more complex analyses are necessary
as the structures do not correspond to real minima.
This approach combining (but not coupling) MD
simulations and averaged QM calculations on se-
lected clusters is not limited to the examples cited
above, and they do not necessarily include an exter-
nal continuum; indeed, it has been used in many
different studies as an alternative to continuum
models. The main difference is that without including
the bulk effects through a continuum dielectric, very
large clusters are required to get the complete picture
of solvation, and thus the QM level of the theory is
usually limited to semiempirical approaches or low-
level ab initio methods.854-862
7.3.2. QM/MM/Continuum: ASC Version
Another example involving discrete and continuum
approaches describes the solute at QM level, the first
solvation shells with MM solvents, and the rest with
a dielectric continuum. In such a way, both the first
solvation shell effects (such as short-range repulsion
and dispersion interactions and solvent entropic
terms related to hydrophobic interactions) and long-
range electrostatic effects in the bulk can be taken
into account with modest computational cost. The
small number of explicit solvent molecules makes the
calculations of quantities that require configurational
samplings much faster compared to a full-scale
microscopic simulation.
In principle, all of the continuum models we have
described in section 2.3 can be coupled with a QM/
MM system. Here, however, we shall limit the
exposition to those belonging to the ASC class of
continuum models.
An example of this approach is the model developed
by Cui,863 in which the solute and a number of solvent
molecules (“solute-solvent cluster”) are described
with QM/MM and the bulk solvent is treated with
the IEFPCM model (see section 2.3.1.3). In this
model, both the MM atoms, represented by fixed
partial charges, and the QM atoms are contained in
a generalized cavity embedded in the continuum.
Following the same formalism used in section 2.4,
the total electrostatic free energy of the “solvated”
QM/MM system thus becomes
where h÷ ) h + J + Y, G÷ ) G + X, and J, Y, and X
are the PCM matrices described in section 2.4.3 (see
eqs 64-67). The terms JQM/MM, YQM/MM, UNN
QM/MM,
UNN
MM/QM, and UNN
MM represent the additional electro-
static interactions between the QM/MM part and the
induced surface charges. There are four types of
interactions, namely, between QM electrons and MM-
induced surface charges, between MM and QM
electron-induced surface charges, between QM nuclei
Gel ) tr P[h÷ + hQM/MM + 12(JQM/MM + YQM/MM)] +
1
2
tr PG÷ + (V÷ NN + EelMM + VNNQM/MM + 12 UNNQM/MM +
1
2
UNN
MM/MM +
1
2
UNN
MM) (207)
3076 Chemical Reviews, 2005, Vol. 105, No. 8 Tomasi et al.

and MM-induced surface charges, between MM and
QM nuclei-induced surface charges, and between MM
and MM-induced surface charges.
This QM/MM/PCM method has been implemented
into the simulation package CHARMM,864 which was
previously interfaced with GAMESS for ab initio QM/
MM calculations,865 and it can be used for single-point
solvation free energy calculations, Monte Carlo simu-
lations, and geometry optimization as well as MD
calculations for reactions in solution.
Another QM/MM/Cont approach is that of coupling
continuum models with the QM/MM model developed
by Gordon and co-workers, known as the effective
fragment potential (EFP) method.426,866,867 The basic
idea behind this method is to replace the chemically
inert part of a system by EFPs while performing a
regular ab initio calculation on the chemically active
part. A simple example of an active region might be
a solute molecule, with a surrounding spectator
region of solvent molecules represented by fragments.
The charge distribution of the fragments is repre-
sented by an arbitrary number of charges, dipoles,
quadrupoles, and octupoles, which interact with the
ab initio Hamiltonian as well as with multipoles on
other fragments. An arbitrary number of dipole
polarizability tensors can be used to calculate the
induced dipole on a fragment due to the electric field
of the ab initio system as well as all other fragments.
These induced dipoles interact with the ab initio
system as well as the other EFPs, in turn changing
their electric fields. All induced dipoles are therefore
iterated to self-consistency. The combination of the
EFP method with continuum solvation initially used
an Onsager-like model868 but subsequently has been
extended to the IEF-PCM model.869 Using the same
formalism presented in section 2.4.2, the electrostatic
component of the free energy of the solute (ab initio
+ EFP) continuum system can be written
where for h÷ , G÷ (P), VNN, an UNN we have used the
same notation as that in eq 207. All of the other
terms are due to the coupling between the EFP and
PCM methods. In particular, (Jstat
EFP)íî ) ∑qstatEFP Víîe
and (Jpol
EFP)íî ) ∑qpolEFPVíîe , where Víîe are the electro-
static potential integrals in the atomic orbital basis
and qstat
EFP) and qpol
EFP are the ASC due to the static
multipoles and the induced dipoles of the fragments,
respectively. In the remaining terms, VEFP denotes
the electrostatic potential induced by the EFP frag-
ments and the subscript indicates the source of such
potential, namely, the static multipoles (stat) and the
induced dipoles (pol).
Subsequently, the EFP/PCM approach has been
reformulated870 using a simultaneous iterative solu-
tion of the QM self-consistent field (SCF) and of the
electrostatic equations (see also section 2.4.3): in this
way, bulk solvation of large solutes can be efficiently
modeled.
7.3.3. ONIOM/Continuum
An alternative formulation of the coupling of QM,
MM, and continuum approaches is given by energy
subtraction methods. The latter are a very general
class of methods in which calculations are done on
various regions of the molecule with various levels
of theory, and the energies obtained at each level are
finally added and subtracted to give suitable correc-
tions. The introduction of continuum solvation into
this kind of method has been realized by coupling the
IEFPCM model we have described in section 2.3.1.3
with the class of subtraction methods developed by
Morokuma and co-workers871-875 and generally known
by the acronym ONIOM. This acronym actually
covers different techniques, combining either two (or
more) different orbital-based techniques or orbital-
based technique with an MM technique.
The concept of the ONIOM methods is extremely
simple. The target calculation is the high level
calculation for a large real system, E(high, real),
which is prohibitively expensive. Instead of doing
such a calculation, an inexpensive low-level calcula-
tion leading to the energy E(low, real) is performed
together with two calculations (one at accurate high
level and the other at the same low level) for a
smaller part of the system, usually indicated as the
model system, leading to the energies E(high, model)
and E(low, model), respectively. Starting from E(low,
model), if one assumes the correction for the high
level, E(high, model) - E(low, model), and the cor-
rection for the real system, E(low, real) - E(low,
model), to be additive, the energy of the real system
at the high level can be estimated from three inde-
pendent calculations as
Clearly this approach is completely general, and it
can be straightforwardly extended to more than two
layers. In addition, it has the advantage of not
requiring a parametrized expression to describe the
interaction of various regions. Any systematic errors
in the way that the lower levels of theory describe
the inner regions will be canceled out.
The combination of ONIOM with IEFPCM has
been achieved in four alternative ways.876 The result-
ing computational schemes differ mainly with respect
to the level of coupling between the solute charge
distribution and the continuum dielectric, which has
important consequences for the computational ef-
ficiency. All of the schemes share the fact that only
one cavity is defined, based on the geometry of the
real system, which is subsequently used for all three
subcalculations. Three of the four schemes can be
placed in a hierarchy. In the highest method in the
class (called ONIOM-PCM/A in the reference paper),
all three subcalculations exploit the reaction field
that is obtained self-consistently by employing the
ONIOM integrated density. Obviously this is the
E ) E(low, real) + E(high, model) -
E(low, model) (209)
G ) [tr Ph÷ + 12PG÷ (P) + VNN + 12 UNN] +
tr P(Jstat
EFP + Jpol
EFP) + 1
2 ∑qstatEFP VstatEFP +
1
2∑qpolEFP VpolEFP + ∑qstatEFP VpolEFP + ∑qNVstatEFP +
∑qNVpolEFP (208)
Quantum Mechanical Continuum Solvation Models Chemical Reviews, 2005, Vol. 105, No. 8 3077

previous versions of the model), other parameters of
the LD solvation model had also to be reoptimized.
These changes include introduction of the field-
dependent grid points and surface-constrained di-
poles on the outer solvent surface, the new relation-
ship between magnitudes of dipoles placed at the
inner and outer grid points, reparametrized van der
Waals and hydrophobic terms, increased extent of
dipole-dipole interactions, and a newly added solute-
polarization term.
The LD model has been extended to study proteins,
and the corresponding version of the model is known
as protein dipoles Langevin dipoles (PDLD).446,893
The computational code of the LD model (called
ChemSol894) has been also implemented in the Qchem
quantum package.894,895
7.4. Other Methods
The electrostatic solvation methods examined in
section 2 are all based on averaged solvent polariza-
tion, expressed in terms of a response function in
which no use is made of the detailed microscopic
description of the solvent. This definition of the
solvent mean field greatly reduces the computational
task but neglects the possible presence of specific
solute-solvent interactions.
We pass now to examine other approaches in which
some elements of the discreteness of solvent mol-
ecules are considered.
We shall limit our attention to two approaches that
still use the mean solvent field approximations, but
from a mean field obtained from a thermally aver-
aged discrete description of the solvent. The two
approaches are based on different formulations of this
averaged description and will be separately treated.
They will be indicated with the acronyms ASEP
(averaged solvent electrostatic potential) and RISM-
SCF (reference site model in a SCF version), but
actually different versions of the two approaches have
been elaborated, and we shall review all of the
published versions to give a fuller idea of the work
in progress in this specific field.
The problem both approaches are addressing is of
finding a procedure able to reduce the strong com-
putational burden due to the exceedingly large
number of solvent freedom degrees a standard dis-
crete description presents.
7.4.1. ASEP-MD
The method we shall consider here has been
elaborated by Olivares del Valle and Aguilar, with
co-workers, in Badajoz.
The starting point is a standard MD calculation:
the solute and solvent coordinates were dumped at
every N (with N properly defined) steps for further
analysis. The resulting configurations are then trans-
lated and rotated in such a way that all of the solvent
coordinates can be referred to a reference system
centered on the solute mass center with the coordi-
nate axis lying along the principal axes of inertia of
the solute. This process is necessary to define a grid
of points equal for all configurations.
For each stored configuration only those solvent
molecules that are within a sphere centered on the
solute molecule and radius equal to half the size of
the simulation box are selected. The solute and the
selected solvent molecules are then included in a
dielectric continuum. The Laplace equation that gives
the polarization of the continuum is solved numeri-
cally, making use of apparent charges (i.e., through
the DPCM approach presented in section 2.3.1.1).
The total electrostatic potential due to the permanent
solvent charges, induced dipoles, and the surface
charge density is calculated at each point of the grid
defined inside the volume occupied by the solute
molecule, and the averaged value of this solvent
electrostatic potential (ASEP) is obtained by averag-
ing over the configurations. The potential obtained
in this way is a numerical potential and hence is not
very suitable for use in quantum calculations; in-
stead, a set of charges chosen in such a way that they
reproduce the value of the ASEP inside the solute
volume is introduced. The charges are usually ar-
ranged on spherical shells centered on the solute
mass center (a solvent diameter apart), but shells
adapted to the molecular shape of the solute and
defined in terms of intersecting spheres centered on
the solute atoms have also been used. Charges placed
on two concentric shells (spherical or molecularly
shaped) are generally sufficient to represent the effect
of the total solvent electrostatic potential. The energy
and the wave function of the solute molecule in the
presence of this averaged solvent configuration are
thus obtained by solving the Schro¨dinger equation
associated with the perturbed Hamiltonian, namely
where HQM is the usual vacuum Hamiltonian of the
solute molecule. The interaction term, HQM-MM, takes
the form
where F is the charge density operator associated
with the solute molecule and 〈V(r)〉 is the value of
ASEP at the point r.
There are different versions of the ASEP procedure.
In the early versions a single QM calculation is
performed with unpolarized MM descriptions of the
solvent molecules896 or with MM molecules including
polarization.897 More recent versions introduce itera-
tions on the QM part, allowing solute polarization.898
These versions have been called by the authors
“coupled methods”. There are two versions: the
partial coupling (PC) in which the MM solvent
molecules are kept unpolarized, and the complete
coupling (CC) in which all of the molecules of the
liquid specimen are polarized at the same level. The
CC version is at present available for pure liquids
only.
The computational efficiency of ASEP has been
checked by comparisons with QM/MM codes;899 the
ASEP method turns out to be more efficient by some
orders of magnitude. The most important advantage
of ASEP-MD when compared with previous QM/MM
calculations is that the number of quantum calcula-
tions is considerably less: around four to eight with
(HQM + HQM-MM)j¾〉 + Ej¾〉 (210)
HQM-MM ) s dr F〈V(r)〉 (211)
Quantum Mechanical Continuum Solvation Models Chemical Reviews, 2005, Vol. 105, No. 8 3079

framework the operator has the form
where the sum runs over the solvent apparent
charges. In the RISM-SCF approach the parallel
operator has the form
where F is the number density, bì is a population
operator, gìR(r) is a pair correlation function between
ì (solute) and R (solvent), and qR is the partial charge
on the solvent site R.
The main difference between the two operators
(besides the obvious differences in the form) is that
the RISM potential (also indicated as microscopic
mean field VRISM) includes molecular level informa-
tion as well as bulk properties of the solvent. To
better understand this point, we briefly summarize
the RISM-SCF main computational strategy, which
can be profitably presented as a sequence of steps.
In advance of the RISM-SCF cycle, the solvent
correlation function is prepared by solving the RISM
integral equations through a standard iterative
scheme with a set of fixed intermolecular potential
functions and with guess partial charges of the solute
molecule as variables. Then using such a correlation
function, the solute-solvent pair correlation func-
tions are obtained. Immediately after the RISM
convergence has been attained, the electrostatic
solvent operator is calculated from eq 214 and
inserted in the solvated Fock operator to calculate
the solute wave function. Then the electrostatic
contributions due to the nucleus and electronic cloud
of the molecule are optimized independently with
those due to partial charges assigned on each inter-
action site, using the least-squares fitting technique.
Normalization constraints to preserve the total charge
for each contribution are imposed during the fitting
procedure. The RISM-SCF iteration is continued
until both the electronic and solvent structures
become self-consistent within given convergence cri-
teria.
The RISM-SCF method has been extended to
include analytical gradients for geometry optimiza-
tions and to the multiconfigurational self-consistent-
field (MCSCF) method,922 which can be used for
exploring the excited state of a molecule in solution.
Sato et al.925 have also reformulated the RISM-
SCF/MCSCF for the 3D-RISM formalism to properly
include the three-dimensional picture of the solvation
structure necessary for complex solutes. The authors
found that this reformulation allows one to reliably
resolve locations and directions of hydrogen bonding
in the hydration shells. At the same time, however,
they also found that the results from the original
RISM-SCF/MCSCF method are reasonably similar to
those following from the 3D-RISM/SCF approach
after reduction of the orientational dependence.
In recent years, applications of RISM-SCF to the
calculation of NMR chemical shifts of solvated mol-
ecules have also been presented926,927 as well as
studies on the solvent reorganization energy associ-
ated with a charging process of organic com-
pounds.928,929 In the latter case both nonequilibrium
solvation and nonelectrostatic effects were included.
In all of the above RISM-SCF versions the short-
range interactions are described by site-site Len-
nard-Jones potentials, neglecting the solute-solvent
short-range QM effects such as the exchange repul-
sion coming from the Pauli exclusion principle and
the charge transfer interaction. To overcome these
limitations, Yoshida and Kato930 proposed an elec-
tronic structure theory in solution using the molec-
ular Ornstein-Zernike equation, which is referred
to as the MOZ-SCF method. In their treatment, the
exchange repulsion/charge transfer interactions are
incorporated for calculating the solute-solvent in-
teractions by introducing an effective potential lo-
cated on solvent molecule. Comparing the MOZ-SCF
results with those from the RISM-SCF and PCM
methods, it was found that the MOZ-SCF exerts an
intermediate effect on the solute electronic state
between the RISM-SCF and PCM.
Subsequently, Sato et al.931 proposed a new ap-
proach to account for the quantum nature in the
short-range interaction between solute and solvent
in solution. Unlike the RISM-SCF theory, this ap-
proach describes not only the electronic structure of
solute but also solute-solvent interactions in terms
of a kind of model-potential method based on the
Hartree-Fock frozen density formulation. In the
treatment, the quantum effect due to solvent, includ-
ing exchange repulsion, is projected onto the solute
Hamiltonian using the spectral representation method,
whereas the integral equation theory of liquids is
employed to calculate the solvent distribution around
solute.
8. Concluding Remarks
To conclude this review we express some consid-
erations about the present state of the art of con-
tinuum models and about perspectives for the future.
The introductory remarks of section 1 represent a
background for the whole review, which is centered
on QM versions of the continuum solvation models.
We shall comment later what this background means
for the perspectives of development.
In section 2 and 3 we have examined, with abun-
dant details, the elaboration of the models for solva-
tion effects in very dilute solutions (discarding non-
linearities and time-dependent effects). Within this
restrictive definition of solvation, we may consider
the elaboration (computational as well as structural)
of models given in these two sections as having
reached maturity. By comparing the present exposi-
tion with that done in the 1994 review,1 one notices
a remarkable increase in the number of methods and
computational codes, now far more efficient and
accurate. Several new original procedures have also
been formulated, almost all completely elaborated
and introduced in distributed codes.
V(rb) ) ∑
i
qi
jrb - sbij
(213)
VRISM(r) ) F ∑
ì,R
bìqR s 4ðr
2
gì,R(r)
r
dr (214)
Quantum Mechanical Continuum Solvation Models Chemical Reviews, 2005, Vol. 105, No. 8 3081

Most of the published applications of these methods
(that we have only partially reported in this review)
are about the energy of the focused subsystem and
related quantities, such as solvent effects on the
geometry. These studies address a large area of
chemical applications: conformational changes, reac-
tion mechanisms, dissociation/association equilibria,
and solvent partition functions. The number of such
applications, generally claiming success, is quite
large; we are not able to give an estimate of their
number, but we know of more than 1000 papers, and
surely there are many more.
In contrast, less attention has been paid to extract-
ing from the QM calculations elements for analysis
and interpretation of solvent effects. Theoretical
chemistry has elaborated a large number of concepts
and computational tools in terms of submolecular
units for the analysis and interpretation of structures
and phenomena for systems in the gas phase; these
tools are easily extended to solvated systems within
a continuum representation. The inclusion of solvent
effects adds a new dimension to the understanding
of chemical structures and processes, which is rich
in possible suggestions. We give an example drawn
from a topic not treated in this review: the study of
reaction mechanisms. Differential solvation forces
related to the specific nature of chemical groups may
accelerate or retard a chemical reaction. This has
been shown in one of the first examples of a detailed
ab initio study of an organic reaction,932 in which the
forces exerted by solvation on the reacting chemical
groups (expressed in terms of localized orbitals) were
explicitly derived and used to give a rationale for
nuclear geometry and electron distribution changes
near the transition state. This example, however, has
not been followed by other researchers, perhaps
because the necessary codes were not collected into
a distributed package. The elaboration and distribu-
tion of easy-to-use interpretation packages is one of
the aspects of conventional energy-based solvation
methods for which a more decisive effort is called for
and on which we encourage other researchers to
intervene.
Section 4 collects under the heading “Nonunifor-
mities in the Continuum Model” topics of different
origins for which there is still the need of further
methodological elaboration. For this reason we have
given bibliographic references also related to the
“history” of the subject. To do further development
in an already considered field, it is in fact advisable
to examine the sources; some ideas given in seminal
papers are often neglected in the immediately fol-
lowing extensions and applications, but they can
become valuable after some years. The spatial non-
uniformities considered in that section are of local
nature, occurring at some boundaries in the system.
These characteristics prompted the extension of the
review from the dielectric nonuniformities around
ions, to which the “historical” exposition mainly
refers, to other boundaries, enlarging the field of
systems where continuum methods may be applied.
Some among them are related to systems to which
continuum methods, mostly in classical versions,
have long been applied such as polyions and molec-
ular polymers in general; others involve heteroge-
neous systems in which the boundary involves a
phase with a structure more stable in time than
normal liquids.
We have given a few examples of successful ap-
plications of QM continuum methods to systems of
this type, but clearly these examples simply indicate
the beginning of methodological studies that still
have to be developed. We have also indicated another
field in which local nonuniformities may play a role,
namely, that of the complex liquids, characterized by
a mesoscopic length scale. This is a novel field to
which the application of continuum models has not
yet attempted, an attempt which seems to us possible
and promising. In the same line of referring to
material systems to which some basic methodological
issues could be applied, we cite here the ionic liquids,
and the isotropic solutions at finite rather than
negligible concentration.
In conclusion, it may be said that for the topics of
section 4 the state of the art indicates a situation in
progress, not yet assessed, but with brilliant perspec-
tives.
Time-dependent (TD) phenomena and response
properties have been considered in sections 5 and 6,
respectively. Remarkable advances have been regis-
tered in both areas in recent years. In particular, for
TD phenomena, continuum methods have started to
represent a valid alternative to simulation methods
(historically the best suited for this kind of studies),
and for response properties they are slowly but
progressively acquiring the reliability of gas phase
calculations. In both cases, an important discriminat-
ing factor has been the easy insertion of high-level
QM descriptions of a molecule in continuum models.
The computational advantage of a continuum ap-
proach with respect to simulation methods must be
emphasized again: simulations have to dedicate a
considerable portion of computational time to de-
scribe solvent degrees of freedom, and this, for
example, makes almost impossible a serious study
of response properties.
Concerning the description of molecular properties
in condensed phases (both in continuum and in
discrete approaches) some problems arise. Quantum
mechanics gives a complete (in principle) description
of the studied system which satisfies well-established
rules: completeness of the set of eigenfunctions,
complete orthogonality, expansion coefficients over
the states directly related to the probability of finding
that state in a measurement experiment, etc. These
formal properties, however, are no longer valid when
the molecule interacts with others. In other words,
the molecular wave function has a different status
in the presence or absence of molecular interactions.
This fact is well-known to people trying to derive the
value of a property of one of the molecular partners
from calculations performed on a supermolecule. To
obtain this value there is the need of additional
assumptions, varying from case to case, because a
universal protocol has never been established. It may
be said that this protocol exists for continuum
models, being the Hartree separation between the
focused molecule and the continuum in general use.
3082 Chemical Reviews, 2005, Vol. 105, No. 8 Tomasi et al.

This model can be good or bad, but in any case it is
a common starting point on which additional refine-
ments can be based.336
This, however, is not the only important issue to
address. In condensed systems exhibiting rapid mo-
tions of the molecular components, such as liquids,
other complications arise. There is, in fact, the need
of an appropriate thermal average, which may change
from property to property. For example, the calcula-
tion of electronic excitation energies in the simple
model of vertical transitions gives poor results if
based on the excited states of a molecule in equilib-
rium with the solvent. The inclusion in the model of
changes in the fast component of the solvent polar-
ization leads to excited states of better quality.
Other examples could be done, all indicating that
the description of properties in condensed systems
is more complex than for isolated molecules. This fact
has been known for a long time, and experimental
values of properties are always published with em-
pirical corrections addressing the effect of the sur-
rounding medium (namely, the local field factors
introduced in section 6.2). These corrections, elabo-
rated many years ago, are all based on quite simple
continuum models. As shown in section 6.2.2, re-
cently more precise expressions of these corrections,
tailored on the single molecule, have been elaborated
with the aid of modern solvation methods.
In parallel to this evolution, the number of proper-
ties for which there are available detailed and tested
computational strategies has greatly increased, as a
comparison of this review with the preceding ones
shows. This rate of increase shows no tendency to
level off, because there are many properties not yet
examined in detail and because their number in-
creases thanks to the ingenuity of experimentalists
who continually devise new ways of probing the
molecular matter.
The expansion of continuum methods in the de-
scriptions of properties has now passed the divide
between simple systems (such as homogeneous solu-
tions) and complex systems, in which the whole
system gives responses that cannot be reduced to
those of the separate components. Even in this new
field the positive results are numerous and give
confidence in further rapid progresses. The realm of
complex systems is quite interesting, being at the
frontier of physical and chemical research and, for
its variety, presenting a challenging variety of fea-
tures each requiring an attentive analysis for its
modeling.
We do not enter into an analysis of the role of the
calculations of properties in the intrinsic context of
theoretical chemistry and of its support to applica-
tions in academic and technological research, because
our views have already been expressed in several
papers.336,933-935 For the scope of this section, it is
sufficient to say that, in our opinion, the state of the
art and the prospects for the topics considered in the
section 6 are decidedly positive. The basic method-
ological elements have been established, the compu-
tational methods (now well documented) have been
validated, and they can be safely employed, even if
there is still a lot of work to refine.
Section 7 is focused on another aspect of the
extension of continuum methods which is quite
important and in progress, namely, the combination
of continuum and discrete representations of the
medium.
The admixture of continuum and discrete models
(using a solute-solvent supermolecule) has been in
use for a long time in different studies including the
description of solute conformations and the analysis
of the direct role of solvent molecules in chemical
reaction mechanisms. We explicitly cite a paper936
published the same year as the first PCM paper8 to
underscore that the use of supermolecules has been
considered in QM continuum methods since their
very beginning, and it has continued to the present.
In that paper, a supermolecule approach plus con-
tinuum was used to compute solvent shifts in elec-
tronic absorption spectra. Nowadays, it is clear that
small solute-solvent clusters can represent an im-
portant evolution of solvation methods to accurately
describe spectroscopic properties (as also discussed
in section 6). There is, howevet, the need to ac-
company this approach with others in which impor-
tant aspects of the material model, such as the
averages on the motions, are better described. There
is a sizable variety of computational strategies ad-
dressing this mixing of continuum and discrete
methods, and some protocols of use are in the process
of refinement. In parallel, discrete simulations, in
their struggle against computational costs, tend to
exploit more extensively the computational advan-
tages of continuum descriptions. This is an example
of admixture of approaches from which our advance-
ment in the description of matter may derive great
benefits.
The perspective of a stronger admixture between
continuum and discrete models is in line with what
was pointed out in section 1.2 about distribution
functions. In the present case, the main continuum
distribution in question is the density distribution of
massive particles, the molecules, interacting among
themselves via the forces described by the intermo-
lecular interaction theories. The dielectric function
is one of the ancillary quantities related to the
molecular density distribution; other functions, of
different physical origin, exist, but they do not play
a role in the present continuum solvation methods.
The discrete molecule distribution can be replaced,
partially or totally with continuum functions, as we
have discussed in several points of this review. This
change of description leads us to neglect a portion of
the information carried by a discrete description. One
needs to evaluate, problem by problem, whether this
reduction in the information compromises the de-
scription of the phenomenon or whether the reduction
in the computational burden is not accompanied by
a loss of important details. In the case of solvation
methods, as we have shown, the advantages of a
continuum description are in general greater than
the disadvantages. The topics considered in section
7 to a large extent address the area in which
continuum methods require discrete corrections.
In section 1.2 we have considered another distribu-
tion, that of electrons in molecules interacting via
Quantum Mechanical Continuum Solvation Models Chemical Reviews, 2005, Vol. 105, No. 8 3083

(406) Abraham, M. H.; Liszi, J. J. Chem. Soc., Faraday 1 1978, 74,
2858.
(407) Abraham, M. H.; Liszi, J. J. Chem. Soc., Faraday 1 1980, 76,
1219.
(408) Abraham, M. H.; Liszi, J.; Papp, E. J. Chem. Soc., Faraday 1
1982, 78, 197.
(409) Abraham, M. H.; Liszi, J.; Kristo´f, E. Aust. J. Chem. 1982, 35,
1273.
(410) Abraham, M. H.; Matteoli, E.; Liszi, J. J. Chem. Soc., Faraday
1 1983, 79, 2781.
(411) Basilevsky, M. V.; Soudackov, A. V.; Voronin, A. I. Chem. Phys.
1998, 235, 281.
(412) Basilevsky, M. V.; Rostov, I. V.; Newton, M. D. Chem. Phys. 1998,
232, 189.
(413) Leontyev, I. V.; Basilevsky, M. V.; Newton, M. D. Theor. Chem.
Acc. 2004, 111, 110.
(414) Buckingham, A. D. Discuss. Faraday Soc. 1957, 21, 151.
(415) Frank, H. S.; Wen, W. Y. Discuss. Faraday Soc. 1957, 24, 133.
(416) Goldman, S.; Bates, R. G. J. Am. Chem. Soc. 1972, 94, 1476
(417) Tremaine, P. R.; Goldman, S. J. Phys. Chem. 1978, 82, 2317.
(418) Pitzer, K. S. J. Phys. Chem. 1983, 87, 1120.
(419) Tanger, J. C.; Pitzer, K. S. J. Phys. Chem. 1989, 93, 4941.
(420) Floris, F.; Persico, M.; Tani, A.; Tomasi, J. Chem. Phys. Lett.
1992, 199, 518.
(421) Floris, F. M.; Persico, M.; Tani, A.; Tomasi, J. Chem. Phys. Lett.
1994, 227, 126.
(422) Floris, F.; Persico, M.; Tani, A.; Tomasi, J. Chem. Phys. 1995,
195, 207.
(423) Floris, F. M.; Martinez, J. M.; Tomasi, J. J. Chem. Phys. 2002,
116, 5448.
(424) Floris, F. M.; Martinez, J. M.; Tomasi, J. J. Chem. Phys. 2002,
116, 5460.
(425) Chillemi, G.; Barone, V.; D’Angelo, P.; Mancini, G.; Persson, I.;
Sanna, N. J. Phys. Chem. B 2005, 109, 9186.
(426) Jensen, J. H.; Day, P. N.; Gordon, M. S.; Basch, H.; Cohen, D.;
Garmer, D. R.; Kraus, M.; Stevens, W. J. Effective fragment
method for modeling intermolecular hydrogen-bonding effects
on quantum-mechanical calculations. In Modeling the Hydrogen
Bond; ACS Symposium Series 569; American Chemical Soci-
ety: Washington, DC, 1994.
(427) Merrill, G. N.; Webb, S. P. J. Phys. Chem. A 2003, 107, 7852.
(428) Merrill, G. N.; Webb, S. P.; Bivin, D. B. J. Phys. Chem. A 2003,
107, 386.
(429) Lee, S.; Kim, J.; Park, J. K.; Kim, K. S. J. Phys. Chem. 1996,
100, 14329.
(430) Katz, A. K.; Bock, C. W.; Glusker, J. P. J. Am. Chem. Soc. 1996,
118, 5752.
(431) Pavlov, M.; Siegbahn, P. E. M.; Sandstro¨m, M. J. Phys. Chem.
A 1998, 102, 219.
(432) Pye, C. C.; Rudolph, W. W. J. Phys. Chem. A 1998, 102, 9933.
(433) Rudolph, W. W.; Mason, R.; Pye, C. C. Phys. Chem. Chem. Phys.
2000, 2, 5030.
(434) Richens, D. T. The Chemistry of Aqua Ions; Wiley: Chichester,
U.K., 1997.
(435) Tunon, I.; Rinaldi, D.; Ruiz-Lopez, M. F.; Rivail, J. L. J. Phys.
Chem. 1995, 99, 3798.
(436) Martinez, J. M.; Pappalardo, R. R.; Sanchez Marcos, E.; Men-
nucci, B.; Tomasi, J. J. Phys. Chem. B 2002, 106, 1118.
(437) Sandberg, L.; Casemyr, R.; Edholm, O. J. Phys. Chem. B 2002,
106, 7889.
(438) Sandberg, L.; Edholm, O. J. Chem. Phys. 2002, 116, 2936.
(439) Jorgensen, W. L.; Maxwell, D. S.; TiradoRives, J. J. Am. Chem.
Soc. 1996, 118, 11225.
(440) Rizzo, R. C.; Jorgensen, W. L. J. Am. Chem. Soc. 1999, 121, 4827.
(441) Sandberg, L.; Edholm, O. J. Phys. Chem. B 2001, 105, 273.
(442) Berendsen, H. J. C.; van der Spoel, D.; van Drunen, R. Comput.
Phys. Commun. 1995, 91, 43.
(443) Gilson, M. K.; Rashin, A.; Fine, R.; Honig, B. J. Mol. Biol. 1985,
183, 503.
(444) Harvey, S. C. Proteins 1989, 5, 78.
(445) Antosiewicz, J.; McCammon, J. A.; Gilson, M. K. J. Mol. Biol.
1994, 238, 415.
(446) Schutz, C. N.; Warshel, A. Proteins 2001, 44, 400.
(447) Bjerrum, N. Z. Phys. Chem. Stoechiom. Verwandschaftsl 1923,
106, 219.
(448) Kirkwood, J. G.; Westheimer, F. H. J. Chem. Phys. 1938, 6, 506.
(449) Simonson, T.; Perahia, D. Comput. Phys. Commun. 1995, 91,
291.
(450) Simonson, T.; Brooks, C. L. J. Am. Chem. Soc. 1996, 118, 8452.
(451) Lamm, G.; Pack, G. R. J. Phys. Chem. B 1997, 101, 959.
(452) Voges, D.; Karshikoff, A. J. Chem. Phys. 1998, 108, 2219.
(453) Jayaram, B.; Liu, Y.; Beveridge, D. L. J. Chem. Phys. 1998, 109,
1465.
(454) Liu, Y.; Ichiye, T. Biophys. Chem. 1999, 78, 97.
(455) Roux, B.; Simonson, T. Biophys. Chem. 1999, 78, 1.
(456) Mallik, B.; Masunov, A.; Lazaridis, T. J. Comput. Chem. 2002,
23, 1090.
(457) Takahashi, T.; Sugiura, J.; Nagayama, K. J. Chem. Phys. 2002,
116, 8232.
(458) Vasilyev, V. J. Comput. Chem. 2002, 23, 1254.
(459) Haberthur, U.; Majeux, N.; Werner, P.; Caflisch, A. J. Comput.
Chem. 2003, 24, 1936.
(460) Kundu, S.; Gupta-Bhaya, P. J. Mol. Struct. (THEOCHEM) 2003,
639, 21.
(461) Simonson, T. Rep. Prog. Phys. 2003, 66, 737.
(462) Patel, S.; Mackerell, A. D.; Brooks, C. L. J. Comput. Chem. 2004,
25, 1504.
(463) Warshel, A. R., S. T.; Churg, A. K. Proc. Natl. Acad. Sci. U.S.A.
1984, 81, 4785.
(464) Mehler, E. L.; Eichele, G. Biochemistry 1984, 23, 3887.
(465) Mehler, E. L. J. Phys. Chem. 1996, 100, 16006.
(466) Hingerty, B. E.; Ritchie, R. H.; Ferrel, T. L.; Turner, J. E.
Biopolymers 1985, 24, 427.
(467) Ramstein, J.; Lavery, R. Proc. Natl. Acad. Sci. U.S.A. 1988, 85,
7231.
(468) Guarnieri, F.; Schmidt, A. B.; Mehler, E. L. Int. J. Quantum
Chem. 1998, 69, 57.
(469) Sandberg, L.; Edholm, O. Proteins 1999, 36, 474.
(470) Hassan, S. A.; Guarnieri, F.; Mehler, E. L. J. Phys. Chem. B
2000, 104, 6478.
(471) Bucher, M.; Porter, T. L. J. Phys. Chem. 1986, 90, 3406.
(472) Ehrenson, S. J. Comput. Chem. 1989, 10, 77.
(473) Mehler, E. L. The Lorenz-Debye-Sack theory and dielectric
screening electrostatic effects in proteins and nucleic acids. In
Molecular Electrostatic Potentials; Concepts and Applications;
Murray, J. S., Sen, K., Eds.; Elsevier: Amsterdam, The Neth-
erlands, 1996.
(474) Gelwart, W. M.; Ben-Shaul, A. J. Phys. Chem. 1996, 100, 13169.
(475) Hoshi, H.; Sakurai, M.; Inoue, Y.; Chujo, R. J. Chem. Phys. 1987,
87, 1107.
(476) Hoshi, H.; Sakurai, M.; Inoue, Y.; Chujo, R. J. Mol. Struct.
(THEOCHEM) 1988, 180, 267.
(477) Furuki, T.; Umeda, A.; Sakurai, M.; Inoue, Y.; Chujo, R.; Harata,
K. J. Comput. Chem. 1994, 15, 90.
(478) Houjou, H.; Inoue, Y.; Sakurai, M. J. Am. Chem. Soc., 1998, 120,
4459.
(479) Furuki, T.; Hosokawa, F.; Sakurai, M.; Inoue, Y.; Chujo, R. J.
Am. Chem. Soc. 1993, 115, 2903.
(480) Furuki, T.; Sakurai, M.; Inoue, Y. J. Phys. Chem. 1995, 99,
12047.
(481) Sakurai, M.; Furuki, T.; Inoue, Y. J. Phys. Chem. 1995, 99,
17789.
(482) Sakurai, M.; Tamagawa, H.; Inoue, Y.; Ariga, K.; Kunitake, T.
J. Phys. Chem. B 1997, 101, 4810.
(483) Tamagawa, H.; Sakurai, M.; Inoue, Y.; Ariga, K.; Kunitake, T.
J. Phys. Chem. B 1997, 101, 4817.
(484) Houjou, H.; Inoue, Y.; Sakurai, M. J. Phys. Chem. B 2001, 105,
867.
(485) Sakurai, M.; Sakata, K.; Saito, S.; Nakajima, S.; Inoue, Y. J.
Am. Chem. Soc. 2003, 125, 3108.
(486) Bonaccorsi, R.; Scrocco, E.; Tomasi, J. Int. J. Quantum Chem.
1986, 29, 717.
(487) Bonaccorsi, R.; Hodoscek, M.; Tomasi, J. THEOCHEM 1988, 41,
105.
(488) Bonaccorsi, R.; Ojalvo, E.; Tomasi, J. Collect. Czech. Chem.
Commun. 1988, 53, 2320.
(489) Bonaccorsi, R.; Floris, F.; Palla, P.; Tomasi, J. Thermochim. Acta
1990, 162, 213.
(490) Bonaccorsi, R.; Ojalvo, E.; Palla, P.; Tomasi, J. Chem. Phys. 1990,
143, 245.
(491) Frediani, L.; Pomelli, C. S.; Tomasi, J. Phys. Chem. Chem. Phys.
2000, 2, 4876.
(492) Frediani, L.; Mennucci, B.; Cammi, R. J. Phys. Chem. B 2004,
108, 13796.
(493) Benjamin, I. Chem. Phys. Lett. 1998, 287, 480.
(494) Benjamin, I. Chem. Rev. 1996, 96, 1449.
(495) Benjamin, I.; Kharkats, Y. I. Electrochim. Acta 1998, 44, 133.
(496) Michael, D.; Benjamin, I. J. Electroanal. Chem. 1998, 450, 335.
(497) Vieceli, J. C.; Benjamin, I. J. Chem. Phys. 2002, 117, 4532.
(498) Benjamin, I. Chem. Phys. Lett. 2004, 393, 453.
(499) Vieceli, J.; Benjamin, I. Chem. Phys. Lett. 2004, 385, 79.
(500) King, G.; Lee, F. S.; Warshel, A. J. Chem. Phys. 1991, 95, 4366.
(501) Simonson, T.; Perahia, D. Proc. Natl. Acad. Sci. U.S.A. 1995,
92, 1082.
(502) Zhou, F.; Schulten, K. J. Phys. Chem. 1995, 99, 2194
(503) Stern, H. A.; Feller, S. E. J. Chem. Phys. 2003, 118, 3401.
(504) Corni, S.; Tomasi, J. Chem. Phys. Lett. 2001, 342, 135.
(505) Corni, S.; Cappelli, C.; Cammi, R.; Tomasi, J. J. Phys. Chem. A
2001, 105, 8310.
(506) Corni, S.; Tomasi, J. J. Chem. Phys. 2002, 117, 7266.
(507) Corni, S.; Tomasi, J. J. Chem. Phys. 2002, 116, 1156.
(508) Corni, S.; Tomasi, J. J. Chem. Phys. 2003, 118, 6481.
(509) Andreussi, O.; Corni, S.; Mennucci, B.; Tomasi, J. J. Chem. Phys.
2004, 121, 10190.
(510) Ashcroft, N. W.; Mermin, N. D. Solid State Physics; Saunders
College Publishing: Fort Worth, TX, 1976.
3088 Chemical Reviews, 2005, Vol. 105, No. 8 Tomasi et al.

(511) Lundqvist, S. In Theory of the Inhomogeneous Electron Gas;
Lundqvist, S., March, N. H., Eds.; Plenum Press: New York,
1983.
(512) Lindhard, J. K. Danske Vidensk. Selsk. Mat.-Fys. Medd. 1954,
28.
(513) Kliewer, K. L.; Fuchs, R. Phys. Rev. 1969, 181, 552.
(514) Rahman, T. S.; Maradudin, A. A. Phys. Rev. B 1980, 21, 504.
(515) Jorgensen, S.; Ratner, M. A.; Mikkelsen, K. V. J. Chem. Phys.
2001, 115, 3792.
(516) Jorgensen, S.; Ratner, M. A.; Mikkelsen, K. V. J. Chem. Phys.
2002, 116, 10902.
(517) Jorgensen, S.; Ratner, M. A.; Mikkelsen, K. V. Chem. Phys. 2002,
278, 53.
(518) Bottcher, C. J. F.; Bordewijk, P. Theory of Electric Polarization,
2nd ed.; Elsevier: Amsterdam, The Netherlands, 1978.
(519) Landau, L. D.; Lifshitz, E. M. Electrodynamics of Continuous
Media; Butterworth-Heinemann: Boston, MA, 1999.
(520) Georgievskii, Y.; Hsu, C.-P.; Marcus, R. A. J. Chem. Phys. 1999,
110, 5307.
(521) Marcus, R. A. J. Chem. Phys. 1956, 24, 966.
(522) Cammi, R.; Tomasi, J. Int. J. Quantum Chem. Symp. 1995, 29,
465.
(523) Jortner, J. Mol. Phys. 1962, 5, 257.
(524) Basilevsky, M. V.; Chudinov, G. E. Chem. Phys. 1991, 157, 345.
(525) Cossi, M.; Barone, V. J. Phys. Chem. A 2000, 104, 10614.
(526) Aguilar, M. A. J. Phys. Chem. A 2001, 105, 10393.
(527) Cammi, R.; Mennucci, B. J. Chem. Phys. 1999, 110, 9877.
(528) Cossi, M.; Barone, V. J. Chem. Phys. 2001, 115, 4708.
(529) Caricato, M.; Mennucci, B.; Tomasi, J. J. Phys. Chem. A 2004,
108, 6248.
(530) Cammi, R.; Corni, S.; Mennucci, B.; Tomasi, J. J. Chem. Phys.
2005, 122, 104513.
(531) Chandler, D.; Schweizer, K. S.; Wolynes, P. G. Phys. Rev. Lett.
1982, 49.
(532) Chen, Y.-C.; Lebowitz, J. L.; Nielaba, P. J. Chem. Phys. 1989,
91, 340.
(533) Wolynes, P. G. J. Chem. Phys. 1987, 86, 5133.
(534) Nichols, I. A. L.; Calef, D. F. J. Chem. Phys. 1988, 89, 3783.
(535) Rips, I.; Klafter, J.; Jortner, J. J. Chem. Phys. 1988, 89, 4288.
(536) Bagchi, B.; Chandra, A. J. Chem. Phys. 1989, 90, 7338.
(537) Fried, L. E.; Mukamel, S. J. Chem. Phys. 1990, 93, 932.
(538) Raineri, F. O.; Zhou, Y. Q.; Friedman, H. L. Chem. Phys. 1991,
152, 201.
(539) Raineri, F. O.; Resat, H.; Perng, B. C.; Hirata, F.; Friedman, H.
L. J. Chem. Phys. 1994, 100, 1477.
(540) Biswas, R.; Bagchi, B. J. Phys. Chem. 1996, 100, 4261.
(541) Hsu, C. P.; Song, X. Y.; Marcus, R. A. J. Phys. Chem. B 1997,
101, 2546.
(542) Song, X. Y.; Chandler, D. J. Chem. Phys. 1998, 108, 2594.
(543) Song, X. Y.; Chandler, D.; Marcus, R. A. J. Phys. Chem. 1996,
100, 11954.
(544) Ingrosso, F.; Mennucci, B.; Tomasi, J. J. Mol. Liq. 2003, 108,
21.
(545) Caricato, M.; Ingrosso, F.; Mennucci, B.; Tomasi, J. J. Chem.
Phys. 2005, 122, 154501.
(546) Mennucci, B. Theor. Chem. Acc. 2005, submitted for publication.
(547) Parsons, D. F.; Vener, M. V.; Basilevsky, M. V. J. Phys. Chem.
A 1999, 103, 1171.
(548) Berg, M. Chem. Phys. Lett. 1994, 228, 317.
(549) Fourkas, J. T.; Benigno, A.; Berg, M. J. Chem. Phys. 1993, 99,
8552.
(550) Fourkas, J. T.; Berg, M. J. Chem. Phys. 1993, 98, 7773.
(551) Gardecki, J.; Horng, M. L.; Papazyan, A.; Maroncelli, M. J. Mol.
Liq. 1995, 65-6, 49.
(552) Reynolds, L.; Gardecki, J. A.; Frankland, S. J. V.; Horng, M. L.;
Maroncelli, M. J. Phys. Chem. 1996, 100, 10337.
(553) Berg, M. A.; Hubble, H. W. Chem. Phys. 1998, 233, 257.
(554) Berg, M. A. J. Chem. Phys. 1999, 110, 8577.
(555) Fleming, G. R.; Cho, M. Annu. Rev. Phys. Chem. 1996, 47, 109.
(556) Yan, Y. J.; Mukamel, S. J. Chem. Phys. 1988, 89, 5160.
(557) Mukamel, S. Principles of Nonlinear Optical Spectroscopy;
Oxford University Press: New York, 1995.
(558) Li, B. L.; Johnson, A. E.; Mukamel, S.; Myers, A. B. J. Am. Chem.
Soc. 1994, 116, 11039.
(559) Kelley, A. M. J. Phys. Chem. A 1999, 103, 6891.
(560) Marcus, R. A. J. Chem. Phys. 1965, 43, 679.
(561) Hush, N. Prog. Inorg. Chem. 1967, 8, 391.
(562) Marcus, R. A. J. Phys. Chem. 1989, 93, 3078.
(563) Marcus, R. A. Rev. Mod. Phys. 1993, 65, 599.
(564) Newton, M. D. Adv. Chem. Phys. 1999, 106, 303.
(565) Liu, Y.-P.; Newton, M. D. J. Phys. Chem. 1995, 99, 9.
(566) Bursulaya, B. D.; Zichi, D. A.; Kim, H. J. J. Phys. Chem. 1995,
99, 10069.
(567) Kim, H. J. J. Chem. Phys. 1996, 105, 6818.
(568) Schmitt, U. W.; Voth, G. A. J. Phys. Chem. B 1998, 102, 5547.
(569) Matyushov, D. V.; Voth, G. A. J. Phys. Chem. A 1999, 103, 10981.
(570) Matyushov, D. V.; Voth, G. A. J. Phys. Chem. A 2000, 104, 6470.
(571) Matyushov, D. V.; Voth, G. A. J. Phys. Chem. A 2000, 104, 6485.
(572) Naka, K.; Morita, A.; Kato, S. J. Chem. Phys. 1999, 110, 3484.
(573) Kim, H. J. J. Chem. Phys. 1996, 105, 6833.
(574) Scholes, G. D. Annu. Rev. Phys. Chem. 2003, 54, 57.
(575) Forster, T. Ann. Phys. 1948, 2, 55.
(576) Dexter, D. L. J. Chem. Phys. 1953, 21, 836.
(577) Knox, R. S.; van Amerongen, H. J. Phys. Chem. B 2002, 106,
5289.
(578) Juzeliunas, G.; Andrews, D. L. J. Chem. Phys. B 1994, 106, 5289.
(579) Agranovich, V. M.; Galanin, M. D. Electronic Excitation Energy
Transfer in Condensed Matter; North-Holland: Amsterdam, The
Netherlands, 1982.
(580) Tretiak, S.; Middleton, C.; Chernyak, V.; Mukamel, S. J. Phys.
Chem. B 2000, 104, 4519.
(581) Tretiak, S.; Middleton, C.; Chernyak, V.; Mukamel, S. J. Phys.
Chem. B 2000, 104, 9540.
(582) Hsu, C.-P.; Fleming, G. R.; Head-Gordon, M.; Head-Gordon, T.
J. Chem. Phys. 2001, 114, 3065.
(583) Iozzi, M. F.; Mennucci, B.; Tomasi, J.; Cammi, R. J. Chem. Phys.
2004, 120, 7029.
(584) Bauernschmitt, R.; Haser, M.; Treutler, O.; Ahlrichs, R. Chem.
Phys. Lett. 1997, 264, 573.
(585) Casida, M. E.; Jamorski, C.; Casida, K. C.; Salahub, D. R. J.
Chem. Phys. 1998, 108, 4439.
(586) Stratmann, R. E.; Scuseria, G. E.; Frisch, M. J. J. Chem. Phys.
1998, 109, 8218.
(587) van Gisbergen, S. J. A.; Kootstra, F.; Schipper, P. R. T.;
Gritsenko, O. V.; Snijders, J. G.; Baerends, E. J. Phys. Rev. A
1998, 57, 2556.
(588) Tozer, D. J.; Handy, N. C. Phys. Chem. Chem. Phys. 2000, 2,
2117.
(589) Marques, M. A. L.; Gross, E. K. U. Annu. Rev. Phys. Chem. 2004,
55, 427.
(590) Davydov, A. S. Theory of Molecular Excitons; McGraw Hill: New
York, 1962.
(591) Mennucci, B.; Tomasi, J.; Cammi, R. Phys. Rev. B 2004, 70,
205212.
(592) Gersten, J.; Nitzan, A. J. Chem. Phys. 1981, 75, 1139.
(593) Eckardt, W.; Penzar, Z. Phys. Rev. B 1986, 34, 8444.
(594) Leung, P. T.; Hider, M. H. J. Chem. Phys. 1993, 98, 5019.
(595) Girard, C.; Martin, O. J. F.; Dereux, A. Phys. Rev. Lett.. 1995,
75, 17.
(596) Basilevsky, M. V.; Parsons, D. F.; Vener, M. V. J. Chem. Phys.
1998, 108, 1103.
(597) Cammi, R.; Tomasi, J. Int. J. Quantum Chem. 1996, 60, 297.
(598) Nielsen, C. B.; Sauer, S. P. A.; Mikkelsen, K. V. J. Chem. Phys.
2003, 119, 3849.
(599) Cammi, R.; Cossi, M.; Mennucci, B.; Tomasi, J. J. Chem. Phys.
1996, 105, 10556.
(600) Mennucci, B.; Tomasi, J.; Cammi, R.; Cheeseman, J. R.; Frisch,
M. J.; Devlin, F. J.; Gabriel, S.; Stephens, P. J. J. Phys. Chem.
A 2002, 106, 6102.
(601) Kramer, P.; Saraceno, M. Geometry of the Time-Dependent
Variational Principle in Quantum Mechanics; Springer-Verlag:
Berlin, Germany, 1981.
(602) Moccia, R. Int. J. Quantum Chem. 1973, 7, 779.
(603) Moccia, R. Int. J. Quantum Chem. 1974, 8, 293.
(604) Arickx, F.; Broeckhove, J.; Kesteloot, E.; Lathouwers, L.; Van
Leuven, P. Chem. Phys. Lett. 1986, 128, 310.
(605) Broeckhove, J.; Lathouwers, L.; Kesteloot, E.; Van Leuven, P.
Chem. Phys. Lett. 1988, 149, 547.
(606) Broeckhove, J.; Lathouwers, L.; van Leuven, P. J. Mol. Struct.
(THEOCHEM) 1989, 199, 245.
(607) McLachlan, A. D.; Ball, M. A. Rev. Mod. Phys. 1964, 36, 844.
(608) Dalgarno, A. In Perturbation Theory and its Application in
Quantum Mechanics; Wilcox, C. H., Ed.; Wiley: New York, 1967.
(609) Langhoff, P. W.; Epstein, S. T.; Karplus, M. Rev. Mod. Phys.
1972, 44, 602.
(610) Lowdin, P. O.; Mukerjee, P. K. Chem. Phys. Lett. 1972, 14, 1.
(611) Olsen, J.; Jorgensen, P., Time-dependent response theory with
applications to SCF and MCSCF wavefunctions. In Modern
Electronic Structure Theory; Yarkony, D. R., Ed.; World Scien-
tific: Singapore, 1995; Vol. 2.
(612) Lavenda, B. H. Nuovo Cimento Soc. Ital. Fis. B 2003, 118, 143.
(613) Christiansen, O.; Jorgensen, P.; Hattig, C. Int. J. Quantum
Chem. 1998, 68, 1.
(614) Cammi, R.; Mennucci, B.; Tomasi, J. J. Phys. Chem. A 2000,
104, 5631.
(615) Casida, M. E. In Recent Advances in Computational Chemistry;
Chong, D. P., Ed.; World Scientific: Singapore, 1995; Vol. 1.
(616) Bauernschmitt, R.; Ahlrichs, R. Chem. Phys. Lett. 1996, 256, 454.
(617) Hirata, S.; Head-Gordon, M. Chem. Phys. Lett. 1999, 314, 291.
(618) Fonseca, T.; Kim, H. J.; Hynes, J. T. J. Photochem. Photobiol.,
A: Chem. 1994, 82, 67.
(619) Gedeck, P.; Schneider, S. J. Photochem. Photobiol., A: Chem.
1997, 105, 165.
(620) Kim, H. J.; Hynes, J. T. J. Photochem. Photobiol., A: Chem. 1997,
105, 337.
(621) Gedeck, P.; Schneider, S. J. Photochem. Photobiol., A: Chem.
1999, 121, 7.
Quantum Mechanical Continuum Solvation Models Chemical Reviews, 2005, Vol. 105, No. 8 3089

(725) Rega, N.; Cossi, M.; Barone, V. J. Am. Chem. Soc. 1998, 120,
5723.
(726) Jolibois, F.; Cadet, J.; Grand, A.; Subra, R.; Rega, N.; Barone,
V. J. Am. Chem. Soc. 1998, 120, 1864.
(727) Adamo, C.; Heitzmann, M.; Meilleur, F.; Rega, N.; Scalmani, G.;
Grand, A.; Cadet, J.; Barone, V. J. Am. Chem. Soc. 2001, 123,
7113.
(728) Saracino, G. A. A.; Tedeschi, A.; D’Errico, G.; Improta, R.; Franco,
L.; Ruzzi, M.; Corvaia, C.; Barone, V. J. Phys. Chem. A 2002,
106, 10700.
(729) Cossi, M.; Crescenzi, O. J. Chem. Phys. 2003, 118, 8863.
(730) D’Amore, M.; Improta, R.; Barone, V. J. Phys. Chem. A 2003,
107, 6264.
(731) Pavone, M.; Benzi, C.; De Angelis, F.; Barone, V. Chem. Phys.
Lett. 2004, 395, 120.
(732) Cossi, M.; Crescenzi, O. Theor. Chem. Acc. 2004, 111, 162.
(733) Crescenzi, O.; Pavone, M.; De Angelis, F.; Barone, V. J. Phys.
Chem. B 2005, 109, 445.
(734) Benzi, C.; Cossi, M.; Barone, V.; Tarroni, R.; Zannoni, C. J. Phys.
Chem. B 2005, 109, 2584.
(735) Manalo, M. N.; de Dios, A. C.; Cammi, R. J. Phys. Chem. A 2000,
104, 9600.
(736) Manalo, M. N.; de Dios, A. C. Magn. Reson. Chem. 2002, 40,
781.
(737) Wu, A. A.; Cremer, D.; Gauss, J. J. Phys. Chem. A 2003, 107,
8737.
(738) Tuttle, T.; Kraka, E.; Wu, A.; Cremer, D. J. Am. Chem. Soc. 2004,
126, 5093.
(739) Zaccari, D. G.; Snyder, J. P.; Peralta, J. E.; Taurian, O. E.;
Contreras, R. H.; Barone, V. Mol. Phys. 2002, 100, 705.
(740) Zaccari, D. G.; Barone, V.; Peralta, J. E.; Contreras, R. H.;
Taurian, O. E.; Diez, E.; Esteban, A. Int. J. Mol. Sci. 2003, 4,
93.
(741) Ciofini, I.; Reviakine, R.; Arbuznikov, A.; Kaupp, M. Theor.
Chem. Acc. 2004, 111, 132.
(742) Rinkevicius, Z.; Telyatnyk, L.; Vahtras, O.; Ruud, K. J. Chem.
Phys. 2004, 121, 5051.
(743) Buckingham, A. D.; Schafer, T.; Schneider, W. G. J. Chem. Phys.
1960, 32, 1227.
(744) Kamlet, M. J.; Abboud, J. L. M.; Taft, R. W. Prog. Phys. Org.
Chem. 1980, 13, 485.
(745) Witanowski, M.; Biedrzycka, Z.; Sicinska, W.; Grabowski, Z. J.
Magn. Reson. 1998, 131, 54.
(746) Witanowski, M.; Sicinska, W.; Biedrzycka, Z.; Webb, G. A. J.
Mol. Struct. 1999, 476, 133.
(747) Witanowski, M.; Biedrzycka, Z.; Sicinska, W.; Webb, G. A. J.
Mol. Struct. 2000, 516, 107.
(748) Witanowski, M.; Biedrzycka, Z.; Sicinska, W.; Grabowski, Z. J.
Mol. Struct. 2002, 602, 199.
(749) Witanowski, M.; Biedrzycka, Z.; Sicinska, W.; Grabowski, Z. J.
Magn. Reson. 2003, 164, 212.
(750) Ramsey, N. F. Phys. Rev. 1950, 78, 699.
(751) Ramsey, N. F. Phys. Rev. 1953, 91, 303.
(752) Ditchfield, R. Mol. Phys. 1974, 27, 789.
(753) Wolinski, K.; Hinton, J. F.; Pulay, P. J. Am. Chem. Soc. 1990,
112, 8251.
(754) Gauss, J. J. Chem. Phys. 1993, 99, 3629.
(755) Gauss, J.; Stanton, J. F. J. Chem. Phys. 1995, 102, 251.
(756) Gauss, J.; Stanton, J. F. J. Chem. Phys. 1996, 104, 2574.
(757) Cheeseman, J. R.; Trucks, G. W.; Keith, T. A.; Frisch, M. J. J.
Chem. Phys. 1996, 104, 5497.
(758) Helgaker, T.; Wilson, P. J.; Amos, R. D.; Handy, N. C. J. Chem.
Phys. 2000, 113, 2983.
(759) Ho¨ller, R. L. L. Chem. Phys. Lett. 1981, 84, 94.
(760) Hinton, J. F.; Guthrie, P.; Pulay, P.; Wolinski, K. J. Am. Chem.
Soc. 1992, 114, 1604.
(761) Chesnut, D. B.; Rusidoski, B. E. J. Mol. Struct. (THEOCHEM)
1994, 314, 19.
(762) Pecul, M. J. S. Chem. Phys. 1998, 234, 111.
(763) Benedict, H.; Shenderovich, I. G.; Malkina, O. L.; Malkin, V. G.;
Denisov, G. S.; Golubev, N. S.; Limbach, H. H. J. Am. Chem.
Soc. 2000, 122, 1979.
(764) Pecul, M.; Leszczynski, J.; Sadlej, J. J. Phys. Chem. A 2000, 104,
8105.
(765) Pecul, M.; Leszczynski, J.; Sadlej, J. J. Chem. Phys. 2000, 112,
7930.
(766) Pecul, M.; Sadlej, J.; Leszczynski, J. J. Chem. Phys. 2001, 115,
5498.
(767) Del Bene, J. E.; Jordan, M. J. T. J. Phys. Chem. A 2002, 106,
5385.
(768) Del Bene, J. E.; Elguero, J.; Alkorta, I.; Yanez, M.; Mo, O. J.
Chem. Phys. 2004, 120, 3237.
(769) Helgaker, T.; Watson, M.; Handy, N. C. J. Chem. Phys. 2000,
113, 9402.
(770) Engels, B.; Eriksson, L. A.; Lunell, S. Adv. Quantum Chem. 1996,
27, 297.
(771) Improta, R.; Barone, V. Chem. Rev. 2004, 104, 1231.
(772) Schreckenbach, G.; Ziegler, T. J. Phys. Chem. A 1997, 101, 3388.
(773) vanLenthe, E.; Wormer, P. E. S.; vanderAvoird, A. J. Chem.
Phys. 1997, 107, 2488.
(774) Kaupp, M.; Remenyi, C.; Vaara, J.; Malkina, O. L.; Malkin, V.
G. J. Am. Chem. Soc. 2002, 124, 2709.
(775) Kaupp, M.; Reviakine, R.; Malkina, O. L.; Arbuznikov, A.;
Schimmelpfennig, B.; Malkin, V. G. J. Comput. Chem. 2002, 23,
794.
(776) Rosenfeld, L. Z. Phys. 1928, 52, 161.
(777) Hansen, A. E.; Bouman, T. D. Adv. Chem. Phys. 1980, 44, 545.
(778) Ruiz Lopez, M. F.; Rinaldi, D. J. Mol. Struct. (THEOCHEM)
1983, 93, 277.
(779) Furche, F.; Ahlrichs, R.; Wachsmann, C.; Weber, E.; Sobanski,
A.; Vogtle, F.; Grimme, S. J. Am. Chem. Soc. 2000, 122, 1717.
(780) Autschbach, J.; Ziegler, T.; Gisbergen, S. J. A. v.; Baerends, E.
J. J. Chem. Phys. 2002, 116, 6930.
(781) Diedrich, C.; Grimme, S. J. Phys. Chem. A 2003, 107, 2524.
(782) Kongsted, J.; Hansen, A. E.; Pedersen, T. B.; Osted, A.; Mikkels-
en, K. V.; Christiansen, O. Chem. Phys. Lett. 2004, 391, 259.
(783) Kongsted, J.; Pedersen, T. B.; Osted, A.; Hansen, A. E.; Mikkels-
en, K. V.; Christiansen, O. J. Phys. Chem. A 2004, 108, 3632.
(784) Pecul, M.; Marchesan, D.; Ruud, K.; Coriani, S. J. Chem. Phys.
2005, 122, 024106.
(785) Moffitt, W.; Moscowitz, A. J. Chem. Phys. 1959, 30, 648.
(786) Condon, E. U. Rev. Mod. Phys. 1937, 9, 432.
(787) Eyring, H.; Walter, J.; Kimball, G. E. Quantum Chemistry;
Wiley: New York, 1944.
(788) Amos, R. D. Chem. Phys. Lett. 1982, 87, 23.
(789) Bak, K. L.; Jorgensen, P.; Helgaker, T.; Ruud, K. Faraday
Discuss. 1994, 121.
(790) Cheeseman, J. R.; Frisch, M. J.; Devlin, F. J.; Stephens, P. J. J.
Phys. Chem. A 2000, 104, 1039.
(791) Stephens, P. J.; Devlin, F. J.; Cheeseman, J. R.; Frisch, M. J. J.
Phys. Chem. A 2001, 105, 5356.
(792) Polavarapu, P. L. Chirality 2002, 14, 768-781.
(793) Stephens, P. J.; Devlin, F. J.; Cheeseman, J. R.; Frisch, M. J.;
Mennucci, B.; Tomasi, J. Tetrahedron-Asymmetry 2000, 11,
2443.
(794) da Silva, C. O.; Mennucci, B.; Vreven, T. J. Org. Chem. 2004,
69, 8161.
(795) Marchesan, D.; Coriani, S.; Forzato, C.; Nitti, P.; Pitacco, G.;
Ruud, K. J. Phys. Chem. A 2005, 109, 1449.
(796) Kongsted, J.; Pedersen, T. B.; Strange, M.; Osted, A.; Hansen,
A. E.; Mikkelsen, K. V.; Pawlowski, F.; Jorgensen, P.; Hattig,
C. Chem. Phys. Lett. 2005, 401, 385.
(797) Holzwarth, G.; Hsu, E. C.; Mosher, H. S.; Faulkner, T. R.;
Moscowitz, A. J. Am. Chem. Soc. 1974, 96, 251.
(798) Nafie, L. A.; Cheng, J. C.; Stephens, P. J. J. Am. Chem. Soc.
1975, 97, 3842.
(799) Barron, L. D.; Bogaard, M. P.; Buckingham, A. D. J. Am. Chem.
Soc. 1973, 95, 603.
(800) Hug, W.; Kint, S.; Bailey, G. F.; Scherer, J. R. J. Am. Chem.
Soc. 1975, 97, 5589.
(801) Stephens, P. J.; Ashvar, C. S.; Devlin, F. J.; Cheeseman, J. R.;
Frisch, M. J. Mol. Phys. 1996, 89, 579.
(802) Devlin, F. J.; Stephens, P. J.; Cheeseman, J. R.; Frisch, M. J. J.
Phys. Chem. A 1997, 101, 9912.
(803) Stephens, P. J.; Chabalowski, C. F.; Devlin, F. J.; Jalkanen, K.
J. Chem. Phys. Lett. 1994, 225, 247.
(804) Stephens, P. J.; Devlin, F. J. Chirality 2000, 12, 172.
(805) Cappelli, C.; Corni, S.; Mennucci, B.; Cammi, R.; Tomasi, J. J.
Phys. Chem. A 2002, 106, 12331.
(806) Frimand, K.; Bohr, H.; Jalkanen, K. J.; Suhai, S. Chem. Phys.
2000, 255, 165.
(807) Abdali, S.; Jalkanen, K. J.; Bohr, H.; Suhai, S.; Nieminen, R.
M. Chem. Phys. 2002, 282, 219.
(808) Stephens, P. J.; Devlin, F. J. Chirality 2000, 12, 172.
(809) Schaefer, M.; Bartels, C.; Karplus, M. J. Mol. Biol. 1998, 284,
835.
(810) Hassan, S. A.; Guarnieri, F.; Mehler, E. L. J. Phys. Chem. B
2000, 104, 6490.
(811) Hassan, S. A.; Mehler, E. L.; Zhang, D. Q.; Weinstein, H. Proteins
2003, 51, 109.
(812) Hockney, R. W.; Eastwood, J. W. Computer Simulation Using
Particles, 2nd ed.; IOP: Bristol, U.K., 1988.
(813) Essmann, U.; Perera, L.; Berkowitz, M. L.; Darden, T.; Lee, H.;
Pedersen, L. G. J. Chem. Phys. 1995, 103, 8577.
(814) Allen, M. P.; Tildesley, D. J. Computer Simulation of Liquids;
Oxford University Press: New York, 1989.
(815) van Gunsteren, W. F.; Berendsen, H. J. C. Angew. Chem.-Int.
Ed. 1990, 29, 992.
(816) Steinbach, P. J.; Brooks, B. R. J. Comput. Chem. 1994, 15, 667.
(817) Hummer, G.; Soumpasis, D. M.; Neumann, M. J. Phys.-Condens.
Matter 1994, 6, A141.
(818) Hu¨nenberger, P. H.; van Gunsteren, W. F. J. Chem. Phys. 1998,
108, 6117.
(819) Alper, H.; Levy, R. M. J. Chem. Phys. 1993, 99, 9847.
(820) Wang, L.; Hermans, J. J. Phys. Chem. 1995, 99, 12001.
(821) Barker, J. A. Mol. Phys. 1994, 83, 1057.
Quantum Mechanical Continuum Solvation Models Chemical Reviews, 2005, Vol. 105, No. 8 3091