Nonbond Cutoffs

Contents

Introduction

The energy expression of Eq. 2-6, which is representative of current forcefields, is computationally tractable only for systems with relatively small numbers of atoms. The number of internal coordinates grows linearly with the size of a molecule, so the computational work involved in the first nine terms in Eq. 2-6 also grows linearly. However, inspection of the final summation, which represents the nonbond interactions, reveals a quadratic dependence on the number of atoms in the system. If the system of interest has 1000 atoms, the nonbond summation has about 500,000 terms. If it has 10,000 atoms, the summation has 50,000,000! Therefore, it is common to introduce a nonbond cutoff, i.e., to neglect the nonbond interactions for pairs of atoms separated by distances greater than a cutoff value. Clearly, this has the potential for being a rather drastic simplification. A review of several cutoff methods was published by Brooks et al. (1985).

To appreciate the impact of cutoffs on computational efficiency, consider a protein-substrate-solvent system with 5000 total atoms. A system of this size would be typical of a small protein (100-150 residues) surrounded by 1-2 layers of water.

Figure 2-6 shows how the number of nonbond interactions increases with the cutoff distance. This calculation would run at least 10 times faster with an 8.0 Å cutoff than with no cutoff (assuming that the nonbond term is rate limiting, which it usually is). The trade-off is, of course, that interactions beyond the cutoff distance are not accounted for. Both van der Waals interactions and electrostatic interactions are significant up to 15 Å or more. For example, in a recent calculation of the energy as a function of cutoff distance in the [Ala-Pro-D-Phe]2 crystal, Kitson and Hagler showed that the nonbond energy accounted for changes from 63 to 97% of the asymptotic value as the cutoff distance was increased from 8 to 15 Å (Kitson and Hagler 1988).

Figure 2-7 shows how the van der Waals component of the nonbond energy varies as a function of cutoff distance for the [Ala-Pro-D-Phe]2 crystal. The van der Waals energy changes by 40% when the cutoff distance is increased from 8 to 15 Å. The exact dependence of the energy on the cutoff distance depends on the system itself and should be calibrated for each new system.

Nonbond Cutoff Algorithm

Simply cutting the nonbond interactions off at a given distance leads to discontinuities in the energy and its derivatives. Thus, the Discover program uses a quintic spline function to smoothly turn off the interactions over a range of distance in order to mitigate these effects. This switching function is illustrated in Figure 2-8, which shows the unmodified van der Waals potential as a black curve; the switching function S(r) as a black striped curve; and the resulting, switched potential in gray. The switching function is defined by two variables: the point where the function reaches zero and the range over which the function decreases from one to zero. (The variable names differ between Discover versions 2.9 and 95.0, as illustrated in the figure.) The Discover program uses a fifth-order polynomial for the switching function. It is arranged so that the first and second derivatives are zero at both the inner and outer range of the switching region. Thus, the interaction energy and its first and second derivatives are continuous, although higher derivatives are not.

For efficiency, the Discover program creates a neighbor list containing all pairs of atoms to be used during calculation of the nonbond interactions. Since this list is not updated at every step of the calculation, it includes atoms in a buffer region that might move close enough to contribute before the next update of the neighbor list. In the Discover 2.9 program the buffer region is defined by the difference between CUTDIS and CUTOFF, and in the Discover 95.0 program the buffer width is explicitly specified, along with the cutoff. To ensure that no atoms outside the buffer region can move close enough to interact, the list is automatically updated whenever any atom moves more than one-half the buffer width. Thus, the width of this buffer region, coupled with the velocity with which atoms move, determine the maximum time before the neighbor list is updated.

Charge Groups

The use of cutoffs with the van der Waals interaction potential is quite reasonable. The potential is relatively short range and dies out as 1/r6. By 8-10 Å the energy and forces are quite small. Thus, switching the van der Waals potential to zero at about 10 Å is a reasonable approximation. The Coulombic interactions, on the other hand, die off as 1/r, so even at considerable distances the energy of interaction is not negligible. However, except for a few formally charged groups, most molecules are composed of neutral fragments with dipoles and quadrupoles. In reality, then, the leading term in the electrostatic interaction between molecules or parts of molecules is a dipole-dipole interaction, which falls off as 1/r3.

