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