Generation of Non-bonded Interactions Nonbonded interactions (frequently abreviated "nbond") refer to van der Waals terms and the electrostatic terms between all atom pairs that are not specifically excluded from nonbond calculations as for example are directly bonded atoms *note nbx: (struct.doc)nbx. These terms are defined on atom pairs and to a first aproximation would require the number of atoms squared calculations. To avoid this burden various truncation and approximation schemes can be employed in the program, breaking the nonbonded calculation into two parts, initialization and actual energy calculation. The method of approximation, cutoffs, and other relevant parameters can be entered any time the nbond specification parser is invoked. See the syntax section for a list of all commands that invoke this parser. * Menu: * Syntax:: Syntax of the nonbond specification * Defaults:: Defaults used in the nonbond specification * Function:: Description of the options * Tables:: Using nonbond lookup tables in place of analytic potential energy functions * Cubes:: Alternative way to compute the nonbonded ontact list * Cluster:: Cube-Cluster nonbonded list generation method

[SYNTAX NBONDs] { NBONds } { [INBFrq integer] nonbond-spec } { UPDAte ... } { } { ENERgy ... } { } { MINImize ... } { } { DYNAmics ... } { } { VIBRAN ... } { } { CORRel ... } { } NOTE: The INBFrq value is remembered. If its value is zero, no interpretation of [nonbond-spec] will be made, as well as no modifications of the nonbond lists. It's default value is -1 . In all cases as many keywords and values as desired may be specified. The keywords are: nonbond-spec::= [method-spec] [distances-spec] [misc-specs] [INIT] [RESET] method-spec::= [ ELEC electrostatics-spec ] [ VDW vdw-spec ] [ BYCUbes ] [ NOELectrostatics ] [ NOVDwaals ] [ BYGRoup ] [ GRAPe ] [ LIST ] [ BYCBim ] [ NOGRape ] [ NOLIst ] [ BYCC ] electrostatics-spec::= [ EWALD ewald-spec ] [ FMA fma-spec ] elec-opt-spec [ NOEWald ] [ NOFMA ] [ ACE ace-spec ] elec-opt-spec::= [ ATOM ] [ CDIElec ] [ SHIFted ] [ GROUp [ EXTEnded [GRADients] [QUADrip] ] ] [ RDIElec ] [ SWITched ] [ [ [NOGRad ] [NOQUads] ] ] [ FSWItch ] [ [ NOEXtended ] ] [ FSHIft ] [ MSHIft ] ! MMFF [ TRUNcate ] ! MMFF vdw-spec::= [ VGROUP [ VSWITched ] ] [ ] [ VATOM [ VSHIfted ] ] [ [ VSWItched ] ] [ [ VFSWitch ] ] [ [ VTRUnc [ CTVT real ] ] ! MMFF only distances-spec::= [general-dist] [warning-dist] general-dist::= [ CUTNB real ] [ CTONNB real ] [ CTOFNB real ] [ CTEXNB real ] SOFT [EMAX real ] [MINE real] [MAXE real] VDWE ELEE ! SOFTVDW warning-dist::= [ WMIN real ] [ WRNMXD real ] misc-specs::= [ EPS real ] [ E14Factor real ] [ NBXM integer] - [ NBSCale real] [IMSCale real] [EXOForce] [ NORXN ] [ RXNFLD rxnfld-spec ] [ RXNNB rxnfld-spec ] rxnfld-spec::= [ EPSEXT real ] [ ORDER integer ] [ SHELL real ] fma-spec::= [LEVEL int] [TERMS int] ewald-spec::= (see *Note ewald::(ewald.doc)syntax ) ace-spec::= (see *Note ace::(ace.doc)syntax )

