CHARMM c30b1 umbrel.doc

Charmm Element doc/umbrel.doc $Revision: 1.2 $

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

                Umbrella Sampling along a Reaction Coordinate

                        By J. Kottalam, December 1990

        This module of charmm is used for defining a reaction
coordinate for any molecule based on its structure and impose an
umbrella potential along that reaction coordinate  (i.e., to run
activated dynamics along this coordinate) in order to trace out the
free energy profile during the structural change along the coordinate.

* Menu:

* RXNCOR::       Specifying a reaction coordinate in CHARMM
* Dynamics::     Running dynamics under an umbrella potential
* Statistics::   Extracting results from umbrella dynamics
* Results::      Interpreting and using the results

File: Umbrel ]-[ Node: RXNCOR
Up: Top -=- Next: Dynamics -=- Previous: Top\n

     I. Specifying a reaction coordinate in CHARMM

            A reaction coordinate is either a distance or an angle.
        In some cases it can be a linear combination of distances or
        angles.  For example, the interconversion between the chair 
        and boat forms of cyclohexane can be described in terms of 
        three rotation angles.  These angles rotate carbons 2, 4 and 6
        respectively about the plane of carbons 1, 3 and 5.  In some
        other cases the reaction coordinate can be a ratio of two
        distances.  In proton transfer between atoms A and B the
        ratio of A - proton distance to the A - B distance is a good
        measure of the extent of reaction.

            A distance can be between two points or the perpendicular
        distance of a point from a plane or a line.  An angle can be
        between two (intersecting or non-intersecting) lines and so on.
        Planes and lines are ultimately defined in terms of points
        and points can be defined as the centres of geometry (or mass)
        of a groups of atoms in the macromolecule.

            Thus we have in charmm a general algorithm to define
        any reaction coordinate (of the above form) in terms of the
        cartesian coordinates of arbitrary sets of atoms.

        SYNTAX:  All the commands invoking this module are prefixed by the
        keyword RXNCord.  The geometrical elements such as points, lines and
        planes are referenced by names selected by the user.  The command
        syntax for the definition is as follows:

        (Blanks, except the rightmost, stand for words directly above)

    RXNCord: DEFIne name POINt [MASS] atom_selection
                         DIREction    point1       point2
                                      direction1   direction2
                         LINE         THROugh point1 THROugh point2
                                                     PARAllel_to line1
                                                     PERPend line1 PERP line2
                                                     NORMal_to   plane1
                         PLANe        THROugh point1 THRO point2 THRO point3
                                                     CONTaining line1
                                                     PERPendic  line1
                                                     PARAllel   line1 PARA line2
                         DISTance     point1       point2
                         ANGLe        direction1   direction2   direction3
                                      line1        line2        direction3
                                      line1        plane2       direction3
                                      plane1       plane2       direction3
                         RATIo        distance1    distance2
                         COMBination  name1 weight1 name2 weight2 ...
                         SCOMbination name1 mult1 name2 mult2 ...
    RXNCord: SET name

        The name immediately following the DEFIne keyword is the name of 
        the object being defined.  Names referenced towards the right 
        side are names of objects already defined.
        If the MASS keyword appears in the POINt definition, then the point
        is the centre of mass, otherwise it is the centre of geometry.
        The order of points in DIREction definition is important, unlike
        in DISTance definition.  When two directions are used to define a
        third direction, the third direction is perpendicular to the other
        two.  The ANGLe is the one subtended by the first two arguments
        measured anticlockwise when viewed along direction3. The angle
        defined with a line and plane is the one between the line and the
        plane normal, NOT the plane itself.  RATIo is distance1/distance2.
        COMBination is the weighted mean of all names; SCOMbination is the
        simple sum mult1*name1+mult2*name2...

        The SET command sets a particular geometric element as the reaction
        coordinate.  This element must be a scalar quantity and must be
        already defined. 

        If a name is referenced and not defined, the program will stop
        pointing this out.  If a name is defined but not referenced, no
        message is issued, but this will result in some unnecessary
        computations without affecting the results.

        EXAMPLE:  Let us number the carbon atoms in cyclohexane as C1,...,C6.
        Consider the plane of atoms c1, c3 and c5 and the plane of atoms
        c1, c2 and c3.  Let us call the angle between these two planes as a
        rotation angle.  There are two other rotation angles about the
        c1-c3-c5 central plane.  The mean of these three rotation angles 
        serves as a reaction coordinate for the chair-boat conversion.

        rxncor: define c1 point select atom cycl 1 c1 end
        rxncor: define c2 point select atom cycl 1 c2 end
        rxncor: define c3 point select atom cycl 1 c3 end
        rxncor: define c4 point select atom cycl 1 c4 end
        rxncor: define c5 point select atom cycl 1 c5 end
        rxncor: define c6 point select atom cycl 1 c6 end
        rxncor: define d13 direction c1 c3
        rxncor: define d35 direction c3 c5
        rxncor: define d51 direction c5 c1
        rxncor: define d12 direction c1 c2
        rxncor: define d34 direction c3 c4
        rxncor: define d56 direction c5 c6
        rxncor: define norc direction d35 d13
        rxncor: define nor1 direction d13 d12
        rxncor: define nor2 direction d35 d34
        rxncor: define nor3 direction d51 d56
        rxncor: define alf1 angle norc nor1 d13
        rxncor: define alf2 angle norc nor2 d35
        rxncor: define alf3 angle norc nor3 d51
        rxncor: define mean combi alf1 1.0 alf2 1.0 alf3 1.0
        rxncor: set mean

File: Umbrel ]-[ Node: Dynamics
Up: Top -=- Next: Statistics -=- Previous: RXNCOR\n

    II. Running dynamics under an umbrella potential and collecting
        statistics for free energy

     RXNCord: UMBRella FORM form KUMB ku DEL0 del0

     RXNCord: STATistics LOWDelta lowdel HIDElta hidel DELDelta deldel
                 START start

        The UMBR commad specifies the form and parameters for an umbrella
        potential.  The forms implemented so far:

              form          functional form of potential
              ----          ----------------------------
               1               ku*(delta-del0)**2
               5               ku*(delta-del0)**2 + Ubias(delta)

        FORM 5:

        Ubias is a biasing potential, used in addition to the harmonic
        restraining potential, to truncate high barriers of activated
        processes.  Ideally, Ubias(Rc) would be the negative of the 
        potential of mean force g(Rc), where Rc is the reaction coordinate.
        Ubias is implemented as a cubic spline function based on a
        tabulated data to be read prior to the call of RXNC.  An example
        of the use of FORM 5, and the biasing potential file is shown

     rxncor: define RC  distance PCC PCO
     rxncor: set RC
     open read unit 33 form name bias.pot
     rxncor: bias unit 33
     close unit 33
     rxncor: umbrella kumb  20. del0 1.50 form 5   

        WHERE bias.pot contains:

     * Trial Biasing Potential 
     * example, 1200
         5                   ! number of points (up to 25)
     1.0 -28.0               ! Rc  Ubias(Rc)
     1.3 -18.0
     1.4 -8.0
     1.5  0.0
     1.6  -2.0

        END description of FORM 5                             JG 12/00

        The STAT command specifies when to collect statistics and in what 
        range of the reaction coordinate.  Starting from the 'start'-th step
        of dynamics, the number of occurences of delta (the value of the
        reaction coordinate) will be counted in each interval 'deldel' long.
        This counting will be done in the range 'lowdel' to 'hidel'.

        Dynamics is then run by invoking the DYNAmics command.

File: Umbrel ]-[ Node: Statistics
Up: Top -=- Next: Results -=- Previous: Dynamics\n

   III. Extracting results from umbrella dynamics

              The statistics collected by the RXNCOR STAT command is
        printed out at the end of dynamics by the command

     RXNCord: WRITe [ UNIT unit ]     default is outu

        This will print out a table without any header.  The header
        is left out on purpose to enable this file to be read by any
        plotting program.  The meanings of the numbers are therefore
        explained here.  Recall that the RXNC STAT command essentially
        set up a range of the reaction coordinate (delta) and divided
        this range into small pieces.  For each piece the following 
        information is printed out in one line
           the midpoint of the delta interval
           the free energy at this midpoint (after subtracting the umbrella
           the number of observations for which delta fell in this
                         interval.  (One observation is made at every
                         step of the dynamics)

        Apart from this cumulative statistics, it may be useful to have
        a print out of the delta values versus time.  In fact the trajectory
        of any quantity defined by the RXNC DEFI command can be printed by
        using the RCNC TRACe command.  This command should appear before 
        the dynamics command.

        RXNCord: TRACe name [ UNIT unit ]    (default is outu)

        EXAMPLE: (continuing from the previous example)

        rxncor: trace alf1 unit 22
        rxncor: trace alf2 unit 23
        rxncor: trace alf3 unit 24
        rxncor: trace mean unit 25
        rxncor: umbrella kumb 15.0 form 1 del0 -0.4
        rxncor: statistics lowdelta -0.6 hidelta -0.4 deldel 0.002 -
                     start 5000
        dynamics rest nstep 15000 firstt 300.0 finalt 300.0 ihtfrq 0 -
              teminc 0.0 iunwrite 19 kunit 20 iunread 17
        rxncor: write unit 21
        close unit 21

File: Umbrel ]-[ Node: Results
Up: Top -=- Previous: Statistics -=- Next: Top\n

     IV.  Interpreting and using the results

             The essential results are the first and second columns of
        the output from the RXNC WRIT command.  This essentially is a
        plot of free energy vs. the reaction coordinate in the covered
        range of the reaction coordinate values (delta values).  In order
        to cover the full range of delta the procedure is repeated by
        constraining delta around a particular value by using an umbrella
        potential centered at that value.

             The details are omitted here.  For theory and practice
        please see Kottalam and Case in Journal of American Chemical
        Society, 110, 7690 (1988)

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