Schlick T. Molecular Modeling and Simulation: An Interdisciplinary Guide

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

310 10. Nonbonded Computations

(see Fig. 10.2). This requires alteration of the nonbonded function over the larger

region r ≤ b (rather than 0 <a≤ r ≤ b).

Shift functions include the following two formulations:

(r)=[1− (r/b)

]

, for r ≤ b, (10.10)

(r)=[1− r/b ]

, for r ≤ b. (10.11)

Both decrease monotonically from 1 to 0 as r increases from 0 to b,buttheir

curvature is different (Fig. 10.3, part (c)). These two functions work well for

Coulomb interactions.

In fact, S

(r) above augments the true Coulombic force by a constant (1/b

)

for r ≤ b. This can be seen by writing the modiﬁed energy



E(r) of the scalar

Coulomb potential E(r)=1/r as



E(r) ≡ S

(r) · E(r)=



b − r



(b −r)

. (10.12)

The associated derivative is:





(r)=

− b

−

= E



(r)+

, (10.13)

showing a derivative augmentation by the constant 1/b

The potential shift approach by S

(r) is often termed force shift for this reason.

Note, however, that this name is misleading since S

is applied to the potential

and not to the force.

CHARMM, for example, employs S

(r) above for group-based potential shift

and S

(r) for atom-based shifts. For the van der Waals interactions, the potential

shift function is additive rather than multiplicative, applied as an auxiliary term

so as to dampen the force monotonically to zero:

(X)=





i,j



−A

+ C

+ D



r ≤ b

0 r>b

. (10.14)

Here, C and D are chosen so that both the van der Waals potential and force

functions are zero at r = b. This leads to the expressions in terms of A

and

as:

= −A

−12

+2B

−18

(10.15)

=2A

−6

− 3 B

−12

. (10.16)

See Figure 10.2 (lower panels) for an illustration of this shift function

approximation.

10.4. The Ewald Method 311

10.4 The Ewald Method

In this section, we only sketch the efﬁcient particle-mesh Ewald (PME) method,

popular for biomolecular dynamics. See [284,369,1413] for technical details, and

[71, 283, 1081], for example, for works which emphasize the impact of PME on

nucleic acid simulations. An outline of the Ewald summation technique can also

be found in [428].

Throughout this section, we use notation consistent with other works in this

subject (coming from crystallography), though different from notation used else-

where in this text. The symbol X for the collective Cartesian vector used

elsewhere is identical to the collective position vector x used here. The position

vector of atom i is x

, and the interparticle distance vector from atom j to i is

≡ x

− x

; we express the magnitude of this distance as the Euclidean norm

of the vector, or |r

10.4.1 Periodic Boundary Conditions

The Ewald method is a technique for calculating the electrostatic energy of a sys-

tem on a lattice with periodic boundary conditions. By periodic conditions, we

mean that the modeled system (biomolecules and solvent molecules) is placed in

the unit cell and considered to have inﬁnitely many images in space. This repli-

cation forms an inﬁnite lattice in 3D space. Figure 10.4 shows a standard square

lattice in 2D of dimension L; this lattice generalizes to a cubic geometry in 3D.

Space-Filling Polyhedra

Non-cubic periodic domains can also be modeled in 3D, such as the truncated

octahedron, hexagonal prism, and face-centered cube, as shown in Figure 10.5.

These lattice geometries reﬂect variations in the relationships among the three side

lengths (a, b, c) and the three angles between these sides (α, β, γ) of the domain.

For some commonly used space-ﬁlling lattices we have the following relations:

1. Cubic: a = b = c and α = β = γ =90

2. Hexagonal: (e.g., hexagonal prism): a = b and α = β =90

and

γ = 120

3. Rhombic Dodecahedron (e.g., face-centered cube): a = b = c and α =

γ =60

and β =90

More generally, other space-ﬁlling polyhedra are shown in Figure 10.6.

Minimum-Image Convention

