significant long - range 
 electrostatic interactions arise in many 
 biomolecular systems , such as negatively charged dna and rna , polar 
 or charged membranes , ion channels , and electrostatic steering of 
 protein  protein and enzyme  substrate association . accordingly
the computational cost increases in principle 
 as n , where n is the 
 number of charged particles in the system .
the evaluation of 
 electrostatic interactions can , through the use 
 of the fast fourier transform ( fft ) , be approximated with controlled 
 accuracy by two finite range calculations , one in real space and one 
 in reciprocal space . however , these fft - based approximation methods , 
 including particle  mesh ewald ( pme ) and 
 particle  particle particle  mesh ( pppm ) , generally require simulations to describe infinite three - dimensional 
 lattices where each lattice cell is filled with a copy of the simulated 
 system . as a result ,
to prevent artifacts due to copy  copy interactions , 
 the biologically relevant components of each copy must be surrounded 
 by ample solvent to guarantee enough spacing between copies , unfortunately 
 increasing system size and , thereby , computational cost . moreover , 
 the communication cost of calculating in parallel the two 3-d ffts 
 required by pme and pppm outpaces the computational cost as the number 
 of atoms and number of processors increase , due to the necessity of 
 exchanging data between all of the processors involved in the fft 
 calculation . on massively parallel computers routinely employed today 
 for large system molecular dynamics ( md ) simulations ,
the fft communication 
 becomes the main performance bottleneck , and efforts to maintain scalability 
 for pme have inspired various strategies , such as reserving a subset 
 of processors for the fft calculation .
alternative methods for calculating electrostatics , including the 
 real - space convolution method implemented in the special - purpose anton 
 2 supercomputer , have been developed to 
 avoid the fft communication bottleneck .
in particular , the multilevel 
 summation method ( msm ) effectively replaces the nonscalable fft communication 
 by more densely localized convolution calculations that permit better 
 parallel scaling and utilization of modern vector computational hardware 
 units .
an example where msm is methodologically more suitable 
 than pme 
 arises in studies of the asymmetric environment across a membrane 
 bilayer , where one might want to employ periodic boundaries along 
 the surface but not along the membrane normal .
if pme is employed 
 and fully periodic boundary conditions are required , one needs to 
 use a dual - membrane  dual - volume strategy , where two spatially 
 separated membranes are included and the 
 computational cost is doubled .
interest in supporting systems 
 having planar geometry has led to 
 the development of a 2-d ewald summation method , although this method has not yet been implemented in any 
 mainstream md software .
efforts have also been made to circumvent 
 the enforced periodicity from pme by creating a nonperiodic 
energy 
 step  at the edge of a periodic simulation cell to maintain 
 asymmetric ion concentrations across a membrane , but this approach can not maintain a charge gradient across 
 a membrane .
non - ewald methods , such as 
 the wolf method and the isotropic periodic 
 sum method with discrete fast fourier transform ( 3d - ips / dfft ) ( the latter available in the amber md package ) can also be applied to nonfully periodic systems , 
 but good accuracy for heterogeneous systems requires a large long - range 
 cutoff value .
msm has the advantage of 
 providing a unified methodology that can treat periodic , semiperiodic , 
 and nonperiodic boundary conditions with a single algorithm that , 
 furthermore , offers good parallel scalability .
msm , employing 
 so - called nested interpolation of softened pair 
 potentials in real space , follows a strategy for long - range force 
 calculation that differs from the one employed by pme .
msm was initially 
 introduced for solving integral equations and later applied to long - range electrostatic interactions in 2-d .
the method was then extended to calculate continuous 
 forces in 3-d suitable for md simulation in case of nonperiodic boundary 
 conditions , intended originally to replace 
 the functionality of the fast multipole method ( fmm ) .
the use of multiple spatial scales in msm makes 
 it better suited than pme to multiple time stepping that is employed in the time integration of md simulations .
the difficulties with the parallel scaling of pme and other fft - based 
 methods for ever greater numbers of processors have generated renewed 
 interest in fmm for improving the scalability of md , despite known stability problems caused by fmm producing discontinuous 
 potentials and forces that require the use of computationally costly 
 high order approximations .
the calculation 
 of highly accurate electrostatic forces provides no discernible benefit 
 to md , which instead requires continuous forces for the energy conserving 
 integration methods used for simulating long time scales .
there is 
 also no added benefit from using adaptive solvers , such as fmm or 
 other related oct - tree methods , due to the general uniformity of the 
 particle density arising in applications of md simulation to biomolecular 
 systems .
msm has a hierarchical structure similar to that of fmm , 
 allowing good parallel scalability while also producing continuous 
 forces that provide stable dynamics without the high computational 
 cost .
a gpu - accelerated implementation 
 of msm has been developed already to calculate maps of the coulombic 
 potential .
the methodology of msm can also be applied to other 
 pairwise potentials without truncation , most significantly to dispersion forces ( the long - range part of van der waals forces ) , which have critical 
 effects on membrane properties . in the present 
 study , we compare msm with pme , both as implemented 
 in version 2.10 of the md program namd , for accuracy 
 and efficiency . the parallelization of msm in namd is distinguished 
 from other recent parallelization efforts by making use of a combined domain decomposition and force decomposition 
 approach to provide scaling to large numbers of processors , together 
 with fine - grained simd ( single instruction , multiple data ) parallelism 
 employed to greatly enhance the performance of the local grid calculations . 
 the practical accuracy produced by msm , although shown to be less 
 than that typical for pme , is sufficient to reproduce the pme calculations 
 of various structural , dynamic , dielectric , and mechanical properties 
 of water . in particular , excellent agreement between msm and pme is 
 found for interface potentials at air  water and membrane  water 
 interfaces , for both periodic and semiperiodic simulations , indicating 
 that long - range coulomb interactions are well represented with msm 
 regardless of boundary condition .
applications demonstrate also the 
 suitability of msm for systems with semiperiodic and nonperiodic boundaries . 

semiperiodic simulations have been performed for membrane electroporation 
 induced by an imbalance of charge across a membrane and also for ion 
 conduction through a graphene nanopore with an ion gradient .
