spectroscopy, and ’21.0 kJ/mol from accurate quantum calculations. This

is also apparent from the second virial coe¬cient, which deviates from ex-

perimental values in the sense that the dimer attraction is too

large (Fig. 6.9).20 As noted by Berendsen et al. (1987), there is an incon-

sistency in the use of e¬ective pair potentials when these incorporate average

induced dipoles: the full electrostatic interaction between two molecules is

taken into account, while the polarization energy of induced dipoles equals

only half the electrostatic energy. Half of the electrostatic energy is “used”

for the formation of the induced dipoles (see below for the correct equa-

tions). While polarizable models take this factor of two correctly into ac-

count, e¬ective pair potentials do not. The heat of vaporization should be

corrected with the self-energy of the induced dipoles. It was found that a

water model with such corrections (the extended simple point charge model,

SPC/E) gives better values for density, radial distribution function, di¬usion

constant and dielectric constant than the same model (SPC) without correc-

tions. On the other hand, the heat of vaporization, and also the free energy

change from liquid to gas, is too large, as the molecule retains its enhanced

dipole moment also in the gas phase. As seen in Fig. 6.9, the discrepancy

from gas phase dimer interaction is even larger than for the “classical” ef-

fective pair models. The boiling point should therefore be higher than the

experimental value and free energies of solvation of water into non-polar

environments should also be too large.21

There is no remedy to these artifacts other than replacing e¬ective poten-

tials with polarizable ones.

20 It is through the pioneering work of A. Rahman and F. H. Stillinger in the early 1970s (Rah-

man and Stillinger, 1971; Stillinger and Rahman, 1972, 1974) that the importance of e¬ective

potentials became clear. Their ¬rst simulation of liquid water used the Ben-Naim-Stillinger

(BNS) model that had been derived on the basis of both gas-phase data (second virial coe¬-

cient related to the pure pair potentials) and condensed-phase data (ice). This pair potential

appeared too weak for the liquid phase and could be improved by a simple scaling of energy.

When a modi¬ed version, the ST2 potential (Stillinger and Rahman, 1974), was devised, the

notion of an e¬ective potential was already developed (Berendsen et al., 1987).

21 The SPC dipole moment of 2.274 D is enhanced to 2.351 D in SPC/E (compare the gas phase

value of 1.85 D). The heat of vaporization at 300 K (41.4 kJ/mol) increases to 46.6 kJ/mol,

which decreases the vapor pressure tenfold and increases the boiling point by several tens of

degrees. Amazingly, according to Guissani and Guillot (1993), the SPC/E model follows the

liquid-vapor coexistence curve quite accurately with critical parameters (Tc = 640 K, ρc =

0.29 g/cm’3 , Pc = 160 bar) close to those of real water (Tc = 647.13 K, ρc = 0.322 g/cm’3 ,

Pc = 220.55 bar). It does better than SPC with Tc = 587 K and ρc = 0.27 g/cm’3 (de Pablo

et al., 1990).

176 Molecular dynamics

6.3.7 Choices for polarizability

The ¬rst choice to be made is the representation of the induced dipoles in

the model. We assume that the ¬xed sources in the non-polarized model are

charges qi only.22 One may:

• include induced dipoles (dipolar model ),

• represent the induced dipoles by a positive and negative charge, connected

with a harmonic spring (shell model ),23 or

• modify charges at given positions to include induced dipoles (¬‚uctuating

charge model ).24

The dipole, or in the second case one of the charges, can be placed on atoms

or on virtual points that are de¬ned in terms of atom positions. For the shell

model, the other charge is in principle a massless interaction site. In the

third case there must be a su¬cient number of properly placed charges (four

for every polarizability). It is highly recommended to use isotropic polariz-

abilities on each site; in that case the induced dipole points in the direction

of the electric ¬eld and “ if iterated to self-consistency “ there is no torque

μ — E acting on the dipole or the spring. With anisotropic polarizabilities

there are torques, also in the shell model where the spring constant would

be an anisotropic tensor, which considerably complicates the equations of

motion. In the isotropic dipolar model there is a force only, given by the

tensor product of dipole moment and ¬eld gradient: one needs to evaluate

