Charmm Element doc/umbrel.doc $Revision: 1.1 $
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
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 ...
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
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
below:
...
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.
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
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
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)
NIH/DCRT/Laboratory for Structural Biology
FDA/CBER/OVRR Biophysics Laboratory
Modified, updated and generalized by C.L. Brooks, III
The Scripps Research Institute