standard 
 use of pme for either system would require the aforementioned dual - membrane  dual - volume 
 strategy , doubling the system size .
nonperiodic simulations have also 
 been performed for two well - studied proteins , each inside a spherical 
 water droplet , to further validate the use of msm .
generally speaking , 
 the additional modeling options offered by msm to representative systems 
 makes msm more efficient than pme by reducing the number of atoms 
 needed for simulation .
parallel scalability results are also presented 
 that show the ability of msm to outperform pme when scaling a system 
 of modest size ( less than 100 k atoms ) to over a thousand processors , 
 demonstrating the suitability of msm for large - scale parallel simulation . 

we conclude that the improved modeling offered by msm , in combination 
 with sufficient accuracy , better parallel scaling than pme , and availability 
 in a mainstream md program , makes msm a compelling alternative to 
 pme .
for the 
 convenience of the reader , we provide below a brief overview 
 of the multilevel summation method ( msm ) algorithm . detailed discussions 
 including the error analysis employed are available elsewhere .
readers interested only in the use of msm , not in its algorithmic 
 underpinnings , may proceed to section 3 .
msm approximates the coulombic 
 potential energy for a system of n particles , having 
 position ri and charge qi,1where (i ) denotes 
 for each i the set of atom indices to be excluded 
 from the summation , which includes j = i and , typically for the simulation of biomolecules , also those atoms j that are either covalently bonded to i or to another atom that is covalently bonded to i. the interaction kernel k is split into the sum 
 of a short - range part k0 smoothly truncated 
 at cutoff distance a and slowly varying parts k1,k2, ... ,kl1,kl smoothly truncated at cutoff distances 
 2a , 4a, ...
,2a ,  , respectively , the 
 splitting can be defined systematically 
 in terms of a single unparameterized softening function (r ) that softens 1/r for r  1,2 the slowly varying parts are interpolated from grids of spacing h , 2h, ...
,2h , 2h , respectively , for which interpolation operators   are defined in terms of interactions between 
 grid point positions rml = ( xml , yml , zml ) and their 
 corresponding nodal basis functions ml,3where m and n index the grid points .
the nodal 
 basis functions can be defined 
 in terms of a single dimensionless basis function ,interpolation by piecewise polynomials of 
 degree p provides local support for  with 
 stencil size p + 1 .
continuous forces , deemed important 
 for energy conservation when simulating over long time scales , are produced from a continuously differentiable  and sufficient 
 continuity from . details regarding the softening and interpolation 
 basis functions are provided in section 2.2 .
the approximation is made efficient by nesting the interpolation 
 between levels,4 the 
 charmm force field prescribes a cutoff distance of a = 12  for calculating van der waals forces , which we adopt , for the sake of efficiency , 
 in calculating all short - range nonbonded interactions . the choice 
 of grid spacing h = 2.5  , which is slightly 
 larger than the interatomic spacing , works effectively in practice
. 
 the computational work per grid point is bounded by a constant , and 
 the number of grid points is reduced by about a factor of 1/8 at each 
 successive level .
the computational cost of msm is o((p + ( a / h))n ) , and analysis of the error shows an asymptotic bound of the formexperiments in the present study test c cubic interpolation ( p = 
 3 ) with c splitting , for which the functional forms of  and  , 
 given in section 2.2 , are among the lowest 
 order choices that produce continuous forces .
the coulombic 
 potential energy is approximated by substituting 
 eqs 4 and 2 into eq 1,5where the first summation is the exact short - range 
 part , the second summation is the approximate long - range part , and 
 the final summation removes the excluded interactions from the long - range 
 part . the eilong long - range electrostatic potential 
 is6 algorithmic steps for msm .
the short - range part of eq 5 simply 
 evaluates 
 all particle pair interactions within the cutoff distance a. the algorithm for eilong in eq 6 is made efficient by factoring the interpolation 
 operator in eq 3 , calculating at each grid point 
 intermediate charge qml and 
 potential eml .
the order of 
 the calculation is the following : anterpolation 
 7restriction 
 8grid cutoff calculation 
 9top level grid calculation 
 10prolongation 
 11interpolation 
 12figure 1 shows a diagram 
 depicting the computational steps and their dependencies .
the bottom  left 
 shows the particle positions and charges , and the calculated potentials 
 and forces are on the bottom 
the gridded charges are 
 calculated for each level moving up the left side , and the gridded 
 potentials are calculated moving down the right side .
the horizontal 
 arrows represent the calculation of the short - range part along the 
 bottom and the grid cutoff calculations for each grid level , all of 
 which collectively require the majority of the computational effort . 

forces are obtained analytically by taking the negative gradient of 
 eq 5,13 the interpolation 
 of atomic potentials from the finest grid in 
 eq 12 evaluates for each atom position the nearby 
 nodal basis functions requiring o(p ) operations per atom .
the anterpolation ( the adjoint 
 of interpolation ) in eq 7 performs similar work to spread atomic charges to the grid . 

all of the other functions ml and kl are evaluated at grid points , so can 
 be calculated a priori .
the restrictions in eq 8 and prolongations in eq 11 are gridded versions 
 of anterpolation and interpolation , each requiring just o(p ) operations per grid point when exploiting the 
 regularity of the nodal basis function stencil by factoring the sums 
 along each dimension . after a sufficient number of restrictions , the 
 top level grid calculation in eq 10
most of the 
 computational work for the long - range part is due to the grid cutoff 
 calculation in eq 9 , which evaluates pairs within 
 a cutoff distance of 2a / h grid points , 
 requiring o((a / h ) ) operations per grid point .
the partial potential eml , cutoff is the convolution 
 of grid charges qnl with the stencil kl ( rml , rnl ) = kl(0,rm  nl ) for |m 
an 
 alternative approximation scheme is available through hermite 
 interpolation , an approach that exactly reproduces function values 
 and derivatives at grid points , discussed in more detail in section 2.2 .
when one extends hermite interpolation to three 
 dimensions , there arise 2 = 8 values per grid point that 
 can account for all selections of zero or one derivatives independently 
 in each dimension .
