# CHARMM c28a3 umbrel.doc

```Charmm Element doc/umbrel.doc \$Revision: 1.1 \$

```

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

```

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.

* 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

```
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
plane1
DISTance     point1       point2
line
plane
ANGLe        direction1   direction2   direction3
line1        line2        direction3
line1        plane2       direction3
plane1       plane2       direction3
RATIo        distance1    distance2
COMBination  name1 weight1 name2 weight2 ...

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.

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

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

```
II. Running dynamics under an umbrella potential and collecting

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

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

```
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

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

```
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