Schlick T. Molecular Modeling and Simulation: An Interdisciplinary Guide

Подождите немного. Документ загружается.

290 9. Force Fields

The CHARMM and AMBER force ﬁelds employ this energy-minimum / dis-

tance procedure outlined above, using the following deﬁnitions for V

and r

For each atom type, the pair {

} is speciﬁed as a result of the parameteriza-

tion procedure (crystal ﬁtting or liquid simulations). The parameters for a given ij

interaction are then obtained by setting V

and r

as the geometric and arithmetic

(see below) means, respectively:

√



, (9.36)

=(r

+ r

) , (9.37)

and then using eqs. (9.33)and(9.34) to deﬁne the attractive and repulsive coefﬁ-

cients, respectively. The presence/absence of the

factor in the eq. (9.37) depends

on the original deﬁnition of {r

} (i.e., given radii may have already been divided

by 2).

Slater-Kirkwood Procedure (A

)

The second parameterization procedure is based on the Slater-Kirkwood equation.

The attractive coefﬁcient A

is determined on the basis of the atomic properties

of the interacting atom pair [701,898,973]:

365 α

(α

)

1/2

+(α

)

1/2

, (9.38)

where α represents the experimentally determined atomic polarizability, and N

denotes the number of outer-shell electrons. The numerical factor in the Slater-

Kirkwood equation (that is, 365) is derived from universal constants so that the

energies are produced in kcal/mol. This factor is computed from the relation

3e/2m

,wheree = electron charge,  = Planck’s constant divided by 2π,and

= electron rest mass.

Slater and Kirkwood have shown, based on experiments for noble gases, that

this derivation for A

produces a London attraction potential in good agreement

with experiment.

When A

is obtained from the Slater-Kirkwood equation (9.38), the coefﬁcient

of the repulsive term can then be obtained from A

using equation (9.35).

This leads to a combined van der Waals potential that produces an energy

minimum at r

The MMFF force ﬁeld van der Waals parameterization employs the Slater-

Kirkwood equation to deﬁne V

via V

= A

/[2(r

)

] (from eq. (9.33)) where

the minimum-energy separation r

is determined from special combination rules.

In addition, a “buffered 7/14 attraction/repulsion” van der Waals term is used,

found to better ﬁt rare gas interactions [498]. The form of this potential is:

MMFF

−A

+ γ

)

+ γ

)

+ γ

)

]

, (9.39)

where the ij-dependent A and B values are functions of V

and r

,andγ

and

are constants.

9.7. The Coulomb Potential 291

Such more complex nonbonded terms (also the 6/exponential combination

used by the MM2/3/4 force ﬁeld) represent a tradeoff between the accuracy

of ﬁtting and the computational expense involved in potential evaluation and

differentiation.

9.7 The Coulomb Potential

9.7.1 Coulomb’s Law: Slowly Decaying Potential

Ionic interactions between fully or partially charged groups can be approximated

by Coulomb’s law

for each atom pair {i, j}:

coul

) ∝ q

, (9.40)

where F is the force and q

is the effective charge on atom i. The force is positive

if particles have the same charge (repulsion) and negative if they have opposite

signs (attraction). This leads to the Coulomb potential of form

coul

= K

coul

r

, (9.41)

2 3 4 5 6 7 8 9

−0.2

−0.1

0.1

0.2

0.3

0.4

0.5

LJ Energy [kcal/mol]

C...C

O...H

r (distance) [A]

Figure 9.12. Van der Waals nonbonded potentials for C...C and water O...H distances

using the CHARMM parameters −0.055, −0.1521,and−0.046 for 

and 2.1750,

1.7682,and0.2245 for the distances r

. The coefﬁcients A and B are computed as:

=2(r

+ r

)

√



, B

=(r

+ r

)

√



, where coefﬁcients produce energies

in kcal/mol.

Charles Augustin de Coulomb (1736–1806) was a French physicist who formulated around 1785

the famous inverse square law, now named after him. This result was anticipated 20 years earlier by

Joseph Priestly, one of the discoverers of oxygen and author of a comprehensive book on electricity.

292 9. Force Fields

where  is the dielectric constant and K

coul

is the conversion factor needed to

obtain energies in units of kcal/mol with the charge units used (see Box 9.1).

Unlike the van der Waals potential, the Coulomb interactions decay slowly with

distance (see Figure 9.13). In fact, electrostatic interactions are important for sta-

bilizing biomolecular conformations in solvent and associating distant residues

in the primary sequence closer in the folded structure. However, because of the

O(N

) complexity in evaluating all pairwise terms directly, prior calculations