the calculation of the grid point interactions 
 is then expressed as 8  8 matrix  vector products , which 
 have a straightforward mapping to cpu vector instructions , as discussed 
 in section 2.5 .
the increase in density of 
 the grids is offset by doubling the finest level grid spacing , so 
 that roughly the same amount of arithmetic operations are required 
 for hermite as for cubic interpolation .
the accuracy of hermite interpolation 
 with doubled grid spacing is shown to be between that of cubic and 
 quintic interpolation .
periodicity 
 along a dimension is accomplished by including the 
 periodic images of atoms contained within the cutoff distance into 
 the short - range part of eq 5 and wrapping around 
 the respective edges of the grid when calculating the long - range part . 

the top level is reduced , in this case , to a single grid point , for 
 which the charge will be zero , if the system of atoms is neutrally 
 charged , or will be set to zero , which effectively 
 acts as a neutralizing background potential .
the use of periodic boundary conditions imposes additional constraints 
 on the grid spacing and the number of grid points along each periodic 
 dimension .
the grid spacing must exactly divide the periodic cell 
 length , with the number of grid points chosen to be a power of 2 to 
 maintain the doubling of the cutoff distance and grid spacing at each 
 successive level . by permitting the number of grid points to also 
 have up to one factor of 3 ,
the grid spacings x , y , z can always 
 be chosen so that 2   x , y , z < 3  .
this strategy 
 also benefits simd implementation discussed in section 2.5 , where an odd number of grid points along a dimension 
 will need to be processed only once , namely on the top level .
the function 
 (r ) that softens 1/r for r  1 is defined in terms of a truncated 
 taylor series expansion of s about s = 1 and substituting r = s. for example , the softening used 
 in this study with cubic interpolation isgiving c continuity 
 at the unit sphere and c continuity 
 everywhere else .
additional advantages of an even - powered softening 
 are that ((1/a)(x + y + z ) ) has bounded derivatives and can be calculated without 
 square roots . a piecewise interpolating polynomial of odd degree p having c continuity , as needed 
 for producing continuous forces , can be constructed as a linear blending 
 of the two degree p  1 interpolating polynomials 
 centered on consecutive nodes . the dimensionless basis function ( ) 
 for cubic interpolation used in this study ishigher order interpolation 
 with quintic ( p = 5 ) , 
 septic ( p = 7 ) , and nonic ( p = 9 ) 
 polynomials provides improved accuracy , but at greater 
 computational cost .
empirical results show for degree p interpolation the optimality of using c softening . an alternative approach for interpolating with c continuity
is to construct the hermite interpolant that 
 reproduces the function values and first derivatives at nodes .
hermite 
 interpolation along one dimension requires two basis functions , ( ) for the function values and ( ) for the first derivatives : performing hermite interpolation in 3-d requires storing 2 = 8 values per grid point , to account for 0 or 1 partial 
 derivatives independently along the coordinate dimensions . the 8-element 
 vector of nodal basis functions is similarly constructed by all products 
 of  and  independently for 
 each of the three coordinates .
the anterpolation in eq 7 becomes multiplication 
 of an 8-vector by a scalar , and the interpolation in eq 12 becomes the inner product of two 8-vectors .
the restriction 
 in eq 8 and prolongation in eq 11 become 8  8 matrix  vector products , where the 
 restriction matrices are transposes of the prolongation matrices . 

the grid cutoff calculations in eq 9 and eq 10 are also 8  8 matrix  vector products . 

the grid spacing for hermite interpolation is doubled from that of 
 cubic interpolation to compensate for the  density  
 of each grid point increasing by a factor of 8 .
both 
 msm and pme perform a splitting of the interaction potential , pme 
 into two parts and msm into l + 1 parts .
both methods 
 perform a charge - spreading ( anterpolation ) step to a grid of charges 
 and an interpolation step from a grid of potentials .
the major difference 
 between the two methods is that pme transforms the grid of charges 
 using a discrete 3-d fast fourier transform ( fft ) , followed by a simple 
 scaling of the transformed charges in fourier space and a second fft 
 to transform back to the grid potentials , while msm performs nested 
 interpolation entirely in real space .
the use of the fft for pme is 
 predicated on the assumption that the molecular system is periodic 
 in all three dimensions , whereas msm can be used flexibly with or 
 without periodic boundary conditions .
although the asymptotic 
 complexity of pme is o(n log n ) due to the use of the fft , the overall operation count 
 of pme is relatively small for moderate system sizes , as the fft by 
 itself is economical in the number of arithmetic operations required . 

however , computing the fft requires large strided memory accesses 
 that are inefficient for modern computer architectures .
moreover , 
 the fft does not provide good parallel scaling , as the communication 
 pattern for implementing a 3-d fft is many - to - many , essentially the 
 same as a matrix transpose .
msm uses a hierarchical algorithm with 
 close neighbor communication and arithmetically dense localized 3-d 
 convolutions that map well to multicore processors with vectorized 
 instructions . in spite of having a higher operation count for moderate 
 system sizes ,
the use in msm of a more efficient communication pattern 
 combined with more effective use of computer processing resources 
 offers better parallel scaling and performance compared to pme .
comparison 
 of the force accuracy of msm and pme in namd . shown 
 is the relative error in mass - weighted force versus cutoff distance 
 for simulations of a 50  cube of 4142 water molecules . for msm
, 
 the cutoff distance is used to control accuracy for a given interpolation 
 and splitting .
although msm provides lower numerical accuracy than 
 pme ( specifically , there is about one digit of precision difference 
 at the 12  cutoff distance between cubic msm and quartic pme 
 with 1.2  grid spacing ) , force errors within 1% ( i.e. , less 
 than 0.01 ) are deemed sufficient for application of md to biomolecular 
 systems . the numerical accuracy of msm 
 compared to pme is addressed in figures 2 and 3 for simulations of 
 an equilibrated 50  cube of 4142 water molecules .
the figure 
 compares namd s msm using cubic , hermite , and quintic interpolation , 
 discussed above , to namd s smooth pme method using quartic 
 and sexic interpolation of the structure factors .