To understand the implications, note that the interaction of two monopoles, each of one electronic unit of charge, is about 33 kcal mol-1 at 10 Å, while the interaction energy of two dipoles formed from unit monopoles is no more than about 0.3 kcal mol-1. It is clear that ignoring the monopole-monopole interaction would be grossly misleading, whereas ignoring the dipole-dipole interaction would only be a modest approximation. If the nonbond cutoffs were applied on an atom-by-atom basis, they could artificially split dipoles by having one atom inside the cutoff and one atom outside. Instead of ignoring a relatively small dipole-dipole interaction, this would have the effect of artificially introducing a large, monopole-monopole interaction. To avoid these artifacts, the Discover program applies cutoffs over charge groups.

A charge group is a small group of atoms close to each other which have a net charge of zero or almost zero. In many cases, charge groups are identical to common chemical functional groups. Thus, a carbonyl group, methyl group, or carboxylic group would also be an approximately neutral charge group. The Discover program designates one atom from each charge group as the switching atom. The Discover program generates the neighbor list by considering the distance between the switching atoms of two groups. If the distance is less than the cutoff distance, then all the atom pairs in the two groups are included. If the distance is greater than the cutoff, they are all excluded. Similarly, when calculating the actual interaction energy, the Discover program switches off the interaction between all atom-atom pairs in the two groups based only on the distance between the two switching atoms. This procedure prevents artifacts due to splitting of dipoles.

The size of a charge group, as defined by the greatest distance from the switching atom to another atom in the group, must be significantly smaller than the cutoff distances. Otherwise, the interaction between two atoms close to each other might be ignored because the switching atoms of the two groups are farther apart than the cutoffs. Typical groups are no more than 1-3 Å large, so cutoffs larger than 7-8 Å are reasonable. However, occasionally molecules contain considerably larger groups. The Discover program checks the size of the groups against the cutoff distances, outputs an error message, and terminates if it decides the cutoffs are too short compared to the group size. In this case, you must either increase the cutoffs or redefine the groups to be smaller. The Discover program also warns you about significantly non-neutral groups. Some of these can be expected if the molecule contains formally charged functional groups, such as protonated amines and carboxylates. However, other non-neutral groups usually indicate an error in group definitions.

Double Cutoffs

Double cutoffs are available in the Discover 2.9.5 program only.

The Discover program also incorporates an improvement over a single cutoff distance called double cutoffs, or, as it is sometimes called in dynamics, multiple timesteps. The nonbond interaction potential at a distance is a smooth function that does not vary rapidly. With double cutoffs, two cutoff distances--an inner and outer one--are assigned. The two distances define an inner spherical region and an outer shell around a given atom. Interactions with atoms in the inner region are treated as usual, but interactions with the atoms in the outer shell are updated only at every neighbor list update. It is assumed that their contribution varies only slowly. These double cutoffs are used to make the calculation less expensive by reducing the inner cutoff to a smaller value than could normally be used with a single cutoff. Accuracy is regained at minimal cost by using a large distance for the outer cutoff.

It is important to realize that the effective potential energy surface is not quite continuous when double cutoffs are used. The magnitude of the discontinuities depends on the cutoff distances and the system that is being studied. These discontinuities are only a minor problem for dynamics, where they are manifested as small fluctuations in the total energy. Their effect during minimization depends on the minimizer that is used, because some minimization algorithms, such as conjugate gradients, are quite sensitive to discontinuities in the surface. Other algorithms, such as steepest descents, are relatively robust.

Main access page Theory/Methodology access.

Describe-System access Periodic Boundary Conditions Cell Multipole Method

Copyright Biosym/MSI