The defaults for the nonbond specification reside with the parameter file. The defaults are specified at the begining of the van der Waal section. These defaults are the recommended options. The following command contains all defaults for one of the older protein parameter files, and is equvalent to the command "NBONDS INIT" in it usage when this parameter file is present. NBONDS ELEC ATOM NOEX NOGR NOQU SWIT RDIE VATOM VDW VSWI - CUTNB 8.0 CTEXNB 999.0 CTOFNB 7.5 CTONNB 6.5 WMIN 1.5 WRNMXD 0.5 - EPS 1.0 NORXN EPSEXT 80.0 ORDER 10 SHELL 2.0 CTVTRN 10.0 - E14FAC 1.0 NBSCAL 1.0 IMSCAL 1.0 NBXMOD 5 NOFMA - NOEWALD NOPME KAPPA 1.0 KMAX 5 KMAXSQ 27 ERFMOD -1 MMFF specific defaults: VTRUNC MSHFT E14FAC 0.75 CTVTRN 8.0 SOFTVDW specific defaults: SOFT EMAX 30.0/EPS MINE -300.0/EPS MAXE -2.*MINE ! CDIE SOFT EMAX 15.0/EPS MINE -120.0/EPS MAXE -2.*MINE ! RDIE Values do not change unless explicitly specified, except for the ON/OFF values which cascade when the cutoff values are changed as; CTOFNB=CUTNB-0.5 CTONNB=CTOFNB-1.0 WARNING:: These old defaults have been shown to be detrimental to protein behavior. It is generally better to use the defaults in the parameter sets. The presence of soft core nonbonded terms is reccomdended for calculations in a dense system (docking, loop refinement, NMR refinement). The initial value of RMIN (swithiching distance for the soft core potential) is receommended to be 0.885 RECOMMENDED: Presented here is a suggested list of options. Where specifications are missing, substitute the defaults: For no-cutoff periodic systems: ** ** ** NBONDS ATOM EWALD PMEWALD KAPPA 0.32 FFTX 32 FFTY 32 FFTZ 32 ORDER 6 - CUTNB 12.0 CTOFNB 10.0 VDW VSHIFT (** system size dependent - use about 0.8-1.0 grids per angstrom) For atom based cutoffs: NBONDS ATOM FSHIFT CDIE VDW VSHIFT - CUTNB 13.0 CTOFNB 12.0 CTONNB 8.0 WMIN 1.5 EPS 1.0 or (perhaps better for longer cutoff distances, but more expensive) NBONDS ATOM FSWITCH CDIE VDW VSHIFT - CUTNB 13.0 CTOFNB 12.0 CTONNB 8.0 WMIN 1.5 EPS 1.0 For group based cutoffs (doesn't vectorize well): NBONDS GROUP FSWITCH CDIE VDW VSWITCH - CUTNB 13.0 CTOFNB 12.0 CTONNB 8.0 WMIN 1.5 EPS 1.0 For extended electrostatics : NBONDS GROUP SWITCH CDIE VDW VSWI EXTEND GRAD QUAD - CUTNB 13.0 CTOFNB 12.0 CTONNB 8.0 WMIN 1.5 EPS 1.0 For a better description of these methods and how they perform, see: P.J. Steinbach, B.R. Brooks: "New Spherical-Cutoff Methods for Long-Range Forces in Macromolecular Simulation," J. Comp. Chem. 15, 667-683 (1994). OPTIONS THAT ARE NOT RECOMMENDED (OR REALLY BAD): [ ATOM ] [ CDIElec ] [ SHIFted ] no (obsolete, but used in the past) [ ATOM ] [ CDIElec ] [ SWITched ] NO! Very bad - do not use! [ GROUp ] [ CDIElec ] [ SHIFted ] no (obsolete) [ GROUp ] [ CDIElec ] [ SWITched ] NO! Very bad with non-neutral groups! [ ATOM ] [ RDIElec ] [ SHIFted ] yes, but do you really want RDIE?? [ ATOM ] [ RDIElec ] [ SWITched ] no. switch is bad here.

NBSCale & IMSCale ================= The first time that the primary or image non-bond list is generated, an estimate is made, based on the number of atoms, of how much memory will be needed to store the pair list. If too large an estimate is made, memory will be wasted. If too small an estimate is made, a second (and larger) estimate will be made and the memory allocated on the first attempt is wasted. NBSCale is a correction factor to the intial estimate allowing better control of memory allocation. For example NBSCale 1.5 allocated 50% more memory than CHARMM would usually guess and NBSCale 0.8 allocated 20% less. IMSCale does the same thing when the image pair list is generated. The values of NBSCale and IMSCale must be determined empirically, but they can generate huge memory savings on large systems. These keywords are valid wherever nonbond options may appear, e.g. ENERgy, DYNAmics, MINImiz, and UPDAte. Note that NBSCale must be used in the first statement which generates a nonbond list; an UPDAte without NBSCale followed by DYNAmics with NBSCale is ineffective. For a system of about 17,000 atoms, a value of NBSCALE 1.5 was effective in providing about 25 MB of memory reduction compared to the the default NBSCale value (1.0); while for a system of about 10,000 atoms, the optimum of 1.4 gave a reduction of about 12 MB. For example: mini abnr nstep 500 nprint 10 tolenr 1.e-6 cutnb 14. ctofnb 12. ctonnb 10. - fshift cdie eps 1.0 vshift inbfrq 20 imgfrq 20 cutim 14. nbscale 1.5 Likewise for DNYAmics, ENERgy, or the UPDAte commands; the latter is useful in doing some trial and error probes to determine the optimum NBSCale value, with fixed sizes for CUTNB and CUTIM: update cutnb 14. ctofnb 12. ctonnb 10. fshift cdie eps 1.0 vshift - inbfrq 20 imgfrq 20 cutim 14. nbscale @1 - where the csh or tcsh command line might be something like: % charmm medium 1:1.5 < tst_nbscal.inp >& nbscal_1.5 & The above is based on single processor calculations; the same general idea applies to parallel calculations, but the optimum value for NBSCale will be less than 1.0, perhaps 0.7 to 0.8 for systems in the 10K to 17K atom range. Memory usage can be further reduced using the IMSCale keyword; some experimentation will be required depending on the number of atoms and the cutoffs being used. INBFrq n ======== Update frequency for the non-bonded list. Used in the subroutine ENERGY() to decide whether to update the non-bond list. When set to : 0 --> no updates of the list will be done. +n --> an update is done every time MOD(ECALLS,n).EQ.0 . This is the old frequency-scheme, where an update is done every n steps of dynamics or minimization. -1 --> heuristic testing is performed every time ENERGY() is called and a list update is done if necessary. This is the default, because it is both safer and more economical than frequency-updating. Description of the heuristic testing algorythm : ----------------------------------------------- Every time the energy is called, the distance is computed each atom moved since the last list-update. If any atom moved by more than (CUTNB - CTOFNB)/2 since the last list-update was done, then it is possible that some atom pairs in which the two atoms are now separated by less than CTOFNB are not in the pairs-list. So a list update is done. If all atoms moved by less than (CUTNB - CTOFNB)/2 , then all atom pairs within the CTOFNB distance are already accounted for in the non-bond list and no update is necessary. Description of the code for the heuristic testing : -------------------------------------------------- This section describes how programmers can control the list-updating behavior when their routines call the ENERGY() subroutine. All list-updating decisions, whether they are frequency based or heuristic based, are made in the subroutine UPDECI(ECALLS) , which is called from only one place : at the very beginning of ENERGY(). UPDECI(ECALLS) can be controled through INBFRQ (via the CONTRL.FCM common block) and ECALLS (via the ENERGY.FCM common block) as follows : If INBFRQ = +n --> non-bond list is performed when MOD(ECALLS,n).EQ.0 . Image and H-bond lists are updated according to IMGFRQ and IHBFRQ. If INBFRQ = 0 --> non-bond list update is not performed. Image and H-bond lists are updated according to IMGFRQ and IHBFRQ. If INBFRQ = -1 --> all lists are updated when necessary (heuristic test). (note that ECALLS is incremented by ICALL every time ENERGY(,,,ICAL) is called. In most cases, ICALL=1 ) The current implementation of UPDECI() will work (without modifications) to decide whether the image/crystal non-bond lists need updating, provided the periodicity parameters don't change (i.e. constant Volume). UPDECI() is easily adapted to variable Volume dynamics/minimizations. This is described in comments of the routine itself. Further computational economy in update-testing : ------------------------------------------------- A programmer can sometimes skip the heuristic test itself, making the decision whether to do list-updating even more economical. This option is only available if the size of the step taken since the last call to ENERGY() is known. For an example of usage, see the subroutine ENERG() in TRAVEL. NON-BOND ENERGY TERMS. ====================== The electrostatic options are separate from the van der Waal options, though some keywords are shared between them. The following is a description of all options and keywords. 1) Electrostatics The ELEC keyword (default) invokes electrostatics. The NOELec keyword supresses all electrostatic energy terms and options. There are two basic methods for electrostatics, GROUp and ATOM. A model based on the GROUp method is the extended electrostatics model which approximates the full electrostatic interaction and eliminates the need for a cutoff function. This model is based on the partitioning of the electrostatic term into two contributions. One comes from the interaction between particles which are spatially close and is treated by conventional pairwise summation. The second contribution comes from interactions between particles which are spatially distant from one another and is treated by a multipole moment approximation. [The original model was described in B. R. Brooks, R. E. Bruccoleri, B. D. Olafson, D. J. States, S. Swaminathan, M. Karplus. J. Comp. Chem., 4, 187, (1983) and more recently in R.H. Stote, D.J. States and M. Karplus, J. Chimie Physique (1991)] A) Functional Forms Atom electrostatics indicates that interactions are computed on an atom-atom pair basis. There are two options that specify the radial energy functional form. The keywords CDIE and RDIE select the basic functional form. The SWIT, SHIF, FSWI, and FSHI keywords determine the long-range truncation option. CDIE - Constant dielectric. Energy is proportional to 1/R. RDIE - Distance dielectric. Energy is proportional to 1/(R-squared) SWIT - Switching function used from CTONNB to CTOFNB values. SHIF - Shifted potential acting to CTOFNB and zero beyond. FSWI - Switching function acting on force only. Energy is integral of force. FSHI - Classical force shift method for CDIE (force has a constant offset). B) Atom electrostatics Atom electrostatics indicates that interactions are computed on an atom-atom pair basis. There are two options that specify the radial energy functional form. The keywords CDIE and RDIE select the basic functional form. The SWIT and SHIF keywords determine the long-range truncation option. [ ATOM ] [ CDIElec ] [ SHIFted ] [ RDIElec ] [ SWITched ] [ FSWItch ] [ FSHIft ] CDIE - Constant dielectric. Energy is proportional to 1/R. RDIE - Distance dielectric. Energy is proportional to 1/(R-squared) SWIT - Switching function used from CTONNB to CTOFNB values. SHIF - Shifted potential acting to CTOFNB and zero beyond. FSWI - Switching function acting on force only. Energy is integral of force. FSHI - Classical force shift method for CDIE (force has a constant offset) C) Group Electrostatics electrostatics-spec::= [ GROUp [ EXTEnded [ GRADients ] [ QUADrip ] ] ] [ CDIElec ] [ SWITched ] [ [ [ NOGRad ] [ NOQUads ] ] ] [ RDIElec ] [ FSWItch ] [ [ NOEXtended ] ] SWIT - Switching function used from CTONNB to CTOFNB values. FSWI - Switching function, but QiQj/Rcut is added before switching. (FSWI has no effect on neutral groups). D) Electrostatic Distances electrostatic-dist::= [ CUTNB real ] [ CTEXNB real ] [ CTONNB real ] [ CTOFNB real ] [ EMAX real ] [ MINE real ] [ MAXE real ] ! SOFTVDW only EMAX - Twice the VDW energy from which the soft core potential becomes active (if SOFT key word is employed). It has a linear form for r<rc (E(rc)=EMAX/2) : Esoftl=EMAX/2+alfa*(r-rc) unless the VDWE key word turns on the exponential form : EsoftE=EMAX+alfa*r**beta For the exponential form E(0)=EMAX, so EMAX is also the VDW energy at r=0 for the exponential form MINE - Twice the energy from which the electrostatic attractive soft potential begins. ELEE turns on the exponential form. MAXE - Twice the energy from which the electrostatic repulsive soft potential begins. ELEE turns on the exponential form. The soft core potential is currently implemented only in fast energy routine, so fast option has to be used to allow it. The form of the soft potential for the electrostatics is the same as for the VDW, the defaults are different for VDW, electrostatic repulsion and attraction. The defaults need to be modified for some cases (i.e. for the spc water model) to prevent shifting of energy minima. The SOFT key word turns on both VDW and electrostatic soft core. To turn off the soft core set EMAX = 0 with the SOFT keyword (in energy, minimization or nbonds call). CTEXNB - defines the cutoff distance beyond which interaction pairs are excluded from the Extended Electrostatics calculation. E) Extended (group) Electrostatics electrostatics-spec::= [ ATOM ] [ CDIElec ] [ SHIFted ] [ GROUp [ EXTEnded [ GRADients ] [ QUADrip ] ] ] [ RDIElec ] [ SWITched ] [ [ [ NOGRad ] [ NOQUads ] ] ] [ [ NOEXtended ] ] EXTE - invokes the extended electrostatics command for calculating long range electrostatic interactions. NOEX - suppress the extended calculation. GRAD - keyword flags the inclusion of the field of the extended gradient in calculating the force on each atom,i.e. include first and second derivatives. QUAD - flags the inclusion of the quadrupole in the multipole expansion. F) Reaction Fields misc-specs::= [ EPS real ] [ E14Factor real ] [ NORXN ] [ RXNFLD rxnfld-spec ] [ RXNNB rxnfld-spec ] rxnfld-spec::= [ EPSEXT real ] [ ORDER integer ] [ SHELL real ] 2) Van Der Waal Interactions The VDW keyword (default) invokes the van der waal energy term. To supress this term, the NOVDw keyword may be used. A) Distance specified van der Waal Function [ VGROUP [ VSWITched ] ] [ ] vdw-spec::= [ VATOM [ VSHIfted ] ] [ [ VSWItched ] ] [ [ VFSWitch ] ] 3) Miscellaneous options and keywords A) Dielectric specification misc-specs::= [ EPS real ] [ E14Factor real ] B) Warning Distance Specifications warning-dist::= [ WMIN real ] [ WRNMXD real ] WRNMXD - keyword defines a warning cutoff for maximum atom displacement from the last close contact list update (used in EXTEnded) C) Initialization In all cases as many keywords and values as desired may be specified. The key words, their functions, and defaults are: 1) The method to be used 2) Distance cutoff in generating the list of pairs CUTNB value (default 8.0) 3) Distance cut at which the switching function eliminates all contributions from a pair in calculating energies. Once specified, This value is not reset unless respecified. CTOFNB value (default CUTNB-0.5) 4) Distance cut at which the smoothing function begins to reduce a pair's contribution. This value is not used with SHFT. Once specified, This value is not reset unless respecified. CTONNB value (default CTOFNB-1.0) 6) Dielectric constant for the extened electrostatics routines (RDIE option sets the dielectric equal to r times the EPS value) EPS value (default 1.0 for r dielectric) EPS 0.0 or NOELec (zero elecrostatic energy) 7) Warning cutoff for minimum atom to atom distance. Pairs are checked during close contact list compilation. WMIN value (default 1.5) 8) Warning cutoff for maximum atom displacement from the last close contact list update (used only in EXTEnded) WRNMXD value (default 0.5) 9) The presence of soft core nonbonded interactions (turned on only if SOFT key word is present) D) Exclusion Lists By default, vdw and electrostatic interactions between two bonded 1-2 interactions) and two atoms bonded to a common atom (1-3 interactions) atoms are excluded from the calculation of energy and forces. Also, special vdw parameters and an electrostatic scale factor (E14FAC) can be applied to atom pairs separated by 3 bonds (1-4 interactions). The control of the exclusion list is by the integer variable, NBXMod. NBXMod = 0 Preserve the existing exclusion lists NBXMod = +/- 1 Add nothing to the exclusion list NBXMod = +/- 2 Add only 1-2 (bond) interactions NBXMod = +/- 3 Add 1-2 and 1-3 (bond and angle) NBXMod = +/- 4 Add 1-2 1-3 and 1-4'S NBXMod = +/- 5 Add 1-2 1-3 and special 1-4 interactions A positive NBXMod value causes the explicit exclusions in the PSF (inb array) to be added to the exclusion list. A negative value causes the use of only the bond connectivity data to construct the exclusion list (thus, ignoring the PSF data). ALGORITHMS There are four algorithms used in calulating the nonbonded energies, each making different approximations in an attempt to speed the calulation. 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 caluculations 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 SHIFt option is similar to SWITch except, the potential: E = (QI*QJ/(EPS*R)) * ( 1 - (R/CTOFNB)**2 )**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. CDIE 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. The short range forces are identical to those calculated with the other options, reflecting the decrease in dielectric shielding at short ranges. 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. Description of the Extended Electrostatics method ------------------------------------------------- For the long range forces there is effectively no cutoff in the electrostatic energy when using the Extended Electrostatics model. The Extended Electrostatics model approximates the full electrostatic interaction by partitioning the electric potential and the resulting forces at any point ri into a near and extended contribution. The near contribution arises from the charged particles which are spatially close to ri while the extended contribution arises from the particles which are spatially distant from ri. The total electrostatic potential can be written as a sum of the two. The near region is defined in terms of a radial distance, CUTNB, for each atom. Interactions between atoms separated by a distance greater than CUTNB are calculated using a time saving multipole approximation when the nbond list is updated. These interactions are stored together with their first (NOGRad) or first and second (GRADients) derivatives. Interactions between particles within CUTNB are calculated by the conventional pairwise additive scheme. (For a more complete development of the model, see R.H. Stote, D.J. States and M. Karplus, J. Chimie Physique Vol. 11/12, 1991). 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. In using this model the GROUp method for constructing the nonbond list must be used. The interactions between particles within CUTNB are truncated rather than having a SHIFt or SWITch function applied. Additionally, as one is calculating all electrostatic interactions in the system, the dielectric constant should be set to 1.0. Not Available at this time: 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. IMPLEMENTATION AND DATA STRUCTURES The initialization and list compilation is performed by the subroutine NBONDS which in turn will call a lower level routine that will do all of the work. It functions by guessing how much space will be needed to store the close contact list, allocating that space (and space for electrostatic potentials and gradients if necessary) on the heap, and calling the appropriate subroutine to actually compile the nonbonded list (NBONDG,...). If sufficient space was not available 1.5 times as much is allocated and another attempt is made. ENBOND evaluates the nonbonded energy, calling EEXEL to evaluate the stored electric potentials and fields. Double precision is used for all arithmatic. All of the nonbonded cutoffs and lists are stored on the heap. BNBND is the descriptor array passed through most of the program (in some of the analysis routines an additional array BNBNDC is used for the comparision data structure). BNBND holds heap adresses and LNBND holds the lengths of the elements in the data structure. To actually access the data it is necessary to include INBND.FCM (an index common block) and specify HEAP(BNBND(xxx)) where xxx is the desired element name in INBND.FCM. This is arrangement has the advantage of allowing dynamic storage allocation and easy modification of the types of information passed from routine to routine. The nonbonded data structure is described in: source/fcm/inbnd.fcm FAST VECTOR/SHARED-MEMORY PARALLEL ROUTINES There are 5 sets of standard routines used to compute nonbond energies and forces. The can be summarized: FAST CRAYVEC - optimized vector code for a Cray (array processor) FAST PARVECT - optimized vector code for an SMP machine FAST VECTOR - general optimized vector code FAST SCALAR - general optimized scalar code FAST OFF - the generic - support everything routine These routines are processed in the order listed. The highest gaining priority based on what options have been compiled and what the user requests. If the user does not specify a fast option and all code is compiled, then the CRAYVEC code will probably be used, unless this routine does not support the requested options (in which case, the next routine is tried). For example, the PARVECT code supports PME Ewald, but CRAYVEC does not, so a calculation with PME will run with PARVECT (unless otherwise specified). To determine which routine is actually doing your calculation, use "PRNLEV 6" to list energy routines as they are processed. OTHER The option EXOForce, forces update of exclusion lists. The option SOFT will turn on the soft core nonbonded potential. The option GRAPe will perform all the nonbond interactions in the specialized molecular dynamics hardware, called GRavity PipE. The option NOLIst will perform all the no cuttof nonbond interactions without the use of nonbond list, the same as on GRAPE machine.

