Schlick T. Molecular Modeling and Simulation: An Interdisciplinary Guide

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

320 10. Nonbonded Computations

equation, also written as ∇

u =0or Δu =0.TheLaplacian is a differential operator

from scalar to scalar ﬁelds whereas the divergence is a differential operator from vector to

scalar ﬁelds.

10.4.4 Finite-Dielectric Correction

The decomposition above, while convergent, is strictly correct only for an inﬁnite

dielectric ( = ∞) medium or for unit cells with a zero dipole moment. The cor-

rection term (E

cor,

) for a nonuniform ﬁeld associated with a macroscopic crystal

in a dielectric continuum with external dielectric constant  was only derived 60

years after Ewald’s original derivation [293], and yields:

cor,

= E

coul

( =1)− E

coul

( = ∞)=

2π





j=1



. (10.36)

Though reliable, Ewald’s algorithm as corrected in [293] was still O(N

) in com-

putational complexity. This is because the long-range reciprocal-space Fourier

sum requires O(N

) to be a sufﬁciently accurate approximation for large β (see

below).

10.4.5 Ewald Sum Complexity

The key breakthroughs for reducing the computational complexity of the Ewald

sum came in two steps.

Optimization of β

First, by optimizing the parameter β that controls the width of the screen-

ing Gaussians (see Figure 10.7), the relative convergence rates of the real and

reciprocal-space series are adjusted to optimize the work involved [402]. Namely,

as β increases, the real-space sum converges more rapidly and the reciprocal

sum more slowly. Both sums are truncated in practice (ﬁnite number of terms).

Thus, for example, a sufﬁciently large β yields an accurate direct-space sum

with an appropriate truncation, resulting in O(N ) work for the direct-space sum

rather than O(N

); however, the reciprocal-space sum still requires O(N

) work.

The optimal work balance between the two components can yield an overall

O(N

3/2

) method by adjusting β [991]. This is much better than O(N

) but still

considerably expensive for large biomolecular systems.

Mesh Interpolation

The second breakthrough in the Ewald sum came by noting that the trigonometric-

function values in the Fourier series used to represent the reciprocal-space term

10.4. The Ewald Method 321

can be evaluated through a smooth interpolation of the potential over a regular

grid. The resulting particle-mesh Ewald (PME) method reduces the overall com-

putation to O(N log N). The smoothing can be done by Lagrange interpolation

[284] or by B-spline interpolation [369].

Variations

Credit for this interpolation idea is due to Hockney and Eastwood who devel-

oped several methods in the early 1970s for simulations of hot gas plasmas in

fusion machines. Included is a ‘particle-particle particle-mesh’ (P

M) scheme, as

detailed in their text [554]. The P

M method for evaluating long-range forces in

large systems splits the interparticle force summation into a short-range, rapidly-

varying part and a smooth, slowly-varying remainder. A direct ‘particle-particle’

sum is used to compute the former, while a ‘particle-mesh’ interpolation proce-

dure is used to approximate the latter on a uniform grid. A ‘Q-minimizing’ method

is also used to optimize the scheme’s parameters (given a desired accuracy), such

as mesh size, cutoff radius, and charge assignment scheme.

The P

M method and its cousins are thus closely related to the fast PME

schemes used in biomolecular dynamics. Application of the P

M scheme as de-

scribed in [554] to molecular dynamics has also demonstrated good performance

[792, 793], with guidelines for the scheme’s parameters obtained by optimizing

an auxiliary function.

10.4.6 Resulting Ewald Summation

Following a series of mathematical manipulations (see [1267], for example), the

resulting composition of the Ewald summation for E

coul

of eq. (10.24)at =1

has ﬁve terms:

coul

= E

real

+ E

recip

+ E

cor,self

+ E

cor,ex

+ E

cor,

. (10.37)

The ﬁrst term, E

real

, corresponds to a real-space (or direct) sum of the electro-

static energy due to the point charges screened by oppositely-charged Gaussians.

The second sum, E

recip

,istheassociated canceling term (periodic sum

of Gaussians) summed in reciprocal space using smooth interpolation of

Fourier-series values. The last three terms are correction terms.

