Energy Manipulations: Minimization and Dynamics The main purpose of CHARMM is the evaluation and manipulation of the potential energy of a macromolecular system. In order to compute the energy, several conditions must be met. There are also several support commands which directly relate to energy evaluation. * Menu: * Description:: Description of the energy commands * Skipe:: Selection of particular energy terms * Interaction:: Computation of interaction energies and forces. * Fast:: Requirements for using the fast routines * Needs:: Requirements for all energy evaluations * Optional:: Optional actions to be taken beforehand
Syntax for Energy Commands There are two direct energy evaluation commands. One is parsed through the minimization parser and the other involves a direct call to GETE. See *note Minimiz:(chmdoc/minimiz.do,,) an d*note Gete: (usage.doc)interface. In addition to getting the energy, the forces are also obtained. The ENERgy command. (processed through the minimization parser) [SYNTAX ENERgy] ENERgy [ nonbond-spec ] [ hbond-spec ] [ image-spec ] [ print-spec ] [ COMP ] [ INBFrq 0 ] [ IHBFrq 0 ] [ IMGFrq 0 ] [NOUPdate] hbond-spec *note Hbonds:(hbonds.doc). nonbond-spec *note Nbonds:(nbonds.doc). image-spec *note Images:(images.doc)Update. If the COMP keyword is specified, then the comparison coordinate set is used, but this disables the use of the fast routines. The keyword NOUPdate turns off all update routines, and thus requires all lists to be present already. The GETE command. (a direct call to GETE) [SYNTAX GETEnergy] GETE [ COMP ] [ PRINt [ UNIT int ] ] For this command to work, all list must be set up. This is best done through the UPDAte command. The COMP keyword will cause the comparison coordinate set to be used. The PRINt keyword will result in a subsequent call to PRINTE in order to print the energy. If the PRINt keyword is not specified, then NO indication that the energy has been called will be given. The UPDAte command (sets up required lists for GETE) [SYNTAX UPDAte lists] UPDAte [ nonbond-spec ] [ hbond-spec ] [ image-spec ] [ COMP ] [ INBFrq 0 ] [ IHBFrq 0 ] [ IMGfrq 0 ] [ EXSG {list-of-segment-names} | EXOF ] The update command will set up the codes lists and also create a nonbond list (unless INBFrq is 0) and a new hbond list (unless IHBFrq is 0). If the COMP keyword is specified, then the comparison coordinates will be used in setting up the nonbond and hbond lists. EXSG keword with optional following list of segment names allows to exclude some nonbonded interactions (ELEC & VDW). If list of names is empty ALL INTERsegment nonbonded interactions will be excluded. If list is not empty all INTER and INTRA segment nonbonded interactions for listed segments will be ecluded. EXOF turns off this option. H-bond energies (HBON) are not affected at the moment (Dec 3, 1991).
Skipping selected energy terms There is a facility to skip any desired energy terms during energy evaluation. For each energy term there is associated a logical flag determining whether that energy term is to be computed. Specifications are processed sequentially. The default operation is INCLude which implies that subsequent energy term are to be removed from the energy calculation. NOTE: that EXCLude implies that the energy term is to be computed. If for some reason, the list presented here is out of date, the data in SKIPE(energy.src) and in ENER.FCM of the source should be consulted. Syntax: [SYNTAX SKIP energy terms] [ INCLude ] [ EXCLude ] SKIPe repeat( [ ALL ] ) [ NONE ] [ item ] item::= [ BOND ] [ ANGL ] [ UREY ] [ DIHE ] [ IMPR ] [ VDW ] [ ELEC ] [ HBON ] [ USER ] [ HARM ] [ CDIH ] [ CIC ] [ CDRO ] [ NOE ] [ SBOU ] [ IMNB ] [ IMEL ] [ IMHB ] [ XTLV ] [ XTLE ] [ EXTE ] [ RXNF ] [ ST2 ] [ IMST ] [ TSM ] [ QMEL ] [ QMVDW] [ ASP ] [ EHARM] [ GEO ] [ MDIP ] [ STRB ] [ VATT ] [ VREP ] [ IMVREP ] [IMVATT] [ OOPL ] description: BOND - bond energy ANGL - angle energy UREY - Urey-Bradley energy term DIHE - dihedral energy IMPR - improper dihedral energy VDW - van der Waal energy ELEC - electrostatic energy HBON - hydrogen bond energy USER - user supplied energy (USERLINK) HARM - harmonic positional constraint energy CDIH - constrained dihedral energy CIC - internal coordinate constraint energy CDRO - quartic droplet potential energy NOE - NOE general distance restraints SBOU - solvent boundary energy IMNB - image van der Waal energy IMEL - image electrostatic energy IMHB - image hydrogen bond energy XTLV - crystal van der Waal energy XTLE - crystal electrostatic energy EXTE - extended electrostatic energy RXNF - reaction field energy ST2 - ST2 water-water energy IMST - image ST2 water-water energy TSM - TMS free energy term. QMEL - energy for the quantum mechanical atoms and their electrostatic interactions with the MM atoms using the AM1 or MNDO semi-empirical approximations QMVDW - van der Waals energy between the quantum mechanical and molecular mechanical atoms ASP - solvation free energy term based on Wesson and Eisenberg surface area method EHARM - second harmonic restraint term (for implicit Euler integration) GEO - Mean-Field-Potential energy MDIP - MDIPole mean fields constraints STRB - strech-bend interaction (MMFF) VATT - VdW attraction (MMFF) VREP - VdW repulsion (MMFF) IMVREP - image VdW repulsion (MMFF) IMVATT - image VdW attraction (MMFF) OOPL - out-of-plane (MMFF) Examples; SKIP ALL EXCL BOND - do just bond energy SKIP EXCL ALL - return flags to default state SKIP ELEC VDW - throw out electrostatics and van der Waals energy
Interaction energies and forces The INTEraction command computes the energy and forces between any two selections of atoms. [SYNTAX INTEraction energy] INTEraction [ COMP ] [ NOPRint ] 2x(atom-selection) [UNIT int] If only one atom selection is given, then a self energy will be computed. This routine is quite efficient and may be used within a Charmm loop without too much overhead, though there are some restrictions. The COMP keyword causes the comparicon coordinates to be used. The NOPRint keyword will prevent the results from being printed. This routine works in the same manner as the GETE command in that all of the lists (CODES, nonbond, and Hbond) must be specified before invoking this command. One difference is that SHAKE will not be respected with this command (i.e. if the coordinates don't satisfy the constraints, neither will the energy). The following energy terms may be computed by this routine (unless supressed with the SKIP command); Bond - Energy defined by the two atoms involved. Angles - Energy allocated to the central atom (auto energy only). Dihedral - Energy defined between central two atoms Improper - Energy defined by first atom (auto energy only) van der Waal - ATOM option only. Energy defined by two atoms involved. Electrostatic - ATOM option only. Energy defined by two atoms involved. Hbond - Energy defined by heavy atom donor and acceptor atom. Harmonic cons - Energy allocated to central atom (auto energy only). Dihedral cons - Energy defined by central two atoms. User energy - Atom selections may be passed to USERE in the selection common (DEFIne command). Fill forces and energies as desired. All other energy terms will be zeroed. For terms listed "auto energy only", the corresponding atom must be present in both atom selections. For the remaining terms, one atom of the pair must be present in each of the atom selections. The energy division matches the method used in the analysis facility. This command will not work with the selection of images atoms, or the selection of ST2 waters. All energy terms not listed above will not be computed. The nonbond list must be generated with the ATOM and VATOM options. The individual energy terms are stored in the energy common and are available in commands and titles via the "?energy-term" substitution. The forces for all kept energy terms will be returned in the force arrays. Note, that it is possible for atoms to have a force that were not selected in either selection specification. This may happen for angle or dihedral terms on the first and last atoms. It may also happen in a similar manner for improper dihedrals, hydrogen bonding terms, and dihedral constraints.
[SYNTAX FASTer ] FASTer {integer} {OFF } {ON } {DEFAult} {SCALar } ! for testing only {VECTor } ! for testing only Instead of using an integer value, FASTer command can be issued with one of the following keywords. Keyword Equivalent integer ---------------- ---------- FASTer OFF -1 DEFAult 0 ON 1 SCALar 2 VECTor 3 The FASTer keyword or integer defines which versions of the energy routines to be used. FASTer -1 : Always use slow routines FASTer 0 : Use fast routine if possible, no error if cannot (default) FASTer 1 : Use best optimized routine for the current machine (Error message if cannot) FASTer 2 : Use fast scalar routine (Error message if cannot) FASTer 3 : Use fast vector routine (Error message if cannot) There exist a general and a fast version of the internal energy routines (bond, angle, dihedral, and improper dihedral). The is also a fast version of nonbond energy evaluation (roughly 30-50% faster). These routines were designed for long minimization or dynamics calculations. To request the FAST routine, the FASTer command should be used with a positive integer or an appropriate keyword. A negative integer will disable the fast energy routines. If the fast routines are requested and it is not possible to use the fast routines, a warning will be issued, and the general routines will be used in their place. The fast routines are more efficient in several ways; (1) arrays are included in common files rather than passed (2) second derivatives have been removed (3) analysis and print options have been removed The restrictions are that; (1) the MAIN coordinate set must be used in the energy evaluations (2) second derivatives may not be requested (3) The PSF, parameter, and codes arrays must be used (from the common files) (4) a limited set of nonbond options must be used. The current nonbond options supported by the fast nonbond routine are as follows. ATOM [CDIE] [SHIFt ] VATOM [VSHIft ] [RDIE] [SWITch ] [VSWItch ] [FSWItch] [VFSWitch] [FSHIft ] GROUP [CDIE] [SWITch ] VGROUP [VSWItch ] [RDIE] [FSWItch]
Requirements before energy manipulations can take place Before the energy of a system can be evaluated and manipulated, a number of data structures must be present. First, a PSF must be present. Second, a parameter set must be present. It must contain all parameters which are required by the PSF being used. Third, coordinates must be defined for every atom in the system. An undefined coordinate has a particular value, and if two coordinates have the same value, division by zero will occur in the evaluation of the energy. If the positions of hydrogens are required, the hydrogen bond generation routine, see *note Hbond: (hbonds.doc), must be called before the energy is evaluated. Fourth, provisions must be made for having a hydrogen bond list and a non-bonded interaction list. Having non-zero frequencies for updating this lists is one way, one can also read these lists in, see *note read:(io.doc)read, or generate them with separate commands, see *note HBgen:(hbonds.doc), or *note NBgen:(nbonds.doc).
Optional actions you can take to modify the energy manipulations There exist several commands which can modify the way the potential energy is calculated or can affect the way energy manipulations are performed. The Constraint command, see *note Cons:(cons.doc), can be used to constraints of various kinds. First, it can be used to set flags for particular atoms which will prevent them from being moved during minimization or dynamics. Second, it can be used to add positional constraint term to the potential energy. This term will be harmonic about some reference position. The user is free to set the force constant. Third, the user can place a harmonic constraint on the value of particular torsion angles in an attempt to force the geometry of a molecule. Other constraints are also available. The SHAKe command, see *note shake:(cons.doc)SHAKE, is used to set constraints on bond lengths and also bond angles during dynamics. It is very valuable in that it permits a larger step size to be used during dynamics. This is vital for dynamics where hydrogens are explicitly represented as the low mass and high force constant of bonds involving hydrogen require a ridiculously small step size. The user interface commands can be used to modify the calculation of the potential and to add another term to the potential energy. See *note Modify:(usage.doc)interface for details.
The following command line substitution values may be included in any command or title. To get the total energy, the syntax; ...... ?TOTE ..... should be used. Energy related properties: 'TOTE' - total energy 'TOTK' - total kinetic energy 'ENER' - total potential energy 'TEMP' - temperature (from KE) 'GRMS' - rms gradient 'BPRE' - boundary pressure applied 'VTOT' - total verlet energy (no HFC) 'VKIN' - total verlet kinetic energy (no HFC) 'EHFC' - high frequency correction energy 'EHYS' - slow growth hysteresis energy correction 'VOLU' - the volume of the primitive unit cell = A.(B x C)/XNSYMM. Defined only if images are present. 'PRSE' - the pressure calculated from the external virial. 'PRSI' - the pressure calculated from the internal virial. 'VIRE' - the external virial. 'VIRI' - the internal virial. 'VIRK' - the virial "kinetic energy". Energy term names: 'BOND' - bond (1-2) energy 'ANGL' - angle (1-3) energy 'UREY' - additional 1-3 urey bradley energy 'DIHE' - dihedral 1-4 energy 'IMPR' - improper planar of chiral energy 'VDW ' - van der waal energy 'ELEC' - electrostatic energy 'HBON' - hydrogen bonding energy 'HARM' - harmonic positional restraint energy 'CDIH' - dihedral restraint energy 'CIC ' - internal coordinate restraint energy 'CDRO' - droplet restraint energy (approx const press) 'USER' - user supplied energy term 'IMNB' - primary-image van der waal energy 'IMEL' - primary-image electrostatic energy 'IMHB' - primary-image hydrogen bond energy 'SBOU' - solvent boundary lookup table energy 'NOE' - general distance restraint energy (for NOE) 'XTLV' - crystal van der waal energy 'XTLE' - crystal electrostatic energy 'EXTE' - extended electrostatic energy 'RXNF' - reaction field electrostatic energy 'ST2' - ST2 water-water energy 'IMST' - primary-image ST2 water-water energy 'TSM' - TMS free energy term 'QMEL' - QM/MM electrostatic energy term 'QMVD' - QM/MM van der Waals energy term Energy Pressure/Virial Terms: 'VEXX' - External Virial 'VEXY' - 'VEXZ' - 'VEYX' - 'VEYY' - 'VEYZ' - 'VEZX' - 'VEZY' - 'VEZZ' - 'VIXX' - Internal Virial 'VIXY' - 'VIXZ' - 'VIYX' - 'VIYY' - 'VIYZ' - 'VIZX' - 'VIZY' - 'VIZZ' - 'PEXX' - External Pressure 'PEXY' - 'PEXZ' - 'PEYX' - 'PEYY' - 'PEYZ' - 'PEZX' - 'PEZY' - 'PEZZ' - 'PIXX' - Internal Pressure 'PIXY' - 'PIXZ' - 'PIYX' - 'PIYY' - 'PIYZ' - 'PIZX' - 'PIZY' - 'PIZZ' - Examples: 1. Save the structure with a lower NOE restraint energy. READ COOR CARD UNIT 1 ! Read the first structure READ COOR CARD COMP UNIT 2 ! Read the second structure ENERGY ! Compute energy of first structure SET 1 ?NOE ! save the NOE energy value ENERGY COMP ! Compute the energy of the second structure IF ?NOE LT @1 COOR COPY ! replace first structure if second has ! a lower energy. 2. Write some energy values when saving coordinates .... COOR ORIE RMS MASS ENERGY OPEN WRITE CARD UNIT 22 NAME RESULT.CRD WRITE COOR CARD UNIT 22 * Final coordinates * energy=?ENER and electrostatic energy=?ELEC * mass weighted rms deviation from xray structure is ?RMS *
NIH/DCRT/Laboratory for Structural Biology
FDA/CBER/OVRR Biophysics Laboratory