Nonbond Lookup Tables The nonbond energy terms may be specified with a user supplied binary lookup table. The command; READ TABLE UNIT int will invoke this feature and disable all other energy term options. The nonbond list specifiers will still be used (cutoff distances...). This feature is not designed for casual users, and is not supported with test cases. Also, in version17, there is an uncorrected bug in the second derivative determination. To use this feature, first read the common file ETABLE.FCM for a description of variables, and then create a file the the routine REATBL (consult the source) can read. The sources for this option are contained in the file ETABLE.FLX.

SUMMARY The purpose of using finite cutoffs in energy calculations is to reduce, from O(N*N) to O(N), the number of nonbonded interaction terms that actually need to be computed. CHARMM has three ways of building the nonbonded interaction list: BYGRoups, BYCUbes, and 'By-Clusters-in-Cubes,' or BYCC. The BYCC method is a combination of the earlier BYGR and BYCU methods. For a given set of atoms and a given cut-off distance, all three algorithms should generate the same non-bonded list. BYCBim extends the BYCUbes method to systems with images or periodic boundaries, but is more restricted in other options that are compatible with it. It will not work with replica, extended electrostatics, or constant pressure dynamics since no group list is generated. However, BYCBim does work in parallel mode; BYCUbes does not. The differences between BYCUbes, BYGRoups, and BYCC are in speed and memory requirements. This section gives a description of the algorithms, followed by a description of some of the unique aspects of BYCC. Finally, a few general guidelines for choosing between the earlier BYCUbes and BYGRoups methods are presented. ALGORITHMS Atom-based calculations: The basis of the efficiency of the BYGRoups algorithm over a brute-force comparison of each atomic position with all others in the system is that it clusters atoms into chemical groups, initially ignoring the individual atoms, themselves. This significantly reduces the number of pairs of particles that need to be examined. Effectively, BYGR speeds up the calculation by reducing the particle density, which it does by simply redefining the particles. Once a list of group-group pairs satisfying the initial distance criterion is made, only atoms from this relatively short list are then considered for further atom- atom distance testing. In contrast, the efficiency of the BYCUbes algorithm is based on the partitioning of the system into small cubical regions. A synopsis of how it works: (1) Find a rectangular parallelepiped that bounds the system (with margins), is aligned with the Cartesian axes, and has sides that are integral multiples of the cutoff distance. Divide this parallelepiped, or box, into cubes whose sides are equal in length to the cutoff distance. (2) Identify the cube that contains each atom. (3) For each cube C, loop over each atom A contained in C. For each such A, loop over each atom A' contained in the 27-cube surrounding region, which is the (3 cube x 3 cube x 3 cube) region that contains the central cube. If the A--A' pair falls within the cutoff distance, check various other inclusion criteria; e.g. that the pair is not a 1-2 or 1-3 excluded pair. Hence, BYCU speeds up the generation of the non-bonded list by distance- testing only pairs of atoms that are in nearby regions. The 'By-Clusters-in-Cubes', or BYCC, algorithm incorporates both the partitioning technique of BYCUbes and also the atomic grouping technique of BYGRoups. BYCC divides the system into a cubical grid and compares particles only in adjacent cubes, as BYCU does. However, the particles it compares are clusters of atoms, in the spirit of BYGR, as opposed to individual atoms. The use of both techniques in BYCC allows for greater efficiency than is possible using either technique alone, and for this reason BYCC is generally faster than the other algorithms. In addition, because it does the final atom-atom distance calculations, the exclusions, and the formation of the non-bonded list all in one final loop, the routine needs to store only a cluster-cluster pair- list as a work array. Thus, the memory requirements are reduced relative to BYCU, which stores a much longer atom-atom pairlist (essentially a second non-bonded list) internally. Like BYCU, the computational time for BYCC increases linearly with the number of atoms and sigmoidally as a function of cut-off distance. When (cut-off distance) << (radius of system) the time dependence on cut-off distance is essentially third order, but depends on the system, the machine type, and the actual cutoff value. The computational time for BYGR increases quadratically with the size of the system. NB: The energy calculations are always O(N) if a non- bonded list is being used, regardless of which method is used to generate the list. However, generation of the list itself can be either O(N*N) (e.g. BYGR) or O(N) (e.g. BYCU or BYCC). Group-based calculations: The BYCC option, like the others, supports calculations based on chemical groups instead of atoms, and an additional option for extended electrostatics. However, clusters are not used in group-based calculations, since their role is in effect subsumed by the groups. The relative speed advantage of BYCC over BYCU is therefore diminished in group-based computations. For the BYCC or BYCU options, when extended electrostatics are requested, CPU time for nonbonded list generation and calculation of extended electrostatics will depend on the extended electrostatics cut-off distance, CTEXNB. Hence, smaller CTEXNB values will significantly speed up the calculations. For BYGR, CPU time is independent of CTEXNB. Note that for group-based calculations with either BYCC, BYCU, or BYGR, the groups are treated as point-like. This means that, unlike the case for cluster size in atom-based calculations, there is no dependence of CPU time on group size in the generation of the non-bonded list. BYCC: CLUSTERS With the BYCC option, a requirement for the atom-based calculations is that clusters of atoms need to be created. A cluster is a set of atoms that have mutually close connectivity relationships. Generally, small, dense, uniformly-sized clusters yield the most efficient non-bonded list generation. The use of a separate clustering scheme instead of the usual chemical grouping arrangement used elsewhere in CHARMM (e.g. in the BYGRoups algorithm) provides a separate spatial organization of the system that can be manipulated and optimized largely independently of the other CHARMM functions. Because the criteria for the grouping of atoms are chemical, whereas those for the clustering of atoms are spatial, the optimal arrangement of atoms for CHARMM groups is not necessarily (and, indeed, not usually) the optimal arrangement for clusters. The command: MKCLuster [CADL] [CMARgin] [NOBRidge] generates clusters from the structural information in the .psf file. This command MUST be executed prior to using the BYCC option in NBONd. 1) Keyword CADL Consecutive atom distance limit -- distance between consecutive atoms (as defined by order in RTF/PSF) in a cluster beyond which the cluster will be split. Helps prevent the generation of overly large clusters. Default : 4.0 Angstroms. 2) Keyword CMAR Limiting cluster margin width. The cluster margin is the calculated distance that is added to the non-bonded cutoff distance, when partitioning the clusters in the system, to ensure that all atom-atom pairs are in fact included in the distance testing. This margin depends on the largest cluster size and will affect CPU time (larger margin means slower non-bonded list generation). CMAR sets the upper limit for this added margin. Default: 5.0 Angstroms. Note that while the calculated cluster margin is the width that GUARANTEES all-atom testing, in practice CMAR can often be set below the calculated value (speeding up the calculations) without changing the results. This is particularly true when there is a large spread in the cluster sizes or when there are statistical outliers. If CMAR < calculated cluster margin, a warning is issued. The MKCL command also causes the generation of the table of exclusions. Since the BYCC non-bonded list-generating algorithm uses only this table for the exclusions, and since the table does not change in a CHARMM run unless MKCL is executed again, for maximal efficiency it is recommended that NBXMod be set to 0 (see nbonds.doc) when BYCC is used. That is, changes in the exclusion lists will only be reflected in the energies when MKCL is re-executed before non-bonded list generation. Clusters should be re-generated if there is significant distortion in connectivity during a CHARMM run. If the clusters initially generated become too large, the algorithm will slow significantly unless MKCL is re-invoked. BYCC: ACTIVE ATOMS An additional feature of the BYCC algorithm is its ability to handle "active" atom selections, which are specified with the command NBACtive [atom selection] The purpose of this feature is to allow the user to focus calculations on regions of interest without having to alter the psf or coordinate files. It allows the non-bonded list generation routine (BYCC) to ignore completely atoms that are not defined as active, so that only active atoms appear in the non-bonded list. This speeds up non-bonded updates and energy calculations, and, in addition, it allows for more selectivity in energy calculations-- in dynamics and minimization, particularly. An example of the ideal use of this command is the study of a single subunit in a protein containing several. Another is the study of a single side chain, or a group of side chains, on a fixed protein backbone. If the energy calculations are to be consistent with the inactive portion of the system, a buffer region (ideally of width CUTNB) is required surrounding the region of interest. This buffer region should be active and fixed. Active clusters and groups are defined automatically on the basis of the active atom selection. If no active atom selection is given, it is assumed that all atoms are active. In its current implementation, the active atom selection does not affect the bonded energies (note to developers: this should be rectified in the future). This has at least two implications: 1) It is currently recommended that in conjunction with this command, 'inactive' atoms be fixed with the CONS FIX command, since otherwise they will contribute to the bonded terms but not to the non-bonded terms. 2) The connectivity between inactive and active regions, if it exists in the original system, is not broken by the NBACtive command. This means that regions defined as active must remain 'tethered' to the inactive regions. BYCC: EXAMPLE OF ADDITIONAL COMMANDS The additional commands in a typical CHARMM script using the BYCC option for non-bonded list generation would be as follows: The MKCLusters command and the BYCC keyword in the NBOND command are necessary for using the BYCC option. The other commands are optional. ! make clusters MKCL !necessary ! define active atom set and buffer subset (optional) DEFINE ACTIVE SELE (active atom selection) END DEFINE BUFFER SELE (buffer atom selection) END !subset of ACTIVE ! "activate" the selected atoms (optional) NBACtive SELE ACTIVE END !default is all atoms ! fix inactive atoms (optional) CONS FIX SELE ((.NOT.ACTIVE) .OR. BUFFER) END ! call non-bonded list generation routine NBONd - BYCC !necessary PAYOFF THRESHOLD: Comparison of BYCUbes and BYGroups Given any pair of functions F1(N) and F2(N), which are O(N) and O(N*N) respectively, there will always be some constant N0 such that F2(N) > F1(N) for all N > N0. This constant N0 may be referred to as the "payoff threshold," since it is the system size above which the BYCUbes algorithm will be faster than BYGRoups for a given cut-off distance, particle density, machine type, and set of options. The following are some properties of the payoff threshold that can be used as guidelines for choosing between BYGRoups and BYCUbes: (1) The payoff threshold is smaller on a vector machine than on a scalar machine. That is, BYCUbes is more vectorizable than BYGRoups. (2) The shape of the system does not have a big impact on the payoff threshold. BYCUbes operates by drawing a rectangular box around the system and dividing it into cubes, but on a serial machine the empty cubes take little time. (3) The payoff threshold usually grows steeply with cutoff distance. This is because the speed advantage of the BYCUbes algorithm is based upon the compartmentalization of the system. The smaller the compartments (cubes), the more advantage. Since the compartment size is based on the cutoff distance, as the cutoff distance increases, there is less advantage for BYCUbes over BYGRoups. (4) BYCUbes trades memory for time. Its memory requirements are significantly higher than those of BYGR and they may be prohibitive for large systems.

Information and HTML Formatting Courtesy of:

NIH/DCRT/Laboratory for Structural Biology

FDA/CBER/OVRR Biophysics Laboratory

Modified, updated and generalized by C.L. Brooks, III

The Scripps Research Institute