CHARMM c30b1 energy.doc



File: Energy ]-[ Node: Top
Up: (commands.doc) -=- Next: Description\n

             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
* Substitution::        Command line energy substitution parameters
* Running Average::     ESTATS command usage
* SPASIBA::             SPASIBA spectoscopic force field


File: Energy ]-[ Node: Description
Up: Top -=- Next: Skipe -=- Previous: Top\n

                        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,,) and 
*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 ] ]
               [ NOPRint            ]

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).


File: Energy ]-[ Node: Skipe
Up: Top -=- Next: Interaction -=- Previous: Description\n

                      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



File: Energy ]-[ Node: Interaction
Up: Top -=- Next: Fast -=- Previous: Skipe\n

                    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. [T.Lazaridis, July 1999: Now INTE can work with the GROUP option]

        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.


File: Energy ]-[ Node: Fast
Up: Top -=- Next: Needs -=- Previous: Interaction\n

[SYNTAX FASTer ]

FASTer {integer}
       {OFF    }
       {ON     }
       {DEFAult}
       {SCALar } ! for testing only
       {VECTor } ! for testing only
       {CRAYvec } ! Use parallel code designed for a CRAY
       {PARVec  } ! Use parallel/vector code best SMP machines and Convex

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]    
                    


File: Energy ]-[ Node: Needs
Up: Top -=- Next: Optional -=- Previous: Fast\n

        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).


