CHARMM c30b1 cons.doc

File: Cons ]-[ Node: Top
Up: (commands.doc) -=- Next: Harmonic Atom\n


        The following forms of constraints are available in CHARMM:

* Menu:                          command

* Harmonic Atom::       "CONS HARM" Hold atoms in place
* Dihedral::            "CONS DIHE" Hold dihedrals near selected values
* Internal Coord::      "CONS IC"   Holds bonds, angles and
                                    dihedrals near table values
* Quartic Droplet::     "CONS DROP" Puts the entire molecule in a cage
                                    about the center of mass
* RMSD restraints:      "RMSD"      Holds atoms in place relative to
                                    reference structure/structures
* Fixed Atom::          "CONS FIX"  Fix atoms rigidly (sets the IMOVE array)
* Center of Mass::      "CONS HMCM" Constrain center of mass of selected atoms
* SHAKE::               "SHAKE"     Fix bond lengths during dynamics.
* NOE::                 "NOE"       Impose distance restraints from NOE data
* Restrained Distances:: "RESD"     Impose general distance restraints
* External Forces::     "PULL"      Impose externally applied (pulling) force
* Rg/RMSD restraint::   "RGYR"      Impose radius of gyration or rmsd restraint
* Distance Matrix restraint:: "DMCO" Impose a distance matrix restraint
* Sbound: (sbound.doc).      Solvent boundary potential

File: Cons ]-[ Node: Harmonic Atom
Up: Top -=- Next: Dihedral -=- Previous: Top\n

                        Holding atoms in place



CONS HARMonic {[ABSOlute] absolute-specs force-const-spec   coordinate-spec }
              { BESTfit   bestfit-specs  force-const-spec   coordinate-spec }
              { RELAtive  bestfit-specs  force-const-spec 2nd-atom-selection}
              { CLEAr                                                       }

force-const-spec ::=  { FORCE real } atom-selection [MASS]
                      { WEIGhting  }

absolute-specs ::=  [EXPOnent int] [XSCAle real] [YSCAle real] [ZSCAle real]

bestfit-specs  ::=  [ NOROtation ] [ NOTRanslation ]

coordinate-spec::= { [MAIN] } 
                   {  COMP  } 
                   {  KEEP  } 

        The potential energy has a harmonic restraint term which allows
one to prevent large motions of individual atoms.  There are three forms
for this restraint, ABSOLUTE, BESTFIT, and RELATIVE.  It is possible to
combine multiple restraints in one energy calculation, but no atom may
participate in more than one harmonic restraint set.

ABSOLUTE positional restraints.

Absolute positional restraints specify a location a cartesian space where
an atom must remain proximal.  This is the original positional restraint
in CHARMM.  The form for this potential is as follows for coordinates:

  EC = sum over selected atoms of k(i)* [mass(i)] * (x(i)-refx(i))**exponent

where refx is a reference set of coordinates.  If MASS is specified in
the command line, then k is multiplied by the mass of the atom
resulting in a natural frequency of oscillation for the restraint 
of sqrt(k) in AKMA units.   An atom restrained with MASS FORCE 1.0
will oscillate at 8 cycles/picosecond if free of other interactions.
For most operations involving harmonic restraints, mass weighting is
recommended. There are three reasons for this. First, the results obtained
will be similar regardless of what atom representation is used
(extended vs. explicit) for hydrogen atoms. Second, Hydrogen atoms
are allowed greater relative freedom if present. And third, The character
of the normal modes of a molecule are unperturbed with mass weighting
(essential if normal modes or low frequency motions are of interest).

        Note, there is no longer a prefactor of 0.5 on the force
constant specification. This is appropriate in that exponent values
other than "2" are allowed. This differs from the earlier versions of
CHARMM (up to version 16).

       The restraint force constant can be set to any positive
value (specified by the FORCE keyword followed by the desired value).
The force constants may also be obtained from the weight array, in
which case the FORCe keyword is not read.  When using this option,
a negative values may be used for some atoms, however, the total weight
must be positive. 

      The reference coordinates can be the current set at the point when
restraints are specified (the default) or a set can be the comparison
set (COMP keyword).  When multiple CONS HARM commands are used, the
KEEP option preserves the reference coordinates from the previous
restraints.  This is useful in cases where the force constant is
to be modified, but no other changes are desired.

      The variables XSCAle, YSCAle, and ZSCAle are global scale factors
for ABSOLUTE harmonic restraint terms.  The default scale factor is 1.0 for
all terms.  If multiple harmonic restraint sets are used, they may have
different scale factors.  The RELATIVE and BESTFIT types do not allow a
scale factor at present.

Example:   CONS HARM FORCE 1.0 MASS SELE atom * * CA  END  COMP

This command harmonically restrains all alpha-carbons to the current
positions in the comparison coordinates with a force constant of
24 Kcals/mol/A**2 (assuming a mass of 12).

BESTFIT positional restraints

      This restraint is similar to the absolute restraints except that
the reference coordinates are (implicitly) rotated and translated so
as to bestfit the selected atoms.  This best fit is done in a manner that
minimizes the restraint energy.  Due to the nature of the best fit, this
restraint term does not add any net force or torque to the system.

Note 1: An exponent may not be specified (it is set to 2)
Note 2: Global scale factors do not apply
Note 3: At present, there is no Hessian code for this restrain
Note 4: There is no energy partition (ANAL command) for this restraint


These commands will restrain segments A and B to the geometries they have
in the current comparison coordinates.  This restraint will not add a net
force or torque to the system (unlike the ABSOLUTE restraint type).  Segments
A and B can move (rotate/translate) independently with no change in the
restraint energy.  See also CONS RMSD below.

RELATIVE positional restraints.

      The relative positional restraints are similar to the bestfit
restraint except there is no reference coordinates.  In this case,
one part of the system is restrained to have the same shape as another
part of the system.  The two restraint sets are implicitly bestfit
by an optimal rotation/translation (minimizing the restraint energy).
Both sets of atoms 

Note 1: An exponent may not be specified (it is set to 2)
Note 2: Global scale factors do not apply
Note 3: At present, there is no Hessian code for this restrain
Note 4: There is no energy partition (ANAL command) for this restraint
Note 5: The atoms of the two sets are matched on-to-one in sequential order.
Note 6: If the two sets do not have the same number of atoms, an
            error will be issued and the set lists will be truncated.
Note 7: Both sets must be specified and must not use set number 1.


This command will force two replicas (A1 and A2) to have the same coordinates
based on the values in the weighting array (as best fit weights).


This command will force two segments (A and B) to have the same shape, but
they may have very different locations and orientations.  Atoms are matched
one-to-one by selected atom number.


      It is important to understand some aspects of how the restraints are
set in order to get the most flexibility out of this command. When CHARMM
is loaded, each atom has associated with it a harmonic force constant
initially set to zero. Each call to the CONS HARM command changes the value
of this constant for only those atoms specified.  When this command is
invoked with an atom selection (and KEPP is not specified), only the
reference coordinates (XREF,YREF,ZREF) for selected atoms are modified.

Each atom may participate in AT MOST one harmonic restraint term.
This is a coding limitation designed to maximize compatibility with older
CHARMM scripts (i.e. doing a series of minimizations with a decreasing
series of force constants).  This could be easily modified with a bit of
work to increase the capability (at the expense of script compatibility).

      When multiple restraint sets are used, it is important to note
that all selections should be exclusive.  When they are not exclusive,
then atoms will be assigned to the restraint of the most recent
CONS HARM command which selected that atom.  In other words, the restraint
set number is an atomic property.  If restraint sets are broken up, then an
error message will be issued.  If an entire set is replaced, then the new
restraint replaces the old one (without a warning message).


        The harmonic restraints may no longer be read and written to files.
The PRINT command still functions for harmonic restraints for information.
To examine or modify the internal harmonic restrain data, the SCALar
command (arrays: CONStraints,XREF,YREF, and ZREF) may be used
(see *Note scalar::(chmdoc/scalar) ).   In addition, one may look at
the contributions to the energy in detail using the ANALysis command,
see *note anal:(chmdoc/analys).


File: Cons ]-[ Node: Dihedral
Up: Top -=- Next: Internal Coord -=- Previous: Harmonic Atom\n

                  Holding dihedrals near selected values

        Using this form of the CONS command, one may put restraints on
the dihedral angles formed by sets of any four atoms. The improper
torsion potential is used to maintain said angles.

        The command for setting the dihedral restraints is as follows:



CONS DIHEdral [BYNUM int int int int] [FORCe real] [MIN real] [PERIod int] 
              [   atom-selection    ]              [ COMP   ]   [WIDTh real]
              [     4X(atom-spec)   ]              [ MAIN   ]


Syntactic ordering:  DIHE or CLDH must follow CONS, and FORCE, MIN and
PERIod  must follow DIHE.

where:      atom-spec ::= { segid resid iupac }
                          { resnumber iupac   }

        DIHEdral adds a torsion angle to the list of restrained angles
using the specified atoms, force constant, minimum and periodicity.
If an atom selection is used, then the first 4 selected atoms (in order)
will define the dihedral angle.  If either MAIN or COMP is specified and
[MIN real] is not, then the minimum angle value will be determined by
the current dihedral angle value in the corresponding coordinate set.

If the PERIodicity is zero (improper type), then the force constant
has units of kcal/mol/radian/radian, else it has units of kcal/mol.

   Ecdih = FORCE * max(0, abs( phi - MIN*pi/180) - WIDTH)**2  [ PERIod = 0 ]
   Ecdih = FORCE * (1-cos( PERIod* (phi - MIN*pi/180 )) )     [ PERIod > 0 ]

CLDH clears the list of restrained dihedrals so that different angles
or new restraint parameters can be specified.

Other commands:

        The PRINT CONS command, see *note print:(io.doc)print,
will work for restraints.

File: Cons ]-[ Node: Internal Coord
Up: Top -=- Next: Quartic Droplet -=- Previous: Dihedral\n

                Holding Internal Coordinates near selected values



CONStraint IC  [BOND real [EXPOnent integer] [UPPEr]]
                   [ANGLe real] [DIHEdral real] [IMPRoper real]

        Using this form of the CONS command, one may put restraints on
any internal coordinate. For this energy term, the IC table is
used.  At each energy call, the reference (zero-force) value of each IC
is set to the value currently in the IC table.
All nonzero bond entries are restrained with the bond constant,
using the optional EXPOnent (default 2) in the potential K*(S-S0)**EXPOnent.
Second derivatives are currently supported only with EXPOnent=2.
The angle, dihedral, and improper terms are only harmonic.
The DIHEdral term only applies to IC's of normal type, and the
IMPRoper term only applies to the improper IC type (those with a "*")

If UPPEr is specified the reference bond length is taken as an upper
limit and the restraint potential is applied only if S>S0; this is
intended for use with distance restraints from NMR NOE data.
All nonzero angle entries are restrained with the angle constant. All
dihedrals are restrained with the dihedral constant using the improper
dihedral energy potential. If any IC entry contains an undefined atom
(zeroes), then the associated bonds,angles, and dihedral will not be

      The force constant has units of kcal/mol/radian/radian for both
angle and dihedral restraints.  The bond force constant has units of

        This restraint term is very flexible in that the user may
chose which bonds... to restrain by editing an IC table. The major
drawback is that all bonds must have the same force constant. The same is
true for angles and dihedrals. By listing some IC's several times, the
effective force constant is increased. Also, if only angle restraints are
desired, then the bond and dihedral constants can be set to zero eliminating
their contribution.

File: Cons ]-[ Node: Quartic Droplet
Up: Top -=- Next: Fixed Atom -=- Previous: Internal Coord\n

                        The Quartic Droplet Potential


CONStraint DROPlet [FORCe real] [EXPOnent integer] [NOMAss]

        This restraint term is designed to put the entire molecule
in a cage. Is is based on the center of mass (or center of geometry if
NOMAss is specified) so that no net force or torque is introduced by this
restraint term. The potential function is;

        Edroplet= FORC* sum over atoms (( r-rcm )**EXPO )*mass(i))