Only coordinates of the unit cell need to be recorded and propagated. As an atom

leaves the unit cell by crossing the boundary, an image enters to replace it, and

hence the total number of particles is conserved. In biomolecular simulations,

the unit cell is usually set to be large enough to limit such occurrences and

312 10. Nonbonded Computations

60 A

Figure 10.4. Periodic domain in 2D, showing the unit cell (center) and its images (sur-

rounding replicas), and in 3D (rectangular), as used for a solvated protein (BPTI) system

[1089].

to avoid artifacts that can be caused by artiﬁcial periodic boundary conditions

[510,580,581,1203]. Such artifacts may be more pronounced in a solvent of low

dielectric permittivity (not water), for a solute having a large overall charge, and

when the solute cavity is relatively large compared to the size of the unit cell (due

to artifactually small solvation of extended solute). Speciﬁcally, it has been found

that artiﬁcial periodicity perturbs the potentials of mean force associated with

conformational equilibria of solvated biomolecules and can tend to overstabilize

compact folded forms. A reasonable cell size involves a solvated biomolecule

with at least a 10

A layer of water molecules [291].

In practical terms, a minimum-image convention is typically used so that each

atom i interacts only with the closest periodic image of the other N − 1 parti-

cles. In addition, a spherical cutoff (as described in the last section) is applied to

restrict this number of calculated interactions further. For a consistent combina-

tion of spherical cutoffs and the minimum-image convention, the cutoff distance

(b of last section) is at most L/2,whereL is the dimension of the side of the

box. Note from Table 10.1 that applying periodic boundary conditions only about

doubles the computational work with respect to using 12

A cutoffs (for a ﬁnite

system). Further, implementing the particle-mesh Ewald fast summation scheme

(to approximate all nonbonded interactions in the periodic model) requires only

slightly more work than the periodic summation with this cutoff.

Choice of Geometry

The geometry of the periodic domain used also affects the total number of pair-

wise interactions considered. Note from Figure 10.5, for example, that a truncated

octahedron may be more compact than a cubic lattice for a solvated protein, while

a hexagonal prism is appropriate for solvated DNA systems. Programs to opti-

mize such models, that is, to orient and solvate the solute macromolecule in a

given domain subject to a minimal water layer thickness (e.g., 10

A) so as to

10.4. The Ewald Method 313

77A

52A

65A

Truncated

Octahedron (24)

Hexagonal

Prism (12)

Rhombic

Dodecahedron (14)

Figure 10.5. Examples of nonrectangular periodic domains in 3D used for biomolecular

simulations: truncated octahedron, containing a solvated protein (villin) [338]; hexag-

onal prism, containing a solvated DNA dodecamer [1230]; and rhombic dodecahedron

(face-centered cube), containing a polymerase/DNA complex [1134]; water molecules (not

shown) ﬁll the domain. A side and top view is shown for each geometry, as well as a

space-ﬁlling rendering (at top) indicating the number of vertices in each case.

Squashed

Dodecahedron (18)

Elongated

Dodecahedron (18)

Triangular Prism

(6)

Figure 10.6. Three additional space-ﬁlling polyhedra are shown, with the number of

vertices for each domain indicated.

314 10. Nonbonded Computations

yield the lowest number of water molecules, were pioneered by Mihaly Mezei

[857] (see references and details for program Simulaid on the author’s website)

and recently generalized, with enhanced efﬁciency, for additional space-ﬁlling

polyhedra in the program PBCAID [1027].

10.4.2 Ewald Sum and Crystallography

The Ewald construct has its roots in the crystallographic community. There are

three speciﬁc commonalities between crystallography and the Ewald summation

for fast electrostatics:

1. Reciprocal Term. The notion of the reciprocal space is used in addi-

tion to the direct space. In crystallography, an Ewald sphere describes the

diffraction conditions in terms of the reciprocal rather than direct lattice.

The reciprocal lattice is an orthogonal system related to the orthogonal

