6.3 Algorithms for Generating Non-bonded Interactions

There are four algorithms used in calculating the nonbonded energies, each making different approximations in an attempt to speed the calculation. Electrostatic interactions are the most difficult to deal with for two reasons. They do not fall off quickly with distance (so it is inappropriate to simply ignore all interactions beyond some cutoff), and they depend on odd powers of r necessitating expensive square root calculations for each pair evaluated. The approximations used to make the electrostatics calculation more tractable are setting the dielectric constant equal to r or using a constant dielectric but only calculating distant interactions periodically (and storing the value in between).

Setting the dielectric constant equal to the atom atom distance times a constant factor (determined by the EPS keyword value) makes the computation easier by eliminating the need to calculate square roots and by making the calculated contribution fall off more quickly. It also introduces problems. The force calculated using an r dependent dielectric will be larger than the force from a constant dielectric at short distances (5.0 angstroms or less by comparsion to a constant dielectric of 2.5). In addition, the electrostatic contribution still falls off relatively slowly and large distance cutoffs are needed. As the number of atom pairs included will be proportional to the cutoff cubed, this is a significant disadvantage.

The RESI options use an r dependent dielectric (as do all of Bruce Gelin's protein calculations). In the RESI searches, the program constructs a rectangular box around every residue, and examines every residue–residue pair that can yield atom atom pairs within the cutoff. Every atom–atom pair within those residue–residue pairs is added to the non-bonded list. When switching functions are used, RESI and ATOM will yield the same energy, but RESI will generate a longer list.

The SHFT option is similar to ATOM except, the potential:

E= (QI*QJ/EPS)*( 1.0/R**2 + R**4/(2.0*CTOFNB**6) - 1.5/CTOFNB**2 ) is used when ( R < CTOFNB ) and zero otherwise. This potential and it first derivative approach zero as R becomes CTOFNB, without the messy computation of switching functions and steep forces at large R.

The EXEL and EXFL routines use a constant dielectric at large distances joined to the short distance function. The long distance terms are approximated using interpolation and periodic updating. At CUTNB the constant dielectric potential is joined to a dielectric equal to r potential with a shift and scale that preserve the continuity of the energy and forces. CONS uses a constant dielectric everywhere. This requires a square root to be calculated in the inner loop of ENBOND, slowing things down a bit, but it is physically more reasonable and widely employed by other groups doing empirical energy modelling (ex. ST2 water).

This form allows a small CUTNB (5.0 angstroms with EPS=2.5) even though the electrostatic terms are still varying rapidly at that distance. The short range forces are identical to those calculated with the other options, reflecting the decrease in dielectric shielding at short ranges. For the long range forces there is effectively no cutoff in the electrostatic energy. Interactions between atoms beyond the CUTNB are calculated when the list is updated and stored together with their first (EXFL,CONS) or first and second (EXEL) derivatives. The energy is calculated by explicitly evaluating pairs in the list and using the stored potentials, fields, and gradients to approximate the distant pairs. In essence the routines assume that for distant pairs the atom movements will be small enough that the changes in their electrostatic interactions can be accurately calculated using local expansions.

The constant dielectric routines compile the close contact list using the same two stage minimum rectangle box search that is described above. In this way the efficiency of a residue by residue search is exploited while being certain that all necessary pairs are included. For close residue pairs an atom by atom search is then performed. Atom pairs are either included in the list of close contacts or their electrostatic interactions are calculated and stored.

An option is offered to increase the accuracy of residue residue interactions by using a multipole expansion of one residue evaluated for each atom of the other. This cutoff for this treatment is CUTMP. For residue pairs outside of CUTMP only a single multipole evaluation is made and second order polynomial expansion is used to extrapolate to each atom. Ordinarily this is sufficient and CUTMP is set to 0.0.

The DUMMY option results in no non-bonded interactions whatsoever. The creation of a dummy list can be done very quickly. This is useful when one is analyzing only the other terms in the energy.