File:Cons ]-[ Node: Fixed Atom
Up:Top -=- Next: Center of Mass -=- Previous: Quartic Droplet\n

                    How to fix atoms rigidly in place


Syntax: CONS FIX atom-selection-spec { [PURG]                     }
                                     { [BOND] [THET] [PHI] [IMPH] }

        This command will fix all selected atoms and unfix all non-selected
atoms.  For example, the command; CONS FIX SELE NONE END 
will remove all fixing of atoms (except for lonepairs).

        This command fixes atoms in place by setting flags in an array
(IMOVE) which tells the minimization and dynamics alogrithms which atoms
are free to move. If atoms are fixed, it is possible to save
computer time by not calculating energy terms which involve only fixed
atoms. The nonbond and hydrogen bond algorithms in CHARMM check IMOVE
and delete pairs of atoms that are fixed in place from the nbond and
hbond lists respectively. In addition the PURG or individual energy term
options specified with the CONS FIX command allow all or some of the
internal coordinate energies associated with fixed atoms to be deleted.
Interactions between fixed and moving atoms are maintained.

*** NOTE *** because some energy terms are deleted from fixed systems,
the total energy calculated with fixed atoms will be different from the
total energy of the same system with all atoms free. The forces on the
moveable atoms will however be identical.  The purpose of this feature is
to remove the computational cost of energy terms that do not change for
simulations where a large fraction of the atoms are fixed.  It is not
recommended for any other purpose.

        The way CHARMM keeps track of fixed atoms is by the IMOVE array
in the PSF. The IMOVE array is 0 if the atom is free to move, and has
the value 1 if it is fixed.  A value of -1 indicates that this atom 
is a lonepair. 

***** WARNING ***** The purge options modify the PSF. The effects of
this command cannot be undone by the subsequent releasing of atoms.

***** WARNING ***** The fixing of atoms does not work for constant
pressure simulations.

File:Cons ]-[ Node: Center of Mass
Up:Top -=- Next: SHAKE -=- Previous: Fixed Atom\n

                  Constrain centers of mass for selected atoms


Syntax: CONS HMCM FORCe real [WEIGhting] reference-spec atom-selection
        reference-spec ::= REFX real REFY real REFZ real
        This command will harmonically restrain centers of mass from the
selected atoms to the absolute reference point specified with REFX, REFY and
REFZ. The force constant of the harmonic potential is set with the FORCe 
keyword. Mass weighting is switched off by default but can be selected
by using the WEIG key. 

        The primary use of this command is during the reconstruction of
all-atom representations from low resolution models with virtual particles
at side chain centers.


        CONS HMCM FORCE 50.0 WEIG REFX 10.4 REFY 12.1 REFZ 1.3 -
             SELECT RESID 21 .AND. .NOT. -
             ( TYPE H* .OR. TYPE N .OR. TYPE C .OR. TYPE O ) -

        This will create a harmonic restraint with a force constant
of 50 kcal/mol that holds the side chain center of mass at residue 21
of a protein near (10.4, 12.1, 1.3).

File: Cons ]-[ Node: SHAKE
Up: Top -=- Next: NOE -=- Previous: Center of Mass\n

                Fixing bond lengths or angles during dynamics.

        SHAKE is a method of fixing bond lengths and, optionally, bond
angles during dynamics, minimization (not ABNR and Newton-Raphson methods),
coordinate modification (COOR SHAKe command), and vibrational analysis
(explore command). The method was brought to CHARMM by Wilfred Van
Gunsteren (WFVG), and is referenced in J. Comp. Phys. 23:327 (1977).
When hydrogens are present in a structure, it will allow a two-fold
increase in the dynamics step size if SHAKE is used on the bonds.

        To use SHAKE, one specifies the SHAKE command before any
SHAKE constraints usage. The SHAKE command has the following syntax:

[SYNTAX SHAKe constraints]

SHAKE { OFF                                               } 
      { shake-opt  fast-opt  2x(atom-selection) [NOREset] }

shake-opt:== [BONH]  { [MAIN]     } [TOL real] [MXITer integer]
             [BOND]  { COMP       }
             [ANGH]  { PARAmeters }        [SHKScale real] 

fast-opt:== { [ FAST  [ WATEr water-resn ] ] }
            {   NOFAst                       }

        BONH specifies that all bonds involving hydrogens are to be
fixed. BOND specifies all bonds. ANGH specifies that all angles
involving hydrogen must be fixed. ANGL specifies that all angles must be
shaken. BOND is implied if any angles are fixed, otherwise, only the 1-3
distances would be fixed. Coordinates must be read in before the SHAKE
command is issued, unless the PARAmeter option is specified.

      SHAKE constraints are applied only for atom pairs where one atom
is in the first atom selection and one atom in the second atom selection.
The default atom selection is ALL for both sets.  

        TOL specifies the allowed relative deviations from the reference
values (default: 10**-10). MXITer is the maximum number of iterations
SHAKE tries before giving up (default: 500).

        When the SHAKE command is used, it will check that there are
degrees of freedom available for all atoms to satisfy all their
constraints. Angles cannot be fixed with SHAKE if one has explicit
hydrogen arginines in the structure as the CZ carbon has too many
constraints. This is a general problem for any structure which has too
many branches close together.

        SHAKE is not recommended for fixing angles. The algorithm
converges very slowly in the case where one has three angles centered on
a tetravalent atom and the constraints are satisfiable only using out of
plane motions.

        The use of SHAKE modifies the output of the dynamics command.
The number appearing to the right of the step number is the number of
iterations SHAKE required to satisfy all the constraints. This number
should generally be small.

        When ST2's are present, SHAKE constraints are automatically
applied for the O-H bonds and H-O-H angles.

        There is a PARAmeter option for the SHAKe command. This option
causes the shake bond distances to be found from the parameter table
rather than from the current set of coordinates. This option is
NOT compatible with the use on angle SHAKE constraints, and it will
give an error if this is tried.

        With these commands, the bond energy may be zeroed without
any minimization with the command sequence;

[SYNTAX SHAKe FAST constraints]
SHAKe FAST [WATEr SELEct water_selection END] [OLDWatershake]
         [ MXITer <iterations> TOL <tolerance> ] [PARAmeter] [COMP]
   This command specifies the use of the new vector/parallel and analagous
   scalar fast SHAKE constraint routines (implemented Aug 2000). 
   Certain assumptions are made when
   this command is issued: The only bonds involved are between heavy atoms
   and hydrogens, except for water molecules included in the 
   WATEr selection ... end sub-command.
   This selection is used to indicate the water molecules
   that have an H-H bond. It is assumed that the selection will include
   all atoms in the water molecule and that said molecule contains exactly
   two X-H bonds and one H-H bond where X is any heavy atom.  Testing
   for "hydrogen-ness" is done via the CHARMm hydrog() function which
   makes it's choice based on atomic mass.  The prefered selection is
   through the use of the RESNAME selection specifier, eg:


   By default, water molecules selected with the WATEr sub-command will
   be constrained via the use of a special water-SHAKE routine which
   uses the direct inversion method. This algorithm is from 25 to 30 %
   faster than the normal iterative, scalar SHAKE routine. For the rest
   of the heavy atom -hydrogen bonds, a vector/parallel version of
   the original SHAKE routine is used.  This is about 5X the scalar SHAKE.
   If the optional keyword OLDWatershake is used, the vector/parallel
   (not the watershake) routines are used.

   The rest of the keywords are the same as in the original SHAKE command.
   Note: that FAST has to be the second word in command line.

File: Cons ]-[ Node: NOE
Up: Top -=- Previous: SHAKE -=- Next: Restrained Distances\n

[SYNTAX NOE constraints]

            Invoke the module

            Reset all NOE restraint lists. This command clears all
            existing NOE restraints. Resets scale factor to 1.0

            Turn on the restraint between a given atom specified
            by ASSIgn and a point specified by CNOX, CNOY and CNOZ
            intead of a restraint between two atoms.
            The use of this restraint is desirable for docking, 
            and loop refinements.

 ASSIgn [KMIN real] [RMIN real] [KMAX real] [RMAX real] [FMAX real]
                 [TCON real] [REXP real] {2X(atom_selection)}
                {[CNOX real] [CNOY real] [CNOZ real] 1X(atom selection) } 

            Assign a restraining potential between the atoms the
            first selection and the atoms of the second selection.

            Where multiple atoms are selected, 
                     R = [ average( Rij**(1/REXP) ) ]**REXP
            where (i) runs over the first atom selection and (j)
            runs over the second atom selection.
            The default REXP value is 1.0 (a simple average).
            An REXP value of 3.0 may be optimal for NOE averaging.

                   /  0.5*KMIN*(RAVE-RMIN)**2    R<RMIN
                 /    0.0                        RMIN<R<RMAX
                 \    0.5*KMAX*(RAVE-RMAX)**2    RMAX<RAVE<RLIM
                   \  FMAX*(RAVE-(RLIM+RMAX)/2)  RAVE>RLIM

                     RAVE=R                     TCON=0

                     RAVE=RRAVE**(-1/3)         TCON>0

            for initial conditions, RRAVE=RMAX**(-3)
            DELTA is the integration time step.  For minimization,
            the value is either 0.001ps or the previous simulation value.

            Where: RLIM = RMAX+FMAX/KMAX (the value of RAVE where the
                                          force equals FMAX)

            Defaults for each entry: KMIN=0.0, RMIN=0.0,
                                     KMAX=0.0, RMAX=9999.0, FMAX=9999.0
                                     TCON=0.0, REXP=1.0

Also, the old sytax is supported:

 ASSIgn rminvalue  minvariance  maxvariance  2X(atom_selection) 

For this format, KMAX=0.5*Kb*TEMP/(maxvariance**2)

   READ UNIT <integer>
            Reads restraint data structure from card
            file previously written.
   WRITe UNIT <integer> [ANAL]
            Writes out the restraint data in card format to a file on the
            specified unit. A CHARMM title should follow the command.
            SCALE are saved together with the lists in the NOE common block.
            The ANAL option will print out the distances and energy data
            computed with the current main coordinates.
   PRINT [ANAL [CUT real]]
            Same as the WRITe command except to the output file and slightly
            more user friendly form. A positive CUT value will list only
            those that have a distance that exceeds RMAX by more than DCUT.
   SCALe [real]
            Set the scale factor for the NOE energy and forces.
            Default value: 1.0
   TEMPerature real
            Specify the temperature for the old format.
            Return to main command parser.

No other commands (I/O or loops) are supported inside the NOE module.
Looping can be performed outside if necessary.  The units are Kcal/mol/A/A
for force constants and Angstroms for all distances.

EXAMPLE. Set up some NOE restraints for one strand of a DNA-hexamer
in a file to be streamed to from CHARMM.

! First clear the lists
! Since there are many identical atom pairs we use a loop
set 1 1
label loop
!   Sugar protons, same in all six sugars (don't pay any attention to
!     the numeric values)
    ASSIgn  SELE ATOM A @1 H1' END SELE ATOM A @1 H2'' END -
            KMIN 1.0 RMIN 2.7 KMAX 1.0 RMAX 3.0 FMAX 2.0
    ASSIgn  SELE ATOM A @1 H3' END SELE ATOM A @1 H2'' END -
            KMIN 1.0 RMIN 2.7 KMAX 1.0 RMAX 3.0 FMAX 2.0
incr 1 by 1
if 1 le 6 goto loop
! Now do some more specific things

   SCALE 3.0  ! Multiply all energies and forces by 3
   PRINT ANAL  ! See what we have so far
   PRINT ANAL CUT 2.0 ! list 

File: Cons ]-[ Node: Restrained Distances
Up: Top -=- Previous: NOE -=- Next: External\n

      Apply general restrained distances allowing multiple distances to
be specified.  This restraint term has been added to allow for facile
searching of a reaction coordinate, where the reaction coordinate is
estimated to be a linear combination of several distances.

      By Bernard R. Brooks - NIH - March, 1995

[SYNTAX Restrained Distances]

RESDistance [ RESEt ] [ SCALE real ] [ KVAL real  RVAL real [EVAL integer] -
       [ POSItive ] [ IVAL integer ]  repeat( real first-atom second-atom ) ]
       [ NEGAtive ]

   E = 1/EVAL *  Kval * Dref**EVAL


   Dref =  K1*R1**Ival + K2*R2**Ival + ... + Kn*Rn**Ival - Rval

   Where K1,K2,...Kn are the real values in the repeat section and
   R1,R2,...Rn are the distances between specified pair of atoms.

   RESEt       Reset the restraint lists. This command clears the
               existing restraints. Resets the scale factor to 1.0

   SCALe real  Set the scale factor for the energy and forces.
               Default value: 1.0

   POSITIVE    Include this restraint only when Dref is positive.
   NEGATIVE    Include this restraint only when Dref is negative.

If anything else is on the command line then a new restraint is added to the
list of distance restraints.
      KVAL real     The force harmonic constant
      RVAL real     The target distance
      IVAL integer  The exponent for individual distances.
      EVAL integer  The exponent (default 2). EVAL must be positive.
      repeat( real first-atom second-atom )
            The real value is a scale factor for the distance between
            the first and second specified atoms in the pair.

1. Create a reaction coordinate for QM/MM 
2. Set up some restraints to force three atoms to make an equilateral triangle.

!!! 1 !!!  Create a reaction coordinate for QM/MM 

OPEN WRITE FILE UNIT 22 name reaction.path
* trajectory of a minimized reaction path


SET 1 1
SET V -5.0

SKIP NONE             ! make sure all energy terms are enabled
RESDistance  RESET KVAL 2000.0  RVAL @v - 
   1.0   @atom1  @atom2    -1.0   @atom2  @atom3 

PRINT RESDistances    ! print a check of distances
TRAJ WRITE            ! write out the new minimized frame

SKIP RESD             ! turn off the restraint energy term
ENERGY                ! recompute the energy without restraints
WRITE TITLE UNIT 21   ! write out the current restraint distance and energy
* @V ?ENER

INCR 1 BY 1           ! increment the step counter
INCR V BY 0.25        ! increment the restraint value


!!! 2 !!! Make a water nearly an equilateral triangle

set atom1  WAT 1 O
set atom2  WAT 1 H1
set atom3  WAT 1 H2

RESDistance  RESEt

RESDistance  KVAL 1000.0  RVAL 0.0 - 
        1.0   @atom1  @atom2  -
        1.0   @atom1  @atom3  -
       -2.0   @atom2  @atom3
RESDistance  KVAL 1000.0  RVAL 0.0 - 
        1.0   @atom1  @atom2  -
       -2.0   @atom1  @atom3  -
        1.0   @atom2  @atom3
RESDistance  KVAL 1000.0  RVAL 0.0 - 
       -2.0   @atom1  @atom2  -
        1.0   @atom1  @atom3  -
        1.0   @atom2  @atom3

print resdistances
mini abnr nstep 200 nprint 10
print resdistances

!!! 3 !!! Prevent an atom from moving more than 20A from the others,
! but have no restraint energy when no distance is large.

set atom1 SOLV 1 OH2 
set atom2 SOLV 2 OH2  
set atom3 SOLV 3 OH2  
set atom4 SOLV 4 OH2  
set atom5 SOLV 5 OH2  

RESDistance  RESEt

RESDistance  KVAL 1.5E-12  RVAL 6.4E7 IVAL 6 POSITIVE - 
        1.0   @atom1  @atom2  -
        1.0   @atom1  @atom3  -
        1.0   @atom1  @atom4  -
        1.0   @atom1  @atom5  -
        1.0   @atom2  @atom3  -
        1.0   @atom2  @atom4  -
        1.0   @atom2  @atom5  -
        1.0   @atom3  @atom4  -
        1.0   @atom3  @atom5  -
        1.0   @atom4  @atom5   

print resdistances
mini abnr nstep 200 nprint 10
print resdistances

File: Cons ]-[ Node: External Forces
Up: Top -=- Previous: Restrained Distances -=- Next: Rg and RMSD restraint\n

[SYNTAX External Forces]

PULL { FORCe  <real>  }  XDIR <real> YDIR <real> ZDIR <real> [PERIod <real>]
     { EFIEld <real>  }
     { OFF            }
     { LIST           }
                         [WEIGht] atom-selection  

  A force will be applied in the specified direction on the selected atoms
  either as a constant:
  FORCe <value> specified in picoNewtons (pN)
  or oscillating in time: FORCe*COS(TWOPI*TIME/PERIod), FORCe <pN> PERIod <ps>
    time is counted from the start of the dynamcis run.
  The force due to an electrical EFIEld (V/m) (possibly also oscillating)
  may also be specified, in which case partial charges are taken from the psf
  and used to calculate the force.
  If WEIGht is specified the forces are multiplied by the wmain array.

  Each invocation of this command adds a set of forces to the previously
  defined set. 

  PULL OFF turns off all these forces.
  PULL LIST produces a listing. NB! Forces defined by PULL will move
  atoms in the specified direction, which is opposite to that listed
  by the forces from the COOR FORCE command.

File: Cons ]-[ Node: CONS RMSD
Up: Top -=- Previous: External Forces -=- Next: Rg/RMSD restraint\n

The RMSD restraint is useful to manipulate and control macromolecular
conformations.  The restraint is related to the CONS HARM BestFit, which sets
up harmonic restraints with respect of a reference structure.  However,
because all the reference data structure is stored in XREF, YREF, ZREF,
this command allows only a single bestfit restraint.  In addition, it allows
only a restraint to zero value of RMSD.  It is useful to allow multiple such
bestfit RMSD restraint to progress from one conformation to a second
conformation of a molecular system.  The new command CONS RMSD allows such
multiple bestfit restraint.   In fact, that is principally the advantage over
the BESTfit method (only the data structure is changed, the energy subroutines
themselves are the same).  The method can also be used to performed targetted

The form of the new restraint energy is:

                 E   =   Sum_i   KFORCE_i * [RMSD  - OFFSET_i]**2

Where RMSD is the (possibly mass-weighted) root-mean-square-deviation (RMSD)
of the current coordinates with respect to a reference structure, KFORCE_i
is a force constant, and OFFSET_i is a constant value setting a relative
distance with respect to the RMSD of the structure. The restraint energy is
equivalent the normal BestFit energy.  The forces have been checked with the
TEST FIRST command.

All the data structure is stored dynamically on the HEAP and thus, no extra
permanent (static) storage is introduced.  The size of initial HEAP storage
is set by the MAXRmsd integer the first time that the command is issued
(by default this is set to the number of atoms if nothing is specified).

By specifying the RELATIVE keyword, it is possible to impose a 1-D
constraint to simultaneously constrain a given structure with respect to
two end-point structures. This is achieved by constraining the difference
(RMSD2 - RMSD1) instead of just the RMSD1 or RMSD2 values individually, where
RMSD1 and RMSD2 are the RMSD values of given structure from endpoint structure
1 and 2 respectively. This allow full freedom of movements orthogonal to the
relative axis.  For the relative RMSD constraint, the form of the
new restraint energy is:

                 E   =   Sum_i   KFORCE_i * [(RMSD2-RMSD1) - OFFSET_i]**2

The syntax is very similar to all current restraints in CHARMM:

CONS RMSD {RELAtive} {MAXRmsd integer} orient-specs  force-const-spec -
                                       coordinate-spec atom-selection

force-const-spec ::=  { FORCE real } [MASS]  {OFFSet real}

orient-specs     ::=  [ NOROtation ] [ NOTRanslation ]

coordinate-spec::= { [MAIN] } 
                   {  COMP  } 


removes all multiple RMSD restraints


prints all current RMSD restraints with all parameters.

File: Cons ]-[ Node: Rg/RMSD restraint
Up: Top -=- Previous: CONS RMSD -=- Next: Distant Matrix restraint\n

[SYNTAX Rg/RMSD resteaint]

RGYRation { FORCe <real> } REFErence <real> [RMSD] [COMParison] [ORIEnt] 
                         OUTPut_unit <integer> NSAVe_output <integer>
                         SELEction <atom selection> END
          { RESEt      }

  This restraint force restraints the central moment of the selected
  atoms about 1) the center of geometry of the selected atoms (Rg
  restraint) or 2) a specified reference structure structure (RMSD).
  The form of the restraint term is:

                                0  2
       E= 1/2 * CONST * (R   - R  )
                          GY    GY
       2                      2
      R  = 1/N SUM ( r  - R  )       (Rg restraint)
       GY       i     i    CG
      R  = 1/N SUM ( r )
       CG       i     i
       2                   ref 2
      R  = 1/N SUM ( r  - R   )       (RMSD restraint)
       GY       i     i    i

  The specific terms are:
  RGYR             -  Invoke the Rg/RMSD restraint term parser in CHARMM
  FORCe <real>     -  Use a restraint force constant (CONST above) of <real>
  REFErence <real> -  Use a target Rg/RMSD value of <real>, in A.
  RMSD             -  Employ RMSD based restraint instead of Rg restraint.
  COMP             -  Use comparison cordinate set as reference 
                      (default uses main).
  ORIEnt           -  Do a coor orie on coordinates before computing RMSD.
  OUTPut <integer> -  During dynamics write Rg/RMSD history to unit <integer>.
  NSAVe  <integer> -  Save Rg/RMSD history every <integer> steps.
  RESEt            -  CLear Rg/RMSD restraint flags, release memory.
  The following examples illustrate the use of this restraint term.
  1) Add an Rg restraint potential to dynamics run, save trajectory
  information to unit.  Base Rg calculation on Ca positions only.
  open unit 12 write form name traj.rgd
  RGYRestraint Force 50 Reference 12.9 -
      output 12 nsave 50 select type ca end
  2) Add an RMSD based restraint to target a conformational change in a loop.
  open unit 1 read form name open.crd
  read coor card unit 1
  close unit 1
  open unit 1 read form name closed.crd
  read coor card compare unit 1
  close unit 1
  coor orie rms select type ca end
  coor rms select ( ires 12:24 .and. .not. hydrogen ) end
  Calc drms = ?rms / 6
  Calc rmsd = ?rms - @drms
  set count 1
  label domini
      rgyr force 50 reference @rmsd rmsd comp -
           select ( ires 12:24 .and. .not. hydrogen ) end
      mini conj nstep 400 tolenr 0.00001 cutnb 12 ctofnb 10 ctonnb 8 -
           inbfrq -1 atom cdie vatom vswitch fshift bycb
      rgyr reset
      coor rms select ( ires 12:24 .and. .not. hydrogen ) end
      open unit 1 write form name o2c_@count.pdb
      write coor pdb unit 1
  *  Coordinates from frame @count of open to closed path.
  *  Current loop rmsd (residues 12-24 ca only) is ?rms A.
      incr count by 1
      Calc rmsd = @rmsd - @drms
  if rmsd gt 0.1  goto domini
  3) Use RMSD restraint to unfold a helical peptide with RMSD computed
  based minimum RMSD at each step (Oriented).
  set rmsd 2
  label unfold
      rgyr force 100 reference @rmsd rmsd comp orient select type ca end
      mini conj nstep 400 tolenr 0.00001
      rgyr reset
      coor orie rms select type ca end
      open unit 1 write form name frame@rmsd.pdb
      write coor pdb unit 1
  *  Coordiantes from frame with reference rmsd = @rmsd, current rmsd= ?rms
     incr rmsd by 2
  if rmsd le 10 goto unfold

