1.5
1 2 3 3 4 5
1
±_ „ ±iso

0.5
1.4 1.5 1.6 1.7 1.8 1.9
d / ±1/3
Figure 6.10 Inadequacy of the interacting dipole model for polarizability of di
atomics: The parallel and perpendicular components, relative to the isotropic value
of the polarizability are drawn as a function of the internuclear distance, the latter
expressed relative to the 1/3 power of the atomic polarizability that reproduces the
isotropic molecular polarizability. The symbols indicate the values of the parallel
(uptriangles), isotropic (circles) and perpendicular (downtriangles) polarizabili
ties for diatomic molecules. Filled symbols represent theoretical values from recent
NIST tables; open symbols are from experimental values as quoted by Thole (1981).
It should be noted that shell models with interacting shell distributions
will be equally capable of producing correct polarizabilities. Let us discuss
the question of which description should be preferred in polarizable force
¬elds: induced dipoles, harmonic shells or ¬‚uctuating charges? Compu
tational considerations were given above, and combining a preference for
both simplicity and speed, they favor the shell model. More important are
scienti¬c considerations that favor the shell model as well. We follow the
argument of Jordan et al. (1995). Consider two neon atoms: if described
by a repulsion/dispersion interaction and dipolar polarizabilities centered
on the nuclei, there are no electric e¬ects when the atoms approach each
other. However, we know that collisioninduced infrared absorption can be
observed, and the (computed) electrical quadrupole moment of a neon pair
is not zero. The logical explanation is that when the electron clouds attract
6.3 Force ¬eld descriptions 179
each other at larger distances by dispersion, they are pulled in the direction
of the other atom and produce a positive quadrupole moment; when they
repel each other at shorter distances, they are pushed away and produce
a negative quadrupole moment. The quadrupole moment, computed by a
highlevel ab initio method, does indeed follow the interatomic force. A shell
model with the interatomic interaction acting between the shells reproduces
such results. Pursuing this idea, Jordan et al. (1995) succeeded in devis
ing a shell model for the nitrogen molecule which reproduces experimental
static and dynamic properties in the solid, liquid and gas phases, including
a subtle pressureinduced phase transition between a cubic and a hexagonal
solid phase. The model contains three shells, one bound to both nitrogens
(a “bond”) and the other two bound to one nitrogen (“lone pairs”).
Models of this type could be devised for other molecules, but they tend to
become complex and computationally intensive. The number of interaction
sites would be quite large: one shell per bond and at least one additional
shell for each nonhydrogen atom; thus far, complete force ¬elds along these
lines have not been constructed. A somewhat similar and quite successful
shell model for water (the MCDHO model) was published by SaintMartin
et al. (2000).26 The model includes a single shell with an exponential charge
distribution, which is connected to the oxygen nucleus by a spring but also
interacts with the hydrogen charges in the molecule. Parameters were ¬tted
to ab initio calculations of dimer and oligomers. The fact that this model
reproduces experimental properties of the gas, liquid and solid phases quite
well holds the promise that general force ¬elds with transferable terms may
be derived from ab initio calculations on small molecules, possibly with small
empirical corrections, if electronic polarization is properly included.
More simple polarizable models have been quite successful as well. Most
work has been done on models for water, with the aim to construct models
that yield accurate thermodynamic and dynamic properties for a wide range
of phases and conditions. While the development of polarizable models is
still proceeding, we shall review only the features from which basic principles
can be learned. A comprehensive review of models for simulation of water
by Guillot (2002) is available.
26 MCDHO: Mobile Charge Densities in Harmonic Oscillators. The article includes Monte Carlo
simulations of the liquid state. See also Hess et al. (2002), who performed molecular dynamics
on this model.
180 Molecular dynamics
The successful pairadditive simple point charge models27 have been modi
¬ed to include polarization, both with induced dipoles and with shell models.
The earliest attempt of this type was a modi¬cation of the SPC model by
Ahlstr¨m et al. (1989), who added a polarizable point dipole on the oxygen
o
atom while reducing the charges to obtain the gas phase dipole moment, and
others of this type followed with ¬netuning of parameters.28 Van Maaren
and van der Spoel (2001) investigated the properties of a shell model, with
the moving charge attached with a spring to a virtual position on the sym
metry axis of the molecule about 0.014 nm from the oxygen. They retained
the gas phase structure, dipole and quadrupole moment and polarizability,
while optimizing the Lennard“Jones interaction. Several water models have
a similar position for the negative charge for the simple reason that three
charge models cannot satisfy the experimental quadrupole moment without
displacing the negative charge. Yu et al. (2003) developed a simple model
with three atoms plus moving charge bound with a spring to the oxygen,
which they named the “chargeonspring” model. The model was intended
for computational e¬ciency; it uses a large moving charge and is in fact a
shellimplementation of an induced point dipole model. Polarizable mod
els of this kind, whether they use a moving charge or a point dipole, are
moderately successful. They repair the main de¬ciencies of e¬ective pair
potentials, but do not show the accurate allphase behavior that one should
wish. For example, the interaction in the critical region “ and thereby the
27 SPC (Berendsen et al., 1981) uses the three atoms as interaction site with partial charges
on oxygen (’0.82 e) and hydrogens (+0.41 e). The geometry is rigid, rOH = 0.1 nm; HOH
angle = 109.47—¦ . There is a Lennard“Jones interaction on the oxygens: C6 = 2.6169 — 10’3
kJ mol’1 nm6 ; C12 = 2.6332—10’6 kJ mol’1 nm12 . The similar TIP3 model (Jorgensen, 1981)
has the rigid experimental geometry: OH distance of 0.09572 nm, HOH angle of 104.52—¦ , with
hydrogen charge of 0.40 e; C6 = 2.1966 — 10’3 ; C12 = 2.4267 — 10’6 (units as above). This
model did not show a secondneighbor peak in the radial distribution function and was modi¬ed
to TIP3P (Jorgensen et al., 1983), with hydrogen charge of 0.417 e; C6 = 2.4895 — 10’3 and
C12 = 2.4351 — 10’6 . In a subsequent foursite model (TIPS2: Jorgensen, 1982) the negative
charge was displaced to a point M on the bisectrix in the direction of the hydrogens at 0.015
nm from the oxygen position, while the Lennard“Jones interaction remained on the oxygen,
and the parameters were improved in the TIP4P model (Jorgensen et al., 1983). The hydrogen
charge is 0.52 e; C6 = 2.5522 — 10’3 and C12 = 2.5104 — 10’6 . Finally, the SPC/E model
(Berendsen et al., 1987), which includes an energy correction for average polarization, is like
the SPC model but with oxygen charge 0.4238. Van der Spoel et al. (1998) evaluated these
models and optimized them for use with a reaction ¬eld. They conclude that the SPC/E model
is superior in predicting properties of liquid water.
28 These include models by Cieplack et al. (1990) with polarizabilities on O and H and extra
repulsion/dispersion terms; Caldwell et al. (1990) with SPC/E modi¬ed by O and H polariz
abilities and a threebody repulsion term; Dang (1992) with an improvement on the the latter
model; Wallqvist and Berne (1993) with a polarizable and a nonpolarizable model with extra
terms; Chialvo and Cummings (1996) with an evaluation of displacement of negative charge
and a point polarizability on oxygen; Svishchev et al. (1996) with the PPC model which has a
displaced (by 0.011 nm) negative charge and polarizability only in the molecular plane caused
by displacement of the negative charge (this is a partial shell model with enhanced permanent
dipole and reduced polarizability), Dang and Chang (1997) with a revised foursite model with
the negative charge and dipolar polarizability displaced by 0.0215 nm from the oxygen;
6.3 Force ¬eld descriptions 181
critical temperature “ is underestimated,29 and the dielectric constant is
often too large. Both e¬ects are a result of the relatively (too) large polar
ization interaction at short distances, which is then compensated by (too
weak) attraction at longer distances. A shell model developed by Lam
oureux et al. (2003) could be made to ¬t many water properties including
the dielectric constant, but only with reduced polarizability (1.04 instead of
1.44 ˚3 ). A study of ion hydration (Sp˚
A angberg and Hermansson, 2004) with
various water models showed too large solvation enthalpies for polarizable
models. Giese and York (2004) ¬nd an overpolarizability for chains of water
molecules using screened Coulomb interactions; they suggest that exchange
overlap reduces shortrange polarization.
The answer seems to lie in a proper damping of the shortrange Coulombic
¬elds and polarization. This can be accomplished by Tholetype smearing of