10464
J. Phys. Chem. 1996, 100, 10464-10468
A Very Fast Molecular Dynamics Method To Simulate Biomolecular Systems with Realistic Electrostatic Interactions Piero Procacci Centre Europe´ en de Calcul Atomique et Moleculaire (CECAM), Ecole Normale Superieure de Lyon, 46 Alle´ e d’Italie, 69364 Lyon, France
Tom Darden National Institute of EnVironmental Health, Sciences Research Triangle Park, North Carolina 27709
Massimo Marchi* Section de Biophysique des Prote´ ines et des Membranes, DBCM, DSV, CEA, Centre d’EÄ tudes, Saclay, 91191 Gif-sur-YVette Cedex, France ReceiVed: January 26, 1996; In Final Form: April 18, 1996X
In this paper we describe the implementation of a very fast molecular dynamic method to realistically handle electrostatic interactions in simulations of solvated proteins. Our scheme is based on a recently proposed reversible multiple time step (r-RESPA) algorithm and a new modification of the particle mesh Ewald method. While the latter technique provides a fast and accurate representation of the Coulombic interactions for infinite systems, the r-RESPA algorithm exploits a separation of the particle force into components with increasingly longer time scales corresponding to contributions from short-, medium-, and long-range radial shells. By combining the two techniques we are able to reduce considerably the computational cost of molecular dynamics simulation of large biomolecular system without affecting energy conservation and dynamical properties. With respect to single time step simulations employing standard Ewald summation and rigid bond constraints, we obtain a speed-up of about 1 order of magnitude. Finally, our method is about 2.5 times as fast as simulations making use of spherical cutoffs. Since the majority of today biomolecular simulations use spherical cutoffs, we expect that our algorithm will find general applications in the field.
I. Introduction The dielectric properties of the aqueous solvent are of crucial importance to many biological processes such as folding of proteins, electron/proton transfer reactions, and phase transitions. Thus, any theoretical investigation of such phenomena must deal with the problem of how to correctly handle long-range electrostatic interactions. In the past, computer simulations of solvated biomolecules have routinely used a spherical cutoff for Coulombic forces. Although more accurate techniques have been available for a long time1 and routinely applied in simulation of molecular liquids, their computational cost for large systems such as solvated biomolecules has been until recent times too high to represent a valid alternative to spherical cutoff. Moreover, using conventional algorithms, simulations can be run at a reasonable computational cost only by constraining a large fraction of the degrees of freedom in the system. This approximation perturbs the density of states of the system2-4 and therefore all the related thermodynamic properties. Thus, new and more effective algorithms are needed to make it possible to simulate solvated biomolecules with realistic potentials. Recently, two approaches have been successful in improving the computational efficiency of Ewald summations for large periodic systems. First, a new scheme for Ewald summation called smooth particle mesh Ewald (SPME) method which scales with N log N, N being the number of particles, has been proposed.5 This method is based on the interpolation * To whom correspondence should be addressed. X Abstract published in AdVance ACS Abstracts, June 1, 1996.
S0022-3654(96)00295-X CCC: $12.00
of the reciprocal space sums by means of B-spline and evaluations of the resulting convolution using fast Fourier transforms. The SPME algorithm allows the use of Ewald summation at a computational cost comparable to that of spherical cutoff of 10 Å. Second, a method has been presented4 which combines the Ewald summation technique with a reversible multiple time step (r-RESPA) integration algorithm.6,7 In this approach, no degrees of freedom of the system are constrained and the force associated with each particle of the system is partitioned into four components which evolve in time with distinct and increasingly longer time scales. A suitable time scale separation is achieved by subdividing the direct space nonbonded interactions, inclusive of Coulombic and van der Waals contributions, into short-, medium-, and long-range shells. The fastest bonded forces are associated to the component with the shortest time scale, while the reciprocal space nonbonded interaction is included with the medium-range direct space forces. In combination with standard Ewald summation technique, this r-RESPA algorithm allows up to a 4-fold reduction in computational cost for system as large as 20 000 atoms with respect to a standard algorithm involving constraints.4 In this paper, we present a molecular dynamics method which combines the two schemes achieving a significant improvement in performance with respect the SPME and r-RESPA algorithms alone. The use of the SPME technique as opposed to standard Ewald summation within the framework of the multiple time step algorithm has implied modifications of the force partitioning. Such changes have taken into account the diminished weight of the reciprocal space sum contribution to the calculation © 1996 American Chemical Society
Method To Simulate Biomolecular Systems
J. Phys. Chem., Vol. 100, No. 24, 1996 10465
of the total force. The very fast molecular dynamics method presented here makes use of a reduced direct lattice cutoff and consequently smaller size shells for the force subdivision. The SPME contribution was included with the medium range shell as done previously in ref 4. Anticipating our results, we find that for a system composed of 20 627 atoms our technique is able to achieve speed-ups of up to an order of magnitude with respect to standard constraint algorithms for the same level of accuracy in the Coulombic potential. This paper is organized as follows. In section II we describe the interaction potential for solvated proteins. While a brief description of the SPME method is presented in section III, we provide the details of the r-RESPA implementation for solvated proteins in section IV. In section V we discuss the performances of our integration algorithms against samples of hydrated Cphycocyanin and hydrated reaction center containing 7040 and 20 627 atoms, respectively. The paper ends with a Conclusion. II. Interaction Potentials In the past, many force field have been developed for proteins.8-10 Typically, the total potential energy for solvated proteins can be written as
Vtot ) Vnb + Vb
(II.1)
Vnb ) Vvdw + Vqd + Vqr
(II.2)
Vb ) Vstretch + Vbend + Vtor
(II.3)
The first term on the right-hand side of eq II.1 is the nonbonded potential energy which includes the van der Waals, the electrostatic direct space, and the reciprocal space contributions, i.e. Vvdw, Vqd, and Vqr in eq II.2, respectively. The bonded potential energy Vb in eq II.3 is composed of Vstretch, Vbend, and Vtor due to stretching, bending, and torsional interactions, respectively. We do not distinguish between protein and solvent energy terms as their corresponding time scales overlap one with the other. In all the tests presented in this paper we have employed a popular united atom model for the proteins8,11 and a simple point charge (SPC) model for water12 modified to include internal degrees of freedom. Details of the implementation of the force field in our calculations can be found in ref 4. In the protein force field, the potential functions for the bond stretchings and angle bendings are harmonic function of bond length and angle bending, respectively. Torsional potentials instead are of two types: the so-called improper torsions and the normal or proper torsions.8 III. Electrostatic Interactions and the Smooth Particle-Mesh Ewald Method The PME algorithm,13-15 as well as similar fast Fourier transform (FFT) based approaches to Ewald summation,16,17 efficiently approximates the reciprocal Ewald sum, thus allowing larger values of the Ewald convergence constant R than are used in conventional Ewald summation, and consequently shorter cutoffs in the direct space Ewald sum which is evaluated in its standard form. The reciprocal space Ewald sum can be written
exp(-π2m2/R2) S(m)S(-m) Vqr ) ∑ 2πVm*0 m2 1
(III.4)
where, given a reciprocal vector m, the structure factor S(m) is defined in terms of the point charges qi and their Cartesian coordinates ri (and finally their fractional coor-
dinates s1i, s2i, s3i) by N
S(m) ) ∑qj exp(2πim‚rj) ) j)1
N
qj exp(2πi(m1s1j + m2s2j + m3s3j)) ∑ j)1
(III.5)
From the above equation, it is clear that the set of structure factors can be looked at as a discrete Fourier transform of the set of charges irregularly placed within the unit cell. In case the charges were placed on a regular grid within the unit cell, the structure factors could be easily computed by FFT. All of the FFT based approaches involve in some sense a smearing of the charges over nearby grid points to produce a regularly gridded charge distribution. The PME method accomplishes this by interpolation. The complex exponential at a point charge position is rewritten as a sum of interpolation coefficients times their values at the nearby grid points. A typical example would be fourth-order interpolation, in which the exponential is expanded over four neighboring grid points in each dimension or 64 neighboring grid points in three dimensions. Rearrangement of the sum appearing in eq III.5 then expresses the structure factor in terms of the discrete Fourier transform of an approximating gridded charge array.13,15 Details of the algorithm for the computation of electrostatic energies and forces, which require a backward FFT, together with a rigorous error analysis and a comparison with regular Ewald summation are given by Petersen.14 The original PME method13 used Lagrangian interpolation for which the expansion coefficients could be interpreted as “charge smearing” weights. In the more recent smooth PME method15 (SPME), which is employed in this study, B-spline interpolation is adopted in alternative. Since B-splines can be analytically differentiated, the forces are the analytical derivative of the corresponding energy and can be obtained somewhat more conveniently than in the original method. In turn, the interpolation B-spline coefficients cannot be directly interpreted as charge smearing as in the Lagrangian interpolation. The reader will find all technical details about the SPME method in ref 15. IV. Integration r-RESPA Algorithm In a previous paper4 we proposed a MD integration scheme based on reversible and symplectic multiple time step (r-RESPA) algorithms developed by Berne and co-workers6,18-20 in the past years and applied most recently to quantum simulations.7,21 This class of algorithms can be easily derived by application to the phase point Γ of an appropriately partitioned and discretized classical propagator, eiLt. Here L is the Liouvillian of the system defined as
iL )
∑ i)1,N
[
∂ ∂V ∂ r3 i ∂ri ∂ri ∂pi
]
(IV.6)
where ri and pi are respectively the coordinates and the momenta of the particles in the system, while r3 i is the corresponding velocity and V is the interaction potential. As previously, we have partitioned the potential of the system appearing in eq IV.6 into four components of increasingly longer time scale, namely
V ) V(0) + V(1) + V(2) + V(3) The four contributions are respectively
(IV.7)
10466 J. Phys. Chem., Vol. 100, No. 24, 1996
Procacci et al.
V(0) ) Vstretch + Vbend
(IV.8)
(1) V(1) + Vvdw + V(1) qd + Vtors
(IV.9)
(2) V(2) ) Vvdw + V(2) qd + Vqr
(IV.10)
(3) V(3) ) Vvdw + V(3) qd
(IV.11)
The above potential subdivision is identical to that used in ref 4. The major difference between the new r-RESPA-SPME algorithm (called R-SPME hereafter) and the old one lies in the use of larger R’s and smaller cutoffs allowed by SPME at a minimal cost. Consequently, all the parameters of the potential subdivision have to be recalibrated to obtain maximum performance. Here, the 0-th reference system includes only the faster bonded interactions. Indeed, the torsional contributions are slowly varying (compared to the other bonded interactions) and can be safely assigned to the order 1 reference system. Also, nonbonded interaction, typically in a range between 0 and 4.35.3 Å, are included in this reference system. To V(2) we assigned the medium-range van der Waals and Coulombic interactions in direct space within a typical range of 4.3-5.3 to 6.9-7.3 Å and the whole reciprocal space term Vqr. Finally, the order 3 reference system which is the most slowly varying contains the remaining direct space interactions from 6.9-7.3 Å to cutoff distance. Since the direct space interactions are typically truncated at Rc ) 10 Å, a number of spherical shells larger than three is unnecessary. The multiple time step propagator can now be easily derived.4 We finally obtain
{ [ { [ { [
eiLt ) exp -
]
(3) ∂(Vvdw + V(3) qd ) ∂ ∆t3 × ∂ri ∂pi 2
exp -
] ]
(2) ∂(Vvdw + V(2) qd + Vqr) ∂ ∆t2 × ∂ri ∂pi 2
+ Vtors) ∂ ∆t1 × ∂ri ∂pi 2 ∂(Vstretch + Vbend) ∂ ∆t0 ∂ exp exp r3 i ∆t0 × ∂ri ∂pi 2 ∂ri ∂(Vstretch + Vbend) ∂ ∆t0 P0 exp × ∂ri ∂pi 2 exp -
{ [
[ [
[
exp exp -
(1) ∂(Vvdw
+
V(1) qd
] [ ] ]}
]} ]}
(1) ∂(Vvdw + V(1) qd + Vtors) ∂ ∆t1 ∂ri ∂pi 2 (2) ∂(Vvdw + V(2) qd + Vqr) ∂ ∆t2 ∂ri ∂pi 2
[
exp -
P1
×
P2
×
]}
(3) ∂(Vvdw + Vqd)(3) ∂ ∆t3 ∂ri ∂pi 2
P3
(IV.12)
where the sums on the 3N degrees of freedom have been left implicit and P0, P1, P2, and P3 are integers. The integration algorithm for such a subdivision uses therefore four time steps, i.e.
∆t0, ∆t1 ) P0∆t0, ∆t2 ) P1∆t1, ∆t3 ) P2∆t2
(IV.13)
The multiple time step R-SPME algorithm described above has been implemented in the molecular dynamics program ORAC.22
V. Results A. Algorithms and Parameters. As shown in the previous sections, the SPME algorithm can reduce considerably the computational weight of the reciprocal space contribution to electrostatic forces. Compared to standard Ewald summation, it becomes now possible to use a larger R parameter in the direct lattice sum and consequently a smaller cutoff radius. In a previous paper,4 in order to obtain an acceptable accuracy with minimal computational effort we used a cutoff in the direct space Rc ) 12.5 Å with R ) 0.18 Å-1. In this study, we use an R in a range between 0.35 and 0.45 Å-1 and cut the direct lattice interaction at 10 Å. The corresponding reciprocal space contribution is computed by SPME with grid point spacing below 1 Å and a fourth-order B-spline. In all the numerical tests presented here the relative accuracy of the electrostatic sum is kept on average of the order of 10-4. The same level of accuracy is required for the van der Waals energy term. This is achieved by cutoffs of 10 Å and larger. We calibrated the R-SPME algorithms with respect to energy conservation and performance by carrying out short simulations on an equilibrated system composed of 7070 atoms which contained a C-phycocyanin protein molecule and 1335 water molecules. The scaling of the algorithms was investigated by carrying out simulations on a solvated reaction center protein composed of 20 627 atoms. More details on the two systems are given in ref 4. Although the proposed algorithms were tuned for best performance and accuracy using a popular united atom force field,8,11 application to other force fields or all atom models is straightforward. We contrast the performance of the R-SPME algorithms with respect to standard Ewald and SPME single time step constrained algorithms called from now on CN-EW and CN-SPME, respectively. We are also interested in comparing our results with those obtained by one of our previous most efficient r-RESPA Ewald algorithms, RESPA4 in ref 4. Moreover, since the overall majority of the MD applications to biomolecular systems use a spherical cutoff we provide timing for spherical cutoff at 10 Å, named CN-CUT. In all single time step algorithms, bonds involving hydrogen atoms were constrained by using SHAKE.23 All the simulations were run for at most 4 ps at 300 K starting from the same equilibrated configuration. As done in past studies,4,24 the accuracy of each integrator with respect to energy conservation is calculated by the ratio R ) ∆E/∆K, where ∆E and ∆K are the standard deviations of the total and kinetic energy, respectively. B. Numerical Results. In Table 1 we present accuracy and performance results obtained from numerical simulations of solvated C-phycocyanin. In this table and in Table 2 each R-SPME algorithm is uniquely identified by eight parameters: the short-, intermediate-, and long-range cutoff in direct space (r1, r2, r3, respectively), ∆t, P0, P1, P2 which identify the four time steps given in eq IV.13, and the value of R. The first four rows in Table 1 show the changes on the energy conservation ratio R due to varying R between 0.35 and 0.35 Å-1 for a given choice of potential subdivision and time steps. The very visible effects on R can be understood by noticing that in the R-SPME scheme the total electrostatic potential acting on a given atom is scattered among different spatial shells. Thus, the time scales of the different contributions are now dependent on R and affect energy conservation. We notice that while at the extreme R-SPME1 (R ) 0.35 Å-1) is not stable yielding a marked drift in the total energy, R-SPME2 to R-SPME4 (R ) 0.40, 0.42 and 0.45 Å-1, respectively) are stable algorithms. Thus, for the potential subdivision used R is at a minimum when R is between 0.40
Method To Simulate Biomolecular Systems
J. Phys. Chem., Vol. 100, No. 24, 1996 10467
TABLE 1: Numerical Tests for the Integration Algorithmsa ∆t
P0 P1 P2
R-SPME1 8.0 8 R-SPME2 8.0 8 R-SPME3 8.0 8 R-SPME4 8.0 8 R-SPME5 12.0 8 R-SPME6 12.0 6 RESPA4 6.0 4 CN-SPME 1.0 1 CN-EW 1.0 1 CN-CUT 1.0 1
2 2 2 2 2 2 3 1 1 1
r1
r2
r3
2 5.0 6.9 10.0 2 5.0 6.9 10.0 2 5.0 6.9 10.0 2 5.0 6.9 10.0 3 5.0 7.3 10.0 4 4.3 6.9 10.0 2 6.5 11.0 13.0 1 10.0 10.0 10.0 1 13.0 13.0 13.0 1 10.0 10.0 10.0
R 0.35 0.40 0.42 0.45 0.40 0.40 0.18 0.30 0.18
R
CPU
0.262b 2.59 0.088 2.44 0.077 2.40 0.112 2.37 0.079 2.32 0.036 2.59 0.015 7.37 0.018 9.65 0.018 19.28 0.018 6.06
a R-SPME’s are the reversible multiple time step algorithms described in section IV of this paper. RESPA4 is our previous most efficient r-RESPA Ewald algorithms;4 CN-EW and CN-SPME are single time step algorithms with standard Ewald and SPME, respectively, for electrostatics; CN-CUT is a single time step algorithm using spherical cutoffs. The numerical tests refer to 4.0 ps of simulation starting from an equilibrated sample of solvated C-phycocyanin at 300 K. The simulation box was nonorthogonal with cell parameters, a ) 72.95 Å, b ) 52.74 Å, c ) 28.36 Å, R ) 100.61°, β ) 82.41°, γ ) 126.3°, and contained a total of 7070 atoms. CPU times were obtained on a DEC-alpha 3000/800s workstation. ∆t for single time step integrators (entry CN-EW, CN-SPME, and CN-CUT) is the integration time step in femtoseconds. For the multiple time step entries ∆t is defined as ∆t ) P0P1P2∆t0 in femtoseconds. R is the direct space convergence parameter. r1, r2, and r3 are the short-, intermediate-, and long-range cutoff in direct space. Finally, R is the energy conservation ratio and the CPU time (CPU) is given in 103 s per ps of simulation. For further details see text. b Drifting algorithm.
TABLE 2: Performance of the Integrators for a Solvated Reaction Center of Rhodobacter sphateroidesa R-SPME5 R-SPME6 RESPA4 CN-SPME CN-EW CN-CUT
∆t
P0
P1
P2
r1
r2
r3
R
CPU
12.0 12.0 6.0 1.0 1.0 1.0
8 6 4 1 1 1
2 2 3 1 1 1
3 4 2 1 1 1
5.0 4.3 6.5 10.0 13.0 10.0
7.3 6.9 11.0 10.0 13.0 10.0
10.0 10.0 13.0 10.0 13.0 10.0
0.40 0.40 0.18 0.30 0.18
7.7 8.6 24.4 28.7 69.1 18.5
a The simulation box was orthogonal of dimensions 75 × 58 × 65 Å3 and contained 8321 protein atoms, 3 Na+ atoms, and 4101 SPC waters, for a total of 20 627 atoms. For further details Table 1.
and at 0.42 Å-1. In other words, with the parameter R in this range, the time scales of the potential derivatives in the various shells become compatible with the selected time steps and the resulting algorithms are accurate and stable. For R larger than 0.42, R worsens due to the increasing short-ranged nature of the reciprocal lattice forces. In algorithms R-SPME5 and R-SPME6 was used an R of 0.40 Å-1 and varied the short and intermediate range cutoffs (r1 and r2, respectively) to obtain the best performance. The conservation ratio R can also be improved if the correction of the intermediate shell reference system is done more frequently. We recall that the reciprocal lattice SPME sum is always included in the intermediate shell in all R-SPME algorithms of Table 1. As shown by R-SPME6, if we change the time step for the intermediate shell from 4.0 to 3.0 fs, an improvement in energy conservation is obtained. Although in Table 1 R-SPME5 is the most efficient algorithm with respect to CPU time, R-SPME6 best combines accuracy and performance. In the latter algorithm we were able, without appreciably affecting R, to reduce the cutoff of the first shell and to increase the time step of the fourth shell. The resulting algorithm has an excellent conservation R with little loss of computational efficiency as compared to R-SPME4 or RSPME5. Remarkably, R-SPME6 is able to run 30 ps per day on a desktop workstation for a system of 7000 atoms computing the exact Coulomb potential and integrating all the degrees of
freedom of the system including the fast stretchings vibrations involving hydrogens with excellent energy conservation. When compared to single time step algorithms R-SPME6 for a 7000 atoms system is found approximately 7-8 and 3-4 times faster than CN-EW and CN-SPME, respectively. Also, with respect to the most effective r-RESPA algorithms using standard Ewald we are able to gain a factor of 2-3 on the performance. Finally, we stress that R-SPME6 is about 2.5 times faster than the method of constraints with spherical cutoff. In Table 2 we show the performance of R-SPME5 and R-SPME6 for a larger system. Although the SPME algorithm scales approximately linearly for systems up to 20 000 atoms,5 the scaling factor S ) 3.3 of the R-SPME algorithms is slightly larger than the expected S ) 20 627/6070 = 2.9. This difference is totally due to the calculation of the Verlet neighbor list used by our MD program which scales as N2. The use of the Verlet neighbor lists becomes the limiting factor as the number of particle grows: For a system of 100 000 atoms simulated with R-SPME approximately half of CPU the time will be spent in the calculation of the Verlet neighbor list. This shortcoming can be easily overcome by using, for instance, linked cells neighbor list methods which scale linearly with the number of particles in the system.25-27 C. Accuracy. So far we have calculated the accuracy of our r-RESPA integrators according to the relative fluctuations of the total with the kinetic energy, i.e., according to R. This criterion is, however, misleading when used to compare r-RESPA with single time step algorithms. Indeed, for single time step integrators R depends quadratically on the step size, ∆t. Thus, a more accurate integrator, i.e., with smaller ∆t, has always a smaller R. In contrast, for multiple time step r-RESPA integrators such as ours the fluctuations of the total energy depend on the potential subdivision and on the time step sizes in a complicated manner. It has been found in the past3,4,24 that r-RESPA invariably computes the thermodynamic and dynamical properties of the system more accurately than simulations carried out with single time step Verlet algorithms of comparable and smaller R’s. A possible explanation of this fortunate behavior lays in the flexibility of the multiple time step scheme which allows the integration of forces of separate time scales with the most appropriate time step and, thus, with similar accuracy. In order to verify the accuracy of our most efficient integrators, we have computed instantaneous energy differences between the “exact” trajectory and two trajectories generated by R-SPME5 (R ) 0.079) and R-SPME6 (R ) 0.036). All the trajectories began from the same phase space point and the “exact” trajectory was generated by using an accurate (R ) 0.012) double time step integrator with time step 0.25 fs for the fast forces (stretching and bending) and 1 fs for all other forces. In Figure 1 we plot the torsional and bending energy differences, ∆UT(t) and ∆UB(t), respectively. These two energies are representative of time scales around 100 and 400 cm-1, respectively. In the same picture we show the corresponding energy difference between a reference trajectory for the CNCUT algorithm obtained using a time step of 0.5 fs (R ) 0.005) and two other trajectories obtained with a time step of 1.0 fs (R ) 0.018) and 1.5 fs (R ) 0.041). For both multiple time step algorithms, we observe in Figure 1 significant deviations from the reference trajectory after 400 fs. After 200 fs more the deviations are on the same order of magnitude as the respective fluctuations. On the other hand, despite an R value 2-4 times smaller the trajectory obtained for the CN-CUT algorithm with a time step of 1.0 fs behaves similarly. Worse is the behavior of the other CN-CUT algorithm
10468 J. Phys. Chem., Vol. 100, No. 24, 1996
Procacci et al. we obtain is of a size that makes us expect that our method will find general application in the field. Acknowledgment. P.P. acknowledges funding from the research network vs. ERBCHRXCT930351, “Molecular Dynamics and Monte Carlo simulations of quantum and classical systems”, of the European Union Human Capital and Mobility Programme. References and Notes
Figure 1. Time record of the energy difference in kJ/mol of the torsional (top) and bending (bottom) contributions in hydrated Cphycocyanin at 300 K with respect to the corresponding “exact” values (see text) for R-SPME (left) and CN-CUT (right) algorithms.
with a larger R which begins to deviate from the reference trajectories well before 400 fs. VI. Conclusions In this paper we have proposed a molecular dynamics method to simulate biomolecular systems with realistic interactions at minimal computational cost. Our technique makes use of the SPME technique to represent the electrostatic forces among the particles. The equations of motion are solved by a r-RESPA algorithm which exploits the time scale separation of the interactions. The combination of the two recently developed techniques has allowed us to reduce considerably the computational cost of simulating infinite biomolecular systems with accurate electrostatic interactions and including all degrees of freedom. Our technique is independent of the force field used to model the nonbonded interactions and can be easily applied, with some minor modifications, to any of the popular biomolecular interaction potentials. With respect to single time step simulations employing standard Ewald sum and rigid bond constraints, we obtain a speed-up of about 1 order of magnitude. When compared to simulations making use of constraints and a spherical cutoff at 10 Å (the radius of truncation of the direct space sum) our molecular dynamics method is about 2.5 times as fast. In view of the fact that the majority of today biomolecular simulations use spherical cutoffs at around 10 Å and larger, the speed-up
(1) de Leeuw, S. W.; Perram, J. W.; Smith, E. R. Proc. R. Soc. London 1980, A373, 27. (2) Watanabe, M.; Karplus, M. J. Phys. Chem. 1995, 99, 5680. (3) Procacci, P.; Berne, B. J. J. Chem. Phys. 1994, 101, 2421. (4) Procacci, P.; Marchi, M. J. Chem. Phys. 1996, 104, 3003. (5) Essmann, U.; Perera, L.; Berkowitz, M. L.; Darden, T.; Lee, H.; Pedersen, L. G. J. Chem. Phys. 1995, 103, 8577. (6) Tuckerman, M. E.; Berne, B.; Martyna, G. J. Chem. Phys. 1992, 97, 1990. (7) Tuckerman, M. E.; Parrinello, M. J. Chem. Phys. 1994, 101, 1302. (8) Brooks, B. R.; Bruccoleri, R. E.; Olafson, B. D.; States, D.; Swaminanthan, S.; Karplus, M. J. Comput. Chem. 1983, 4, 187. (9) Wiener, S. J.; Kollmann, P. A.; Nguyen, D. T.; Case, D. A. J. Comput. Chem. 1986, 7, 230. (10) van Gunsteren, W.; Berendsen, H. J. C. Groningen Molecular Simulation GROMOS Library Manual; Biomos, Groningen, 1987. (11) Parameter and Topology files for CHARMm, Version 22. Polygen Corp., Copyright 1986, Release 1992. (12) Berendsen, H. J. C.; Postma, J. P. M.; van Gunsteren, W. F.; Hermans, J. In Intermolecular Forces; Pullman, B., Ed.; Reidel: Dordrecht, 1981. (13) Darden, T. A.; York, D. M.; Pedersen, L. G. J. Chem. Phys. 1993, 98, 10089. (14) Petersen, H. J. Chem. Phys. 1995, 103, 3668. (15) Essmann, U.; Perera, L.; Berkowitz, M.; Darden, T.; Lee, H.; Pedersen, L. G. J. Chem. Phys. 1995, 103. (16) York, D.; Yang, W. J. Chem. Phys. 1994, 101, 3298. (17) Luty, B. A.; Tironi, I. G.; van Gunsteren, W. G. J. Chem. Phys. 1995, 103, 3014. (18) Tuckerman, M. E.; Berne, B. J.; Rossi, A. J. Chem. Phys. 1990, 94, 1465. (19) Tuckerman, M. E.; Martyna, G. J.; Berne, B. J. J. Chem. Phys. 1991, 94, 6811. (20) Tuckerman, M. E.; Berne, B. J. J. Chem. Phys. 1991, 95, 8362. (21) Tuckerman, M. E.; Parrinello, M. J. Chem. Phys. 1994, 101, 1316. (22) Marchi, M.; Paci, E.; Procacci, P. ORAC a Program to Simulate Complex Molecular Systems with Realistic Interaction. A paper about the program will be submitted to J. Comput. Chem. Shortly after publication, ORAC will become available upon request to the authors. (23) Ryckaert, J.-P.; Ciccotti, G.; Berendsen, H. J. Comput. Phys. 1977, 23, 327. (24) Humphreys, D. D.; Friesner, R. A.; Berne, B. J. J. Phys. Chem. 1994, 98, 6885. (25) Quentrec, C. B. J. Comput. Phys. 1975, 13, 430. (26) Hockney, R. W. Computer Simulation Using Particles; McGrawHill: New York, 1989. (27) Allen, M. P.; Tildesley, D. J. Computer Simulation of Liquids; Oxford University Press: Walton Street, Oxford OX2 6DP, U.K., 1989.
JP960295W