File: Cons ]-[ Node: Distance Matrix restraint
Up: Top -=- Previous: Rg/RMSD restraint -=- Next: top\n

[SYNTAX Distance Matrix restraint]

DMCOnstrain FORCe real REFErence real OUTPut_unit integer 
            NSAVe_output integer CUTOff real NCONtact integer
            {SELE {atom selection} WEIGt real}(ncontact times) 

         THIS COMMAND ADDS A Quadratic POTENTIAL to restrain the
         reaction coordinate. The reaction coordinate is defined
         as a weighted sum of contacts.
         E= 1/2 * CONST * (RHO - DMC0)
     RHO  = SUM (Weigt * (1-STATE );
             i        i          i
     STATE = ----------------------------
          i   1 + EXP(20*(DIST - (CUTOff+0.25)))
     DIST - distance between centers of geometry of residues 
            forming contact i 

  The following examples illustrate the use of this restraint term.
  1) Add an distance matrix restraint potential to dynamics run, save 
  trajectory information to unit.  This example applies a distance 
  restraint of the form given above between pairs of side chain 
  centers-of-mass for a set of 54 contacts and restrains the system
  to a fractional value of the overall reaction coordinate of 0.625.
  Each restraint term is given a weight based on the amount of time
  the given contact was formed in the native state simulation.  For
  details of the method used here, the reader is referred to:
  [Sheinerman and Brooks, JMB, 278, 439 (1998).]

  open unit 25 write form name "dmc/GB1H.rho"       

  define bb -
  sele segid agb1 .and. -
  (type ca .or. type c .or. type n .or. type o ) end
  define sd -
  sele segid agb1 .and. .not. (bb .or. hydrogen) end
  set dmforce 2000.
  set dmref  0.625
  set dmsave 100
  DMCO FORCe @dmforce     REFE @dmref OUTPut 25 NSAVe @dmsave -
  CUTOff 6.5 NCONtact  54 
  sele sd .and. (resi   2  .or. resi   4 ) end WEIGht  0.8158139980007818  
  sele sd .and. (resi   4  .or. resi  21 ) end WEIGht  0.9665094391591674 
  sele sd .and. (resi   4  .or. resi  24 ) end WEIGht  0.9879192119842544 
  sele sd .and. (resi   4  .or. resi  27 ) end WEIGht  0.9895819946415852 
  sele sd .and. (resi   4  .or. resi  51 ) end WEIGht  0.9408263780596178 
  sele sd .and. (resi   5  .or. resi  18 ) end WEIGht  0.9415959785662404 
  sele sd .and. (resi   6  .or. resi   8 ) end WEIGht  0.9999909781703169  
  sele sd .and. (resi   6  .or. resi  17 ) end WEIGht  0.9995401666453647 
  sele sd .and. (resi   6  .or. resi  31 ) end WEIGht  0.9999999999999796
  sele sd .and. (resi   7  .or. resi  16 ) end WEIGht  0.9796513289404690 
  sele sd .and. (resi   7  .or. resi  52 ) end WEIGht  0.7780795300785477 
  sele sd .and. (resi   7  .or. resi  54 ) end WEIGht  0.6905994027412474 
  sele sd .and. (resi   8  .or. resi  17 ) end WEIGht  0.8323463051907218 
  sele sd .and. (resi   8  .or. resi  35 ) end WEIGht  0.9731361885884368 
  sele sd .and. (resi   8  .or. resi  38 ) end WEIGht  0.9174238684004473 
  sele sd .and. (resi   8  .or. resi  40 ) end WEIGht  0.9754236727233426 
  sele sd .and. (resi   8  .or. resi  55 ) end WEIGht  0.9151151195555089 
  sele sd .and. (resi   9  .or. resi  14 ) end WEIGht  0.6645390074837504 
  sele sd .and. (resi   9  .or. resi  56 ) end WEIGht  0.7285333256999120 
  sele sd .and. (resi  17  .or. resi  19 ) end WEIGht  0.8912444191627017 
  sele sd .and. (resi  17  .or. resi  34 ) end WEIGht  0.8920438869599703 
  sele sd .and. (resi  19  .or. resi  21 ) end WEIGht  0.7259406753340431 
  sele sd .and. (resi  19  .or. resi  30 ) end WEIGht  0.6711322572172727 
  sele sd .and. (resi  19  .or. resi  34 ) end WEIGht  0.7844587325555364 
  sele sd .and. (resi  21  .or. resi  27 ) end WEIGht  0.9999999531742481 
  sele sd .and. (resi  23  .or. resi  25 ) end WEIGht  0.9999999994693732 
  sele sd .and. (resi  23  .or. resi  26 ) end WEIGht  0.9999999999991671 
  sele sd .and. (resi  24  .or. resi  27 ) end WEIGht  0.9979691873576392 
  sele sd .and. (resi  24  .or. resi  51 ) end WEIGht  0.6964187937040895 
  sele sd .and. (resi  25  .or. resi  28 ) end WEIGht  0.8927269972864454 
  sele sd .and. (resi  25  .or. resi  29 ) end WEIGht  0.8210307146883757 
  sele sd .and. (resi  26  .or. resi  29 ) end WEIGht  0.8905543107011445 
  sele sd .and. (resi  27  .or. resi  30 ) end WEIGht  0.9612292542758489 
  sele sd .and. (resi  27  .or. resi  31 ) end WEIGht  0.9999999983165117 
  sele sd .and. (resi  28  .or. resi  53 ) end WEIGht  0.9290808245024628 
  sele sd .and. (resi  30  .or. resi  33 ) end WEIGht  0.9273633852785603 
  sele sd .and. (resi  30  .or. resi  34 ) end WEIGht  0.9852961901326586 
  sele sd .and. (resi  31  .or. resi  53 ) end WEIGht  0.9960510555412353 
  sele sd .and. (resi  32  .or. resi  44 ) end WEIGht  0.9811553525890055 
  sele sd .and. (resi  35  .or. resi  40 ) end WEIGht  0.9996382371796980 
  sele sd .and. (resi  35  .or. resi  44 ) end WEIGht  0.9529057898756480 
  sele sd .and. (resi  35  .or. resi  55 ) end WEIGht  0.9921250365257495 
  sele sd .and. (resi  38  .or. resi  40 ) end WEIGht  0.9469272633832431 
  sele sd .and. (resi  40  .or. resi  55 ) end WEIGht  0.9855702226874301 
  sele sd .and. (resi  40  .or. resi  57 ) end WEIGht  0.7830957174831968 
  sele sd .and. (resi  44  .or. resi  55 ) end WEIGht  0.9999950887381377 
  sele sd .and. (resi  45  .or. resi  54 ) end WEIGht  0.9994979068986675 
  sele sd .and. (resi  47  .or. resi  49 ) end WEIGht  0.9878049213251836 
  sele sd .and. (resi  47  .or. resi  50 ) end WEIGht  0.9999916748639699 
  sele sd .and. (resi  47  .or. resi  52 ) end WEIGht  0.9881469299868514 
  sele sd .and. (resi  50  .or. resi  52 ) end WEIGht  0.9973918765679025 
  sele sd .and. (resi  51  .or. resi  53 ) end WEIGht  0.8168211946297691 
  sele sd .and. (resi  52  .or. resi  54 ) end WEIGht  0.8871232193909784 
  sele sd .and. (resi  54  .or. resi  56 ) end WEIGht  0.8700639041799975 

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