have often only counted interactions within a cutoff radius (e.g., 10

A).

To account for all long-range electrostatic interactions in macromolecules, im-

plementation of fast-electrostatics algorithms, such as the multipole method and

Ewald methods, are needed (see next chapter). Such methods can reduce the com-

putational complexity to O(N). The Ewald and Particle Mesh Ewald treatments

have been found to be more amenable to biomolecular dynamics simulations (e.g.,

[751,1133]), especially for parallel platforms; this trend may change as hardware

and software evolve. See next chapter for details.

For details of many concepts in electrostatic energies of macromolecules,

see Warshel’s textbook [1342]; also notable is the work [1344] which intro-

duces consistent electrostatic calculations in proteins as well a quantum/molecular

mechanics approach, and the comprehensive review [1346], which relates micro-

scopic and macroscopic aspects of biomolecules, such as dielectric constants and

dipole moments.

9.7.2 Dielectric Function

The reduction of the force or energy by a dimensionless factor of 1/ is appro-

priate if the charged particles are immersed in any medium other than a vacuum.

That is, a weaker interaction occurs in a polarizable medium (e.g., water) than in

vacuum, since then charges become “screened”. The distance dependent dielectric

function was ﬁrst introduced in the late 1970s [1344]. Very simple approximations

to the dielectric function, such as (r)=r or (r)=

Dr where

D is a constant,

have been used in ﬁrst-generation force ﬁelds when water molecules were not

represented explicitly in the model [1359].

Sigmoidal Function

More sophisticated formulations rely on a large body of theoretical and exper-

imental data, which suggest screening functions in the general sigmoidal form

[528,548,854,1420] that can also account for different ionic concentrations. Such

a sigmoidal function uses a screening parameter κ related to the ionic strength of

the medium [1215]:

(r)=

D exp(κr) ,

where κ has dimensions of inverse length. Other forms employ distance-

dependent dielectric functions like

(r)=(D

+ D

)/(1 + k exp[−κ(D

+ D

)r]) − D

9.7. The Coulomb Potential 293

where D

is the permittivity of water, D

and κ are parameters, and k is a con-

stant [855]. These screened Coulomb potentials simultaneously model two effects

between two charges in a dielectric medium, like water: a substantially-damped

electrostatic ﬁeld when the charges are separated by moderate distances (e.g.,

several

Angstroms apart), and the diminished dielectric screening, approaching

toward vacuum values, when the charges are in close proximity.

Such screened Coulomb potentials have been used for molecular dynamics sim-

ulations and have been incorporated into a procedure to calculate pH-dependent

properties in proteins (pK

shifts, i.e., shifts in hydrogen dissociation constants)

with excellent accuracy [855, 1346].

The sigmoidal functions are more suitable than the simpler forms above for

implicit solvation models of biomolecules. Still, for detailed structural and ther-

modynamic studies of biomolecules, the trend has been to model water molecules

explicitly (i.e., unit value ). Nonetheless, screened Coulomb potentials are of

great usage in other applications, such as macroscopic simulations of long DNA in

a given monovalent ionic concentration [1125, for example]. In such applications,

the coefﬁcient

D is a function of the salt concentration through a salt-dependent

linear charge density for DNA, and the screening parameter κ is the inverse of the

Debye length [197].

See end of Chapter 10 (section on continuum solvation) and

[107], for example, for further details.

Box 9.1: Coulomb Potential Constant (K

coul

)

Chemists typically use an electrostatic charge unit (esu) to deﬁne the unit of charge. This

unit is deﬁned as the charge repelled by a force of 1 dyne when two equal charges are

separated by 1 cm. In these units, the electron charge is 4.80325 ×10

−10

esu.

In cgs units, the constant of proportionality in the Coulomb energy is unity; in the SI inter-

national system of units, the unit of charge is the coulomb (C) = 1 Ampere second (As).

In coulomb units, the electron charge is 1.6022 ×10

−19

C. The corresponding constant of

proportionality in Coulomb’s law is:

coul

=1/(4π

) , (9.42)

where the permittivity of a vacuum 



=8.8542 × 10

−12

−1

−3

=8.8542 × 10

−12

−1

The electrostatic energy is typically expressed as the sum over pairwise interactions between hy-

drodynamic DNA beads separated by segments of length l

as: [(νl

)

/]



i<j