system associated with the atoms in the unit cell (termed the real-space

lattice). This lattice is used to express the scattering vector in crystals

(see Box 10.1). In the Ewald summation, two terms arise for the compu-

tation of the electrostatic energy, for atom pairs in the direct lattice and

for atoms pairs corresponding to interactions with images of the unit-cell

atoms. The former is evaluated by direct calculation, and the latter (smooth

and long-range) by Fourier transforms.

2. Fast Fourier transforms. FFTs are key ingredients in both applications.

In crystallography, FFTs are used to express the amplitude of diffraction

for the whole crystal from the electron density function associated with

one unit cell of the crystal (see Box 10.1). In the Ewald sum, FFTs are

used to evaluate the electrostatic energy corresponding to the reciprocal

lattice.

3. Electron Density Distributions. The notion of electron density distri-

bution is key in both applications. Crystallographers compute the total

electron density to obtain a spatial arrangement of the atoms in space in

terms of FFTs of the structure factors, which describe the scattering pattern

of each atom with respect to an electron at the cell origin (see Box 10.1).

Modelers use an analogous concept to describe the electrostatic energy

in terms of charge densities rather than point charges to make the inﬁnite

sums convergent. Terms analogous to the crystallographic structure factors

also arise in the Ewald sum FFTs.

Box 10.1: Some Tricks from X-Ray Crystallography

A crystal is a 3D periodic arrangement of atoms (see [1047], for example, for an in-

troduction to crystallography, and Chapter 1 for an introduction to the phase problem).

10.4. The Ewald Method 315

The structural subunit of crystals is the unit cell. The X-ray diffraction pattern by the elec-

trons of a molecular crystalline system is interpreted in terms of both the direct lattice,the

orthogonal system associated with the unit-cell atoms, and the reciprocal lattice, an aux-

iliary orthogonal system used to express the scattering vector in crystals. Speciﬁcally, the

scattering of each atom j with respect to an electron at the origin of the cell it occupies is:

= f

exp (2πix

· s), (10.17)

where i =

√

−1, f

is an experimentally-measured amplitude, x

= {x

}

denotes the Cartesian position vector of atom j in the direct lattice, and s is the scattering

vector corresponding to a triplet of indices (indices of reﬂection), such as s = {h, k, l},

in the reciprocal lattice. Thus for N atoms in the unit cell, we write the complex-valued

structure factor F

as:



j=1

exp (2πix

· s)=A

+ iB

, (10.18)



j=1

cos (2πix

· s) ,B



j=1

sin (2πix

· s) . (10.19)

From these quantities, the phase angles associated with F

are computed as φ

tan

−1

). The intensity of the diffracted X-ray beam is proportional to the square

of the structure factor’s amplitude and depends only on the interatomic vectors, not on the

actual, origin-dependent atomic coordinates:

= A

+ B



i=1



j=1

cos [2π (x

− x

) · s] . (10.20)

The relationship between crystallography and the electrostatic potential arises by relating

the diffracted amplitude to a continuous distribution of electron density, ρ(x), expressed

as electrons per unit volume. The scattered amplitude for the whole crystal, F (s), from a

small volume element dv, results from an effective point charge of ρ(x) dv electrons:

F (s)=



ρ(x)exp(2πis · x) dv , (10.21)

where the integration is taken over the volume V of the unit cell. This is the Fourier trans-

form of the electron density in one unit cell of the crystal, sampled at points s on the

reciprocal lattice. The inverse transformation yields the electron density via integration

over the entire volume of the reciprocal space in which s is deﬁned:

ρ(x)=



∗

F (s)exp(−2πis · x) dv

∗

. (10.22)

Crystallographers compute a discrete analog of this electron density to obtain the atomic

arrangement in the crystal in which each point of the reciprocal lattice (where F (s) is

deﬁned) has an associated weight of F

/V where F

is the structure factor deﬁned in

eq. (10.18):

ρ(x)=

∞



h,k,l=−∞

exp (−2πis · x) . (10.23)