The ﬁrst correction (which is position independent) subtracts the self-

interaction term (each point charge and its corresponding Gaussian charge cloud)

which is included in the ﬁrst two terms. The second correction subtracts the

Coulomb contribution from the nonbonded pairs excluded from the Coulomb

energy (denoted as pairs i, j ∈ Ex below) since they are separately accounted for

in bonded, bond-angle, and possibly dihedral-angle terms. The third correction

accounts for the non-inﬁnite dielectric medium (eq. (10.36)).

322 10. Nonbonded Computations

Following algebraic manipulations, the resulting Ewald sum can be expressed

as follows:

real



i,j=1



|n|

erfc (β |r

+ n|)

+ n|

; (10.38)

recip

2πL



|m|=0

exp(−π

|m|

/β

)

|m|

S(m) S(−m) , (10.39)

where S(m)=



j=1

exp [2πim · x

]; (10.40)

cor,self

−β

√



j=1

; (10.41)

cor,ex

= −



i,j∈Ex

erf (β |r

; (10.42)

cor,

2π

(1 + 2) L





j=1



. (10.43)

The erfc function decays to zero with increasing values of the independent

variable. If β is sufﬁciently large, the only term that contributes to the real-space

sum is for |n| =0, and it can be computed in practice with a spherical cutoff of

order b =10

The reciprocal space or Fourier term corresponds to a summation over the re-

ciprocal vectors m , where the Fourier terms S(m) are charge-weighted structure

factors (see Box 10.1). Eqs. (10.39)and(10.40) can also be written as:

recip

2πL



i,j=1



|m|=0

exp(−π

|m|

/β

)

|m|

exp[2πim · (x

− x

)] .

(10.44)

When the charges are interpolated over a uniform set of grid points, these structure

factors are easily computed by FFTs.

10.4.7 Practical Implementation: Parameters, Accuracy,

and Optimization

Gaussian Width

In practice, the Gaussian width parameter β is determined so that the real-space

term achieves a desired accuracy tolerance. In typical solvated proteins or DNA

simulations, β has an order of magnitude of roughly 10/L (L ranges roughly

10.4. The Ewald Method 323

from 60 to 100

A; values exceeding 70

A have been implicated with smaller

artifacts from the enforced periodicity [580]). The effective cutoff used for the

real-space interactions is around 10 to 12

A. When multiple-timestep schemes

are implemented for molecular dynamics, the parameter β (or the cutoff for the

direct-space term) may be further optimized to distribute the work for the real

and reciprocal terms appropriately, so as to yield the greatest overall speedup; see

[97], for example.

Grid Size and Accuracy

The reciprocal-space uses multidimensional piecewise interpolation (e.g.,

B-splines) for evaluating the Fourier terms, with grid size and number of terms

chosen to achieve the desired accuracy. For instance, moderate accuracy (e.g.,

−4

relative force error) might be achieved with a coarse interpolation grid

(1 to 2

A). Very high accuracy (e.g., 10

−10

relative force error) might be obtained

with a ﬁner grid (∼0.5

A). The reciprocal-space energy and force terms are

expressed as convolutions and can thus be evaluated quickly using FFTs.

Work

has also shown that the ﬁnite number of wave vectors used in the discrete ap-

proximation gives rise to truncation errors due to the exclusion of intramolecular

interactions [1019]; this, and the related existence of fast terms in the reciprocal-

force component [97, 98, 1025, 1236], create a problem for development of

efﬁcient multiple-timestep integrators for molecular dynamics simulations using

PME approximations [89].

Computer Architecture Considerations

In the implementation of the PME method on multiprocessors, the cutoff ra-

dius used for the direct sum may be increased to balance the work between

the real-space and reciprocal-space components. This accelerates the computa-

tions because 3D FFTs are challenging to parallelize well. In contrast, the direct

sum parallelizes easily by spatial decomposition. Using a larger cutoff in the

real-space sum reduces the number of lattice vectors (Fourier terms) needed

for the same accuracy in the reciprocal sum and therefore improves the overall

performance.

Experience to date shows that PME implementations for biomolecules —

especially when tightly adapted to the computer architecture — are very fast; the

best implementations require about the same work as needed to evaluate the same

periodic version of the potential but with cutoffs in the range of 10–12

A[1133].

Figure 10.8 shows performance times for an efﬁciently distributed PME code by

John Board and co-workers at Duke University for a huge water system.

For two functions of time f (t) and g(t) and corresponding Fourier transforms F (f) and F (g),

we deﬁne the convolution of the two original functions f and g, f ∗ g,as:f ∗g =



∞

−∞

f(τ) g(t −

τ)dτ. It can be shown that F (f ∗ g)=F (f) F (g). That is, the Fourier transform of the convolution

of two functions (f ∗ g) is just the product of the individual Fourier transforms of those functions.

324 10. Nonbonded Computations

0 1 2 3 4 5 6 7 8 9

All−pairs electrostatics for 23,832 water molecules (71,496 atoms)

processors (450 MHz Pentium II)

time per timestep, sec

DPMTA

DPME

Figure 10.8. CPU time for evaluating the electrostatic energy for 23,832 water molecules

(71,496 atoms) by codes developed at Duke by John Board and co-workers for the dis-

tributed PME (DPME) and distributed parallel multipole tree algorithm (DPMTA), run at

modest accuracy of 10

−3

relative force accuracy, on a tightly coupled network of 450 MHz

Pentium II processors.

See also a review on modeling electrostatic effects in proteins [1347], which

discusses stability and accuracy issues with traditional PME, or non-spherical

Ewald methods as described above. These problems can be alleviated by using

spherical boundary conditions in combination with a local ﬁeld approach. The

above review and others argue that only spherical Ewald methods rigorously give

correct results for charged systems.

10.5 The Multipole Method

In this section, we present a brief introduction to the efﬁcient fast multipole

technique. For details, see [479,481] and other references cited below.

10.5.1 Basic Hierarchical Strategy

Fast multipole techniques are powerful alternatives to the Ewald schemes intro-

duced above for evaluating the pairwise interactions in large molecular systems.

They are widely used for many important problems in applied mathematics,

engineering, physics, chemistry, and biology [478,480]. Examples in astrophysics

include evaluation of gravitational potentials, and examples in chemistry in-

clude electrostatic potentials in molecular dynamics and quantum mechanics

(Hartree-Fock calculations).

10.5. The Multipole Method 325

Series Expansion

Multipole schemes rely on a power-series expansion that describes the interac-

tion between groups of particles (charged bodies in molecular dynamics) [1226].

The term multipole refers to the moments m



i=1

that appear in the

expansion, such as dipole and quadrupole. To see this, consider the electrostatic

potential Φ at x for unit dielectric written in terms of the charge distribution about

each atom j as (see also eqs. (10.24), (10.25) for the periodic version):

Φ(x)=Φ(x

, x

,...,x



{i,j},i=j

/|r

| =



Φ(x

) ,

(10.45)

where

Φ(x



i=j

/|r

|. (10.46)

An expansion of Φ(x

) yields (see Box 10.3):

Φ(x

)=Φ(0)+∇Φ

x +

Hx + O(|x|

) (10.47)

where the derivatives are evaluated at x =0. This expansion produces expressions

for Φ(0) based on the zeroth moment (Φ(0) = m

/|x

|), for ∇Φ(x) in terms of

the dipole moment (vector μ =[m

]), and for the second-derivative term

based on the quadrupole moment (elements m

, m

, ···,m

). See

Box 10.3 for details.

Note that multipole expansions can be expressed using Cartesian coordinations,

as in Boxes 10.3 and 10.4, or spherical coordinates, as in Subsection 10.5.3.Such

power series can save computational work since the target function can be written

as a linear combination of moments (see Box 10.4 for a simple example).

Domain Decomposition

Accelerated multipole algorithms use a hierarchy of approximations to repre-

sent these interactions on the basis of a spatial particle partitioning in a tree-like

structure (typically oct-tree): the original domain (level 0) is subdivided into 8

level-1 domains, leading to 64 regions in level 2, and so on (see also illustration

in [316]). Thus, each successive reﬁnement of the computational domain pro-

duces ‘offspring’ corresponding to the ‘parents’ at the prior level. This recursive

partitioning, like FFTs, works together with a systematic formulation and manip-

ulation of power series expansions of the target potential to reduce the evaluation

computational complexity from O(N

) to the more modest O(N log N ) and

O(N ), depending on the implementation (see Figure 10.9).

Multipole expansions represent one suitable choice of power series, since they

converge well when groups of particles (charges) are well-separated, allowing a

small number of coefﬁcients in the expansions while maintaining good accuracy;

326 10. Nonbonded Computations

other expansions are possible [45,316,317,460,1329]. These power series expan-

sions are used to compute interactions between well-separated pairs (contained

in well-separated clusters). Interactions between nearby particles are computed

directly.

Summation Protocol

The various algorithms differ by how the tree structure is generated and by the

protocol used to determine which power-series approximation to invoke at ev-

ery stage (on the basis of distance separations monitored by interaction lists).

Figure 10.9, for example, illustrates the deﬁnition of interaction lists for multi-

pole expansions by boxes beyond ﬁrst neighbors with respect to a given particle

(shown in red); starting from second neighbors is an alternative deﬁnition.

Generally, good efﬁciency — linear complexity for moderate-sized systems —

for 3D implementations requires meticulous programming and algorithmic

structure [481]. This programming sophistication is especially important for

applications where the particles are distributed heterogeneously through the com-

putational domain (usually not the case for molecular dynamics applications). For

heterogeneous cases, adaptive schemes are needed to distribute the computations

in a balanced and efﬁcient manner since regular subdivisions may generate some

empty cells and empty cells need not be divided.

level 0: 1 box level 1: 4 boxes level 2: 16 boxes

level 3: 64 boxes level 4: 256 boxes

no well separated boxes no well separated boxes

(only near neighbors)

well separated boxes (yellow)

contribute to multipole

interaction list

previous-level multipole interactions (gray) are excluded

from multipole interaction list (yellow)

near neighbors

multipole interactions

previous interactions

Figure 10.9. Hierarchical domain partitioning approach for fast multipole schemes. The

partitioning is illustrated with respect to a central particle (shown in red).

10.5. The Multipole Method 327

The seasoned programming required for fast multiple algorithms partly ex-

plains why Ewald codes have been more popular in the molecular dynamics

community over the last few years. It is possible that as larger systems (of order

atoms) become more common in macromolecular simulations, the differ-

ence between the O(N log N) (best Ewald implementation) and O(N) (best

multipole implementation) may become signiﬁcant, and the effort required for

multipole codes may well be worth the programming investment.

Besides Coulomb interactions, multipole (and Ewald) methods apply more gen-

erally to any function that can be approximated by a converging power-series

expansion, such as functions of 1/r

or exponentials as in screened Coulomb

(Debye-H¨uckel) expressions (exp(−κr)/r)[154,385](seebelow).

Box 10.3: Multipole Expansion

Consider a charge distribution about each atom j (in a system of N atoms) as deﬁned in

eq. (10.46), where r

= x

− x

=[x

− x

], r

= |r

|,and

=(r

)

. An expansion at position x about the origin can be written as

Φ(x



i=1

[(x

− x

)

+(x

− x

)

+(x

− x

)

]

1/2

=Φ(0)+



i=1



k=1

∂Φ

∂x



i=1



l=1



k=1

∂

∂x

higher

order

terms

≡ Φ(0) + ∇Φ

x +

Hx + O(|x|

)

where derivatives are evaluated at x

= x

=0. The gradient vector ∇Φ

and Hessian matrix H above contain partial derivatives corresponding to the 3N com-

ponents of the vector x. Namely, ∇Φ has the 3N components {∂Φ/∂x

, ∂Φ/∂x

, ···,

∂Φ/∂x

}, and the Hessian is:

H =

⎛

⎜

⎝

∂

Φ/∂x

∂

Φ/∂x

∂x

...... ∂

Φ/∂x

∂x

∂

Φ/∂x

∂x

∂

Φ/∂x

...... ∂

Φ/∂x

∂x

∂

Φ/∂x

∂x

∂

Φ/∂x

∂x

...... ∂

Φ/∂x

⎞

⎟

⎠

The moments {m

} for integers k corresponding to the partial charges {q

} are deﬁned

by m



i=1

where the summation extends over all components of x.Theﬁrst

term in the above expansion is written in terms of the zeroth moment m

Φ(0) =



i=1

+ x

]

1/2



i=1

328 10. Nonbonded Computations

The gradient, or the dipole potential, has 3 components: m

. They represent the

dipole moment μ. The expression for the gradient in terms of the dipole moments can be

derived following the relation for the derivatives:

∂R

−1/2

∂x

= −

−3/2

∂R

∂x

= −r

−3

− x

)



= x

/ |x

We then have:



)=∇Φ

x =



i=1

∂R

−1/2

∂x

∂R

−1/2

∂x

∂R

−1/2

∂x



i=1



+ q





+ m



(μ

· x

)

The second-derivative term, the quadrupole potential, has six terms (because of symme-

try), again with all partial derivatives evaluated at x

= x

=0:



Hx =



i=1

∂

−1/2

∂x

∂

−1/2

∂x

∂

−1/2

∂x

∂

−1/2

∂x

∂

−1/2

∂x

∂

−1/2

∂x

∂

−1/2

∂x

+ m

∂

−1/2

∂x

+ ... +2m

∂

−1/2

∂x

Box 10.4: Example of Work Reduction by Multipole Expansion

Let our potential be deﬁned as:

Φ(x)=



j=1

Φ(x



j=1



i=1

− y

)

That is, Φ(x) is the sum of the following components:

{Φ(x

)} q

− y

)