the ¬eld gradient on each dipolar site. In the isotropic shell model there are

only charges and no ¬eld gradients are needed; however, the number of par-

ticles doubles per polarizable site and the number of Coulomb interactions

quadruples. The ¬‚uctuating charge model can accommodate anisotropic

polarizabilities, but requires an even larger number of interactions.

The question arises if the choice of isotropic polarizabilities is adequate.

After all, many molecules possess an anisotropic polarizability, and this can

never be generated by simple addition of isotropic atomic polarizabilities.

However, the mutual interaction of induced dipoles saves the day: neighbor-

ing induced dipoles enhance an external electric ¬eld when they are lined up

22 The inclusion of ¬xed dipoles is a trivial complication, but one should be careful to obtain forces

from the correct derivatives of the dipolar ¬elds. Forces on dipoles include torques, which are

to be transferred to the rigid molecular frame in which the dipoles are de¬ned.

23 The shell model idea originates from solid state physics, where an ion is represented by a

positive core and a negative shell which carries the interactions of electronic nature. Such

models correctly predict phonon spectra. It was ¬rst used by Dick and Overhauser (1958) on

alkali halide crystals.

24 Although used by Berendsen and van der Velde (1972), the ¬‚uctuating charge model was ¬rst

introduced into the literature by Zhu et al. (1991) and by Rick et al. (1994). See also Stern et

al. (2001).

6.3 Force ¬eld descriptions 177

in the ¬eld direction and depress the ¬eld when they are situated in lateral

directions, thus producing an anisotropic overall polarization. In general

it is very well possible to reproduce experimental molecular polarizability

tensors based on interacting isotropic atomic polarizabilities.25 Consider as

example a homonuclear diatomic molecule, with two points with isotropic

polarizability ± at the nuclear positions, separated by a distance d, in an

external ¬eld E. The ¬eld at each point includes the dipolar ¬eld from the

other point. For the induced dipole per point in the parallel and perpendic-

ular directions to the molecular axis we obtain, denoting ± = ±/(4πµ0 ):

2μ ±E

’ μ=

μ = ± E+ , (6.62)

1 ’ 2± /d3

4πµ0 d3

μ⊥ ±E

= ± E’ ’ μ⊥ =

μ⊥ , (6.63)

4πµ0 d3 1 + ± /d3

resulting in a molecular anisotropic polarizability:

2±

±= (6.64)

1 ’ 2± /d3

2±

±⊥ = (6.65)

1 + ± /d3

±iso = (± + 2±⊥ )/3. (6.66)

The parallel component diverges for a distance of (2± )1/3 . This is a com-

pletely unphysical behavior, leading to far too large anisotropies for di-

atomics (see Fig. 6.10). The remedy (Thole, 1981) is to introduce a damping

at short distances by considering the induced dipole not as a point, but as a

distribution. Both an exponential decay and a linear decay yield satisfactory

results, with a single polarizability for each atom irrespective of its chemical

environment, and with a single screening length, if distances between two

atoms are scaled by the inverse sixth power of the product of the two atom

polarizabilities (van Duijnen and Swart, 1998).

Another empirical approach to predict isotropic polarizabilities of organic

molecules, based on additivity, uses the number of electrons, their e¬ective

quantum numbers and e¬ective nuclear shielding (Glen, 1994).

25 Applequist et al. (1972) included interactions between point dipoles and lists atomic polariz-

abilities both for a poor additive model and a more successful interaction model. To amend

the unrealistic approach to diverging behavior at short distances, Thole (1981) proposed to

insert a damping at small distances, as can be realized by considering the polarizable sites

not as points, but as smeared-out distributions. A simple triangular distribution with a width

scaled by the polarizabilities appeared to yield good results with a single isotropic polarizabil-

ity for each atom type (H,C,N,O), irrespective of its chemical nature. Van Duijnen and Swart

(1998) extended Thole™s model for a wider set of polarizabilities, including ab initio data for

70 molecules and including sulphur and the halogens in addition to Thole™s atom set. They

compared a linear as well as an exponential distribution, with a slight preference for the latter.

178 Molecular dynamics

polarizability relative to iso

2.5

1: H2

2: N2

2 3: O2

±// / ±iso 4: Cl2

5: F2