the employed pme grid spacings of 1.2 and 1.5  are 
 representative of the method s use in practice ; for instance , 
 gromacs uses 1.2  as its default 
 spacing and namd uses 1.5  as its default maximum spacing to 
 provide a safety check to the user . the relative error in average 
 mass - weighted force plotted in figure 2
the exact 
 forces fi have been calculated from a highly 
 accurate pme force calculation using octic ( eighth order ) interpolation 
 with the cutoff between short- and long - range parts set to 40  
 and the grid spacing set to 0.1  , making use of the fact that 
 pme converges to the ewald summation in the limit as the cutoff goes 
 to infinity and the grid spacing goes to zero .
the error in potential 
 energy presented in figure 3 has been averaged 
 over  = 1000 frames of a 2 ns trajectory , measured relative to a 
 highly accurate pme 
 calculation of the potential un .
msm provides lower numerical accuracy than pme , but force 
 errors within 1% are deemed sufficient for md applications , so that 
 the force error is not greater than the error incurred from time stepping 
 employed in the integration of the newtonian equations of motion . 

energy conservation for cubic msm with the 12  cutoff between 
 short- and long - range parts adopted here has been verified by simulating 
 the cube of water for 1 ns in the nve ensemble with 1 fs time stepping 
 and observing that the standard deviation of total energy e is 0.3% of either the standard deviation 
 of the kinetic energy or the standard deviation of the potential energy , 
 which is within the 20% criterion considered for energy conservation 
 in an md simulation .
comparison of the potential energy accuracy of msm and 
 pme in namd . 
 shown for the 50  cube of 4124 water molecules is the relative 
 error in potential energy , averaged over 1000 frames of a 2 ns simulation 
 trajectory , plotted versus cutoff distance .
the parallelization of msm is implemented using charm++ , the parallel 
 language extension of c++ used by namd that provides a message - driven 
 programming paradigm expressed as asynchronously executable objects .
namd employs a hybrid data and force decomposition for scalable 
 parallelism , decomposing the atoms spatially into uniform patches and also decomposing into objects the work of calculating 
 the short - range nonbonded forces between nearest neighbors .
msm and pme both make use of the existing infrastructure in namd 
 for calculating their short - range parts , modifying the functional 
 form of the interaction to fit each respective method .
for the long - range 
 part of msm , the grids are decomposed into blocks of grid points , 
 analogous to the namd patches . as with the short - range force calculation
, 
 sufficient work for scalable parallelization is available by carrying 
 out the grid cutoff calculations involving interactions between neighboring 
 msm grid blocks through separately schedulable work objects .
the restriction 
 and prolongation calculations introduced above require a much smaller 
 amount of work , so they are best kept on the same processor as their 
 corresponding grid blocks , as the overhead of communicating data to 
 a remote processor exceeds the runtime of the actual calculation step . 

similarly , the anterpolation and interpolation calculations involving 
 atom coordinates are kept processor - wise with the patches .
each 
 msm grid block receives charges from below , either from anterpolation 
 at the patch level or from restriction of a lower grid level .
after 
 all expected sub - blocks of charge are received and summed , the grid 
 block , assuming that it is not the top - level block , performs a restriction 
 to a local buffer and sends sub - blocks of the restricted charge up 
 to the upper level grid blocks . after sending the restricted charges
, 
 the grid block sends its charges to all of the block interaction work 
 objects to which it contributes .
a block interaction work object receives 
 a block of charges and performs its part of the grid cutoff calculation 
 to a block of potentials , then sends that block of potentials to its 
 designated grid block .
after all expected contributions to its block 
 of potentials are received and summed by the grid block , including 
 contributions from block interaction work objects and from prolongation 
 of a higher grid level , the grid block will either send the potentials 
 back to the patch level for interpolation , if it is a first level 
 grid block , or it will perform a prolongation to a local buffer and 
 send sub - blocks of the prolongated potential down to the lower level 
 grid blocks .
performance is improved through static ( i.e. , not 
 measurement - based ) 
 load balancing of the work objects and through charm++ message - driven 
 task prioritization .
the msm grid blocks are assigned to the multiprocessor 
 nodes in a round - robin manner ; for large processor counts , the processors 
 will outnumber the grid blocks .
the msm block interaction work objects 
 are similarly assigned so as to achieve an even node - level load distribution , 
 with preference given to nodes that hold either the source charge 
 block or the target potential block for a given block interaction , 
 in order to reduce the amount of internode communication required . 

once the node assignments have been made , the grid blocks and block 
 interaction work objects for each node are equitably divided among 
 its processors .
message priorities are assigned to give highest priority 
 first to sending restricted charges from the lowest grid levels up 
 to the highest , followed by the block interactions and prolongated 
 potentials from the highest grid levels down to the lowest .
the idea 
 is to prioritize the critical path of communication that goes through 
 the top grid level and back , in an effort to reduce the latencies 
 involved with receiving all contributions to a grid block of potentials 
 from prolongation and grid cutoff calculation before further prolongation 
 or interpolation can occur .
msm makes use of the existing namd reduction 
 infrastructure for summing the long - range potential energy .
nearly 
 all contemporary microprocessors offer some degree of support for 
 fine - grained parallel computing using so - called single - instruction 
 multiple - data ( simd ) processing units that execute , in lock - step , 
 the same arithmetic operation in parallel on several independent data 
 values .
cpu and gpu simd hardware is programmed using single - program 
 multiple - data ( spmd ) programming languages such as cuda , opencl , ispc , and openmp ( with simd 
 directives ) , or at the machine instruction level through the use of 
 vendor - specific assembly language or so - called compiler intrinsic 
 functions .
the msm algorithm is well suited to fine - grained parallel 
 execution on such cpu and gpu architectures , due to its predictably 
 and uniformly strided memory access patterns , and because its innermost 
 loops over grids can be decomposed into operations on vectors that 
 match the size and stride of the simd processing units provided by 
 the underlying computing hardware .
msm with hermite interpolation 
 is capable of most naturally utilizing the simd processing units of 
 cpus , available through the intel x86 compiler intrinsics for sse2 
 ( up to 4-element vector instructions ) or avx ( up to 8-element vector 
 instructions ) .