+ q

− y

)

+ ... + q

− y

)

{Φ(x

)} q

− y

)

+ q

− y

)

+ ... + q

− y

)

{Φ(x

)} q

− y

)

+ q

− y

)

+ ... + q

− y

)

Evaluation by straightforward summation requires O(N

) work, but our function of x

and y simpliﬁes as:

(x − y)

= x

− 3x

y +3xy

− y

10.5. The Multipole Method 329

This is a ﬁnite power series, with degree p =3and moments m



i=1

. Therefore

we can rewrite each Φ(x

) as a linear combination of the moments {m

} for k =0,1,2,

Φ(x

)=x

− 3x

+3x

− m

The work can thus be substantially reduced: once the moments are computed in O(Np)

work, evaluating each Φ(x

) requires only O(p) work, where p is the power of the

expansion, 3 here. This reduction O(N

) to O(3N) is signiﬁcant!

10.5.2 Historical Perspective

Hierarchical Reﬁnements

Coming from the astrophysics community, Appel ﬁrst introduced the multi-

pole approach for solving the N-body problem by a hierarchical, power-series

approach [53]. Barnes and Hut accelerated the association scheme between par-

ticles at successive reﬁnement levels to produce a method that is asymptotically

O(N log N) [90]. To see this let D

N×N