File: Energy ]-[ Node: Optional
Up: Top -=- Previous: Needs -=- Next: Substitution\n

        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.



File: Energy ]-[ Node: Substitution
Up: Top -=- Previous: Optional -=- Next: Running Average\n

      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,
             or unless specified with the VOLUme keyword.
 '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
 'STRB'  - Strech-Bend coupling energy (MMFF)
 'OOPL'  - Out-off-plane energy (MMFF)
 'VDW '  - van der waal energy
 'ELEC'  - electrostatic energy
 'HBON'  - hydrogen bonding energy
 'USER'  - user supplied energy term
 'HARM'  - harmonic positional restraint energy
 'CDIH'  - dihedral restraint energy
 'CIC '  - internal coordinate restraint energy
 'CDRO'  - droplet restraint energy (approx const press)
 'NOE'   - general distance restraint energy (for NOE)
 'SBOU'  - solvent boundary lookup table energy
 'IMNB'  - primary-image van der waal energy
 'IMEL'  - primary-image electrostatic energy
 'IMHB'  - primary-image hydrogen bond energy
 'EXTE'  - extended electrostatic energy
 'EWKS'  - Ewald k-space sum energy term
 'EWSE'  - Ewald self energy term
 'RXNF'  - reaction field electrostatic energy
 'ST2'   - ST2 water-water energy
 'IMST'  - primary-image ST2 water-water energy
 'TSM'   - TMS free energy term
 'QMEL'  - Quantum (QM) energy with QM/MM electrostatics
 'QMVD'  - Quantum (QM/MM) van der Waal term
 'ASP'   - Atomic solvation parameter (surface) energy
 'EHAR'  - Restraint term for Implicit Euler integration
 'GEO '  - Mean-Field-Potential energy term
 'MDIP'  - Dipole Mean-Field-Potential energy term
 'PRMS'  - Replica/Path RMS deviation energy 
 'PANG'  - Replica/Path RMS angle deviation energy 
 'SSBP'  - ???????  (undocumented)
 'BK4D'  - 4-D energy
 'SHEL'  - ???????  (undocumented)
 'RESD'  - Restrained Distance energy
 'SHAP'  - Shape restraint energy
 'PULL'  - Pulling force energy
 'POLA'  - Polarizable water energy
 'DMC '  - Distance map restraint energy
 'RGY '  - Radius of Gyration restraint energy
 'EWEX'  - Ewald exclusion correction energy
 'EWQC'  - Ewald total charge correction energy
 'EWUT'  - Ewald utility energy term (for misc. corrections)

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
*



File: Energy ]-[ Node: Running Average
Up: Top -=- Next: SPASIBA -=- Previous: Substitution\n

                    Running Energy Averages (ESTATS)

     The ESTATS command is a basic statistical facility that allows the 
calculation and manipulation of the mean and variance of the potential
energy and its components over a number of potential energy calculations,
without the need for writing out trajectories or coordinate files--i.e.
the calculations are done "on the fly." ESTATS can be used in dynamics runs 
or in other sampling procedures that result in serial calls to the ENERGY 
subroutine.  The facility will collect data points at specified sampling 
intervals along a collection run for a specified step length and calculate
the running statistics.  An initial portion of the collection run may be skipped
(e.g. for eliminating the equilibration period from the statistics during 
dynamics). Results may be written either to standard output or to a file.
The facility will, if requested, "variable-ize" the calculated averages, i.e. 
allow assignment of the values to CHARMM script variables.  The facility can
also write the individual potential energy values to a file.
     
Syntax:
     
ESTAts   [LENGTH <integer>] [SKIP <integer>] [IPRF <integer>]
         [NPRI <integer>] [IUNW <integer>] [NEPR <integer>] 
         [IUPE <integer>]
         [UPLM <real>] [LOLM <real>] [FRPI] [VARI]

         [BOND] [ANGLe] [UREY-Bradley] [DIHEdral] [IMPRoper]
         [VDWaals] [ELECtrostatics] [HBONding] [USER] 
         [SBOUnd]  [ASP]
         
         LENGth  length of trajectory (number of total energy calculations)
            from which sampling is to take place  (default 0).
         SKIP  specifies a length of energy data points (calls to ENERGY) 
            after which the data collection is to begin (default 0).
         IPRFreq  specifies the frequency with which data points will be 
            collected (i.e. every IPRFrequency energy calculations).
         [BOND], [ANGLE], etc.
            the energy term keywords specify which components of the potential
            energy are to have their statistics calculated.  HBONding is
            the hydrogen bonding energy; USER is the user-defined energy;
            SBOUnd is the solvent boundary potential; ASP is the implicit
            solvation energy (e.g. from eef1).  Statistics on the total
            potential energy are always calculated.
         IUNWrite  fortran unit onto which statistics are to be written 
            (default is no printing)
         NPRInt  period for writing the energy statistics to standard output
            (default is no printing)
         IUPE  fortran unit onto which the potential energies are to be 
            written (default is no printing)
         NEPRint period for writing potential energies
         UPLM  limit above which an energy value will be discarded from the 
            statistics (default 99999999).
         LOLM  limit below which an energy value will be discarded from the
            statistics (default -99999999).
         FPRI keyword specifying the final statistics are to be written to 
            standard output at the end of the collection
         STOP  stops the data collection and prints current statistics
         VARI keyword specifying that values of averages and variances will
            be assignable to CHARMM script variables.
            The values can be accessed as follows:
            ?AENE,?VENE  mean and variance for potential energy
            ?ABON,?VBON  mean and variance for bonds
            ?AANG,?VANG  for angles
            ?AURE,?VURE  for Urey-Bradley terms
            ?ADIH,?VDIH  for dihedral terms
            ?AIMP,?VIMP  for improper dihedral terms
            ?AVDW,?VVDW  for van der Waals 
            ?AELE,?VELE  for electrostatics
            ?AHBO,?VHBO  for hydrogen bond terms
            ?AUSE,?VUSE  for user energy
            ?ASBO,?VSBO  for solvent boundary potential (sbound)
            ?AASP,?VASP  for solvation term
            
Note that ALL component energy terms for which statistics are being calculated
must be in the proper range (>LOLM and <UPLM) in order for a given data point
to be included.  Discarded data points will result in statistics that are
based on less than LENGTH data points.  The number of discarded data points
will be printed to standard output with the final statistics at the end of
the collection.

EXAMPLE:

ESTAts LENGTH 1000000 SKIP 100000 IPRFreq 5 NPRINT -1 FPRInt -
VDW ELEC BOND ANGL IMPR SBOU DIHE -
IUNWrite 11 NUPRint 10000 NEPR 1000 IUPE 10

This specifies that the statistics are to be done on 180,000 data points
(1,000,000 - 100,000)/5.  Statistics will be done on the specified
energy terms in addition to the potential energy.  Statistics will be
written every 10,000 steps to unit 11 and the potential energies
will be written every 1000 steps to unit 10.  No printing will be 
done to standard output (NPRINT -1) except for the final statistics
(FPRInt).

This statistics file is written out according to the following format:

EPOT             10    -1912.29336237      226.63620520
BOND             10      212.58550818       91.35922427
ANGL             10      299.99516787       95.65303762
UREY             10       39.09373234       18.64669506

The first column indicates the energy term.
The second column indicates the number of data points included in the
  calculations (number of values over which statistics are taken).
The third column gives the average and the fourth column gives the 
  fluctuation (standard deviation).

(Differences between ESTATS and statistics calculated in standard 
dynamics:
1) In ESTATS, the denominator in the standard deviation formula is
 (N-1)^(1/2), where N is the number of data points.  In dynamics, 
 the denominator is N^(1/2).
2) In dynamics, the initial energy terms are considered step "0" and 
  not included in the statistics; hence for a direct comparison, it
  is necessary to specify SKIP 1 in the ESTATS command and increase
  the LENGTH of the collection by 1. 