hermite interpolation directly computes 8-element grid 
 points , with grid point interactions expressed as 8  8 matrix  vector 
 products .
a single precision implementation of hermite interpolation 
 can map each grid point directly to two sse2 registers or to a single 
 avx register . with the matrices stored in row
a similar optimization is available to the other interpolation schemes 
 that employ single - element grid points .
the optimization involves 
 clustering the grid points into 2-point cubes . with this data reorganization
, the clustered grid point 
 interactions can be expressed as 8  8 matrix  vector products 
 using the same vector instructions as used for hermite interpolation . 

the grid point clustering might require up to one extra layer of grid 
 point padding in each of the x- , y- , and z - directions .
figure 4 illustrates the grid cutoff calculation from eq 9 modified to operate on a sequence of vectors , where a vector 
 of potentials for a clustered grid point is calculated by summing 
 over the  sphere  of clustered grid points comprising 
 vectors of charge .
for each cluster 
 of 
 grid point potentials , the algorithm loops over the sphere of charge 
 clusters and sums their contributions , calculating the sum of matrix 
the 
 grid stencil matrix elements are the precomputed kl(rml , rnl ) values corresponding to 
 the positions of the individual grid points .
many - core processors such as gpus and the intel xeon phi 
 natively 
 operate on larger vector lengths .
xeon phi processors support vector 
 lengths of up to 16 single precision floating point values , and contemporary 
 gpus natively operate on even larger vectors or arrays containing 
 64 to more than 512 elements at a time . in the case of xeon phi ,
it 
 should be possible to organize the msm data structures and use vector 
 mask operations to process two pairs of 8-element vectors per instruction . 

a gpu implementation of msm can be constructed by clustering grid 
 points into much larger tiles , and using a 3-d work 
 decomposition , with each gpu thread block or workgroup assigned to different tiles .
the most critical 
 issue for gpu performance is the use of memory layouts for the clustered 
 grid point tiles that enable so - called coalesced memory 
 access patterns that minimize the number of hardware memory cycles . 

the best parallel work decomposition strategy for gpus hinges upon 
 the amount of parallelism available in the work assigned to one gpu . 
 for msm workloads that result in a sufficiently large number of grid 
 points , a spatially oriented decomposition of grid point tiles and 
 their associated interactions allows the use of so - called gather style algorithms that perform potential accumulation 
 in on - chip registers . for msm workloads that present the gpu with 
 insufficient work for the spatial approach to be profitable
, a parallel 
 decomposition over interactions between msm grid blocks will provide 
 increased parallelism using a scatter type algorithm , 
 at the cost of an extra parallel reduction step required for summing 
 the partial potentials for each msm grid block prior to transferring 
 results back to the host cpu .
all simulations were carried out with a developmental 
 version of 
 namd implementing the parallelized msm 
 algorithm that has been subsequently released in namd version 2.10 . 

the r - respa multiple time - step integrator was applied with time steps of 2 and 4 fs for 
 short - range nonbonded and long - range electrostatic interactions , respectively . 

the settle algorithm maintained rigid 
 geometry for water molecules while rattle constrained the length of covalent hydrogen bonds .
the calculation 
 of nonbonded interactions excluded pairs of atoms covalently bonded 
 to each other or to a common atom .
all semiperiodic systems were 
 initially equilibrated for 500 ps with full periodic boundary conditions 
 in the npt ensemble , where pressure was kept constant at 1 atm by 
 the langevin piston method , before running 
 production simulations in the nvt ensemble using semiperiodic boundary 
 conditions .
electrostatic interactions are treated either using 
 msm or pme . the namd cutoff distance , 
 set to 12  , defines 
 both the van der waals truncation distance and the splitting between 
 short- and long - range parts of the electrostatic interactions for 
 both methods .
the msm simulations use the default cubic interpolation 
 with c splitting and the default grid 
 spacing of 2.5  , except for the performance comparison , which 
 uses hermite interpolation with c splitting 
 and a grid spacing of 5  .
the pme simulations use the default 
 quartic interpolation with a grid spacing of 1.2  . below 
 we explain how semiperiodic and completely nonperiodic systems 
 are simulated in the case of msm electrostatic force evaluation . in 
 the case of fully periodic systems ,
pme and msm simulations are specified 
 identically in namd , except for the choice of electrostatic algorithm . for semiperiodic systems ,
the periodicities 
 are established by defining just one or two of the three cell basis 
 vectors .
msm does not require any special conditions for the cell 
 basis vectors , except for linear independence .
however , any constraints 
 for containing atoms within the boundaries along nonperiodic dimensions 
 are best described by establishing orthogonal periodic basis vectors 
 aligned to the axes of the x , y , z - coordinate system chosen accordingly . to demonstrate the 
 use of msm
, we consider a membrane bilayer with asymmetric content 
 across the membrane in the form of a difference in ion concentration 
 above and below the membrane .
for this purpose , we utilize a simulation 
 cell that is periodic in the x , y - plane only . along the z - axis , we keep the simulated 
 system constrained by a harmonic restoring potential to a finite interval 
 with a lower and upper boundary at z =  a. the boundary is realized through a containing force along 
 the z - direction of k ( z  a ) , for z > a and + k ( z + a ) , for z < a , where k = 3 kcal/(mol ) ; 
 the condition is implemented using the tclbc scriptable boundary condition 
 feature of namd .
to demonstrate 
 a simulation of a 
 completely nonperiodic system , we consider a protein solvated together 
 with ions in a spherical water droplet .
the components of the system , 
 in particular the water molecules , have a strong tendency to remain 
 together at room temperature when placed in vacuum , though once in 
 a while a water molecule or even more rarely an ion or the whole protein 
 will evaporate .
to prevent such evaporation from 
 happening , we surround a water droplet by a repulsive boundary surface 
 of a radius a chosen 5  larger than the radius 
 of the droplet . for this purpose
we introduce a repulsive radial force 
 directed toward the center with magnitude k ( r  a ) , for r > a ; here , r is the radial position of any 
 atom in the system and the force constant is k = 
 10 kcal/(mol ) , implemented using the spherical 
 boundary condition feature of namd .
