Fast Electrostatic Force Computations

From macromolecular structures to aqueous biological systems, accurate computation of electrostatic and van der Waals interactions is one of the most difficult tasks in computer modeling. Simulations of peptides [131] and membranes [6] as well as of ions in aqueous solutions [10,12,45,119,126] have provided clear-cut evidence of artifacts due to the use of cutoffs. MD simulations involving explicit solvent molecules have usually been performed using periodic boundary conditions (PBC) to reduce surface tension effects [24]. Therefore, it is important to develop fast force computation algorithms to simulate the N-body system with PBC at a cost that is not much greater than that of cutoff scheme but with better accuracy.

Ewald summation [36,48,50,52,60,83,87,117,130,142,141] is such a method to compute the long-range Coulomb forces with PBC. The Coulomb potential energy function per cubic cell under PBC is a conditionally convergent sum. Abel's theorem shows that a conditionally convergent sum can be turned into an absolutely and uniformly convergent sum by using a sequence of monotonic and uniformly bounded functions [90]. In the standard Ewald summation a Gaussian function is chosen as the convergence factor. With the optimal parameter the computational complexity of the standard Ewald sum is O(N^3/2). In Fourier-based Ewald sum [46,98,122,158], Fast Fourier Transforms (FFT) are used to compute the reciprocal sum over reciprocal vectors more efficiently. However, the Fourier-based Ewald does not parallelize well due to heavy interprocessor communication. Note that the choice of Gaussian function as the convergence factor is arbitrary. It was shown that the convergence characteristics cound be changed by varying the convergence factor [65,143]. The objective of this subproject is to find such a convergence factor in order to develop an Ewald sum that will be efficiently parallelized. We will use multilevel summation [13,26,27]: the force will be split in ``fast'' and ``smooth'' parts. The ``fast'' parts will be computed using cutoff computations in the particles. The ``smooth'' parts will be solved recursively in a hierarchy of grids. For periodic boundary conditions we will use two new methods: (i) Compute a direct multilevel sum of replicated periodic cells where the coarsest grid has zero charge (and thus there is nothing else to do); (ii) Compute the Ewald-sum using wavelet transforms [38,43,101,111] of compact support on a grid level, for which fast transforms exist or will be devised. Important computational tradeoffs of these methods will be evaluated, namely, larger size of (i) versus more complicated communication patterns of (ii).