[exp(−κr

)/r

where ν is an effective linear charge density of the DNA,  is the dielectric constant of water, and

1/κ is the Debye length at the monovalent salt concentration c

, given in Molar units (κ ≈ 0.33

√

inverse

Angstrom units at room temperature for 1:1 electrolyte solutions).

294 9. Force Fields

0 0.1 0.2 0.3 0.4 0.5

−5000

−4000

−3000

−2000

−1000

1000

2000

3000

4000

5000

Coulombic energy [kcal/mol]

r (distance) [A]

C...C

O...H

Figure 9.13. Coulomb potentials K

coul

for C...C and O...H interactions using the

CHARMM partial charge parameters (in esu units) of −0.1800 (Arg carbon), −0.8340

(water oxygen), and 0.4170 (water hydrogen).

Thus, to obtain energy in kcal/mol using

Angstroms for distance and esu for partial charges,

the conversion constant is

coul

1Jm

4π (8.8542 × 10

−12

(6.0221 × 10

mol

−1

)(10

A) (4184

−1

kcal) (1.6022 × 10

−19

)

4π (8.8542 × 10

−12

)esu

≈ 332

kcal

mol

esu

. (9.43)

9.7.3 Partial Charges

Selection of partial charges to ﬁt atom centers has been a difﬁcult issue in molec-

ular mechanics calculations, since different quantum mechanical approaches

yield signiﬁcantly different values [1195]; earlier, determinations of atomic par-

tial charges based on X-ray diffraction data were suggested [980]. However,

extensive developments in high-level quantum-derived electrostatic potentials

[218, 265, 805] are being applied to determine partial atomic charges for iso-

lated model compounds, with various adjustments to correct for the neglect of

many-body polarization effects in liquid water and other factors (e.g., [1330]).

9.8. Parameterization 295

For reference, the MMFF force ﬁeld uses a simple modiﬁcation of the Coulomb

potential of form [498]

MMFF

coul

)=K

coul

 (r

+ δ)

, (9.44)

where  is the dielectric constant (unity by default), δ =0.05

A(thebuffering

constant), and n is 1 (default) or 2 (distance-dependent model). This form is re-

quired to dampen the attraction between oppositely charged atoms in combination

with the buffered van der Waals term used.

9.8 Parameterization

9.8.1 A Package Deal

The general parameterization process for potential energy functions is a difﬁcult

task. Several important decisions must be made regarding choices for the func-

tional form and numerical values for the parameters. Even if one is given a speciﬁc

energy form and a set of structural and energetic data to reproduce, the combina-

tions of parameters that can be used are endless. Unrealistic choices for one group

of parameters can be compensated for by adjustment of another.

In theory, the energy terms should have clear physical signiﬁcance with pa-

rameters calibrated by empirical ﬁtting of crystal data, rotational barriers of

analogous small molecules, and vibrational frequencies. However, an approxi-

mation is inherent in the extension of data from small to large systems. Moreover,

interaction with solvent and counterions, reﬂected in the experimental data, must

be interpreted and incorporated in the energy model. In summary, much freedom

and manipulation are possible in constructing empirical energy surfaces. Only if

constructed and parameterized correctly will the energy model generate reliable

structural predictions, as reviewed in [620,802].

In the case of nucleic acid sugars, the importance of parameter choices in the

energy function has already been realized. For example, different potential energy

models have produced results that are qualitatively different regardingsugar pseu-

dorotation [522,750,909,937,1359]. This is particularly possible when choosing

appropriate equilibrium values for endocyclic bond angles [937] or torsion-angle

parameters about sugar atoms [218]. Since the unusual puckering geometry and

ring closure constraints produce signiﬁcant deviations from tetrahedral bond an-

gle arrangements, it is not clear what equilibrium values should be used in the

harmonic bending terms, more appropriate for small ﬂuctuations.

9.8.2 Force Field Comparisons

Some force-ﬁeld dependent conformations for DNA have been discussed with re-

spect to the AMBER and CHARMM force ﬁelds [379,381,415,803,804,1330].

296 9. Force Fields

With some earlier versions, average structures tended to be more B-like with

AMBER and intermediate between A and B-DNA with CHARMM in large part

due to differences in sugar conformations; newer force ﬁelds better balance lo-

cal geometry with global helical propensities and can reproduce the equilibrium

between A and B-DNA in solution (both when Ewald and nonbonded cutoffs are

used for electrostatics) [221,804,1330]. Differences in backbone and base geome-

tries, as well as dynamic properties, are also believed to be force-ﬁeld dependent.

Such disparities, rectiﬁable with improved parameters, warrant a particularly cau-

tionary note in interpreting structural transitions (such as between A and B-DNA

[219]) in molecular dynamics simulations.

With recent improvements in molecular dynamics algorithms, enhanced

sampling methods, computational speed, and rapid performance on parallel archi-

tecture, long-time dynamics simulations of proteins and DNA, as well as folding

studies, have become possible. Such simulations have revealed force ﬁeld inac-

curacies not observed in shorter-length simulations. For example, AMBER was

found to over-stabilize α-helices [445, 934], while CHARMM tends to prefer

π-helices [378]. More generally, evidence for conformational discrepancies was

presented over a wide range of protein force ﬁelds [882,1411].

To overcome these problems, peptide backbone parameters were reparameter-

ized for AMBER [568], CHARMM [806, 807], and OPLS-AA [629]. A new

AMBER force ﬁeld was also created for protein simulations [341]. Parameter

improvements for nucleic acids in AMBER [988] and GROMOS [1206]andfor

lipid parameters in CHARMM [556, 652] were also reported.

Still, recent long-time peptide folding studies have raised concerns about con-

formational biases in existing force ﬁelds (e.g., overstabilization of α helices) that

even the most recent force ﬁeld corrections have not resolved [132, 424, 426]. In-

deed, it is impossible for any force ﬁeld to accurately reproduce all the complex

interactions and properties of real systems, and this stems not only from force

ﬁeld approximations but also from hardware and software limitations, time con-

straints, convergence issues, and the limited resolution or possible errors of some

experimental data. Nonetheless, it is interesting that very different computational

approaches (all atom on one hand versus coarse grained with implicit solvation)

lead to similar results (e.g., folded structures with high α-helical content for a

β protein) that depart from experimental predictions [365, 424, 426, 841]; this

can either indicate strong force ﬁeld biases that persist through various levels

of approximation or perhaps physical explanations that modeling can illuminate

and eventually help resolve. Recent modeling of folding processes has certainly

indicated that conformational space may be more heterogeneous that originally

believed [365,427].

Much work continues to further reﬁne and develop force ﬁelds to mirror more

accurately complex systems and produce realistic simulation data for biological

systems over longer time-scales to resolve such issues.

9.8. Parameterization 297

9.8.3 Force Field Performance

Several discussions of force ﬁeld performance [426, 455, 668, 1039, 1068,

1069, 1172] highlight general issues that remain unresolved and in need of

improvement:

• Determination of partial charges;

• Improvement of electrostatic potentials (e.g., use of distributed multipole

analysis);

• Methods for solvent representation;

• Interpretation of results in the absence of solvent;

• The approximation reﬂected by Cartesian vs. torsion space representation;

and

• General interpretation of conﬂicting results by different models and

potentials.

Frequent comparisons among force-ﬁeld results with respect to experimental and

high-accuracy ab initio data, such as done in [105,500], reveal the sizable errors

made by most force ﬁelds and inherent deﬁciencies. Improvements in the elec-

trostatic formulation that reﬁnes the simple atom-centered charged models are

repeatedly urged by computational chemists.

As force ﬁelds improve, some of these issues may be resolved, but force ﬁelds

cannot be said to be “converging”to one another [1123]. Still, modern force ﬁelds

perform comparably in molecular dynamics simulations [1016].

Still, I emphasize that force ﬁelds need not be perfect to be useful! Their overall

utility is in generating qualitative and quantitative insights into structural, ener-

getic, and dynamic properties of complex systems through systematic studies (for

example by varying critical parameters), especially when trends are compared for

a related group of biological systems (like single base-pair variants of DNA).

Ultimately, predictions or interpretations can be tested experimentally.

298 9. Force Fields

3’

5’

Nonbonded Computations

Chapter 10 Notation

YMBOL DEFINITION

Matrices

D diffusion tensor (also deﬁned for Coulomb manipulations)

Hessian of potential function

mass matrix

hydrodynamic tensor

Vectors

m reciprocal lattice vector, components {L/n

,L/n

}

periodic domain vector, components {n

L, n

arbitrary vector (deﬁned for matrix/vector product)

interatomic distance vector

scattering vector, components {h, k, l} (indices of

reﬂection)

v(x)

gravitational ﬁeld vector at point x

position vector of atom j, components x

force term

complex-valued structure factor, F

= A

+ iB

scattering vector at atom j

nonbonded force component, −∇E

R Langevin random force

X, x

collective position vector

dipole moment

∇Φ

gradient of potential function

Scalars & Functions

a, b, c, ˜c distance parameters

hydrodynamic radius

ionic concentration

T. Schlick, Molecular Modeling and Simulation: An Interdisciplinary Guide, 299

Interdisciplinary Applied Mathematics 21, DOI 10.1007/978-1-4419-6351-2

10,

 Springer Science+Business Media, LLC 2010