the center of mass of the protein 
 is constrained to the center of mass of the droplet in order to keep 
 the protein near the droplet center and well solvated .
in the following , we illustrate 
 the capabilities of molecular dynamics 
 simulations with namd stemming from the msm algorithm .
we first demonstrate 
 the accuracy of the msm algorithm and of its implementation in namd 
 through comparison with namd simulations based on the well - tested 
 pme algorithm .
we then present results of simulations for semiperiodic 
 and nonperiodic systems , until now unfeasible in namd having relied 
 on the pme algorithm .
we finally compare the scaling of namd on multiprocessor 
 machines running pme- and msm - based simulations side - by - side .
properties compared are the water oxygen  oxygen radial 
 distribution function ( a ) , a structural property , and the kirkwood 
 g - factor ( b ) , a dielectric property .
water is fundamental 
 to molecular processes in living cells , often due to the electrostatic 
 interactions among water molecules themselves and among water molecules 
 and proteins or ions that give rise to strong dielectric screening . 

accordingly , water is an excellent test bed for the computational 
 treatment of electrostatic forces . indeed
numerous structural , dynamic , 
 dielectric , and mechanical properties of water can serve to illustrate 
 the accuracy of electrostatic force descriptions by means of msm and 
 pme , for example , density , radial distribution function , diffusion 
 constant , dielectric constant , distance - dependent kirkwood factor 
 or surface tension .
to compare simulation results obtained with msm 
 and pme electrostatic descriptions , a box of water was simulated for 
 6 ns in an npt ensemble with a modified tip3p water model in the charmm 
 force field .
for all water properties calculated , agreement 
 is found between msm and pme , indicating that msm furnishes the same 
 accuracy as pme in electrostatic descriptions .
the water self - diffusion 
 constant was obtained from the mean - square displacement using the 
 einstein relation .
the surface tension was calculated by lz((pxx + pyy)/2 + pzz)/2 , where pij is the ij component of the pressure tensor , and lz is the length of the simulation cell in the direction 
 normal to the surface .
the distance - dependent kirkwood g - factor was 
 resolved as contributions from spherical shells of radius r by employingwhere j is the dipole moment of water and  is the dipole magnitude . 

the dielectric constant  was calculated usingwhere 
 the collective rotational dipole moment 
 is md = j  =  1nj .
compared are resulting electrostatic potentials 
 at an air  water interface ( a ) and at a membrane ( popc)water 
 interface ( b ) .
the distances from water layer center and membrane 
 center are measured along the interface normal .
a key electrostatic property 
 of water arises at the water - membrane interface in living cells and 
 at the air 
only 
 reliable electrostatic treatments can predict the correct interface 
 potential . as shown in figure 6
water 
 and membrane  water interfaces , in case of either a fully periodic 
 boundary or a semiperiodic boundary .
the results show that msm provides 
 a reliable description of the long - range electrostatic interaction . 

if truncated electrostatic interactions are applied with a cutoff 
 of 1.2 nm , an interface potential that is too large is found , as seen 
 in figure 6 . only by elongating the electrostatic 
 interactions , for example , by elongating the cutoff to 1.6 nm ,
is 
 the calculated interface potential close to the value determined by 
 means of simulations employing msm or pme that account for the full 
 long - range interaction .
formation of a pore in a membrane induced by 
 an ion concentration 
 difference .
a majority of cations were placed on one side of the membrane 
 and a majority of anions on the other side , creating an electric field 
 across the membrane .
the electric field causes membrane pore formation . 
 for lipids and water molecules , the following
coloring has been chosen : 
 carbon ( green ) ; oxygen ( red ) ; nitrogen ( blue ) ; phosphorus ( tan ) ; hydrogen 
 ( white ) .
for the sake of clarity , lipid 
 tail hydrogen atoms are not shown . all interface potentials (z )
were 
 calculated 
 viawhere (z ) is the charge 
 density as a function of z. cellular membranes surrounding living cells commonly 
 maintain a voltage gradient for battery - like energy storage and for 
 fueling many types of membrane processes such as atp synthesis or 
 ion vectorial transport .
the potential is due to a difference in charge 
 distribution across the membrane . in the case that the charge distribution 
 difference and , along with it , the electrical field become too large , 
 a pore is formed in the membrane in a process termed electroporation .
for example , by having more cations on one side of the membrane 
 and more anions on the other , a strong electric field can be generated 
 across the membrane , thereby inducing electroporation . in simulations employing pme ,
an electric field across a membrane , 
 created by having different numbers of cations and anions on each 
 side of the membrane , requires actually two different water compartments 
 containing positive and negative ions separated by two membranes in 
 order to permit a periodic boundary condition along the membrane normal . 
 in case of simulations employing msm without periodicity in the z - direction
, only a single membrane is required to separate 
 positive and negative ions , roughly halving the computational effort 
 while describing a situation that is more realistic than arising in 
 the case of pme - based simulations .
md simulations with msm were 
 performed for a membrane bilayer made 
 of popc lipids and employing semiperiodic boundary conditions . with 
 an average nacl concentration of 0.15 m , but more na ions 
 on one side of the membrane and more cl on the 
 other , an electric field of 0.7 v / nm was created across the membrane . 

the numbers of cations and anions were equal and the system charge 
 was neutral .
as illustrated through a simulation snapshot in figure 7 , ions diffused through the pore from one side of the membrane 
 to the other .
lipid head groups kept interacting with the ions within 
 the pore and created a toroidal pore structure .
exactly the same behavior 
 resulted in simulations with pme , but requiring two membranes .
bioengineering is today developing 
 nanoscale sensors for medical diagnostics , for example , for cost - effective 
 dna sequencing .
one type of sensor involves nanopores situated in 
 graphene sheets , the sheets acting much like a cellular membrane . 

md simulations offer guidance to nanosensor development by offering 
 engineers microscopic views of the measuring processes involving nanosensors . 

the ability of msm to simulate single membranes , illustrated above , 
 can be put to good use in the case of graphene nanopores .
