Complexity and Convergence of Electrostatic and van der Waals Energies within PME and Cutoff Methods

In this paper, we report a detailed comparison between the popularly used cutoff and Particle Mesh Ewald (PME) methods in terms of the time complexity and the energy convergence of the long-range electrostatic and van der Waals interaction calculations. For the comparison, we performed various calculations on various representative biological molecules, including seven peptides and proteins, eleven oligonucleotides, and three conformations of a nucleotide-sugar. The results provide useful insights into the appropriate choice of the methods (i.e. the cutoff or PME) and that of the cutoff values for the calculations on different kinds of molecules. It has also been demonstrated that for some cases using different cutoff values for calculating the electrostatic and van der Waals interaction energies will be computationally more efficient.


Introduction
Long-range electrostatic and van der Waals interactions have for a long time been in the center of attention in empirical methods of computational chemistry, mainly within molecular mechanics and dynamics.Particularly, electrostatic interactions belong to permanent pitfalls of computational chemistry [1][2][3][4].There are principally two reasons for that: (i) The electrostatic and van der Waals interactions represent the most expensive part of calculations, and the evaluation of them is so time consuming that it has been necessary to use approximations; (ii) The results of the calculations are sometimes very sensitive on the approximation used.
Until recently, the most frequently used method to handle electrostatic and van der Waals interactions was to ignore all interactions between atoms whose internuclear distance is longer than a certain cutoff value.Such an approach is usually called the Cutoff Method (CM) [5][6][7].There are several possible choices concerning how the cutoff can be used.We have chosen the most straightforward definition as a distance border limit.
Since the calculation of van der Waals interactions is more expensive than that of electrostatic interactions, the CM approach can be slightly modified by using two cutoffs, one for van der Waals and the other for electrostatics [8].These two cutoffs may be set up in such a way, that the final accuracy of results is the same for both components.In other words, unnecessarily more precise (and expensive) calculations of one of these two types of interactions are avoided.
As mentioned above, effects of PME method on geometrical (conformational) features of the system have several times been subjected to studies and reported in the literature.However, we have not found any paper where PME and CM energies would be subjected to such a comparison.Therefore, our effort in this work is to compare efficiency of using two cutoff values in calculations and, especially, to elucidate the speed of convergence of CM and PME energy to its final value.
Several systems have been selected for the study, representing various classes of important biomolecules.

Methods
Altogether, a total of 19 representative molecules were considered in our various testing calculations, for one of them we considered three different conformations.Those include seven peptides and proteins, eleven oligonucleotides, and three conformations of a nucleotide-sugar.More detailed structural information is shown in Table 1.
Unless stated otherwise, all calculations were performed by using the SANDER and ANAL of AMBER 5.0 suite of programs.The initial coordinates were taken from structural databases and/or our previous calculations.First, all molecules were minimized using the distance dependent dielectric constant and the following parameters: scnb=2.0,scee=2.0,dx0=0.01,dxm=0.5, dele=0.001,drms=0.3(remaining parameters were assigned by default values).Then, energies were calculated for various cutoff values using both the PME and cutoff methods.
To examine the dependence of the calculated energy on the used cutoff value, the energy calculations were performed using different cutoff values; the used cutoff value was increased with a step of 1 Å until all nonbonding interactions were included (denoted by the final cutoff below).For convenience, the energy calculated at the final cutoff is called the final energy.For the energy calculation using PME, a rectangular box was created of such a size, that the calculated molecule fits the box with positive tolerance of 3 Å.We then started with a cutoff of 1 Å and subsequently increased it by 1 Å step until the value of B/2 was reached (B is the length of the shortest edge of the rectangular box).This approach was used for both explicitly and implicitly solvated systems.We have started with such a small value of cutoff since our goal was to obtain a complete image.
The calculations with explicit solvent were performed with the ANAL program that allows for separation of solute-solute, solute-solvent, and solvent-solvent interactions.The coordinates of the atoms used in these calculations were taken from equilibrated molecular dynamics (MD) snapshots.

Complexity of electrostatics and van der Waals energy calculation
In this section we only discuss the complexity of the interaction energy calculations within a distance which is shorter than the cutoff.This part of calculations is the same for both cutoff and PME method.We can roughly estimate the complexity of the energy calculations by the following theoretical considerations.
Electrostatic energy is calculated using Coulomb's law -eq.( 1), i.e. one multiplication and one division.