316 10. Nonbonded Computations

This computation is performed using Fast Fourier transforms. The amplitudes of these

structure factors are known from the observed intensities of the X-ray reﬂections, but

the phase angles cannot be measured. This is the heart of the problem of deducing structure

from X-ray diffraction patterns.

10.4.3 Mathematical Morphing of a Conditionally

Convergent Sum

Coulomb Energy in Periodic Domains

The following Ewald sum describes the total Coulomb energy corresponding to a

system in an inﬁnite periodic domain; the conversion factor K

coul

(see eq. (9.41)

of Chapter 9) is omitted for simplicity, and the dielectric constant  =1:

coul



i,j=1



images |n|

+ n|



j=1

Φ(x

) , (10.24)

where

Φ(x



i=1



|n|

+ n|

. (10.25)

The absolute value signs above denote vector magnitudes, and the summation

extends over all images of the unit cell. The vector n denotes the vector n =

L, n

L),wheren

,andn

are integers and L is the box size, and

where the triplet of indices (n

) are related to the magnitude of n, namely

|n|,by|n| = |n

|+ |n

Thus, for example, the ‘image’ corresponding to |n| =0has only one triplet:

n =(0, 0, 0); |n| =1has all triplets that can be written as: n =(±L, 0, 0),

(0, ±L, 0), and (0, 0, ±L);and|n| =2is associated with the vectors that

can be expressed as n =(±2 L, 0, 0), (0, ±2 L, 0), (0, 0, ±2 L), (±L, ±L, 0),

(±L, 0, ±L),and(0, ±L, ±L). Similarly, the reciprocal lattice vectors |m|

