The Ewald Summation method Invoking the Ewald summation for calculating the electrostatic interactions can be specified any time the nbond specification parser is invoked. See the syntax section for a list of all commands that invoke this parser. Prerequisite reading: (nbonds.doc) * Menu: * Syntax:: Syntax of the Ewald summation specification * Defaults:: Defaults used in the specification * Function:: Description of the options * Discussion:: More general discussion of the algorithm
[SYNTAX EWALD] { NBONds } { nonbond-spec } { UPDAte } { } { ENERgy } { } { MINImize } { } { DYNAmics } { } The keywords are: nonbond-spec::= [ method-spec ] method=spec::= { EWALd [ewald-spec] } { PMEWald [pmesh-spec] } { [ NOEWald ] } { [ NOPMewald ] } ewald-spec::= KAPPa real KMAX integer KSQMAX integer [erfc-spec] [ KMXX integer KMXY integer KMXZ integer ] pmesh-spec::= KAPPa real FFTX int FFTY int FFTZ int ORDEr integer [erfc-spec] erfc-spec::= { [ ABROmowitz ] } { CHEBychev } { SPLIne } { INTEpolate } { EXACt_high_precision } { LOWPrecision_exact }
The defaults for the ewald summation are set internally and are currently set to NOEWald, KAPPa=1.0, KMAX=5, KSQMax=27, and NOPMewald, KAPPa=1.0, FFTX=FFTY=FFTZ=32, ORDEr=4
i) The EWALD keyword invokes the Ewald summation for calculation of electrostatic interactions in periodic, neutral systems. The formulation of the Ewald summation dictates that the primary system must be neutral. If otherwise, the summation is not formally correct and some convergence problems may result. The NOEWald (default) suppresses the Ewald method for calculating electrostatic interactions. Van der waals options VSHIFT and VSWITCH are supported with ewald. The algorithm currently supports the atom nonbond list and the CRYSTAL facilty must be used. The PMEWald keyword invokes the Particle Mesh Ewald algorithm for the reciprocal space summation. For details on the PME method, see J. Chem. Phys. 103:8577 (1995). The EWALd algorithm is limited to CUBIC, TETRAGONAL, and ORTHORHOMBIC unit cells. The PMEWald algorithm supports all unit cells that are supported by the CRYSTAL facility. ii) The KAPPa keyword, followed by a real number governs the width of the Gaussian distribution central to the Ewald method. An approximate value of kappa can be chosen by taking KAPPa=5/(CTOFNB*2). See discussion section for detail on choosing an optimum value of KAPPa. iii) The KMAX key word is the number of kvectors (or images of the primary unit cell) that will be summed in any direction. It is the radius of the Ewald summation. For orthorombic cells, the value of kmax may be independently specified in the x, y, and z directions with the keywords KMXX, KMXY, and KMXZ. In the PME version, the number of FFT grid points for the charge mesh is specified by FFTX, FFTY, and FFTZ. iv) The KSQMax key word should be chosen between KMAX squared and 3 times KMAX squared. v) An appropriate, although not optimal, set of parameters can be chosen by taking KAPPA=5/CTOFNB and KMAX=KAPPa*boxlength. The actual values should then be performanced optimized for your particular system. For the PME method, FFTX should be approximately the box length in Angstroms. (for efficiency, FFTX should be a multiple of powers of 2,3, and 5) ORDEr specifies the order of the B-spline interpolation, e.g. cubic is order 4 (default), fifth degree is ORDEr 6. The ORDEr must be an even number and at least 4. vi) EWALd runs in parallel on both shared (PARVECT) and distributed memory parallel computers. PMEWald is currently not parallel. vii) several algorithms are available for the calculation of the complimentary error function, erfc(x). EXACt and LOWPrecision use an interative technique described in section 6.2 of Numerical Recipies. ABRO and CHEB are polynomial approximations. A lookup table (filled at the beginning of the simulation using the EXACt method) can be used with either a linear (INTE) of cubic spline (SPLINe) interpolation. SPLIne is recommended.
The Ewald Summation in Molecular Dynamics Simulation The electrostatic energy of a periodic system can be expressed by a lattice sum over all pair interactions and over all lattice vectors excluding the i=j term in the primary box. Summations carried out in this simple way have been shown to be conditionally convergent. The method developed by Ewald, in essence, mathematically transforms this fairly straightforward summation to two more complicated but rapidly convergent sums. One summation is carried out in reciporcal space while the other is carried out in real space. Based on the formulation by Ewald, the simple lattice sum can be reformulated to give absolutely convergent summations which define the principal value of the electrostatic potential, called the intrinsic potential. Given the periodicity present in both crystal calculations and in dynamics simulations using periodic boundary conditions, the Ewald formulation becomes well suited for the calculation of the electrostatic energy and force. If we consider a system of point charges in the unit or primary cell, we can specify its charge density by ro(r) = sum_i [ q_i * delta(r-r_i)] In the Ewald method this distribution is replaced by two other distributions ro_1(r) = sum_i [ q_i ( delta(r-r_i) - f(r-r_i)] and ro_2(r) = sum_i [q_i f(r-r_i) such that the sum of the two recovers the original. The distribution, f(r), is a spherical distribution generally taken to be Gaussian, the width of the gaussian dictated by the parameter, KAPPa. The charge distributions are situated on the ion lattice positions, but integrate to zero. The potential from the distribution ro_1(r) is a short range potential evaluated in a direct real space summation (truncated at CTOFNB). The diffuse charge distribution placed on the lattice sites reduces to the potential of the corresponding point charge at large r. ro_2(r), being a continuous distribution of Gaussians situated on the periodic lattice positions, is a smoothly varying function of r and thus is well approximated by a superposition of continuous functions. This distribution is, therefore, expanded in a Fourier series and the potential is obtained by solving the Poisson equation. The point of splitting the problem into two parts, is that by a suitable choice of the parameter KAPPa we can get very good convergence of both parts of the summation. For the real space part of the energy, we choose kappa so that the complementary error function term, erfc(kappa*r) decreases rapidly enough with r to make it a good approximation to take only nearest images in the sum and neglect the value for which r > CTOFNB. The reciprocal space sums are rapidly convergent and a spherical cutoff in k space is applied so that the sum over k becomes a sum over {l,m,n}, with (l**2+m**2+n**2) < or = to KSQMAX A large value of KAPPa means that the real space sum is more rapidly convergent but the reciprocal space sum is less rapid. In practice one chooses KAPPa to give good convergence at the cutoff radius, CTOFNB. KMAX is then chosen to such that the reciprocal space calculation converges. The equation (KMAX/(box length)=KAPPa may be used as a rough guide. Optimization with respect to the timing trade offs, ie. how much time is spent in real space vs k-space should be performed before a lengthy production run. The CCP5 notes in several articles in 1993 cover some possible optimization strategies and criteria although a simple line search will suffice. Complete optimiztion of the ewald method for a particular application requires optimizing CTOFNB, KAPPa, and KMAX. A discussion of optimization and error analysis can be found in Kolfka and Perram, Molecular Simulation, 9, 351 (1992).
NIH/DCRT/Laboratory for Structural Biology
FDA/CBER/OVRR Biophysics Laboratory