∑∑
This is calculated for n-times and (n-1)-times for the first and the second summation, respectively.For the overall time complexity T el (n) we can write: ) where c mult and c div are complexities of multiplication and division, respectively.We consider r ij as a constant since it is necessary for calculation both electrostatic and van der Waals energy.
The simplest approximation of van der Waals energy is Lennard-Jones potential -eq. ( 3): It is clear, that we can not calculate term A − using less than 4 multiplications, 2 divisions, and 1 substraction.For the overall time complexity T vdw (n) we then have: where c sub is complexity of subtraction.
For current high performance computers division is more time-consuming than multiplication, that is slightly slower than addition and subtraction (or it is of equal speed).It is then clear from equations ( 2) and ( 4) that calculation of E vdw is more time-consuming than that of E el in all cases.
We have written a program in C language to measure and print the ratio of time complexities for the calculations of E vdw and E el (see Table 2).The results listed in Table 2 agree well with the above theoretical considerations.The ratio of time complexities depends on the used hardware, but we can conclude that the calculation of E vdw is roughly about 3 times slower than that of the corresponding E el .

Comparison in time complexity between the cutoff and PME methods
The time complexities of the cutoff and PME energy calculation methods are shown in Table 3.The data in Table 3 reveal the following remarkable trends: (i) the PME method is in all cases slower than the cutoff method; (ii) for large molecules the PME is about 30% slower.

Peptides and proteins
Seven molecules belonging to this class were studied (see Table 4).The results of energy calculations for the protein P7 are depicted in   E -E f )/E f , where E is the energy value calculated for the cutoff and E f is the final energy (see method section).The convergence of E vdw is slower for larger molecules, i.e., they require to use a larger cutoff value.It applies for both the cutoff and PME methods.However, for PME the convergence is considerably slower.The convergence of electrostatic energy is almost independent of the size of the molecule for the considered peptides and proteins (data not shown).In terms of the energy convergence, the cutoff method seems to be better than PME for peptides and proteins.
Fig 3 shows which minimum cutoff value must be used to have the energy error below 1% for each case.Fig's 1 and 3 demonstrate that the convergence of E vdw is markedly slower than that of E el for both the cutoff and PME methods.A larger cutoff value is necessary for E vdw to reach the same accuracy than for E el .The larger the studied molecule, the greater the difference between the cutoff values required for E vdw and E el reaching the same accuracy of 1%.It means that the computations using different cutoff values for E vdw and E el will be more efficient in this case.A quite small cutoff value could be used to calculate E el .However, this will not lead to a dramatic speed-up of the calculation since E el is relatively inexpensive to calculate.See Table 5 for a brief summary of these results.

Cutoff method is slightly better
The results obtained for peptides and proteins in explicit solvent are not dramatically different from what we discussed above.A notable difference is seen for the convergence of E el with the cutoff method where the curve does not oscillate (see Fig 4), whereas it does in implicit solvent.This is attributed to the coupling of the changes of the solute-solute and solute-solvent contributions to the E el value with the change of the corresponding solvent-solvent contribution.When only the solute-solute and solute-solvent contributions are included in the E el value, the curve becomes similar to that in implicit solvent (see Fig 5).Table 6.Final energies (kcal/mol) calculated using cutoff E f (cutoff) and PME E f (PME) methods.The relative difference of absolute energies δE f (%) was calculated using equation:

