Documentation for TRAVEL ver. 3.10 in CHARMM. By Stefan Fischer. June 1-1996.
*********************************************************
* 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-path.
* 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 ] ]
[ SCALe [{ COMP | WMAIn }] ]
{ END | QUIT }
Subcommands (within TRAVEL module) :
------------------------------------
[VERBose int (2) ]
[DISPlay int (80)]
[CHROno { RESEt | PRINt } ]
TRAJectory READ [UNIT int (40)] [REFErence] [ori]
file1.crd
file2.crd
.
.
fileN.crd
DONE
TRAJectory READ NAME file.trj [UNIT int (40)] [REFErence] [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 [ int [ int ...]]
CPR [cpr_exit] [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 }
cpr_exit :== exit-criteria for CPR
[NCYCle int (1)] [NECAlls int (999999)] [HIGHsad]
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 or higher 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
Using TRAVEL with QM/MM (SCALE keyword)
----------------------------------------
TRAVEL version 3.x or higher has been adapted to run efficiently on combined
quantum/empirical energy potentials. To take advantage of this, some default
parameters must be modified with the CPR command (see "Using CPR with
non-analytical energy potentials").
Also, if the reaction involves the mixing of covalent bond breaking/making
(associated with high energy barriers) together with large-scale motions of
MM atoms (associated with low/no barriers), the resulting energy
surface is going to be very anisotropic over the sampled conformational space.
In some cases, it may be necessary to reduce the steepness of the gradients
for the QM atoms by scaling their coordinates relative to the MM atoms
(note that this scaling is available also for MM atoms, if desired).
This is achieved by multiplying each atomic coordinate by a scaling factor
greater or equal to 1 . The coordinates of the desired QM atoms should be
assigned scaling factors in the 2-10 range. Other atoms should be assigned
a scaling factor of 1 . This effectively "stretches" the energy surface
along the dimensions with the large scaling factors, thereby reducing the
anisotropy. Note that SCALing is not required for every kind of QM/MM
reaction and should only be used if necessary (i.e. CPR fails to converge).
Coordinate-scaling is activated by issuing the SCALe keyword when launching
TRAVEL. By default, the scaling-factors are taken for each atom from the
X,Y,Z values stored in the main coordinate-set. Adding the COMP keyword will
result in these values being taken from the comparisson cooridnate set.
Alternatively, if the WMAIn keyword is added, the scaling factors are taken
from the WMAIN vector (in that case, the X,Y and Z dimensions of a given atom
will all three be assigned the same scaling-factor).
Use the SCALAR command (see SCALAR.DOC) to set the appropriate values of the
X,Y,Z coordinate vectors or of the WMAIN vector, before calling TRAVEL.
The scaling-factors can only be set once, upon the first call to TRAVEL.
Upon COPY or TRAJ WRITE, the coordinates of the system are re-scaled back
to normal.
When coordinate-scaling is activated, re-orientation is automatically disabled
during the CPR refinement (as if issuing the NOORient keyword). But it is
still enabled by default during the reading of the initial structures with
"TRAJ READ".
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 file DYNANAL.DOC for a detailed description of BEGIN, SKIP and STOP.
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.
By default and if no fixed atoms have been declared, the points that are
read-in are all aligned onto a reference structure, so as to minimize the
RMS-difference in coordinates between each point and the reference point
(as would "COOR ORIEnt RMS"). This can be disabled by issuing the "NOORient"
keyword on the "TRAJECTORY READ" command line (generally not recommended).
The aligning status (NOORient or ORIEnt) used during READ will stay in effect
during CPR refinement, unless expressly overridden.
The reference structure can be read-in explicitly (recommended), by specifying
the REFERence keyword with the "TRAJ READ" command. In that case, the first
point in the list (or file, if the NAME keyword is used) becomes the
reference point and does NOT belong to the path.
Alternatively, if the REFERence keyword is not given, the first point to
be read-in will automatically be copied and used as the reference point.
It will ALSO serve as the first point of the path (reactant).
The reference point can only be specified once, before any other points
are read-in. The reference point is always centered and oriented according
to its principal axis (as would "COOR ORIEnt").
When saving the trajectory (TRAJ WRITE), the reference point is NOT written
as part of the path and must be respecified if the refinement is to be
continued in a later TRAVEL session.
If some atoms have been fixed, or if there are crystal images, then all
coordinate re-orienting is automatically disabled (the reference point can
still be read-in, but it will be ignored).
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.
Column:
N = The order of the points along the path.
Idx = The identifier number given to each point in the last TRAVEL
calculation.
Length = The total path-length up to that point (in angstroem).
rms(xN-x1) = The RMS-distance to the reactant (in angstroem).
Energy = The energy (in kcal/mol)
rms(Grad) = The RMS-gradient of the energy (in kcal/mol*angstroem).
LINMIN = -100*SADCYC, if a true saddle-point (according to SADGRA & SADCYC
criteria only, of course). To make sure a true
saddle-point is flagged, SADCYC should be 3*N-1
(N = number of moving atoms) and SADGRA smaller than
.05, which is done automatically when using the
SADDLE keyword with CPR.
+100*SADCYC, if the point was declared as a saddle-point with the
SADDLE command and is accepted as such, but another
lower local energy max. is present on the path
segment starting at this point.
+10*SADCYC, if the point was prematurely flagged as a saddle-
point, normally because of numerical convergence
problems (the gradient has reached machine precision,
so that all SADCYC line-minimizations cannot be
performed). Generally a true saddle-point.
0, otherwise.
Angle = The path curvature, as measured by the angle (in degrees) of the
two path segments joining at this point. If zero, the path is
linear; if 180 the path makes a U-turn.
GrdProj = The RMS-gradient, after its projection onto the plane orthogonal
to the path at this point. Should be zero, if the path was
perfectly following the gradient (achievable by running SCM, but
only in reasonable amounts of CPU time for small systems).
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. Use the command "SADDLE", to avoid removing the saddle-points
from a refined path with the TRAJ DECREAse command.
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. It can also
be limited by specifying an approximate number of energy-calls NECAlls, which
should not be exceeded. NECAlls allows better control than NCYCLE to determine
the amount of CPU-time before exit.
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 for a CPR run is the step-size along
the path (STEPsiz, in Angstroem), used to scan the path for local
energy-maxima. It should be as large as possible to avoid unnecessary
energy-calls, 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. If the reactant
and the product are the same structure (for ex. a full group rotation),
then some STEPsiz value MUST be specified on the CPR command-line.
It is better to err on the side of too large a step-size, because CPR will
reduce STEPsiz during the refinement, if needed. This happens every LOOPred
occurences of algorithm-looping (unless LOOPred is set to 0), when the STEPsiz
is divided by SQRT(2). For ex. if LOOP=4, then the reduced STEPsiz will
be 1/2 its original value after 8 occurences of looping.
IMPORTANT :
When refinement is completed, the current value of STEPSIZE is printed
(if it was reduced). It should be used as input value with the STEP keyword
for the next CPR run.
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 SCM method is better suited to further smooth the path at
the bottom of the valley.
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
(must be larger than finite difference step used for 1st deriv).
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 improve the convergence of the highest
path-points towards the accurate and non-degenerate saddle-points.
A number of CPR parameters have slightly different optimal values for
this second refinement phase. For convenience, they are all modified
simultaneously by adding the SADDLE switch to the CPR command. The only
effect of this switch is to modify the default settings 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)]
Specifying any of these keywords on the CPR command-line in addition to
the SADDLE keyword will overide its "SADDLE" default setting.
The time spend on refining individual saddle-points is significant, since
in the order of 3N line-minimizations will be performed for each saddle-point.
See description of SADCYC above.
Using CPR with non-analytical energy potentials
-----------------------------------------------
If the energy and/or its derivative is computed non-analytically
(for. ex. iteratively in the case of quantum potentials or with finite
difference gradient functions), some of CPR's parameters must be
relaxed, to prevent the lack of precision in the energy/gradient to perturb the
convergence of CPR :
DELTA should be increased, so that it is larger than the one-dimensional
finite-difference step used to compute the derivatives.
TOLMax 0.1 TOLStep 1.0e-5 TOLEne 1.0e-5 are recommended values for numerical
evaluations of the energy, such as resulting from SCF in quantum potentials
(in conjunction with DELTA 1.0e-4 ).
Example : CPR DELTA 1.0e-4 TOLMax 0.1 TOLStep 1.0e-5 TOLEne 1.0e-5
If the SADDle keyword is issued with the CPR command (i.e. a saddle-point
is to be refined), relax the saddle-point convergence criterion. A recommended
value is SADGrad 0.01 . This will still result in RMS-gradients at the
saddle-point which are smaller than SADGrad, because of the larger number of
line-minimizations that will be performed when using the SADDle keyword
(see SADCycle).
Example : CPR SADDle SADGrad 0.01
*****************************
* 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.
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!!! AND !!!
- !!! 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
(unless the NECAlls keyword was used).
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.
Definition of the RMS-difference in coordinates
-----------------------------------------------
Normally and in CHARMM, the RMS-difference in the coordinates of two
structures X1 and X2 is defined as follows :
RMS = |X1-X2| / SQRT(N) ,
where |X1-X2| = SQRT( SUM_N(x1-x2) + SUM_N(y1-y2) + SUM_N(z1-z2) )
and N is the number of non-fixed atoms.
In TRAvel, for all input and output involving a coordinate RMS-difference,
the value used is :
RMS'= |X1-X2| / SQRT(3*N) ,
in other words, the normal RMS = SQRT(3) * RMS'
Prematurely reaching MAXPP
--------------------------
On rare occasions (about once every 100000 CPR-cycles), nearly identical points
are added one next to the other between two path-points, so that the maximum
allowed number of path-points MAXPP is quickly reached and the path-refinement
stops. If this happens, simply continue the refinement after having excised
these extraneous points as well as one or a few path-points more on each side
along the path.
* 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
! This energy-call is mandatory !!! :
ENER IHBFR 0
! Ready to start TRAVEL :
TRAVEL
! Read-in the coordinate-set upon which all path-points will be oriented
! (it will not be part of the path) :
TRAJ READ REFErence
coord_file0.crd
DONE
! 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
! Read-in the coordinate-set upon which all path-points will be oriented
! (it will not be part of the path) :
TRAJ READ REFErence
coord_file0.crd
DONE
! 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 the 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
Modified, updated and generalized by C.L. Brooks, III
The Scripps Research Institute