deﬁned below correspond to vectors (L/n

,L/n

The sum over the images extends from |n| =0to ∞. The prime symbol in

the sum (





) indicates that for n =0we omit the i = j interaction (so the

denominator is well deﬁned).

Conditional Convergence

The sum in eq. (10.24) is only conditionally convergent since the terms decay

as 1/|n|, like the harmonic series.

The value of the sum also depends on the

AseriesS =



∞

n=1

is conditionally convergent if the inﬁnite sum



converges but



| diverges. It is also known that the sum for a conditionally convergent series depends on the

10.4. The Ewald Method 317

nature of the surrounding medium (dielectric constant )[293]. For a unit cell

surrounded by vacuum, a dipolar layer exists on the surface and its contribution

must be canceled.

Ewald’s Trick

Ewald noted a trick that can be used to convert this sum into an expression

involving a sum of two absolutely and rapidly convergent series in direct and re-

ciprocal space. This conversion is accomplished by representing each point charge

as a Gaussian charge density, producing an exponentially decaying function.

This Gaussian transformation must be counteracted by an analogous subtraction

to leave the net result of an effective point charge (delta function charge, see

Fig. 10.7). This canceling distribution is summed in the reciprocal space (and

transformed back to real space) because it is a smoothly-varying periodic function

which can be represented by a rapidly convergent Fourier series.

Essentially, the basic idea can be written for the radial function Φ(r) as the

splitting

Φ(r) ≡

=Φ

real

(r)+Φ

recip

(r) , (10.26)

where

real

(r)=

−

erf(βr)

erfc(βr)

, (10.27)

and

recip

(r)=

erf(βr)

. (10.28)

The ﬁrst term (Φ

real

) is short-range with singularity at the origin, and the second

(Φ

recip

) is a smooth and long-range term (which can be Fourier transformed). The

error function erf is deﬁned as

erf(x)=

√



exp(−t

) dt, (10.29)

√



x −

−

+ ···



, (10.30)

and the complementary error function erfc is

erfc(x)=1− erf(x) . (10.31)

order of summation. For a

=(−1)

n+1

/n, the ‘alternating harmonic series’ 1 −

−

+ ···

converges, but the harmonic series 1+

+ ···diverges, though a

→ 0 as n →∞.Tosee

that the sum for ﬁnite n terms can be as large as we please for the harmonic series, note that terms can

be grouped so that each subsum is larger than



∞

n=1

=1+









(next 8 terms) + ···> 1+

+ ···.

318 10. Nonbonded Computations

In practice, erf(x) can also be evaluated from its relationships to other known

functions, like incomplete gamma functions [1015] or the normal probability

function f(t)=exp(− t

/2)/

√

2π,since

erf(x/

√

2) = 2



f(t) dt ,

and this integral (area under the standard normal curve from 0 to x)isoften

tabulated (see CRC Standard Mathematical Tables).

The Screening Gaussian ρ

More generally (for the periodic potential), we perform this splitting using a

screening spherical Gaussian ρ

, centered at each point charge, parameterized

by β. The parameter β controls the width of the distribution (width =

√

2/β)and

also the rate of convergence of the sum (see below). For each point charge, the

surrounding distribution is given by:

(x)=− q



√



exp [−β

|x|

] . (10.32)

Thus, instead of the cumulative point charge density described by a sum of delta

functions:

ρ(x)=



j=1

δ(x − x

) , (10.33)

we use a sum of localized densities based on the Gaussian of eq. (10.32)todeﬁne

the total Gaussian screening charge density ρ

(x):

(x) ≡



j=1

(x)=−



j=1



√



exp [−β

|x − x

] . (10.34)

The real-space sum can be written in terms of the complementary error func-

tion using Poisson’s equation

for the electrostatic potential, due to the Gaussian

charge cloud followed by integration [428]. (See Box 10.2 for a deﬁnition of the

divergence and Laplace operators). Namely, the solution of Poisson’s equation in

the form ∇

(r)=−4πρ

(r) leads to Φ

(r)=

erf(βr) after two steps

of integration. Thus, the sum in eq. (10.24) is converted into the pair of terms:

coul



j=1



|n|

erfc(β|x

+ n|)

+ n|



j=1



|n|

erf(β|x

+ n|)

+ n|

(10.35)

Simeon Denis Poisson (1781–1840) was a brilliant mathematician whose name frequently ap-

pears in text books. The Poisson equation (1812) in potential theory is a result of his discovery that

Laplace’s equation for the gravitational force holds only at points where no mass is located (see also

Box 10.2).

10.4. The Ewald Method 319

= 20 / L

= 10 / L

(x)

Figure 10.7. Schematic illustration of Ewald’s trick for converting a conditionally con-

vergent sum into two terms that are rapidly and absolutely convergent: each point charge

(represented by delta functions, dashed spikes) is masked by a Gaussian function; corre-

spondingly, a canceling Gaussian function must be added (not shown). The charge density

when point charges are screened by Gaussians is shown for four point charges for two

values of β, the parameter controlling the width of the Gaussian distribution.

The ﬁrst term (E

real

) is short-range, and the second (E

recip

) is the smooth long-

range component of the Coulomb potential.

Box 10.2: The Divergence and Laplacian Operators

The divergence symbol, as used in the Poisson and Poisson-Boltzmann equations, is a dif-

ferential operator on a continuously differentiable vector w ∈R

whose components {w

}

are functions of the coordinate variable x with components {x

}∈R

. The divergence

of w is written as div w or ∇·w and deﬁned as:

div w ≡∇·w ≡

∂w

(x)

∂x

∂w

(x)

∂x

+ ···+

∂w

(x)

∂x

For example, for the gravitational ﬁeld vector v(x) at a point x in space whose distance

from the origin is r,wehavev(x)=−(GR

) x where R

is the earth’s radius and G

is the gravitational constant. For r>R

, it can be shown that the divergence of v is zero:

∇·v(x)=−

∇·x + ∇



−GR



· x

= −GR



∇·x +





∇r · x



=0,

since ∇·x =3and ∇r · x =(x · x)/r = r. The corresponding potential function

u(x)=GR

/r,where∇u = v(x), thus satisﬁes ∇·∇u = ∇·v =0or Laplace’s