we illustrate 
 this capability of msm by describing ion conduction through a nanopore 
 embedded into a single graphene sheet .
the sheet separates , just as 
 the membrane in the previous example , ions initially at different 
 concentrations above and below the graphene sheet .
semiperiodic simulation 
 based on the msm algorithm describing ion 
 conduction through graphene nanopores .
graphene carbon atoms are colored gray . for the sake of clarity , 
 water molecules are not shown .
( c ) 
 number of ions passed through the pore at time t. 
 the slope of the curves shown is the translocation rate , which is 
 affected by the nanopore edge charge . in the simulations , we chose the edge of the graphene nanopore 
 to contain either electrically neutral or positively or negatively 
 charged atoms . after carrying out msm - based simulations on the respective 
 semiperiodic systems
, we compared the ion conductance through neutral , 
 positive , and negative graphene nanopores , where 2 m nacl , described 
 by charmm force parameters , was placed 
 on one side of the graphene and none on the other side .
the nanopore 
 measured 6 nm in diameter , with 30 carbon atoms forming the pore edge . 
 to charge the nanopore ,
12 carbon atoms distributed evenly along the 
 pore edge were selected , and unit positive or negative charges were 
 placed on them . in comparison to employing charged chemicals to decorate 
 the pore , nanopores with these charged 
 carbon atoms have the same physical pore size as the neutral pore . 
 to neutralize the overall system charge ,
12 additional na or cl were added to the 2m - concentration side 
 when the nanopore edge was negative or positive , respectively .
right after the simulation started , ions diffused through the nanopore 
 regardless of the pore charges , as shown in figure 8a . in the cases of charged nanopore edges , counterions translocated 
 faster than the co - ions , until the charges on the two sides of graphene 
 became equal .
the simulations revealed that both positive and negative 
 ions diffuse through the nanopore at the same rate , regardless of 
 pore edge charge , as shown by the slope of ion passage in figure 8c .
for the charged nanopore edges , however , the 
 ion translocation rates are smaller than in the case of the neutral 
 nanopore .
the ion translocation rate through a neutral pore is 2.1 
 ns , while it is 1.7 ns for 
 the positively charged pore edge and 1.6 ns for 
 the negatively charged pore edge . for charged nanopore edges ,
counterions 
 tend to localize around the pore region . as shown in figure 8b ,
therefore , the effective nanopore size is reduced for charged nanopores , 
 and thus , ion translocation rates are smaller .
all simulations 
 were terminated at 18 ns , at which point the ion 
 concentration difference between the two sides of the graphene sheet 
 had decreased to 70% of its original value .
msm can be 
 applied not only to semiperiodic systems but also to 
 systems without any periodic boundary .
for example , proteins can be 
 simulated in a large water droplet with msm .
the total number of water 
 molecules in a droplet is smaller than the number of water molecules 
 in a periodic rectangular or dodecahedral box of the same maximum 
 dimension .
moreover , simulations with nonzero net charge can be performed 
 with msm on a nonperiodic droplet , while the net charge must be chosen 
 to vanish for pme - based simulations .
nonperiodic simulations based on the msm 
 algorithm describing a 
 -hairpin protein motif ( a ) and the protein trp - cage ( c ) inside 
 a water droplet .
protein structures are characterized color - wise as 
 follows : -sheet ( yellow ) ; -helix ( purple ) ; ordered turn 
 ( blue ) ; disordered coil ( white ) . the rmsd values for -hairpin 
 and trp - cage , defined here with reference to initial structure ,
one needs to be aware that biomolecules in single droplets 
 are 
 subject to water surface tension created naturally by the water  air 
 boundary .
such tension arises for droplets suspended in air and can 
 be significant for droplets of small size . in the case of bulk water 
 surface tension does not arise , neither in a real system nor in a 
 system simulated with periodic boundary conditions and completely 
 filled elementary cells .
such avoidance of surface tension may be 
 an advantage of fully periodic msm and pme calculations ; future calculations 
 will test if this advantage is really significant . to illustrate 
 the possibility of simulations of a single water 
 droplet , small proteins , -hairpin and tryptophan cage ( pdb : 1l2y ) , were each simulated in a small spherical water droplet . 

each protein was placed at the center of the droplet , which was 20 
  larger than the radius of gyration of the proteins , as shown 
 in figure 9 .
the -hairpin is composed 
 mostly of -sheet structure , whereas the tryptophan cage is 
 composed mostly of -helical structure .
the protein structures 
 were monitored by means of the root - mean - square displacement ( rmsd ) 
 as shown in figure 9 ; both proteins maintain 
 their structures over 1 s of simulation , as indicated by rmsd 
 values of less than 2  for most of the simulation time . comparing the nonperiodic msm approximation 
 to an exact electrostatics 
 calculation over all pairs of interacting atoms , the droplet containing 
 the -hairpin showed a difference of 0.5 kcal / mol from a total 
 energy of 2.3  10 kcal / mol , with a relative 
 error in average mass - weighted force of 8.5  10 , as measured by eq 14
. performance 
 and scaling of msm and pme are comparable , despite the fact that the 
 namd implementation of pme has had many more years of development 
 effort than msm .
msm is shown to perform best with hermite interpolation , 
 an approach that , as mentioned previously , exactly reproduces function 
 values and derivatives at grid points .
vector 
 products resulting from hermite interpolation in three dimensions 
 are calculated very efficiently using cpu vector instructions .
a comparison of the namd implementations of msm and pme was performed 
 on the -hairpin protein solvated in a cubic box of 3757 water 
 molecules for a total of 11.5 k atoms and simulated with multiple 
 time step integration , in which long - range electrostatic forces were 
 evaluated every 4 fs using 16 processors ( 1 node ) in a dell poweredge 
 compute node of the tacc stampede system .
the comparison shows that 
 pme achieves a simulation rate of 18 ns per day , whereas msm using 
 hermite interpolation with c splitting 
 achieves 15 ns per day and using cubic interpolation with c splitting only 11 ns per day .
performance of msm and pme on apoa1 ( 92 k atoms ) 
 with 1 fs time 
 stepping , on a cray xe6 ( 32 cores per node ) .