Oligonucleotides
Eleven molecules belonging to this class were studied (see Table 7).Note that the absolute values of E f vdw in the explicit solvent are similar for the cutoff and PME methods but the values of E f el calculated by using PME are markedly larger than the corresponding values calculated by using the cutoff method (see Table 6).This difference is caused by the absence of explicit water in the grid for PME calculation.The grid includes charges, which simulate long-range interactions in the molecule.If explicit water molecules are present, they shield grid charges and PME provides absolute energies comparable with energies from the cutoff method (see Table 6).The larger the water box is, the closer to each other the values of E f el for the two methods are.The results of energy calculations for the oligonucleotide O11 are depicted in Fig 6, which shows the convergence of energy to its final value.The general tendencies seen in the results calculated for all of the remaining oligonucleotides are nearly the same (data not shown).The speeds of the convergence of E vdw for both the cutoff and PME methods and that of E el for the cutoff method are somewhat dependent on the size, the number of ions, and the 3D structure of the system (see Fig 7).The convergence of E el for PME method is independent of the size of the molecule.The convergence of E el is dramatically slower than that of E vdw for the cutoff method.It means that using different cutoff values for calculating E el and E vdw should be very efficient in this case; using the two different cutoff values should significantly speed up the computation since E vdw is much more expensive to calculate than the corresponding E el .The convergence of E vdw for PME method is slightly slower than that of E el , but the use of different cutoff values for E el and E vdw would not dramatically affect the speed of computation in this case.
Figure 7 clearly shows that PME is much better than the cutoff method for oligonucleotides.From the figure, one can also clearly see why it was difficult to obtain stable molecular dynamics trajectories on this class of molecules using the cutoff approach to handle electrostatics.One would need to use a very large cutoff value in order to achieve the convergence of the energy.See Table 8 for a brief summary of these results.The results obtained for oligonucleotides in explicit solvent are similar to those discussed above.As seen for proteins and peptides, for oligonucleotides the notable difference also exists in the curve of the convergence of E el for the cutoff method (see Fig 8).The reason for this difference is the same as that implicit s olvent explicit s olvent for proteins, i.e. due to the coupling with the change of the additional solvent-solvent contribution to the E el value.When the solvent-solvent contribution is excluded from the E el calculation, the curve becomes similar to that in implicit solvent (see Fig 9).
Figure 9. Dependency of electrostatic energy (calculated using cutoff method) on cutoff for implicit solvent and explicit solvent (without including solvent-solvent interactions) for oligonucleotide O4.
The same as observed for the absolute energies of the peptides and proteins (see last paragraph of the "Proteins and peptides" section) apply also for the absolute energies of the studied oligonucleotides.
We have selected systems with different space locations of the ions, so the overall shape of the molecule itself as well as of the entire system was different.The basic geometrical parameters are shown in Table 9. * Distance between the aromatic moiety and the base.** Overall size of the system calculated as a sum of distances of all atoms in the system from the atom P5 (the P atom, which is nearer to the base).
The energy convergence for all these systems is shown in Fig's 10 and 11.As shown in the figures, the convergence is clearly dependent on the overall shape of the system.The only exception is for van der Waals interactions the speed of the convergence depends also on the size of molecule.The overall results indicate that for these molecules E vdw converges slowly and E el very slowly for the cutoff method.The difference between the convergence of E el and that of E vdw for the nucleotide-sugar is somewhat between the differences for peptides/proteins and for oligonucleotides.The convergence of E el within the cutoff method is slightly slower than that of E vdw .It means that using different cutoff values for E el and E vdw would be efficient for the computation, but the speed-up of the computation would not be dramatic as the difference between the two cutoff values is small.
Within PME method, E vdw converges slowly and E el very quickly, which is similar to the case of oligonucleotides.The convergence of E vdw is markedly slower than that of E el , indicating that the use of two different cutoff values for E el and E vdw would be efficient for the computation, but the speed-up   implicit solvent explicit solvent - Evdw Eel cutoff PME of the computation would be insignificant, because the computation of E el is relatively faster than that of E vdw .Fig 12 clearly suggests that for nucleotide-sugars, in terms of the energy convergence, the performance of PME method is generally better than the cutoff method.However, when a very highaccuracy is desired for the energy evaluation, the cutoff method can also be used and may even be more efficient.Whereas the cutoff method requires a large cutoff value for E el , PME needs a similarly large cutoff value for E vdw whose calculation is more expensive than the E el calculation.See Table 10 for a brief summary of these results.

