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

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)

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