msm using hermite interpolation 
 with c splitting , a 5  grid spacing , 
 and 12  splitting distance achieves 19.6 ns / day with 1536 cores , 
 whereas pme using quartic interpolation with 1.2  grid spacing , 
 and 12  splitting distance tops out at 17.6 ns / day with 1024 
 cores .
however , a comparison for the 
 case of a larger simulation using 
 more processor cores demonstrates the scaling benefits of msm .
the 
 namd benchmark system apoa1 ( 92 k atoms ) simulated with 1 fs single time stepping using varying numbers of 
 cray xe6 nodes ( 32 cores per node ) of the blue waters jyc test computer 
 system shows that , while pme achieves faster simulation rates for 
 smaller core counts , msm achieves comparable scaling results with 
 slightly faster simulation rates for larger core counts .
the performance crossover 
 point occurs between 256 and 512 cores ; the msm performance continues 
 to scale up to 1536 cores , whereas the pme performance reaches a plateau 
 at 1024 cores .
the multilevel 
 summation method ( msm ) has recently been implemented 
 in the molecular dynamics program namd and provides 
 now a successful alternative in computing electrostatic forces arising 
 in molecular dynamics ( md ) simulations of biomolecular systems . the 
 practical accuracy produced by msm
, although less than the numerical 
 accuracy available in typical use of particle  mesh ewald ( pme ) , is sufficient for md , as shown by a comparison of calculated water 
 properties such as density , diffusion constant , dielectric constant , 
 surface tension , radial distribution function , and distance - dependent 
 kirkwood g - factor .
agreement is also found between msm and pme in 
 calculations of the interface electrostatic potential at air 
msm has key advantages 
 over pme . with nested interpolation of softened 
 pair potentials in real space , msm permits simulations without periodic 
 boundary conditions in some or all of the simulated dimensions . instead , 

pme in its common implementations can only be applied with fully periodic 
 boundary conditions .
for example , simulations with msm can be performed 
 with periodic boundary conditions along the x , y - directions but not along the z - direction . 

such simulations allow different ion concentrations to be placed on 
 either side of a membrane bilayer , producing membrane pore formation 
 in msm - based simulations .
a related example for now feasible simulations 
 with namd is ion conduction through neutral and charged graphene nanopores 
 in the presence of a large ion gradient .
msm can be further applied 
 to spherical systems without any periodicity . to illustrate this , 
 two proteins , a -hairpin motif and the protein tryptophan cage , 

have been simulated in spherical water droplets . in this case , the 
 simulation with msm not only may be more realistic but also requires 
 a smaller atom count , roughly described by the ratio ( 4r/3)/8r  
 0.5 ; now with msm available , the argument that periodic systems , even 
 though larger per elementary cell , avoid surface effects such as surface 
 tension , can actually be tested . in any case
, the simulation examples 
 described in the present study demonstrate that msm offers flexibility 
 in regard to boundary conditions and provides more alternatives in 
 simulation design than does pme .
moreover , msm offers improved 
 parallel scalability over pme . unlike 
 the pme algorithm that is based on the fourier transform , requiring 
 communication that poses a limit to scalability , msm instead
has a 
 hierarchical structure with highly localized arithmetic operations 
 that permit effective utilization of modern vector computational hardware 
 units . according to our performance benchmarks ,
the superior scalability 
 of msm makes it competitive with pme and actually faster than pme 
 in case of simulations of large systems .
ongoing is the development 
 of improved interpolation for msm to 
 provide higher accuracy for a given polynomial degree p without increasing the computational cost .
future work includes 
 also the calculation of dispersion forces without truncation with 
 msm - based namd ; these forces , in particular , their long - range contribution , 
 are considered to be important for membrane properties . with support in namd also for long - range dispersion forces , the 
 present charmm - prescribed 12  cutoff / splitting distance can 
 be used as a true control for msm accuracy .
high performance simulations 
 will then be able to achieve practical accuracy with a reduced splitting 
 distance , where a splitting distance of between 8 and 9  is 
 expected to double the overall simulation performance .
other future 
 work includes extending our earlier development efforts on gpu - accelerated 
 msm for calculating electrostatic potential maps to support msm - based namd on large - scale gpu - accelerated 
 parallel computers .
 the 
 multilevel summation method ( msm ) offers an efficient algorithm 
 utilizing convolution for evaluating long - range forces arising in 
 molecular dynamics simulations . shifting the balance of computation 
 and communication , 
 msm provides key advantages over the ubiquitous 
 particle  mesh ewald ( pme ) method , offering better scaling on 
 parallel computers and permitting more modeling flexibility , with 
 support for periodic systems as does pme but also for semiperiodic 
 and nonperiodic systems . 
 the version of msm available in the simulation 
 program namd is described , and its performance and accuracy are compared 
 with the pme method . 
 the accuracy feasible for msm in practical applications 
 reproduces pme results for water property calculations of density , 
 diffusion constant , dielectric constant , surface tension , radial distribution 
 function , and distance - dependent kirkwood factor , even though the 
 numerical accuracy of pme is higher than that of msm 
 . excellent agreement 
 between msm and pme is found also for interface potentials of air  
 water 
 and membrane  water interfaces , where long - range coulombic interactions 
 are crucial . 
 applications demonstrate also the suitability of msm 
 for systems with semiperiodic and nonperiodic boundaries . for this 
 purpose , 
 simulations have been performed with periodic boundaries 
 along directions parallel to a membrane surface but not along the 
 surface normal , yielding membrane pore formation induced by an imbalance 
 of charge across the membrane . using a similar semiperiodic boundary 
 condition , ion conduction through a graphene nanopore driven by an 
 ion gradient 
 has been simulated . 
 furthermore , proteins have been simulated 
 inside a single spherical water droplet . 
 finally , parallel scalability 
 results show the ability of msm to outperform pme when scaling a system 
 of modest size ( less than 100 k atoms ) to over a thousand processors , 
 demonstrating the suitability of msm for large - scale parallel simulation . 
Introduction
Multilevel Summation Method
Simulation
Protocols
Results and Discussion
Conclusion