J. Phys. Chem. 1996, 100, 2581-2587
2581
Calculating Electrostatic Interactions Using the Particle-Particle Particle-Mesh Method with Nonperiodic Long-Range Interactions Brock A. Luty* and Wilfred F. van Gunsteren Laboratorium fu¨ r Physikalische Chemie, ETH Zentrum, CH-8092 Zu¨ rich, Switzerland ReceiVed: July 10, 1995X
The particle-particle particle-mesh (PPPM) method is an accurate and computationally efficient method for calculating interactions in molecular simulations. The PPPM method is based on separating the total interaction between particles into the sum of short-range interactions, which are computed by direct particleparticle summation, and long-range interactions, which are calculated by solving Poisson’s equation using periodic boundary conditions (PBCs). For nonperiodic systems there is a concern that use of PBCs may introduce an artificial periodicity into the system. In this article we propose, and test, the use of boundary conditions for the long-range interactions which mimic surrounding the explicit system by high dielectric material such as water.
1. Introduction Accurate treatment of the physical interactions in a molecular system is an essential requirement for performing a reliable computer simulation. The long-range nature of Coulomb’s law, however, makes evaluation of the electrostatic interactions particularly difficult to handle in a computationally efficient manner. At the present state of the art, the correct treatment of electrostatic interactions, including the efficient handling of the long-range effects, can be considered a bottleneck to further development of the molecular simulation methods.1 For naturally periodic systems, such as crystals, the electrostatic interactions can be accurately and efficiently evaluated using lattice sums. In the lattice sum methods, Coulomb’s law is split into the sum of two functions. The first is a smooth periodic long-range function which is well-suited for solution using Fourier methods. The second is a short-ranged, rapidly varying function, which is best handled by direct summation over neighboring particles. The most popular of the lattice sum methods is the Ewald method.2,3 In previous publications,4,5 we compared the computationally efficiency and underlying theory of the Ewald method to the particle-particle particlemesh (PPPM) method developed by Hockney and Eastwood.6 In contrast to the Ewald method, which sums the Fourier series of the long-range function analytically, the PPPM method employs the computationally efficient fast Fourier transform (FFT). Even with the computational overhead of constructing the discrete grids required for the FFT, the PPPM method is a competitive alternative to the Ewald method for relatively small systems. For larger systems, the PPPM method becomes increasingly rewarding as a result of the N log2(N) scaling of the FFT relative to the N2 scaling of the analytical summation (N is the number of charges in the system). When simulating a nonperiodic system, which is necessarily of a finite extent, there will be a surface or boundary which requires consideration. In the simplest approximation, we can consider the explicit system to be surrounded by a vacuum. Under these circumstances, the atoms in the neighborhood of the surface will feel unbalanced forces which will tend to distort their behavior from that found in the interior of an extended system. If we are interested in, for example, a bulk fluid or a * Author to whom correspondence should be addressed. X Abstract published in AdVance ACS Abstracts, January 15, 1996.
0022-3654/96/20100-2581$12.00/0
macromolecule in solution, then these surface effects are undesirable artifacts which should be minimized. The distortive effects of the vacuum outside the explicitly simulated region can be reduced by surrounding the system with an implicit bath. The influence of the bath can be approximated by altering the dynamics of the outer layers of explicit atoms and by appropriately modifying the expression for the energy and forces on the explicit atoms. Typically, the positions of the outer layers of explicit atoms are restrained and/or their dynamics are coupled to a heat bath to simulate energy exchange with the surroundings.7 As a first order approximation, the longer-range forces of the bath on the explicit system can be estimated by a static equilibrated configuration of bath molecules8 or from the radial distribution functions expected for the bath molecules.9,10 Both of these methods ignore the dielectric response of the bath and therefore correspond to surrounding the explicit system with a hydrophobic wall.11 For a polar system such as water, this causes a nonphysical preferential orientation of the explicit molecules near the boundary with the implicit system.11 As a better approximation for mimicking the longer-range electrostatic effects of a surrounding polar system, one can consider the explicit system as a low dielectric cavity immersed in a high dielectric continuum bath and solve the corresponding macroscopic electrostatic equations. If the geometry of the explicit system can be approximated by a simple sphere, then the reaction of a surrounding dielectric continuum to an explicit charge can be determined analytically. This reaction potential is found by solving Poisson’s equation with a dielectric discontinuity at the surface of the sphere and takes the form of an infinite series. If the surrounding bath has a large dielectric constant (e.g., water) relative to the explicit system (usually taken as vacuum), then Friedman12 suggested that the infinite series can be, to a large extent, approximated by a single term. This single term is nicely visualized as each charge in the explicit system generating an image charge outside the sphere, which in turn reacts back on all the charges in the explicit system. The implementation and testing of this method within a molecular simulation framework has been discussed by Wallqvist.13 Specific attention must be given to the fact that as a charge approaches the dielectric boundary between the explicit and implicit system, the attraction between the charge and its image becomes infinite. To compensate, it is © 1996 American Chemical Society
2582 J. Phys. Chem., Vol. 100, No. 7, 1996 necessary to augment the force field to include a repulsive term between the explicit system and the boundary or, for example, ignore the reaction field generated from the outer most layers of explicit water.14 Using a simple repulsion potential between the dielectric boundary and the explicit system, Wallqvist found that at least a shell of two water layers was necessary to remove any spurious effects of the boundary conditions. If the explicit system does not conform to a spherical geometry, then the newly developing computational continuum electrostatic methods can be used to estimate the reaction potentials. Juffer has examined the situation, where not only is the system of an irregular geometry, but the boundary is allowed to fluctuate during the simulation.15 He extensively investigated the effect of the parameters of the boundary atoms on the structure and dynamics of the explicitly simulated region. The resulting technique, termed dynamic surface boundary condition (DSBC), performed quite well, reproducing many of the properties found in bulk simulations. He did however note that the surface waters displayed a preferential orientation, which was attributed to a lack of specific interactions between boundary atoms and atoms just outside the boundary, and that there was some difficulty tuning the parameters to keep explicit atoms from evaporating into the continuum region. The computational cost of the method is also relatively high, and time-saving techniques, such as the infrequent updating of the reaction field, were required. A very popular method, which avoids the difficulties associated with setting up and determining the appropriate parameters for boundary atoms, is to simply use periodic boundary conditions in a nonperiodic system. As mentioned above, the lattice sum methods can then be used to accurately evaluate electrostatic interactions. However, there is some concern that using full periodic boundary conditions on nonperiodic systems may introduce artificial periodic effects into the system. Partially for this reason, and partially for computational expediency, the electrostatic interactions are typically truncated. In the minimum image truncation method, one only considers the interactions of a given particle with the nearest periodic images of the remaining particles. For systems with only short-range repulsion and dispersion interactions, the minimum-image truncation is a reasonable assumption when the simulated region has a minimum extent greater than 6 times the characteristic size of a particle.3 For systems which include long-range electrostatic interactions, the minimum image approximation causes artifacts. For example, in a cubic system the minimum image approximation will favor the placing of similarly charged species in opposite corners of the cube, thus imposing a structure on what should be an isotropic system.16 In the more common spherical truncation method, the interactions of a given particle are only calculated with neighboring particles within a sphere centered on the particle of interest. Typically, the radius is taken less than half the minimum extent of the system; therefore the minimum image requirement is also satisfied. However, it has been dramatically demonstrated that even this somewhat more isotropic approximation introduces an artifical structure on the system with similarly charged species segregating just outside the chosen spherical cutoff distance.17 For homogeneous systems (e.g., ions dissolved in water) a very simple correction scheme is available that reduces this artifical effect. In the generalized reaction field method,18 each charge is considered individually to be located at the origin of an explicit spherical system. As in Friedman’s method, the macroscopic electrostatic equations are then solved using spherical symmetry to estimate the reaction of a surrounding implicit dielectric continuum (which may include statistically
Luty and van Gunsteren distributed ions) to the explicit system. This long-range interaction is then added to a short-range explicit particleparticle summation over the neighboring atoms within the explicit sphere. When summed over all particles, the result is a very simple modification of the pair forces which entails only a slight increase in computational costs relative to the evaluation of Coulomb’s law. The radius of the explicit interaction sphere, passed which point the continuum approximation must be valid, is not known a priori and must be estimated. Tironi et al.18 present results for radii of 0.9 and 1.2 nm which show that the generalized reaction field gives a vast improvement over the spherical cutoff for a system of ions dissolved in water. But even with the correction, there are still noticeable differences between the two cutoff radii which imply that this distance should be somewhat greater than 0.9 nm. It is not obvious how to extend this method to nonhomogeneous systems (e.g., a macromolecule in solution), where the explicit sphere can be surrounded by portions of the system which may not respond as a uniform dielectric. In this work, we propose a hybrid method which is a combination of the ideas of the PPPM with Friedman’s reaction field method. As in the PPPM method, a short-range interaction between nonbonded atoms, which includes dispersion and repulsion and a short-range component of the electrostatics, is constructed. These interactions will be handled using explicit particle-particle summation and spherical truncation with minimum image boundary conditions. The remaining longrange interactions will not be determined using the usual periodic boundary conditions but will be estimated using boundary conditions to mimic surrounding the explicit system with a high dielectric continuum. A charge near the center of the explicit system will feel interactions very similar to those which would be estimated using Friedman’s method. As a charge approaches the edge of the explicit system, an increasing portion of the long-range force will come from the dielectric continuum and it may feel an unbalanced force relative to the long-range force from the explicit system. The forces in this region will depend on the exact nature of the unphysical splitting of the electrostatic interactions. It is hoped that the short-range forces, which act periodically, will impart to the waters near the dielectric boundary a physically reasonable structure and dynamics and that any disturbance will not propagate deeply into the explicit region. We note that this method avoids the problem associated with determining the parameters for an explicit boundary and, by only applying the reaction field to the smooth long-range portion of the interaction, avoids the singularity associated with a point charge crossing a dielectric boundary. Finally, a major advantage is that, in the limit that the surrounding continuum has a very large dielectric constant relative to the explicit system (i.e., in the same limit as Friedman’s image charge approximation), we can take over the very well developed algorithms of Hockney and Eastwood with only minor modifications.
2. Theory The pairwise nonbonded interaction energy, U, between particles at positions ri and rj, can, without loss of generality, be separated into the sum of two functions:
U(rij) ) Us(rij) + Ul(rij)
(1)
where rij ) |ri - rj| is the distance between particles. We choose the short-range function, Us, to include repulsion and dispersion and arbitrarily, but conventionally, choose the short-
PPM Calculation of Electrostatic Interactions
J. Phys. Chem., Vol. 100, No. 7, 1996 2583
range portion of the electrostatic interaction to be given by the complementary error function divided by the distance:
Us(r) )
Aij
Bij qiqj erfc(Rr) + r r12 r6 4π01
R3 exp(-R2r2) 3/2 π
2
(3)
With a sufficiently large R, the short-range portion of the electrostatic interaction can be adjusted to be significant only over the same range as the dispersion term. The long-range contribution to the interaction energy is found by solving Poisson’s equation for the long-range electrostatic potential, φl, due to the shielding charge distribution F(r) ) ∑jqjγ(r - rj):
-∇2φl(r) ) F(r)/01
G ˆ (k) ) e-k /(4R )/01k2
(2)
where repulsion and dispersion are modeled using a LennardJones 12-6 function with coefficients Aij and Bij, and qi and qj are the charges on the atoms. The electrostatic term is equivalent to the interaction of a point charge, e.g., qi, with a charge qj at rj, minus the interaction of qi with a shielding charge at rj, all in a medium with uniform dielectric constant 1. The shielding charge at rj has a magnitude qj and is distributed as a Gaussian with an extent governed by a constant R:
γ(r) )
Fourier transform of the Gaussian shielding charge distribution divided by the Fourier transform of the negative Laplacian:5,6,20
(4)
If we assume periodic boundary conditions for the differential equation and use Fourier methods, we arrive at the lattice sum methods of Ewald and of Hockney and Eastwood. If, however, we consider the explicit system to be surrounded by a dielectric continuum, we must solve a more general equation which takes into account the difference in the dielectric medium in the regions. Numerical methods such as finite-difference or boundary element methods can be used in these circumstances (see Davis and McCammon19). As mentioned in the Introduction, the dielectric constant of water is large (∼80) relative to that used for the explicit system (typically 1 ) 1). If, following Friedman, we take the limit that the ratio approaches infinity, then we are led to the considerable simplification of replacing the dielectric boundary between the explicit system and surrounding continuum by a perfectly conducting boundary at the surface of the explicit system. For a conducting boundary, the boundary potential will be identically zero, and if we are using Fourier methods in an orthorhombic geometry, the cosine terms will vanish leaving only the sine components. With this approximation, we can use the following algorithm, which is very similar to the original PPPM method. The short-range interactions, including a short-range component of the electrostatic interactions, are evaluated using standard periodic boundary conditions with the minimum image convention and spherical truncation. Any technique, such as the Verlet pair list,3 may be used to speed this portion of the calculation. For an orthorhombic system, the long-range part of the interactions are evaluated using the sine FFT to solve Poisson’s equation with conducting boundary conditions: (1) assign the charges from the atoms to a grid which encompasses the simulation region; (2) sine FFT the gridded charge density to obtain Fˆ(k) on the lattice; (3) multiply Fˆ(k) by the influence function G ˆ (k) at each lattice point to get φˆ l(k); (4) apply the sine FFT again, this time to φˆ l(k), which yields φl(r) at the grid points; (5) interpolate to get the potential at, and the forces on, each of the charges. If the discretization and mapping/ interpolation introduced no error into the solution of Poisson’s equation, then the influence function G ˆ (k) would simply be the
2
(5)
the reciprocal lattice vectors are k ) (π)(l1/L1, l2/L2, l3/L3), where L1, L2, and L3 are the lengths of the box edges and l1, l2, and l3 independently enumerate the grid points 0...N - 1. By taking into account the errors that the mesh solution will introduce into the solution of the Poisson equation, Hockney and Eastwood were able to derive an “optical influence function” which minimizes, in a least squares sense, the difference between the mesh force and the analytical solution.6 The expression for the optimal influence function for a nonperiodic system is presented in the Appendix (eq A1). It is important to note that for a given system the influence function is only calculated once and that the use of the optimal influence function (eq A1) significantly reduces the error relative to the simple estimate (eq 5). 3. Method There are two free parameters which must be estimated for the conducting boundary method. First, in contrast to the periodic boundary condition, the behavior of the system will depend on the exact splitting of the electrostatic interactions. In the particular case of the Gaussian shielding charges used here, this translates to a dependence on R, which governs the width of the Gaussian. The second parameter relates to the exact placement of the conducting box relative to the explicit system. Because of the irregular (corrugated) surface of the explicit system, it may be appropriate to place the conducting surface inside the surface which is used for calculation of the short-range force and coordinate wrapping (here performed on a molecular basis). For a cubic system, this reduces to a single parameter, δ: where 2δ is the difference in the length of the edge of the box which is used for the conducting boundary in the solution of Poisson’s equation minus the length of the edge of the box of the explicit periodic system. A negative δ then corresponds to placing the conducting box inside the explicit box. The effects of conducting boundary conditions, and, in particular, the parameters R and δ, were investigated by performing molecular dynamics simulations on cubic boxes of 512 SPC (simple point charge) waters.21 The equations of motion were integrated using a Verlet leap frog scheme3 with a time step of 0.002 ps. Bond lengths and angles were constrained using the SHAKE algorithm22 with a relative tolerance of 10-5. The temperature of the system was maintained at 300 K by weak coupling to a thermal bath23 using a time constant of 0.1 ps. The short-range interactions, including a portion of the electrostatic interaction, were evaluated using explicit summation over nearest neighbor particles within a spherical cutoff distance of 0.8 nm. A list containing potential neighboring particles was constructed every 20 time steps using a “skin” of 0.1 nm.3 The long-range electrostatics were evaluated at every time step using the PPPM method with either periodic boundary conditions or conducting boundary conditions. The 27-point TSC weighting function6 was used to allocate the charges to, and interpolate the potential from, the grid. A grid size of (32)3 was used, resulting in grid spacings, h, of 0.070.08 nm. An optimal four point central difference6 was used to estimate the forces from the potentials. Solution of Poisson’s equation was accomplished using IBM’s ESSL24 DRCFT3/ DCRFT3 or DSINF routines for performing the periodic or sine FFTs, respectively. We note that using the highly optimized
2584 J. Phys. Chem., Vol. 100, No. 7, 1996
Luty and van Gunsteren
TABLE 1: Boundary Conditions, Half the Difference in Conducting Box Size and Explicit System Periodic Box Size, δ, the Value of the Parameter r, Which Determines the Width of the Gaussian Shielding Charge, and the Short-Range Cutoff Radius Rc, Used in the Simulation of the Model System Consisting of 512 SPC Water Moleculesa boundary condition
δ, nm
R, nm-1
Rc, nm
DT, 10-9 m2/s
τ1, ps
τ2, ps
Q, (kJ/mol)/step
periodic conducting conducting conducting conducting
0.0 0.1 -0.1 -0.1
4.0 4.0 4.0 4.0 3.2
0.8 0.8 0.8 0.8 0.9
3.9 4.4 5.2 4.2 4.2
2.9 2.7 2.5 2.8 3.0
1.07 0.98 0.90 0.99 1.04