be the matrix deﬁned by



/|r

| i = j

0 i = j

(10.48)

and q be the vector of N partial charges. Hence, the potential for each atom j due

to the charges induced by all other atoms is

Φ(x

)=q



i=j

/|r

| = {Dq}

the jth component of the matrix/vector product Dq. Summing up the values of the

components of the resultant product (N potentials evaluated at N points) yields

the desired potential. This is clearly O(N

) computation by a straightforward

matrix/vector multiplication. However, if the potential Φ is sufﬁciently smooth

when distances are sufﬁciently large, the matrix D can be approximated by low-

rank submatrices; in this case, an O(Np) scheme results where p is the rank of the

approximating matrices. See also the simple example in Box 10.4 and Figure 10.9.

Hierarchical Protocol

This O(N log N) method works roughly as follows.

A hierarchy of boxes is introduced in the computational domain so that each

reﬁnement level l is a subdivision (e.g., into 8) of the domain in level l −1.Boxes

at the same reﬁnement level that share a boundary are considered near neighbors,

and those at the same reﬁnement level that do not share a boundary are considered

well-separated.

As the algorithm sweeps through reﬁnement levels, a multipole expansion is

associated with each box at that level to describe the far-ﬁeld potential contribu-

tion from the particles in that box. (The clusters that contribute to this far ﬁeld