can be found by a least squares optimization procedure. This method

su¬ers from some arbitrariness due to the choice of grid points and

their weights. Another, more robust, method is to ¬t the charges to

multipoles derived from accurate quantum calculations.5

With the use of pair-additive Coulomb interactions, the omission

of explicit polarization, and/or the incomplete treatment of long-

range interactions, the empirically optimized partial charges do not

and should not correspond to the ab initio-derived charges. The

Coulomb interactions should then include the average polarization

and the average e¬ects of the omission of polarizing particles in the

environment. Modi¬cation of partial charges cannot achieve the cor-

rect results, however, as it will completely miss the dielectric solva-

5 See Jensen (2006) for a general, and Sigfridsson and Ryde (1998) for a more speci¬c discussion.

The latter authors advocate the multipole-¬tting method, as do Swart et al. (2001), who

use density-functional theory to derive the multipole moments, and list partial charges for 31

molecules and all aminoacids.

6.3 Force ¬eld descriptions 157

tion energy of (partial) charges in electronically polarizable environ-

ments. If average polarization enhances dipole moments, as in water,

the partial charges are enhanced, while reducing the charges may be

appropriate to mimic the interactions in a polarizable environment.

These de¬ciencies are further discussed in Section 6.3.6.

6.3.3 More sophisticated force ¬elds

Force ¬elds that go beyond the simple type described above may include the

following extra or replacing features:

(i) Polarizability This is the single most important improvement, which

is further detailed in Section 6.3.6.

(ii) Virtual interaction sites Several force ¬elds use interaction sites that

do not coincide with atomic positions. For example, one may place a

partial charge at the position of a “bond” midway between two atoms

1 and 2: r = 1 (r 1 + r 2 ). Such sites are always a (vector) function of

2

n atomic positions:

r = r(r 1 , r 2 , . . . r n ), (6.31)

and move with these positions. Virtual sites have no mass and do not

participate directly in the equations of motion; they are reconstructed

after every dynamic step. The potential energy V (r, . . .) depends on

r 1 . . . r n via its dependence on r and the force F acting on a virtual

site is distributed among the atoms on which the site depends (we

write the ±-component):

‚V ‚r ‚r

F1± = ’ · =F · . (6.32)

‚r ‚x1± ‚x1±

For the simple case of the halfway-site F 1 = F 2 = 1 F . Other linear

2

combinations are similarly simple; more complex virtual sites as out-

of-plane constructions are more complicated but follow from the same

equation.

(iii) Dummy particles These are sites that carry mass and participate in

the equations of motion. They replace real atoms and are meant to

simplify rigid-body motions. The replaced atoms are reconstructed

as virtual sites. For example, the 12 atoms with 12 constraints of a

rigid benzene molecule C6 H6 can be replaced by three dummy atoms

with three constraints, having the same total mass and moments of

inertia. All interaction sites (atoms in this case) can be reconstructed

158 Molecular dynamics

by linear combinations from the dummies.6 Feenstra et al. (1999)

have used dummy sites to eliminate fast motions of hydrogen atoms

in proteins, enabling an increase of time step from the usual 2 fs to as

much as 7 fs. Dummy atoms do not really belong in this list because

they are not a part of the force-¬eld description.

(iv) Coupling terms Force ¬elds that aim at accurate reproduction of

vibrational properties include coupling terms between bond, bond-

angle and dihedral displacements.

(v) Flexible constraints Internal vibrations with frequencies higher than

kB T /h exhibit essential quantum behavior. At very high frequencies

the corresponding degrees of freedom are in the ground state and can

be considered static. It seems logical, therefore, to treat such degrees

of freedom as constraints. However, it is quite tricky to separate the

real quantum degree of freedom from the classical degrees of freedom,

as the ¬‚uctuating force on anharmonic bonds shifts the harmonic os-

cillator both in position and in energy. Constraining the quantum

vibration amounts to constraining the bond length to the position

where the net force vanishes. Such “¬‚exible constraints” were ¬rst

proposed by Zhou et al. (2000) and have been implemented in a po-

larizable ab initio water model by Hess et al. (2002) in molecular

dynamics and by Saint-Martin et al. (2005) in Monte Carlo algo-

rithms. The di¬erence with the usual holonomic constraints is not

very large.

(vi) Charge distributions Descriptions of the electronic charge distribu-

tions in terms of point charges is not quite appropriate if accuracy at

short distances between sources is required. In fact, the electron dis-

tributions in atoms have a substantial width and nearby distributions

will interpenetrate to a certain extent. The modi¬ed (damped) short-

range interactions are better represented by charge distributions than

by point charges. Exponential shapes (as from the distribution of

Slater-type orbitals) are the most appropriate.

(vii) Multipoles To increase the accuracy of the representation of charge

distributions while avoiding too many additional virtual sites, dipoles

“ and sometimes quadrupoles “ may be added to the description. The

disadvantage is that the equations of motion become more compli-

cated, as even for dipoles the force requires the computation of electric

¬eld gradients, due to charges, dipoles and quadrupoles. Dipoles and

quadrupoles are subjected to torques, which requires distribution of

6 The term “dummy” is not always reserved for sites as described here, but is often used to

indicate virtual sites as well.

6.3 Force ¬eld descriptions 159

potential force

0.9 1.0 1.1 0.9 1.0 1.1

r/rc r/rc

(a)

0.9 1.0 1.1 0.9 1.0 1.1

r/rc r/rc

(b)

0.9 1.0 1.1 0.9 1.0 1.1

r/rc r/rc

(c)

Figure 6.5 Truncation schemes for an r’6 dispersion potential. Thick lines: poten-

tial (left) and force (right) for (a) truncated potential, (b) truncated force = shifted

potential, (c) shifted force. Dashed lines give the correct potential and force.

forces over particles that de¬ne the multipole axes. Its use is not

recommended.

(viii) QM/MM methods combine energies and forces derived from quan-

tumchemical calculations for selected parts of the system with force

¬elds for the remainder (see Section 6.3.10).

(ix) Ab initio molecular dynamics applies DFT-derived forces during dy-

namics evolution: see Section 6.3.1.

6.3.4 Long-range dispersion interactions

The non-bonded potential and force calculations are usually based on lists

of atom or group pairs that contain only pairs within a given distance. The

reason for this is a computational one: if all pairwise interactions are in-

cluded, the algorithm has an N 2 complexity, which runs out of hand for

large systems. This implies the use of a cut-o¬ distance, beyond which

the interaction is either neglected, or treated in a di¬erent way that is less

computationally demanding than N 2 . For the r’6 dispersion potential such

cut-o¬s can be applied without gross errors, but for the r’1 Coulomb po-

tential simple cut-o¬s give unacceptable errors.

160 Molecular dynamics