This preview shows page 1. Sign up to view the full content.
Unformatted text preview: Nonbonded Interactions Overview Computation of the NonBonded Potential Recall that the nonbonded contribution to the potential function takes the form qi qj Aij Bij Vnb (R) =  6 + 12 + . rij rij rij
i<j i<j Since each of these sums contains O(N 2 ) terms, direct evaluation of the nonbonded potential is impractical in systems containing thousands or more atoms. In such cases, it is usually necessary to estimate Vnb using faster, approximate methods. Jay Taylor (ASU) APM 530  Lecture 5 Fall 2010 1 / 38 Nonbonded Interactions Spherical Cutoffs Cutoff Techniques Cutoff techniques reduce the cost of the nonbonded computations to order O(N) by ignoring nonbonded interactions between nuclei separated by more than some maximum distance, b. Typically, the interaction function is altered multiplicatively: f (r ) is replaced by f (r )S(r ) where S(r ) is zero for r > b. The cutoffs can be applied using distances between atoms or between group centers. Groupbased methods treat atoms that belong to the same residue or base similarly. Cutoffs can be applied directly to the energy and to the force field. Jay Taylor (ASU) APM 530  Lecture 5 Fall 2010 2 / 38 Nonbonded Interactions Spherical Cutoffs Truncation Truncation methods abruptly set the interaction term equal to 0 when r b: 1 if r < b S(r ) = 0 if r b. Advantage: easy to implement. Disadvantage: because the truncated potential function is discontinuous, optimization and MD routines run poorly. Jay Taylor (ASU) APM 530  Lecture 5 Fall 2010 3 / 38 Nonbonded Interactions Spherical Cutoffs Switching Switch functions smoothly reduce the energy to 0 One example is if 1 1 + y (r )2 (2y (r )  3)) if S(r ) = 0 if where y (r ) = (r 2  a2 )/(b 2  a2 ). over an interval [a, b]. r <a ar b r b, S(r ) as defined above is continuously differentiable, so the corresponding force field is continuous. If a is close to b, then the derivative in [a, b] will be large. Jay Taylor (ASU) APM 530  Lecture 5 Fall 2010 4 / 38 Nonbonded Interactions Spherical Cutoffs Shift Functions Shift functions gradually reduce the energy over the entire interval [0, b]. One formulation used in CHARMM is: S(r ) = [1  r /b]2 0 if r < b if r b. There is a tradeoff between smoothness and underestimation of the potential: Switch functions have larger derivatives but only underestimate the energy over the region [a, b]. Shift functions underestimate the energy at all distances, but vary more smoothly. Jay Taylor (ASU) APM 530  Lecture 5 Fall 2010 5 / 38 Nonbonded Interactions Spherical Cutoffs Bookkeeping Calculating the distance rij for every pair of atoms in the molecule requires O(N 2 ) steps. To implement cutoff techniques in order O(N) calculations, each atom is assigned a pairlist of atoms that are within some distance b + c at the start of the calculation. c should be taken large enough that only those atoms contained within the pairlist are likely to come within distance b of the associated nucleus over the course of the calculation. Jay Taylor (ASU) APM 530  Lecture 5 Fall 2010 6 / 38 Nonbonded Interactions Spherical Cutoffs Cutoff Techniques in Practice Truncation is almost never used because of the problems caused by discontinuities. Shift functions are used with van der Waals forces, which decay rapidly over long distances (like r 6 ). Because the Coulomb potential decays slowly (like r 1 ), longrange electrostatic interactions are especially important in charged molecules like DNA.
Many studies have reported that structural modeling using spherical cutoffs of the electrostatic potential gives poor results. Groupbased methods perform especially poorly. Accurate results may be attainable using atombased shifts with cutoffs of 12 (Norberg & Nilsson, 2000). A Jay Taylor (ASU) APM 530  Lecture 5 Fall 2010 7 / 38 Nonbonded Interactions Ewald summation Overview Ewald summation originated in crystallography as a method for calculating slowly convergent electrostatic potentials. We consider a group of point charges in an infinite periodic lattice. Each point charge is screened by a smooth charge distribution of the opposite sign. This leads to a rapidly converging potential (direct sum). The screening charges are compensated by a smooth distribution that leads to slowly converging potential (indirect sum). Because the compensatory distribution is smooth, it can be converted to a rapidly converging Fourier series (Poisson summation). Jay Taylor (ASU) APM 530  Lecture 5 Fall 2010 8 / 38 Nonbonded Interactions Ewald summation Periodic Extension of the Charge Distribution To apply Ewald summation to the electrostatic potential encountered in MM, we imagine that the molecule sits in a cubic cell of side length L that is used to form an infinite periodic lattice. We then seek to calculate
N 1 qi qj Vel = 2 n xij + n i,j=1 where n = L(n1 , n2 , n3 ) denotes a cellcoordinate vector; The over the first sum means that terms with i = j are omitted when n = (0, 0, 0). Ewald summation can also be carried out over other periodic lattices. Jay Taylor (ASU) APM 530  Lecture 5 Fall 2010 9 / 38 Nonbonded Interactions Ewald summation The Screening and Canceling Distributions The charge distribution is decomposed as i qi xi = where is a probability distribution with a smooth density and + xi denotes the translate of to xi .
Jay Taylor (ASU) APM 530  Lecture 5 Fall 2010 10 / 38 qi i  ( + xi ) + qi ( + xi )
i i Nonbonded Interactions Ewald summation The Screening Distribution The screening distribution is usually taken to be a radiallysymmetric threedimensional Gaussian distribution with density: g (r) 3 e 
2r 2 . As increase, converges weakly to a point mass at 0. g (r) has smooth, rapidlydecaying tails. Jay Taylor (ASU) APM 530  Lecture 5 Fall 2010 11 / 38 Nonbonded Interactions Ewald summation Poisson's Equation Given a charge distribution (r), the electrostatic potential (r) can be found by solving Poisson's equation 2 (r) = 4(r). Recall that the solution can be expressed as a convolution (r ) (r) = dr , 3 r  r  R where G (r, r ) =
1 rr  is the fundamental solution in R3 : 2 G (r, r ) = 4(r  r ). Jay Taylor (ASU) APM 530  Lecture 5 Fall 2010 12 / 38 Nonbonded Interactions Ewald summation The Potential of a Gaussian Charge The potential of a Gaussian screening charge located at zero is g (r ) g (s + r) Gauss (r) = dr = ds  s R3 r  r R3 exp  2 s + r2 3 = ds s R3 The last integral can be evaluated by noting that since g (r) is radially symmetric, so is Gauss (r), and therefore we can take r = (0, 0, r ) and write s + r2 = s 2 + 2rs cos() + r 2 , where is the angle between s and the zaxis. Jay Taylor (ASU) APM 530  Lecture 5 Fall 2010 13 / 38 Nonbonded Interactions Ewald summation Transforming to spherical coordinates (s, , ) gives Gauss (r) = 2 3 2 1 2 2 2 = s ds sin()d d e  (s +2rs cos()+r ) s 0 0 0 3 2 2 2 2 2 = 2 e  r se  s ds sin()e 2 rs cos() d 0 0 3 1 2 2 2 2 2 = 2 e  r se  s ds e 2 rsu du 1 0 3 1  2 r 2 2 2 2 2 = e e  s e 2 rs  e 2 rs ds 2 r 0 Jay Taylor (ASU) APM 530  Lecture 5 Fall 2010 14 / 38 Nonbonded Interactions Ewald summation = = = = 1  2 r 2  2 (s 2 +2rs) 2 2 e e  e  (s +2rs) ds r 0 1  2 (s+r )2 2 2 e  e  (sr ) ds r 0 1  2 s 2  2 s 2 e ds  e ds r r r 2 1 r  2 s 2 2 1 r t 2 e ds = e dt. r 0 r 0 Thus, the potential induced by a Gaussian charge distribution is just 1 Gauss (r) = erf (r ), r where erf (x) denotes the error function.
Jay Taylor (ASU) APM 530  Lecture 5 Fall 2010 15 / 38 Nonbonded Interactions Ewald summation The Potential of a Single Screened Charge Since the potential induced at r by a point charge at 0 is 1 , it follows that r the combined potential of a point charge and its screening Gaussian is dir (r) = 1 1 1  erf (r ) erfc(r ) r r r 2 2 where erfc(x) = x e t dt is the complementary error function. In particular, this function decays exponentially as r : dir (r) < r 2 1 2 2 e  r . Jay Taylor (ASU) APM 530  Lecture 5 Fall 2010 16 / 38 Nonbonded Interactions Ewald summation The Direct Potential Energy To compute the total potential energy of the point charges in the presence of the screened point charges, we calculate Vdir = 1 qi qj dir (n + rj  ri ) 2 n
i=1 N N j=1 = 1 2 n i,j=1 N qi qj erfc(rij + n) rij + n . Because erfc(r ) decays exponentially, this infinite sum converges rapidly. Furthermore, by choosing sufficiently large, we can compute the sum to any degree of accuracy using only order O(N) terms. Jay Taylor (ASU) APM 530  Lecture 5 Fall 2010 17 / 38 Nonbonded Interactions Ewald summation The Compensatory Potential Since the compensatory charge distribution consists of a superposition of Gaussian distributions, each of which has potential erf (r )/r , the total potential energy of the point charges caused by the canceling charges is 1 qi qj gauss (ri  (rj + n)) 2 n
i=1 N N Vcomp = j=1 = N 1 qi qj erf (rij + n) . 2 n rij + n i,j=1 However, since erf (r ) tends to 1 as r , this sum is slowly converging. Jay Taylor (ASU) APM 530  Lecture 5 Fall 2010 18 / 38 Nonbonded Interactions Ewald summation Poisson Summation The key step in the Ewald method is that Vcomp can be converted to a rapidly converging sum by using the Poisson summation formula: F (x) where n f (x + n) = 1 ^ f (k)e 2ixk V k F is the periodic summation of a function f : Rd R. is a lattice in Rd with cell volume V and dual lattice . ^ f is the Fourier transform of f : ^(k) = f f (x)e 2ikx dx.
Rd
Jay Taylor (ASU) APM 530  Lecture 5 Fall 2010 19 / 38 Nonbonded Interactions Ewald summation Poisson Summation of the Compensatory Potential To evaluate the compensatory potential, we take f (r) = N N erf (rij + r) 1 1 qi qj qi qj (rij + r) 2 rij + r 2
i,j=1 i,j=1 and set Vrecip = n f (n) = 1 ^ f (k). V k Note: To use Poisson summation, we need to remove the restriction in Vcomp on the sum over i, j when n = 0. We will correct for this later. Jay Taylor (ASU) APM 530  Lecture 5 Fall 2010 20 / 38 Nonbonded Interactions Ewald summation ^ Evaluating f Notice that the Fourier transform of f can be expressed as ^ f (k) = where ^ (k) = 1 erf (r)e 2ikr dr. r
N ^ qi qj (k)e 2k(ri rj ) , i,j=1 R3 Jay Taylor (ASU) APM 530  Lecture 5 Fall 2010 21 / 38 Nonbonded Interactions Ewald summation ^ Although can be calculated directly from the integral, Poisson's equation provides an easier route: 2 (r) = 4g (r). Expanding and g in Fourier series and substituting into Poisson's equation gives: 2 2ikr ^ ^ (k)e = 4 2 k2 (k)e 2ikr
k = 4 k k g (k)e 2ikr ^ which implies that ^ (k) = 1 g (k). ^ k2
Fall 2010 22 / 38 Jay Taylor (ASU) APM 530  Lecture 5 Nonbonded Interactions Ewald summation Since the Fourier transform of g is g (k) = ^ 3 e 
2 r2 R3 e 2ikr dr = e  2 k2 / 2 , the equality on the previous slide gives ^ (k) = and then
2 2 2 N  2 k2 / 2 e  k / 2k(ri rj ) ^(k) = e f qi qj e S(k)2 . k2 2k2 1 2 k2 / 2 e k2 i,j=1 Jay Taylor (ASU) APM 530  Lecture 5 Fall 2010 23 / 38 Nonbonded Interactions Ewald summation The Structure Factor In crystallography, the function S is known as the chargeweighted structure factor: S(k) =
N j=1 qj e 2krj . S is the Fourier transform of the discrete charge density within the cell located at the origin. It takes order O(N) computations to evaluate S(k) at each wavelength k. Jay Taylor (ASU) APM 530  Lecture 5 Fall 2010 24 / 38 Nonbonded Interactions Ewald summation A Rapidly Converging Series for the Reciprocal Sum Substituting the Fourier transform of f back into the Poisson summation formula gives the series Vrecip
2 2 2 1 e  k / = S(k)2 , 2V k2 k which is rapidly converging since S() is bounded and thus the terms decay exponentially as the wave number k increases. The Fourier series converges rapidly because the function (z) = erf (z)/z is smooth (in fact, entire) when < . Jay Taylor (ASU) APM 530  Lecture 5 Fall 2010 25 / 38 Nonbonded Interactions Ewald summation The SelfInteraction Term Recall that the reciprocal sum includes some terms not present in Vcomp :
N 1 qi qj erf (rij + n) 2 n rij + n i,j=1 Vcomp = Vrecip = N 1 qi qj erf (rij + n) . 2 n rij + n i,j=1 Since the extra terms correspond to nonphysical selfinteractions between the point charges and the compensatory Gaussian charge distribution, they need to be subtracted when calculating the total energy. Jay Taylor (ASU) APM 530  Lecture 5 Fall 2010 26 / 38 Nonbonded Interactions Ewald summation The SelfInteraction Term Since 1 2 lim erf (r) = lim r0 r r0 r r 0 2 2 e t dt = , the correction for the selfinteraction energy is
N  2 Vself Vcomp  Vrecip = qi . i=1 Notice that this sum contains O(N) terms and does not depend on the locations of the nuclei. Jay Taylor (ASU) APM 530  Lecture 5 Fall 2010 27 / 38 Nonbonded Interactions Ewald summation The Ewald Representation Putting everything together, we obtain the following decomposition of the lattice electrostatic potential Vel = Vdir + Vrecip  Vself N erfc(rij + n) 1 = qi qj 2 n rij + n
i,j=1
2 2 2 1 e  k / + S(k)2 2V k2  k N i=1 qi2 , where both infinite sums are rapidly converging.
Jay Taylor (ASU) APM 530  Lecture 5 Fall 2010 28 / 38 Nonbonded Interactions Ewald summation FiniteDielectric Correction A fourth term must be added to the electrostatic energy whenever the dipole moment of a unit cell is nonzero: 2 N 2 = 3 qj rj . 3L j=1 Vcell The vector quantity inside the modulus is the dipole moment of a unit cell. This sum can be evaluated in order N steps. Jay Taylor (ASU) APM 530  Lecture 5 Fall 2010 29 / 38 Nonbonded Interactions Ewald summation Optimization of Recall that controls the variance of the Gaussian screening charges. As increases, the direct sum Vdir converges more rapidly and so fewer terms need to be summed to achieve a desired accuracy. On the other hand, increasing causes the reciprocal sum Vrecip to converge more slowly, requiring more terms to be summed. The selfinteraction term contains N terms and is independent of . These observations suggest that should be chosen neither too large nor too small to minimize the total work needed to estimate Vel to some specified degree of accuracy. Jay Taylor (ASU) APM 530  Lecture 5 Fall 2010 30 / 38 Nonbonded Interactions Ewald summation Optimization of Perram et al. (1988) gave the following heuristic for choosing : Let denote the acceptable relative error. If is chosen so that only charges within some distance b are included in the direct sum, then to satisfy the error condition we need:  ln()/b. If td is the time required to compute one term in the direct sum, then the total amount of time required for this sum will be roughly Td Cd td N (b 3 N/V ) C1 N 2 3 . Jay Taylor (ASU) APM 530  Lecture 5 Fall 2010 31 / 38 Nonbonded Interactions Ewald summation Optimization of Similarly, the time required to compute the indirect sum can be estimated. If K denotes the maximum magnitude of the wave numbers that are included in the sum, then to achieve a relative error of , we need K  ln(). If tr is the time required to compute one term in the reciprocal sum, then the total amount of time required for this sum will be roughly Tr Cr tr K 3 N C2 N 3 . Jay Taylor (ASU) APM 530  Lecture 5 Fall 2010 32 / 38 Nonbonded Interactions Ewald summation Optimization of The total time required to sum Vdir and Vrecip to the desired accuracy is then T = Tdir + Trec = C1 N 2 3 + C2 N 3 . To find the value of that minimizes T , we set opt (C1 N/C2 )1/6 and
dT d 3 = 0, which gives Tmin C4 N 2 N 1/2 + C5 NN 1/2 CN 3/2 . Thus, by choosing appropriately, Ewald summation can be completed in order O(N 3/2 ) time. Jay Taylor (ASU) APM 530  Lecture 5 Fall 2010 33 / 38 Nonbonded Interactions Particle Mesh Ewald Methods Fourier Transforms are a bottleneck in Ewald summation. The time required to calculate Vrecip is determined by two factors: Order O( 3 ) terms must be included in this sum to achieve acceptable levels of accuracy. Each evaluation of the structure function S(k) requires order O(N) time. If the Fast Fourier transform (FFT) could be used, then we could reduce this time to order O( 3 log( 3 )). The obstacle is that, apart from crystals, the charges are not located on a lattice within each cell. Jay Taylor (ASU) APM 530  Lecture 5 Fall 2010 34 / 38 Nonbonded Interactions Particle Mesh Ewald Methods Mesh Interpolation ParticleMesh Ewald (PME) methods address this problem by mesh interpolation: Interpolate the charge distribution onto a periodic lattice (particle mesh) containing K points. Use the FFT to calculate the structure function at all K reciprocal wave numbers. This can be done in order O(K log(K )) time. Of course, K must be taken sufficiently large that the interpolated charge distribution is a good approximation to the true distribution. Because the interpolation can be done locally, K O(N) suffices. Jay Taylor (ASU) APM 530  Lecture 5 Fall 2010 35 / 38 Nonbonded Interactions Particle Mesh Ewald Methods ParticleMesh Ewald Methods PME can be implemented in O(N log(N)) time: can be chosen large enough to reduce Vdir to an order O(N) calculation. Order O(N) points are then needed in the mesh. The discrete charge density can be interpolated to the mesh in O(N) time.
The first formulation of PME used Lagrange interpolation (Darden et al., 1993). Essmann et al. (1995) reformulated PME using Bspline interpolation, which leads to improved accuracy and smooth gradients. Jay Taylor (ASU) APM 530  Lecture 5 Fall 2010 36 / 38 Nonbonded Interactions Particle Mesh Ewald Methods Time Comparison The table shows the times required to complete one energy and one gradient calculation for a solvated BPTI molecule (14275 atoms) using different methods to calculate the electrostatic potential. Method 12 switch A PME All nonbonded Source: Schlick (2006). Time (s) 3.1 8.5 54.8 Jay Taylor (ASU) APM 530  Lecture 5 Fall 2010 37 / 38 Nonbonded Interactions References References Darden, T., York, D. and Pedersen, L. (1993) Particle mesh Ewald: An N log(N) method for Ewald sums in large systems. J. Chem. Phys. 98: 1008910092. Essmann, U. et al. (1995) A smooth particle mesh Ewald method. J. Chem. Phys. 103: 85778593. Frenkel, D. and Smit, B. (2002) Understanding Molecular Simulation: From Algorithms to Applications. Second Edition, Academic Press. Perram, J., Petersen, J. and De Leeuw, S. (1988) An algorithm for the simulation of condensed matter which grows as the 3/2 power of the number of particles. Mol. Phys. 65: 875893. Schlick, T. (2006) Molecular Modeling and Simulation. Springer. Jay Taylor (ASU) APM 530  Lecture 5 Fall 2010 38 / 38 ...
View
Full
Document
This note was uploaded on 03/11/2012 for the course APM 530 taught by Professor Staff during the Fall '10 term at ASU.
 Fall '10
 Staff

Click to edit the document details