File: Energy ]-[ Node: SPASIBA
Up: Top -=- Next: Top -=- Previous: Running Average\n


                           SPASIBA Force Field


     The SPASIBA force field is derived from spectroscopic studies.  In the
current implementation, the van der Waals and electrostatic interactions
of the CHARMM force field are combined with the Urey-Bradley Shimanouchi terms
for bond stretching, valence angle bending, torsional and improper torsional
internal terms derived from spectroscopic data such as vibrational wavenumbers,
dipole moments, rotation barriers for the SPASIBA force field.

     See for a complete description of the SPASIBA (modified UBFF) the
following publications:

1. P. Derreumaux and G. Vergoten, "A new spectroscopic molecular mechanics
   force field. Parameters for proteins," J.Chem.Phys. 102(21), 8586-8605
   (1995).

2. P. Lagant, D. Nolde, R. Stote, G. Vergoten and M. Karplus, "Increasing
   normal mode analysis accuracy:  The SPASIBA spectroscopic force field
   and the CHARMM program," in preparation.

[Note] CHARMM must be compiled with the SPAS keyword in the pref.dat file
       to use the SPASIAB force field.

     The parameters for the SPASIBA force field are read in via a modified
parameter file.  The specific parameters are the K (Urey-Bradley) for bonds,
H and F/F' terms, the Kappa (internal tensions for a tetrahedral atom),
L(CH2), trans and gauche force constants between valence angles.
The SPAS keyword in the parameter file signals the use of the SPASIBA force
field in the appropriately compiled CHARMM version.

     Energy, first and second derivatives are calculated by the SPASIBA
force field.  The evaluation of the a,b,c,d coefficients (kappa), L and
Trans-Gauche terms is performed at each energy calculation (as they are
conformational dependent).

     The topology and psf files are unchanged relative to the standard
CHARMM format.  However, the parameter file differs from the standard format
in the following way.  The general parameter file spa_all22_prot.inp is
derived from the CHARMM all22_prot(mod).inp.  We added the F, kappa, LCH2
and TG terms  

     An example for a small molecule (ethanediol) is given below:

             HC  HM   
             |   |
          HC CT- C9 -S -HS
             |   |
             HC  HM

     -----------------------------------------------------------
* - parameter file par.etsh 
* Parameters from SPASIBA force field
* Test 
* To be used with top.etsh topology file
*

BOND
CT  HC     320.   1.11                ! K, r0    
CT  C9     140.0  1.53                !K in kcal mol-1 Angstroms-2 ,r0 in A
C9  S      139.   1.81 
C9  HM     302.5  1.11
S   HS     278.   1.336 

THETAS
HM  C9  HM     30.6   108.5    10.07    0.        !H, teta0, F
HC  CT  HC     29.    107.7    10.5     0.        !H kcalmol-1 rad-2
C9  CT  HC     14.9   109.4    67.9     0.        !F kcalmol-1 A-2
CT  C9  S      24.0   114.7    14.0     0.        !teta0 (deg)
CT  C9  HM     14.1   109.4    69.4     0.
S   C9  HM     13.0   108.9    55.0     0.
C9  S   HS     53.3    96.0     0.      0.

PHI
HM  CT  S   HS      0.185   1      0.
CT  C9  S   HS      0.185   1      0.
HM  C9  S   HS      0.185   1      0.   
X   CT  C9  X       0.16    1      0. 

NONBONDED  NBXMOD 5  ATOM CDIEL SHIFT VATOM VDISTANCE VSHIFT -
     CUTNB 5.5 CTOFNB 4.5 CTONNB 4. EPS 1.0  E14FAC 0.5  WMIN 1.5
!
HC    0.00    -0.01       0.77       0.00    -0.005      0.77
HM    0.00    -0.01       0.77       0.00    -0.005      0.77
CT    0.00    -0.06       0.900      0.00    -0.03       0.900
C9    0.00    -0.06       0.900      0.00    -0.03       0.900
S     0.00    -0.20       1.000      0.      -0.10       1.00
HS    0.00    -0.02       0.50       0.      -0.01       0.5 

KAPPA 
CT  C9  HC  HC  HC     -3.6         !internal tension occuring in the set of 6 
C9  CT  S   HM  HM     -3.6         !valence angles around a tetrahedral atom

LCH2                                !  HM         HM
C9  X   X   HM  HM      2.877       !       C9 
                                    !   X         X
14TG
HC  CT  C9  HM       15.100    -4.0 !   HC-CT-C9-HM 
HM  C9  S   HS        0.1      -0.1

END


The input file remains unchanged relative to standard CHARMM input files.
The Spasiba force field is activated via the SPAS keyword in the parameter file

*  ETHSH  champ de forces   SPASIBA  
*
 
! Open and read topology file
OPEN READ UNIT 11 CARD NAME etsh.top
READ RTF UNIT 11 CARD
CLOSE UNIT 11


! Open and read parameter file
OPEN READ UNIT 12 CARD NAME etsh.par
READ PARAMETERS UNIT 12 CARD ! PRINT
CLOSE UNIT 12

OPEN READ UNIT 16 CARD NAME etsh.psf
READ PSF CARD UNIT 16
CLOSE UNIT 16

OPEN READ UNIT 14 CARD NAME etsh.crd
READ COORDINATES CARD UNIT 14
CLOSE UNIT 14 

ENERGY 
!MINIMIZE ABNR INBFRQ 10 NSTEP 50000 STEP 0.02 TOLG 1.0E-04
MINIMIZE NRAP NSTEP 1500 STEP 0.02 INBFRQ 10 TOLG 1.0E-06

!OPEN WRITE UNIT 13 CARD NAME etshabnr.crd
OPEN WRITE UNIT 13 CARD NAME etshnrap.crd
WRITE COORDINATES CARD UNIT 13
* Output coordinates  (CH3-CH2-SH) 
*

VIBRAN NMOD 27 
DIAG 
PED MODE 1 THRU 27 TEMP 298.0 TOL 0.01

STOP


Bibliography on the SPASIBA Force Field
---------------------------------------

1. Carboxylic acids: (acetic ,pivalic, succinic, adipic ,L-glutamic acids):
   M. Chhiba, P. Derreumaux and G. Vergoten, J .Mol. Struct., 317, 171-184
   (1994)

2. Linear alkenes:
   M. Chhiba and G. Vergoten  J. Mol. Struct., 326, 35-58 (1994)

3. Aliphatic esters: 
   F. Tristram, V. Durier and G. Vergoten, J. Mol. Struct., 377, 47-56 (1996)

4. Aliphatic alcohols: 
   F. Tristram, V. Durier and G. Vergoten, J. Mol. Struct., 378, 249-256 (1996)

5. Esters:
   M. Chhiba, F. Tristram and G. Vergoten, J. Mol. Struct., 405, 113-122 (1997)

6. alpha-D-Glucose: 
   V. Durier, F. Tristram and G. Vergoten, J. Mol. Struct.(Theochem), 395-396,
   81-90 (1997)

CHARMM .doc Homepage


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