CHARMM c28a3 ewald.doc

File: Ewald ]-[ Node: Top
Up: (commands.doc) -=- Next: Syntax -=- Previous: Top

                      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

File: Ewald ]-[ Node: Syntax
Up: Top -=- Next: Defaults -=- Previous: Top


{  NBONds          }        {  nonbond-spec                   }
{  UPDAte          }        {                                 }
{  ENERgy          }        {                                 }
{  MINImize        }        {                                 }
{  DYNAmics        }        {                                 }

The keywords are:

nonbond-spec::= [ method-spec ]

               { [ NOEWald ]                                          }
               {                                                      }
method-spec::= {   EWALd [ewald-spec] { [ NOPMewald [std-ew-spec] ] } }
               {                      {   PMEWald [pmesh-spec]      } } 

ewald-spec::=   KAPPa real  [erfc-spec]

std-ew-spec::= { [ KMAX integer ]                        } KSQMAX integer
               { KMXX integer KMXY integer KMXZ integer  }

pmesh-spec::=   FFTX int FFTY int FFTZ int  ORDEr integer [QCOR real (***) ]

erfc-spec::=    { SPLIne      { [EWMIn real] [EWMAx real] [EWNPts int] } }
                { INTErpolate {                                        } }
                {                                                        }
                { ABROmowitz                                             }
                { CHEBychev                                              }
                { EXACt_high_precision                                   }
                { LOWPrecision_exact                                     }
                { ERFMode int                                            }

File: Ewald ]-[ Node: Defaults
Up: Top -=- Next: Function -=- Previous: Syntax

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, QCOR=1.0

Recommended values for Ewald are:
   EWALD PMEWald KAPPa 0.34 ORDEr 6 - 
   FFTX intboxvx FFTY intboxvy FFTZ intboxvz  -
   CTOFNB 12.0 CUTNB 14.0 QCOR 1.0(***)

Where intboxv* is an integer value similar to or larger than the corresponding
unit cell dimension that has prime factors of 2,3, and 5 only (2,3 preferred).
grid point spacing should be between 0.8 and 1.2 Angstroms.

These recommended values should give relative force errors of roughly 10**-5.
To reduce the total PME cost at the expense of accuracy, decrease the cutoff
distances while increasing KAPPa (keep the product near 4) reduces the real
space cost.  To reduce the K-space cost, either reduce ORDEr from 6 to 4 or
increase the grid spacing up to perhaps 1.5 Angstroms.

(***) The QCOR value should be 1.0 for vacuum, solid, or finite systems.
For periodic systems in solution, it should be reduced (or set to zero) by an
amount that depends on how the net charge is distributed and on the effective
dielectric constant.  For a treatise on this correction term, see:
S. Bogusz, T. Cheatham, and B. Brooks, JCP (1998) 108, 7070-7084 and references
contained therein (esp. Hummer and Levy).

File: Nbonds ]-[ Node: Function
Up: Top -=- Previous: Defaults -=- Next: Discussion

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 and group nonbond lists 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.  This is fairly conservative.
Values of 4/CTOFNB lead to small force errors (roughly 10**-5).  See
discussion section for details 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.  PME runs in parallel on distributed
memory computers.

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.

viii) Ewald with MMFF

A version of EWALD was developed for MMFF.  The usual MMFF electrostatic
term: qq/(r+d)  is split into two terms:  qq/r -  qq*d/(r*(r+d))
The first term is handled by the Ewald method in the usual manner
(real-space and k-space parts) and the second term is truncated
at the cutoff distance using a switching function (from CTONNB to CTOFNB).
Since the second term is quite small at the cutoff distance, the use of a
switching function should not introduce significant artificial forces.

File: Ewald ]-[ Node: Discussion
Up: Top -=- Previous: Function -=- Next: Top

         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)]


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).   For PME, see Feller, Pastor,
Rojnuckarin, Bogusz, and Brooks. J. Phys. Chem., 100, 42, 17011 (1996)
and some of Tom Darden's published work.

CHARMM .doc Homepage

Information and HTML Formatting Courtesy of:

NIH/DCRT/Laboratory for Structural Biology
FDA/CBER/OVRR Biophysics Laboratory
Modified, updated and generalized by C.L. Brooks, III
The Scripps Research Institute