Documentation for TRAVEL ver. 2.x in CHARMM. By Stefan Fischer. Feb.12-1993.
********************************************************* * TRAVel (Trajectory Refinement Algorithms) Command * ********************************************************* This module offers access to the CPR (Conjugate Peak Refinement) algorithm for finding reaction-coordinates (described in S.Fischer and M.Karplus, Chem. Phys. Letters vol.194, p.252, 1992), as well as to several other tools for refining and analyzing a reaction-coordinate. * Menu: * Syntax:: Syntax of the TRAVel command * MainCmd:: TRAVel main command and miscellaneous subcommands * TrajManip:: TRAVel TRAJectory subcommand (input/output/analysis) * CPRcmd:: TRAVel CPR subcommand description * SDPcmd:: TRAVel CROS and SDP subcommand description * SCMcmd:: TRAVel SCM subcommand description * Usage:: TRAVEL Usage Note
************************************* * Syntax for the TRAVEL Command * ************************************* Keywords in [...] are optional. Default values are given in (...). Choose one from list : {...|...|...} Main command (entering and leaving the TRAVEL module) : ------------------------------------------------------- TRAVel [MAXPoints int (100) ] [ { XAXI | YAXI | ZAXI } [ ROTAtion ] ] { END | QUIT } Subcommands (within TRAVEL module) : ------------------------------------ [VERBose int (2) ] [DISPlay int (80)] [CHROno { RESEt | PRINt } ] TRAJectory READ [UNIT int (40)] [ ORIEnt | NOORient ] file1.crd file2.crd . . fileN.crd DONE TRAJectory READ NAME file.trj [UNIT int (40)] [ori] [BEGIn int (1)] [SKIP int (1)] [STOP int (0)] TRAJectory WRITe NAME file.trj [UNIT int (40)] TRAJectory ANALyze [SCAN [STEP real]] TRAJectory { INCRease | DECRease } [ NGRIdpoints int | STEPsize real ] [ori] SADDle int CPR [NCYCle int (1)] [HIGHsad] [SADDle] [NGRIdpoint int (5) | STEPsiz real] [saddledef] [relaxspeed] [oscillation] [miscrefine] [linextrem] [ori] CROS [MINDist real (e-4)] [MINStep int (50)] [FIRStep real] [ANGL real (20)] SDP [SAVDistance real (.05)] [MINGrad real (e-3)] [NREActant int (maxpoints/2)] [NPROduct int (maxpoints/2)] [MODE int (4)] [ANGLe real (30.)] [MINCycle int (10)] [linextrem] SCM [NCYCle int (1)] [PROJtol real (.01)] [ANGLe real (90.)] [MINUpdate int (2)] [MAXUpdate int (20)] [linextrem] [ori] COPY [COMP] { SADDle | ORDEr int | INDEx int } saddledef :== saddle-point definition keywords : [SADGrad real (.05)] [SADCycle int (1)] relaxspeed :== relaxation-speed keywords : [TOL1proj real (1.)] [TOL2proj real (3.)] [LINMin int (3N-1)] oscillation :== oscillation detection and prevention keywords : [LOOPreduc int (4)] [FRAMe int (10)] [TOLOscill real (.15)] [PROJincr real (2.)] miscrefine :== miscellaneous refinement keywords : [REMOvemod int (0)] [NTANgent int (6)] [DELTa real (2e-9)] linextrem :== one-dimensional line-extremization keywords : [BRAKetstep real (.1)] [BRKScale real (2.)] [LXEVal int (50)] [BRKMagnif real (5.)] [FIRStep real (.05)] [EXITmode int (3)] [TOLMax real (1.0e-4)] [TOLGrad real (.05)] [TOLStep real (1.0e-14)] [TOLEne real (1.0e-14)] ori :== coordinate orientation keywords : [ ORIEnt | NOORient ] CHARMM command variables set by the CPR command : ?SADI :== index of the highest fully refined saddle-point ?SADO :== order along the path of that saddle-point ?SADE :== energy of that saddle-point
**************************** * TRAVel Main Commands * **************************** TRAVEL [MAXPoints int] ---------------------- When entering the TRAVEL module from CHARMM, it is possible to specify the maximum number of path-points that will make up a reaction-path during a given CPR or SDP refinement. If this value is reached during refinement and even more path-points are needed, then the refinement will stop gracefully, allowing to save the path with a TRAJECTORY WRITE command and later restart TRAVEL with a larger value for MAXPOINTS. Using TRAVEL with IMAGES ------------------------ TRAVEL version 2.x supports IMAGES. Before starting TRAVEL, - the image centering MUST be turned off : IMAGE FIXED SELE all END - the image updating should be set to automatic : NBOND IMGFRQ -1 [ { XAXI | YAXI | ZAXI } [ ROTAtion ] ] keywords : When the only involved symmetry operation is a translation along a single axis and/or a rotation along a single axis, then this axis and/or that rotation must be declared explicitly to TRAVEL. Currently, only the three main axis are supported. If there are translations in more than one dimension, then these keywords are not necessary and should NOT be used. Examples. 1) A periodic DNA helix along the y-axis requires : TRAVEL YAXI ROT 2) An dimeric protein with C2 symmetry : TRAVEL ZAXI ROT 3) A periodic chain along the x-axis, without rotational symmetry : TRAVEL XAXI 4) A lipid-bilayer with periodic boundaries : TRAVEL 5) A crystal : TRAVEL Miscellaneous ------------- Inside the TRAVEL module, it is possible to use all the CHARMM commands handled by the SUBROUTINE MISCOM(), such as STREAM, GOTO, LABEL, SET, INCR, DECR, etc. See the documentation on "Miscellaneous Commands", for more details. VERBOSE int ----------- Sets the amount of details being printed out. Values 0-2 are useful for production runs, values 3-6 should only be used for debugging. The print-level variable PRNLEV is set to 3 when entering the TRAVEL module, but this can be overridden by issuing the PRNLevel command after the TRAVEL command. END --- This exits the TRAVEL module, back into the CHARMM command processor. But all data structures set up while last in TRAVEL are maintained in memory, allowing to later re-enter the TRAVEL module to continue refinement. Note that when re-entering TRAVEL, the value of MAXPOINTS can not be changed. QUIT ---- Exits the TRAVEL module, erasing and freeing all its data structures. It allows to re-enter TRAVEL as if it had never been called before.
************************************************ * TRAVEL I/O and trajectory manipulations. * ************************************************ The TRAVEL module supports its own I/O commands. TRAJECTORY READ --------------- There are two ways to read coordinate-sets into TRAVEL. 1) From a series of formated CHARMM coordinate files. 2) From an unformated CHARMM dynamics trajectory file. The TRAJECTORY READ command can be issued several times (in either of its two forms). The successive path-points will be appended to constitute the input-path for CPR, CROS, SDP or SCM. When reading from an unformated dynamics trajectory file (for example the output of an earlier TRAVEL refinement), it is possible to restrict the input to subsection of the trajectory file with the BEGIN, SKIP, STOP keywords. See the documentation on the ACCUMULATION command for a detailed description of BEGIN, SKIP and STOP. By default and if no fixed atoms have been set-up, the points that are read-in are aligned onto the first point along the path, so as to minimize the RMS-difference in coordinates between a given point and the first point (as with "COOR ORIEnt RMS"). The first point itself is centered and oriented according to its principal axis (as with "COOR ORIEnt"). This can be disabled by issuing the "NOORient" keyword on the "TRAJECTORY READ" command line. If some atoms have been fixed, then all coordinate re-orienting is automatically disabled. Once the path refinement has started (by issuing the CPR command), no more points can be read-in and the TRAJECTORY READ command is disabled. SADDLE ------ After the input-path has been read-in, it is necessary to flag any points along the path that have already been refined to saddle-points. This is to avoid having CPR re-refine those points again. For example if the saddle-point to be flagged is point number "n" along the read-in path, then the command to issue is "SADDLE n", which can be repeated on successive lines if there are more than one such saddle-point to be flagged. TRAJECTORY WRITE NAME --------------------- When the desired number of refinement cycles have been completed, the resulting reaction-path is written out to a CHARMM unformated dynamics trajectory file. This file can then be read in later with the TRAJECTORY READ NAME command, in order to continue the refinement. TRAJECTORY ANALYZE ------------------ After a reaction-path has been read-in and/or refined, this command will print for every path-point the following values : The total path-length up to that point (in angstroem), the RMS-distance to the reactant (in angstroem), the energy, the RMS-gradient, the number of successive line-minimizations which resulted in a RMS-gradient smaller than SADGRA, the path curvature (degrees) and the projection of the gradient at that point onto the tangent to the path. With the SCAN switch, the analysis will be extended to points interpolating between the path-points. The density of the interpolations is determined by default from the value of STEPsize and can be specified with the STEP keyword as an option (but must be specified if a trajectory analysis is done with the SCAN switch and no refinement was yet performed). TRAJECTORY INCREASE or TRAJECTORY DECREASE -------------------------------------------- This command allows to increase or to decrease the density of points along the path. If increasing, new points will be inserted by interpolation between the existing points, so that the distance between adjacent points corresponds to the specified STEPsize. Conversely, if decreasing, points will be removed if they are closer than STEPsize to the point preceding them along the path. COPY [COMP] ----------- Once a reaction-path has been read-in and/or refined, one point can be selected and copied to the CHARMM main or comparison coordinate sets, where it will be available for further manipulation after exiting the TRAVEL module. The point can be selected according to its order number n along the path (ORDER n) or according to its current index i (INDEX i). It is also possible to select the currently highest saddle-point along the path (SADDLE), if such a point is already refined and is the global maximum along the path.
********************************* * Conjugate Peak Refinement * ********************************* CPR --- This command is issued in order to refine a path as a whole, using the Conjugate Peak Refinement method. If the command is issued without the "SADDLE" argument, then only a moderate amount of time is spend on refining individual saddle-points, which will be improved until they reach a given RMS-gradient (specified by "SADGRAD") for the first time (SADCYC=1). The number of overall refinement cycles is specified with NCYCLE. If the switch HIGHSAD is used, then the refinement will stop as soon as the highest saddle-point has been refined. This allows to get a quick estimate of the barrier height, before proceeding with a full refinement of the saddle-point. A critical value to be set when starting refinement is the step-size along the path (STEPSIZE, in Angstroem). It should be as large as possible, while still small enough to adequately probe for the details of the underlying energy surface. It is set by default to the RMS-distance between reactant and product, divided by NGRIDPOINTS+1 . It is better to err on the side of too large a step-size, because CPR will reduce it during refinement progress, if needed. Starting with too small a step-size will result in a path following very closely the bottom of the energy valley, but at the cost of a significant amount of CPU time that would be better spend on refining the saddle regions of the path. The SDP or the SCM method are better suited to further improve the path in the bottom of the valley. IMPORTANT : When refinement is completed, the current (reduced) value of STEPSIZE is printed. It should be noted and used as input value for later continuing refinements with CPR. Other refinement parameters are : SADGRAD = Desired RMS-Gradient at a saddle. SADCYC = Number of cycles SADGRA must be satisfied at a saddle. TOL1PROJ = Gradient projection tolerance when refining a path-point. TOL2PROJ = Gradient projection tolerance when adding a path-point. LINMIN = Maximum number of line minimizations per cycle. LOOPRED = Reduce STEPSIZE by sqr(2) when MOD(numb. of looping, LOOPRED) = 0 FRAME = Frame-length (in number of cycles) for oscillation-detection (maximum allowed value is 20). TOLOSCIL = Energy and Gradient oscillation tolerance ratio. If 0.0 , then oscillation will not be detected. PROJINCR = Factor of increase in TOLPROJ when oscillation is detected. REMOVEMOD= Remove mode. If > 0 , do n line-minimizations for adjacent minimum upon point removal. If <= 0 , add minimum only to avoid looping. NTANGENT = Number of energy probes on path tangent in search for local max. DELTA = Finite difference step, in Angstroem, for one-dim. 2nd derivative. The line-extremization parameters are (only used for fine-tuning) : BRAKETST = Maximal 1-dim. bracketing step, in Angstroem. FIRSTEP = First braketing step. BRKSCALE = Dynamic bracketing Scaling factor. BRKMAGNI = Bracket magnification-limit factor. TOLMAX = Tolerated gradient, 1-dim. maximizations. TOLGRAD = Tolerated gradient, 1-dim. minimizations. TOLSTEP = Smallest 1-dim. extremization step, in Angstroem. TOLENE = Smallest fractional energy change. XITMOD = Exit-Mode with respect to TOLGRA. LXEVAL = Max. number of energy-evaluations during line-extremizations. Note that all values set with the CPR command are remembered when temporarily exiting TRAVEL with the END command. CPR SADDLE ---------- Once a reaction-path has been fully refined with the default settings for SADGRAD and SADCYC, it is possible to proceed with the complete convergence of the highest points along the path towards the exact saddle-points. This is done by adding the SADDLE switch to the CPR command. The only effect of this switch is to modify the default values for the following keywords : [SADGrad real (e-3)] [SADCycle int (3N-1)] [TOL1proj real (.5)] [TOL2proj real (.5)] [PROJincr real (1.)] [LOOPreduc int (0)] [NTANgent int (10)] With these default settings, the time spend on refining individual saddle-points will increase significantly, since in the order of 3N minimizations will be performed every cycle.
***************************** * Steepest Descent Path * ***************************** SDP --- The SDP command provides a carefully controlled descent along the adiabatic valley, down from a saddle-point that has been fully refined with CPR. Four modes of descent are available : Mode 1 : The size of the step down the gradient is reduced until the angle between the path and the gradient is less than the value specified by ANGLE. This mode is the truest to the definition of the adiabatic path, but it also is very slow. It is recommended only for very small molecules. Mode 2 : The step is taken along the gradient until the new gradient is orthogonal to it (= strict steepest descent). Mode 3 : The step along the gradient is taken as large as possible, as long as the energy keeps decreasing (= loose steepest descent). Mode 4 : The steps are those of a strict conjugate-gradient descent. While not following the exact bottom of the adiabatic valley from step to step, the average path obtained by saving path-points after several steps is not significantly worse than the path obtained by saving after several steps in Mode 1. This mode is the fastest and the one recommended for large molecules. If a very accurate adiabatic path is demanded, the resulting path can be further improved with the SCM method. This mode is the default. The input to SDP is a chain of points which straddles the saddle. Up to NREACTANT and NPRODUCT points are then added to the chain, on the reactant and product side respectively. The total number of points can not exceed MAXP, though, so when using SDP, it is advisable to start TRAVEL with a large enough value for MAXP. A new point is added to the path every time the sum of the steps taken since the last addition reaches SAVDISTANCE. When the path enters a region where the energy-gradient is less than specified with MINGRAD, SDP stops extending the path on that side of the saddle-point, provided that it already did at least MINCYCLE steps (to allow moving away from the saddle region, where the gradient is also vanishingly small). CROSsing -------- This command is useful in preparing a minimal chain for input to the SDP from a path refined with CPR. A fully refined saddle-point and two points lying on each side of it are the required input to CROS. The two surrounding points serve as initial guess for the saddle-point crossing direction. The output is also a chain with three points, the second of which is the unmodified saddle-point, while the two surrounding points are now located closer to the saddle and to the bottom of the adiabatic valley. By reiterating the CROS command, these points are gradually improved and yield the reactive mode at the top of the barrier. CROS is NOT designed to refine a saddle point, which instead must be provided to it. If the path provided to CROS has more than three points, it is assumed that this path is the result of a CPR refinement during the same TRAVEL session and all points are deleted, except for the highest saddle-point and its two direct neighbors (provided that CPR had already fully refined the saddle-point, otherwise the CROS command is ignored).
************************************** * Synchronous Chain Minimization * ************************************** SCM --- Before issuing this command, a whole path partially refined with CPR or SDP must be read-in. Unlike CPR and SDP, SCM does not change the number of points along the path. It smoothes an existing path and brings it closer to the bottom of the adiabatic valley by synchronously minimizing all its points, under the constraint that the points move within hyper-planes orthogonal to the path. These planes are updated no more often than every MINUPDATE cycles of conjugate minimization (per point), but no later than after MAXUPDATE of such cycles. How many conjugate line-minimizations are done between plane-updates at a given time on a given point is controlled by the value of the angle between the two path segments joining at that point. This angle can vary from 0.0 degrees (when the path is linear) to 180.0 degrees (when the path is reversing its direction). As long as this angle is decreasing, successive conjugate line-minimizations will be continued (up to MAXUPDATE). On the other hand, if this angle is increasing, minimization of that point is stopped when the angle exceeds ANGLE (provided at least MINUPDATE conjugate line-minimizations have already been performed). This global minimization stops when the projection of the gradient onto the path is less than PROJTOL at every point. This can be quite time consuming, so it is recommended to force an exit from SCM by specifying the maximum number of global cycles NCYCLE. For every cycle of NCYCLE, there will be between MINUPDATE and MAXUPDATE conjugate minimizations done at each point, so that NCYCLE should be kept small to allow for periodic "TRAJECTORY WRITEs". Points that are declared as saddle-points (see "SADDLE n" above) are kept unchanged during the SCM refinement.
******************** * TRAVEL Usage * ******************** Before invoking TRAVEL for the first time, the following must have been done : - A valid topology and parameter file was read into CHARMM. - The molecular system (can have multiple segments) was GENErated. - One set of coordinates were read into the main CHARMM coordinate-set. - The desired Images were set-up (optional), with Image Centering turned off. - All desired atoms were fixed (optional). - The desired non-bond settings were set, if different from defaults. - The FASTER mode was set to a value compatible with the non-bond settings. - At least one successful ENERgy call was made with all these settings. Before invoking the CPR command for the first time, the following must have been done within the TRAVEL module : - Setting the desired amount of information printed out during refinement. This is done with the VERBOSE command (optional). - Reading in the initial coordinate sets (at least two : the reactant and the product conformations, preferably well minimized). Either from a series of formated CHARMM coordinate files or from an unformated CHARMM dynamics trajectory file. Note that if some atoms have been fixed, then the fixed atoms of the reactant, of the product and of all intermediates must have exactly the same Cartesian coordinates. If not, TRAVEL will set the fixed atom coordinates of all structures to those of the reactant and give a warning ! - Flag those path-points that are already known saddle-points with the SADDLE command, so that they do not get refined again. After the reaction path has undergone a number of refinement cycles, save the whole trajectory to an unformated CHARMM dynamics file. This is done with the TRAJECTORY WRITE command. This saving should be done frequently (every few cycles), since a given CPR refinement cycle can be unpredictably long and result in exceeded CPU time. If the CPR command is re-issued after the path was already fully refined during the same TRAVEL session, it will be ignored. When exiting a series of refinement cycles, CPR will print the last value of STEPsiz (only if that value underwent automatic reduction during refinement). That value of STEPsize is important and should be used, if refinement of the trajectory that was written out is later continued. This is done by giving the STEP keyword in conjunction with the first (and only the first, since the step-size will be dynamically reduced as needed) CPR command. Common problems : ----------------- Make sure that the internal coordinates of atoms that are not significantly involved in the reaction are roughly the same for the reactant and the product. For example, a phenyl side-chain has two rotationnaly equivalent orientations, in which the delta1 and epsilon1 atoms are exchanged with the delta2 and epsilon2 atoms. When the reactant and the product coordinate sets are from different origins (for example from different X-ray structures) or the result from different calculations (for example separated by a long dynamics run), then for every group or side-chain that has nearly equivalent orientations (Val, Leu, Phe, Tyr, Arg, Asp, Glu, -CH3 and -CH2- come to mind) it is worth checking that the corresponding atoms are named consistently in the reactant and product coordinates. Otherwise, CPR will proceed with refining all the transition paths associated with these changes in equivalent positions (rotating the Phe ring in the example above), which is very CPU time consuming and will make it difficult to interpret the actual pathway one is interested in. * A TRAVEL input-file, beginning a CPR refinement : * ================================================= * ! Generate the system : STRE gene.str ! Fix atoms, if desired CONS FIX SELECT nomov END ! Read-in any coordinate set and call the energy : OPEN UNIT 13 READ FORM NAME system.crd READ COOR CARDS UNIT 13 CLOSE UNIT 13 ENER IHBFR 0 ! Ready to start TRAVEL : TRAVEL ! Reading-in the initial path-points : TRAJ READ coord_file1.crd coord_file2.crd coord_file3.crd DONE ! Starting refinement, letting CPR choose a STEPSize : CPR NCYCLE 50 ! Writing out the path to a trajectory file : TRAJ WRITE NAME path.trj ! Continue refining, while frequently saving the path : CPR NCYCLE 50 TRAJ WRITE NAME path.trj . . . CPR NCYCLE 50 TRAJ WRITE NAME path.trj ! Analyzing the path : TRAJ ANAL ! Copying the saddle (if already found) to the main-coordinate set : COPY SADDLE ! Returning to CHARMM : QUIT RETURN * A TRAVEL input-file, continuing a CPR refinement : * ================================================== * ! Generate the system : STRE gene.str ! Fix atoms, if desired CONS FIX SELECT nomov END ! Read-in any coordinate set and call the energy : OPEN UNIT 13 READ FORM NAME system.crd READ COOR CARDS UNIT 13 CLOSE UNIT 13 ENER IHBFR 0 ! Ready to start TRAVEL : TRAVEL ! Reading-in the 10 first points of a previously saved trajectory : TRAJ READ NAME path.trj BEGIN 1 STOP 10 ! Merging some points from another trajectory : TRAJ READ NAME path2.trj BEGIN 20 STOP 40 ! Adding a point from a coordinate file : TRAJ READ coord_file.crd DONE ! Declaring the 10th and the 23rd points as already refined saddle-points : SADDLE 10 SADDLE 23 ! Starting refinement, setting STEPSize to the value found at the end of ! the previous refinement : CPR NCYCLE 50 STEP 0.02 ! Writing out the path to a trajectory file : TRAJ WRITE NAME path.trj ! Continue refining, while frequently saving the path : CPR NCYCLE 50 TRAJ WRITE NAME path.trj . . . CPR NCYCLE 50 TRAJ WRITE NAME path.trj ! Setting more stringent requirements for saddle-points : CPR NCYCLE 50 SADDLE ! Continuing the stringent refinement, while frequently saving the path : CPR NCYCLE 50 TRAJ WRITE NAME path.trj . . . CPR NCYCLE 50 TRAJ WRITE NAME path.trj ! Analyzing the path : TRAJ ANAL ! Copying the saddle (if already found) to the main-coordinate set : COPY SADDLE ! Returning to CHARMM : QUIT RETURN * A TRAVEL input-file, getting a smooth steepest-descent path : * ============================================================= * ! Generate the system : STRE gene.str ! Fix atoms, if desired CONS FIX SELECT nomov END ! Read-in any coordinate set and call the energy : OPEN UNIT 13 READ FORM NAME system.crd READ COOR CARDS UNIT 13 CLOSE UNIT 13 ENER IHBFR 0 ! Start TRAVEL with more space to put new path-points : TRAVEL MAXPoints 200 ! Reading-in the saddle-point and two neighbors from a fully refined CPR path : TRAJ READ NAME cpr_path.trj BEGIN 11 STOP 13 ! Analyzing the path : TRAJ ANAL ! Improving the barrier crossing mode : CROSsmode ! Repeat to improve further : CROSsmode ! Verify that the 2 surrounding points are close to the saddle-point : TRAJ ANAL ! Specify the distance between saved points (in Angstroem) in steepest descent: SDP SAVDist 0.05 ! Writing out the path to a trajectory file : TRAJ WRITE NAME path.trj ! Smooth and further refine the path : SCM NCYCle 5 TRAJ WRITE NAME path.trj ! Analyzing the path : TRAJ ANAL ! Returning to CHARMM : QUIT RETURN
NIH/DCRT/Laboratory for Structural Biology
FDA/CBER/OVRR Biophysics Laboratory