PME method is better
The results obtained for the nucleotide-sugar in explicit solvent well correspond with those obtained for peptides and oligonucleotides (see Fig 's 13 and 14).The same applies for the cutoff and PME absolute energies (data not shown).implicit solvent explicit

Possibilities to improve the convergence performance by changing parameters
The above results were obtained using the version 5.0 of AMBER.There are two problematic issues that can be improved either by changing PME parameter values for the calculations or by using a more recent version of AMBER (6.0 or 7.0) which uses more parameters to control the computation and also different default PME parameter values.The first issue is associated with the observation that the PME energies calculated with very small cutoff values are dramatically different from the corresponding final energies converged at the large cutoff values.However, this problem can easily be solved by increasing the spline order and/or reducing the grid spacing.We observed a much improved convergence performance by just reducing the grid spacing (data not shown).On the other hand, actually one does not need to worry about the dramatically different energies calculated by using very small cutoff values because no one would really use such a small cutoff value in practical calculations.
The second problem, which looks more serious, is associated with the observation that the convergence of vdW energy is different for the calculations with the cutoff and PME methods.In principle, the convergence of vdW energy should be in both cases exactly the same because vdW evaluation is not dependent on using PME.As far as we have known [33], the PME-based calculations implemented in AMBER is associated with an "atom based pairlist", in contrary to the "(neutral charge-) group based pairlist" used in the calculations without using PME in AMBER 5.0 or the earlier versions.It is this difference that causes the notable discrepancy between the vdW energy convergencies.However, when using AMBER 6.0 or 7.0, the additional homogeneous long-range correction to the vdW may optionally be enabled (VDWMETH=1), which includes a bulk density vdW correction.This setup improves the vdW convergence considerably, as shown in Fig. 15. .Dependency of energy error on cutoff (for peptide P2) as calculated by AMBER 5.0 and 6.0 with different setup.Denotation: Series I: AMBER 6.0 -cutoff method (USE_PME=0, EEDMETH=4); Series II: Amber 6.0 -PME, E vdw without correction (USE_PME=1, EEDMETH=1, VDWMETH=0); Series III: Amber 6.0 -PME, Evdw with correction (USE_PME=1, EEDMETH=1, VDWMETH=1); Series IV: Amber 5.0 -cutoff method; Series V: Amber 5.0 -PME method.
A detailed comparison of the results indeed shows the same vdW energy convergences for the PME and cutoff methods implemented in AMBER 6.0 and the PME method implemented in AMBER 5.0.It appears that the calculation with the cutoff method by using AMBER 5.0 without the additional vdW energy correction converges faster than the corresponding calculation by using AMBER 6/7 with the additional vdW energy correction.For PME, on the other hand, AMBER 6/7 versions are better tuned than AMBER 5.0.
Finally, it should be mentioned that we only examined single-point energies in discussion of the convergence of the energy calculations.Obviously, these energetic results depend on conformation.To exclude this drawback, we were interested rather in relative energies, or, better to say, in energy tendencies.Our conclusions are supported by the fact that for all molecules within one general type identical trends were found.We also note here that this kind of inspection is more detailed that, for example, a molecular dynamics (MD) simulation would be.In the case of MD simulation many energy calculations must be performed and the errors may mutually be canceled, so the final observation may somewhat be deformed.Therefore, we believe that our approach to compare cutoff and PME is reasonable.

Conclusions
Long-range electrostatic and van der Waals interactions still belong to the most problematic part of energy calculation within empirical force field approach of computational chemistry, mainly within molecular mechanics and dynamics.Due to time complexity it is usually not possible to include all interactions into calculations and some approximations must be used.Our effort with this paper has been to compare two most frequently used approximations to handle long-range electrostatic and van der Waals interactions, cutoff, and PME methods.Altogether, 19 molecules have been calculated, for one of them we considered three different conformations.These calculated molecules include peptides, proteins, oligonucleotides, and nucleotide-sugar.The main program package used is AMBER, versions 5.0, 6.0, and 7.0.Some results were obtained also by software written by us.
First, we have compared the time complexities of the calculations of the electrostatic and van der Waals energy components.We have shown by theoretical considerations and also by practical computations that the van der Waals energy component is roughly about 3 times more expensive to calculate than the corresponding electrostatic component.Thus, it is computationally more efficient to use different cutoff values for calculating the electrostatic and van der Waals interactions, particularly for the cases when the van der Waals component requires a lower cutoff value.We have also compared the overall speeds of the calculations with the cutoff and PME methods and concluded that the calculation with the PME method is generally slower, and in some cases is substantially slower, which is consistent with previous results in literature.
Then, we have examined the electrostatic and van der Waals energy convergence to their final values for different molecules, for different cutoff values, and for different methods.This has been performed for both the implicitly and explicitly solvated systems; the implicit solvent is accounted for through the use of a distance-dependent dielectric constant.The borderline cutoff value is the maximum cutoff value including all interactions for the cutoff method or one half of the shortest edge of the PME box for the PME method.The qualitative results are basically the same for both the implicitly and explicitly solvated systems.We have also compared the absolute electrostatic energies calculated by using the cutoff and PME methods.While these values are very different in the implicitly solvated systems, they become almost identical in the explicitly solvated systems.The reason is that the charges in PME grid are not shielded by water molecules in the implicit solvent and, therefore, create an artificial energy term that is not used with the explicit solvent.In general, the energy convergence is found to be dependent on the shape of the system.It has also been concluded that van der Waals energy component convergence is strongly dependent on size of the system measured by the number of atoms.The larger the system, the slower the energy convergence.
It has been observed that the energy convergence is dependent on the type of molecule.Specifically, it is different for peptides, oligonucleotides, and the nucleotide-sugar conformations.For example, with both the PME and cutoff methods peptides and proteins exhibit a slower convergence of van der Waals energy than that of E el .For E vdw , the cutoff method exhibits the faster convergence than PME and, therefore, using the cutoff method, rather than the PME, is preferred for peptides and proteins, at least from the energy point of view.Because E el is relatively inexpensive to calculate as compared to E vdw , the use of different cutoff values for calculating different energy components is insignificant for peptides and proteins.
The situation is quite different for oligonucleotides, where the convergence of E el is very slow with the cutoff method and, therefore, we recommend not to use the cutoff method for this class of compounds unless the systems are very small and the used cutoff value can cover the entire system.This is in a good agreement with the literature data leading to a recommendation for using PME exclusively on oligonucleotides.Using the cutoff method on oligonucleotides would require very careful tuning of the cutoff value.On the other hand, if the cutoff method is used, employing two cutoff values (one for E vdw and one for E el ) could significantly save the time of computation.
The situation for the nucleotide-sugar is somewhat in between those for peptides and proteins and for oligonucleotides.The results show that both the PME and cutoff methods may be used for this class of compounds.However, when the cutoff method is used, the cutoff values should carefully be tuned.In this case, using different cutoff values for E vdw and E el could also save some time of computation.But using the PME appears safer for oligonucleotides.
Fig 1, which describes convergence of the energies to their final values.The results calculated for all remaining peptides and proteins are similar to what depicted in Fig 1.The only difference is that the speed of convergence of E vdw (for both the cutoff and PME methods) depends on the size of molecule, as illustrated in Fig 2.

Figure 1 .
Figure 1.Dependency of energy error on cutoff (for protein P7).Energy error (EE) is calculated using the equation: EE = 100.(E-E f)/E f , where E is the energy value calculated for the cutoff and E f is the final energy (see method section).

Figure 2 .
Figure 2. Dependency of Evdw error on cutoff (for both cutoff (left) and PME (right) methods for all studied peptides and proteins, the darker color of the curve, the larger the molecule is).For description of energy error see Figure 1.

Figure 3 .
Figure 3. Smallest values of cutoff for which energy error is smaller then 1% (for peptides and proteins).

Figure 4 .
Figure 4. Dependency of energy error on cutoff for peptide P2 in explicit solvent.For description of energy error see Figure 1.

Figure 5 .
Figure5.Dependency of electrostatic energy (calculated using cutoff method) on cutoff for implicit solvent (default value for dielectric multiplicative constant of 1.0) and explicit solvent (without including solvent-solvent interactions) for peptide P2.

Figure 7 .
Figure 7. Smallest values of cutoff for which energy error is smaller then 1% (for oligonucleotides).

Figure 8 .
Figure 8. Dependency of energy error on cutoff for oligonucleotide O4 in explicit solvent.For description of energy error see Figure 1.

Figure 10 .
Figure 10.Dependency of energy error on cutoff for conformers of nucleotide-sugar.Energy error is

Figure 11 .
Figure 11.Dependency of energy error on cutoff for conformers of nucleotide-sugar.Energy error is calculated for E vdw using cutoff method (Part a)) and PME method (Part b)).

Figure 13 .
Figure 13.Dependency of energy error on cutoff for nucleotide-sugar N3 in explicit solvent

Figure 14 .
Figure14.Dependency of electrostatic energy (calculated using cutoff method) on cutoff for implicit solvent and explicit solvent (without including solvent-solvent interactions) for nucleotide sugar N3.

Table 1 .
General overview of studied molecules.For more details see Tables4 and 7in Results section.

Table 2 .
Ratio of computational complexities of van der Waals (T vdw (n)) and electrostatic (T el (n)) energies calculation.

Table 3 .
Comparison of computational complexity of cutoff and PME method.Calculations were executed on a Sillicon Graphics workstation with the operating system Irix 6.1.

Table 4 .
Studied peptides and proteins.

Table 5 .
Summary of results for peptides and proteins.

Table 8 .
Summary of results for oligonucleotides.

Table 9 .
Selected geometrical parameters of three different conformations of the nucleotide-sugar.

Table 10 .
Summary of results for nucleotide-sugars.