Dynamics: Description and Discussion
There are three separate dynamics integrators available in CHARMM:
Name Keyword Module
Original Verlet ORIG dynamcv.src (default)
Leapfrog Verlet LEAP dynamc.src (default with NIH compile flag)
Velocity Verlet VVER dynamvv.src
4-D L-F Verlet VER4 dynam4.src
All methods are based on the Verlet scheme, and when used without
any special features, provide identical trajectories for short
simulations. All methods allow SHAKE.
The ORIG integrator is a standard 3-step Verlet integrator
with few frills. It allows:
Langevin Dynamics (LANG)
Thermodynamic Simulation Method (TSM)
The LEAP integrator is similar to the ORIG integrator, but does
provide increased accuracy (esp. for single precision version of
CHARMM). It allows:
Langevin dynamics (LANG) (with accurate temperatures printed)
Constant Temperature and Pressure (CPT) (based on Berendsen's method)
Accurate pressures with SHAKE
High frequency correction to the total energy
Parallel code
Free energy equilibration indicator (deltaF*V) (with PERT)
Thermodynamic Simulation Method (TSM)
The VVER integrator also provides increase accuracy. It allows:
Constant Temperature (NOSE) (Nose-Hoover method)
Multiple Time Step (MTS)
The VER4 integrator enables the energy embedding technique that entails
placing a molecule into a higher spatial dimension [Crippen, G. M. &
Havel,T.F. (1990) J.Chem.Inf.Comput.Sci. Vol 30, 222-227].
The possibility of surmounting energy barriers with these added
degrees of freedom may lead to lower energy minima. Here, this is
accomplished by molecular dynamics in four dimensions. Specifically,
another cartesian coordinates was added to the usual X, Y, and Z
coordinates in the LEAPfrog VERLet algorithm.
In order to generate a dynamics trajectory, all requirements
for evaluating the energy must be met.
See *note Energy:(energy.doc)Needs.
* Menu:
* Syntax:: Syntax of the dynamics command
* Description:: Description of the keywords and options
* Recommended:: Recommended input options and values
* Discussion:: Running dynamics
* Output:: Output from a dynamics run
* Trajectory:: Trajectory manipulation and I/O
* Merge:: Merging or breaking up trajectory files into
different size pieces. Resampling at a larger
interval. Least squares fit to a reference.
* Reorient:: Reorienting a coordinate trajectory
* RMSDyn:: Computes the RMS difference between two trajectories
* Format:: formatting and unformatting a dynamics trajectory
* Monitor:(monitor.doc). Monitor dihedral transitions
* CPT dynamics:(pressure.doc). CPT dynamics
* 4-D dynamics:(fourd.doc). 4-D dynamics
* Pressure:(pressure.doc)Pressure. The pressure command
* MTS:(mts.doc). Multiple Time Scales Method
* Nose:(nose.doc). Nose-Hoover Dynamics
Syntax for the Dynamics Command
DYNAmics {[LEAPfrog]} {[VERLet]} {STRT } {[TIMEstp real]} [NSTEp integer] -
{ ORIG } {LANGevin} {STARt } { AKMA real }
{ CPT } {VVER } {RESTart}
{ LEAPfrog } {VER4 } {[FIL4dimension]}
{[SKBOnd]} {[SKANgle]}
{[SKDIhedral]}
{[SKVDerWaals]} {[SKELectrostatics]}
nonbond-spec hbond-spec frequency-spec -
unit-spec temperature-spec options-spec -
cpt-spec four-dimension-spec
hbond-spec::= see *note Hbonds:(hbonds.doc).
nonbond-spec::= See *note Nbonds:(nbonds.doc).
frequency-spec::= [INBFrq integer] [IEQFrq integer] [IHBFrq integer]
[IHTFrq integer] [IPRFrq integer] [NPRInt integer]
[NSAVC integer] [NSAVV integer] [NTRFrq integer]
[ILBFrq integer] [ISVFRQ integer]
unit-spec::= [IUNCrd integer] [IUNRea integer] [IUNVel integer]
[IUNWri integer] [KUNIt integer] [CRAShu integer]
[BACKup integer]
temperature-spec::= [FINAlt real] [FIRStt real] [TEMInc real]
[TSTRuc real] [TWINDH real] [TWINDL real]
[TBATh real]
options-spec::= [IASOrs integer] [IASVel integer] [ICHEcw integer]
[ISCAle integer] [ISCVel integer] [ISEEd integer]
[SCALe real] [NDEGg integer] [RBUFfer real]
[AVERage] [ECHEck real] [TOL real]
cpt-spec::= [ PCONst {[PINTernal]} [COMPressibility real] ]
{ PEXTernal } [PCOUpling real]
[PREFerence real]
[VOLUme real]
[ TCONst [TCOUpling real] [TREFerence real] ]
four-dimension-spec::= [K4DInitial real] [INC4Dforce integer]
[DEC4Dforce integer] [MULTK4di real]
[E4FILLcoordinates real]
Options common to minimization and dynamics
The following table describes the keywords which apply to both
minimization and dynamics.
Keyword Default Purpose
NSTEP 100 The number of steps to be taken. This is the number of
dynamics steps which is also equal to the number of
energy evaluations.
INBFRQ 50 The frequency of regenerating the non-bonded list.
The list is regenerated if the current step number
modulo INBFRQ is zero and if INBFRQ is non-zero.
Specifying zero prevents the non-bonded list from being
regenerated at all.
INBFRQ = -1 --> all lists are updated when necessary (heuristic test).
IHBFRQ 50 The frequency of regenerating the hydrogen bond list.
Analogous to INBFRQ
ILBFRQ 50 The frequency for checking whether an atom is in the
Langevin region, defined by RBUF, or not.
non-bond- The specifications for generating the non-bonded list.
-spec See *note Nbonds:(nbonds.doc).
hbond- The specifications for generating the hydrogen bond list.
-spec See *note Hbonds:(hbonds.doc).
[ STRT ] STRT The dynamics is assumed to start from the input
[ ] coordinates using an assignment of velocities given by
[ ] IASVEL. No restart file is read.
[ REST ] The dynamics is restarted by reading the restart file
from unit IUNREA.
TIMESTP 0.001 Time step for dynamics in picoseconds. The default value
is 0.001 picoseconds.
IUNREA -1 Fortran unit from which the dynamics restart file should
be read. A value of -1 means don't read any file
IUNWRI -1 Fortran unit on which the dynamics restart file for
the present run is to be written. A value of -1 means
don't read any file. Formatted output.
IUNCRD -1 Fortran unit on which the coordinates of the dynamics run
are to be saved. A value of -1 means no coordinates should
be written. Unformatted output.
IUNVEL -1 Fortran unit on which the velocities of the dynamics run
are to be saved. -1 means don't write. Unformatted output.
KUNIT -1 Fortran unit on which the total energy and some of its
components along with the temperature during the run are
written using formatted output.
CRASHU -1 Fortran unit where a single DCL command file will be
written. If the machine crashes before a restart file
is written, this file won't be touched. If the crash
occurs after a restart is written but before the run
completes, this file will contain the line, "$
@CRASH". If the run completes, the file will contain
the line, "$ @COMPLET". This allows for an automatic
recovery system after crashes.
NSAVC 10 The step frequency for writing coordinates.
NSAVV 10 The step frequency for writing velocities.
NPRINT 10 The step frequency for storing on KUNIT as well as printing
on unit 6, the energy data of the dynamics run.
IPRFRQ 100 The step frequency for calculating averages and rms
fluctuations of the major energy values. If this
number is less than NTRFRQ and NTRFRQ is not equal to
0, square root of negative number errors will occur.
ISVFRQ NSTEP The step frequency for writing a restart file.
IHTFRQ 0 The step frequency for heating the molecule in increments
of TEMINC degrees in the heating portion of a dynamics
run. Zero means do no heating.
IEQFRQ 0 The step frequency for assigning or scaling velocities to
FINALT temperature during the equilibration stage of the
dynamics run.
NTRFRQ 0 The step frequency for stopping the rotation and translation
of the molecule during dynamics. This operation is done
automatically after any heating.
FIRSTT 0.0 The initial temperature at which the velocities have to be
assigned at to begin the dynamics run. Important only
for the initial stage of a dynamics run.
FINALT 300.0 The desired final (equilibrium) temperature
for the system. Important for all stages except initiation.
TEMINC 5.0 The temperature increment to be given to the system every
IHTFRQ steps. Important in the heating stage.
TSTRUC -999. The temperature at which the starting structure has been
equilibrated. Used to assign velocities so that equal
partition of energy will yield the correct equilibrated
temperature. -999. is a default which causes the
program to assign velocities at T=1.25*FIRSTT.
TWINDH 10.0 The temperature deviation from FINALT to be allowed on the
high temperature side.(+ve). i.e. high side of the
temperature window. Useful during equilibration.
TWINDL -10.0 The temperature deviation from FINALT to be allowed on the
low temperature side.(-ve). i.e. low side of the
temperature window. Useful during equilibration.
TBATH FINALT The temperature of the heatbath in Langevin dynamics. When
set to zero it allows one to do purely dissipative
(quenched) dynamics.
RBUF 0.0 Inner radius of the buffer, or Langevin, region sphere. All
atoms with radial positions greater than RBUF angstroms are
propagated by Langevin dynamics, if the dynamics keyword
LANGevin has been specified.
IASORS 0 The option for scaling or assigning of velocities during
heating (every IHTFRQ steps) or equilibration
(every IEQFRQ steps). This keyword does not control
the initial assignment of velocities.
.eq. 0 - scale velocities. (use ISCVEL option)
.ne. 0 - assign velocities. (use IASVEL option)
IASVEL 1 The option for the choice of method for the assignment of
velocities during heating and equilibration when IASORS
is nonzero. This option also controls the initial
assignment of velocities (when not RESTart)
regardless of the IASORS value.
.eq. 0 - Use the comparison coordinate values
in AKMA units (sorry) with the STRT option.
If NTRFRQ is positive,
then net trans/rot will be removed first.
This option supresses other assignments of velocity.
.gt. 0 - gaussian distribution of velocity. (+ve)
.lt. 0 - uniform distribution of velocity. (-ve)
kinetic energy of 3N velocity components are same.
ISEED 314159 The seed for the random number generator used for
assigning velocities.
ISCVEL 0 The option for two ways of scaling velocities.
.eq. 0 - single scale factor for all atoms
.ne. 0 - a scale factor for each atom proportional to the
kinetic energy average ratio between the system
and along every degree of freedom for that atom.
ICHECW 1 The option for checking to see if the average temperature
of the system lies within the allotted temperature window
(between FINALT+TWINDH and FINALT+TWINDL) every
IEQFRQ steps.
.eq. 0 - do not check
i.e. assign or scale velocities.
.ne. 0 - check window
i.e. assign or scale velocities only if average
temperature lies outside the window.
ISCALE 0 This option is to allow the user to scale the velocities
by a factor SCALE at the beginning of a restart run.
This may be useful in changing the desired temperature.
.eq. 0 no scaling done (usual input value)
.ne. 0 scale velocities by SCALE.
WARNING:
Please use this option only when you are changing the
temperature of the run.
SCALE 1. Scale factor for the previous option.
NDEGF computed Number of degrees of freedom to use in computing the
temperature. If not specified on any call, the value is
computed. This specification is not remembered between
successive calls to dynamics.
AVERAGE no When saving coordinates every NSAVC steps, this option will
cause the average structure of the last NSAVC dynamics steps
to be written instead of the final snapshot coordinate set.
This option is primarily used for making smooth movies.
ECHECK 20.0 The maximum amount the total energy may change on any step.
TOL 1.0E-10 The shake tolerance (if SHAKE is in use).
PCONst false Flag to indicate that constant pressure code will be used.
PINTernal true Flag to indicate that the internal pressure will be coupled
the reference pressure.
PEXTernal false Flag to indicate that the external pressure will be coupled
to the reference pressure.
PCOUpling 0.0 The coupling decay time in picoseconds for the pressure.
A good value for this is 5 ps.
COMPress 0.0 The compressibility in atm**-1. A good value for proteins
is 4.63e-5
PREFerence 1.0 The reference pressure in atmospheres.
VOLUme computed The volume in Angstroms**3 to use for the pressure
calculation denominator. This value is calculated if
the CRYStal feature is use.
TCONst false Flag to indicate that constant temperature code will be used.
TCOUpling 0.0 The coupling decay time in picoseconds for the temperature.
A good value for this is 5 ps.
TREFerence FINALT The reference temperature for constant temperature
simulations.
Recommended CHARMM input for dynamics.
This section is intended only as a guide in setting up
a dynamics simulation input file. Changes should be made as necessary
according to personal tastes and project requirements.
1) For heating and early equilibration:
DYNAMICS LEAP VERLET RESTART(*) NSTEP 20000 TIMESTEP 0.001(+) -
IPRFRQ 1000 IHTFRQ 1000 IEQFRQ 5000 NTRFRQ 5000 -
IUNREA 30 IUNWRI 31 IUNCRD 50 IUNVEL -1 KUNIT 70 -
NPRINT 100 NSAVC 100 NSAVV 0 INBFRQ 25 -
hbond-spec nonbond-spec -
FIRSTT 100.0 FINALT 300.0 TEMINC 100.0 -
IASORS 1 IASVEL 1 ISCVEL 0 ICHECW 0 TWINDH 10.0 TWINDL -10.0
(*) Except for first run, the use STRT in place of RESTART
(+) If bonds to hydrogen atoms are SHAKEd
2) For late equilibration and analysis runs:
DYNAMICS LEAP VERLET RESTART NSTEP 20000 TIMESTEP 0.001 -
IPRFRQ 1000 IHTFRQ 2000 IEQFRQ 5000(*) NTRFRQ 5000 -
IUNREA 30 IUNWRI 31 IUNCRD 50 IUNVEL -1 KUNIT 70 -
NPRINT 100 NSAVC 100 NSAVV 0 IHBFRQ 0 INBFRQ 25 -
hbond-spec nonbond-spec -
FIRSTT 100.0 FINALT 300.0 TEMINC 100.0 -
IASORS 0 IASVEL 1 ISCVEL 0 ICHECW 1 TWINDH 10.0 TWINDL -10.0
(*) Window checking should be disabled for the analysis run (i.e. IEQFRQ=0)
if you want a real microcanonical ensemble.
3) For heating, equilibration and analysis runs using Langevin dynamics:(+)
DYNA LEAP LANGEVIN STRT(*) NSTEP 20000 TIMESTEP 0.001 -
IPRFRQ 1000 IHTFRQ 0 IEQFRQ 0 NTRFRQ 0 -
IUNREA 30 IUNWRI 31 IUNCRD 50 IUNVEL -1 KUNIT 70 -
NPRINT 100 NSAVC 100 NSAVV 0 IHBFRQ 0 INBFRQ 25 -
ILBFRQ 1000 RBUFFER 0.0 TBATH 300.0 -
hbond-spec nonbond-spec -
FIRSTT 300.0 FINALT 300.0 -
IASORS 0 IASVEL 1 ISCVEL 0 ICHECW 0 TWINDH 0.0 TWINDL 0.0
(+) Note that the friction coefficients, in units of 1/ps, must first
be initialized by filling the array FBETA with the SCALAR command
SCALAR FBETA SET <real> <optional atom selection>
(*) Except for first run, the use STRT in place of RESTART
4) For quenched molecular dynamics:
For the first run (STRT), read velocities into the comparison
coordinate set, or this should directly follow a former dynamics command.
DYNA VERLET STRT(*) NSTEP 10000 TIMESTEP 0.001 -
IPRFRQ 1000 IHTFRQ 200 IEQFRQ 200 NTRFRQ 400 -
IUNREA 30 IUNWRI 31 IUNCRD 50 IUNVEL -1 KUNIT 70 -
NPRINT 50 NSAVC 50 NSAVV 0 IHBFRQ 0 INBFRQ 25 -
hbond-spec nonbond-spec -
TSTRUC 300.0 FIRSTT 300.0 FINALT 0.0 TEMINC -30.0 -
IASORS 0 IASVEL 0 ISCVEL 0 ICHECW 1 TWINDH 0.0
or equivalently with Langevin (dissipative) dynamics
DYNA LANGEVIN STRT(*) NSTEP 10000 TIMESTEP 0.001 -
IPRFRQ 1000 IHTFRQ 0 IEQFRQ 0 NTRFRQ 4000 -
IUNREA 30 IUNWRI 31 IUNCRD 50 IUNVEL -1 KUNIT 70 -
NPRINT 50 NSAVC 50 NSAVV 0 IHBFRQ 0 INBFRQ 25 -
hbond-spec nonbond-spec -
TSTRUC 300.0 FIRSTT 300.0 FINALT 300.0 -
ILBFRQ 1000 RBUFFER 0.0 TBATH 0.0 -
IASORS 1 IASVEL 1 ISCVEL 0 ICHECW 0 TWINDH 0.0
(*) For first run, use RESTART otherwise
The IASVEL 0 option causes the comparison coordinates to be used
for the initial velocities (AKMA units).
For subsequent runs the options IASORS 1 and IASVEL 1 may be used
if random velocities are to be periodically assigned.
5) For constant temperature and/or pressure dynamics
DYNA LEAP VERLET STRT(*) NSTEP 20000 TIMESTEP 0.001 -
IPRFRQ 1000 IHTFRQ 0 IEQFRQ 0 NTRFRQ 0 -
IUNREA 30 IUNWRI 31 IUNCRD 50 IUNVEL -1 KUNIT 70 -
NPRINT 100 NSAVC 100 NSAVV 0 IHBFRQ 0 INBFRQ 25 -
PCONst PINTernal COMPress 4.63e-5 PCOUpling 5.0 PREFerence 1.0 -
TCONst TCOUpling 5.0 TREFerence 300.0 -
hbond-spec nonbond-spec -
FIRSTT 300.0 FINALT 300.0 -
IASORS 0 IASVEL 1 ISCVEL 0 ICHECW 0 TWINDH 0.0 TWINDL 0.0
Running Molecular Dynamics
The theoretical basis for dynamical simulations is elementary
physics. The force on a particle is equal to the negative gradient of
the potential energy of the particle. CHARMM can solve this equation
numerically for all atoms in the molecule. A simple second order predictor
two step method due to Verlet is used for integration.
The choice of the integration step size is very important.
One must weigh the increased accuracy of using a small step size against
the longer real time that can be simulated with a given amount of
execution time when a larger step size is used. The time step may be
entered in picoseconds (using the TIMESTP keyword).
CHARMM provides information on the accuracy of the numerical
solution. Since the system has no external forces, the total energy
should be conserved. Numerical errors will result in some fluctuations
in the total energy so a good test is to compare the fluctuations in
total energy to the fluctuations in kinetic energy as these fluctuations
are proportional to the heat capacity of the system. See the next node
for a description of dynamics output.
Because the force constants for the bonds and bond angles are
fairly large, it is reasonable under certain circumstances to constrain
their values during dynamics. Such constraints are applicable if the
harmonic motions are weakly coupled to other motions. The advantage of
such constraints is that the step size of the numerical integration may
be increased without sacrificing accuracy as these terms have the
largest gradients in macromolecules simulated at physiological
temperatures. We use the SHAKE algorithm for applying the constraints,
see *note shake:(cons.doc)SHAKE. SHAKE can be applied to just
the bonds involved with hydrogens, all bonds, all bonds and the angles
involving hydrogens, or all bonds and angles.
A dynamics run has basically four parts; initialization,
heating, equilibration, and the simulation itself. Initialization means
providing an initial position and velocity for all the atoms. Heating is
the process of increasing the kinetic energy of the system up to a final
temperature at which the simulation will be conducted. Equilibration is
the process where the kinetic energy and the potential energy of the
system evenly distribute themselves throughout the system. Only when the
average temperature of the system stabilizes can one collect the
trajectory information for analysis.
The initial coordinates of a simulation are obtained after
applying the minimization algorithm to a complete coordinate set. One
cannot start with a system with a large potential energy as it will
quickly heat up to unreasonable temperatures. For initializing the
velocities, the user can specify velocities from the comparison coordinates
(IASVEL 0), a uniform distribution of kinetic energy along each coordinate
with random sign of the motion along each axis (IASVEL -1) or a
Gaussian distribution of velocities (IASVEL 1 the default). The temperature
at which velocities are assigned is determined by FIRSTT and TSTRUC
by the algorithm:
Tassign = 2*(FIRSTT-TSTRUC) + TSTRUC.
For a harmonic system equilibrated to TSTRUC equal partition of the
energy will result in an equilibrated temperature of roughly FIRSTT.
If TSTRUC is not specified 1.25*FIRSTT will be used for assignment.
Velocities may also be passed to dynamics in the comparison
coordinate set (as opposed to assignment). This allows the user
considerable flexibility in setting up the initial conditions.
The heating of system is performed gently by increasing the
kinetic energy by a small amount periodically. The number of integration
steps between heating applications, the final temperature, and the
kinetic energy increment are all user specified. In addition, there is a
choice in the method of increasing the kinetic energy of the system. One
may scale existing velocities or reassign them. The velocities can be
scaled by either one scale factor calculated for the kinetic energy of
the system averaged over many time steps or by scale factors established
for each atom base ed on the ratio of its time averaged kinetic energy
with that of the system. If reassignment is chosen, the velocities can
have either a uniform or Gaussian distribution.
To equilibrate the structure, one can specify a window around
the final temperature where velocity adjustments will be made. The
choice of velocity adjustments is the same as described above for
heating.
For the actual run, CHARMM will output the position and
velocities of all atoms at intervals specified by the user. The
temperature window can be set larger so that any gross conformational
changes which result in a different potential energy will cause the
temperature to be maintained.
At any time energy is added to the system, the angular momentum
of the system will be reduced to zero and translational motion will be
stopped. One can also request that these operations be performed at any
time during the dynamics run.
The use of a restart file is essential for running dynamics.
The restart file contains information about the most recent coordinate
sets necessary for the VERLET algorithm. In addition the values of
the energy accumulators are stored. All other information (such as
SHAKE, fixed atoms, harmonic constraints, friction coefficients) has
to be regenerated before invoking a dynamics restart.
When the run is initiated, a restart file must be written using
the IUNWRI keyword. As the dynamics routine complete NCYCLE, see *note
Output::, steps of dynamics, the Fortran unit specified by IUNWRI will
be rewound and a restart file will be written. In case of crashes, one
has restart files corresponding to various points in the run. The CRASHU
variable may prove valuable. Successive runs of CHARMM to continue the
dynamics run must read the previous restart file using the IUNREA
keyword and write it out for the next part of the run.
Restarts may be done to reset various options, or to break up a long run
into several shorter runs. Restart files will only run with the version
of CHARMM they are created with.
There are many numbers giving the frequency of actions to be
taken during dynamics such as updating the non-bonded list, heating the
molecule etc. Some of these numbers are adjusted along with the number
of steps to run so that numbers all have a common divisor. At the
present time, there are combinations which result in errors. At some
point an attempt may be made to catalog all the actions, and check for
erroneous processing.
If one is interested in simulating the motion of part of the
system with the rest of the system remaining fixed, it is possible to
fix atoms in place, see *note fix:(cons.doc)fixed atom. If this
is done, there are several effect on the dynamics. First, since the
system is now anchored in space, the center of mass motion and total
angular velocity is never stopped. Second, the number of degrees of
freedom used for calculating the temperature is set to the number of
free atoms times 3 minus 6. Third, the coordinate and velocity
trajectory files will contain the position of the fixed atoms only
once, and all other records will hold just the moving atoms. This
saves a great deal of disk space.
Trajectory files can be merged, broken in smaller pieces, and
sampled at different intervals. Likewise, said operations can be
performed on coordinate trajectories while rotating the coordinates to
match a given coordinate set.
When the DYNAmics command exits, the main coordinate set
contains the final coordinate positions from the last energy evaluation and
the comparison coordinates will contain the final velocities In AKMA units.
Finally, a brief discussion of the Langevin dynamics algorithm is
presented. The Langevin dynamics algorithm presently in CHARMM was intented
to be used primarily with the "Stochastic Boundary Molecular Dynamics" method
and consequently has been restricted to an algorithm which is valid only for
the case FBETA*TIMESTEP<<1.0. That is for cases where relatively small
friction coefficients are used. Typically values of FBETA*TIMESTEP up to
about 0.3 still produce a stable dynamics which also satisfy the
fluctuation-dissipation theorem.
The algorithm itself reduces to the Verlet algorithm when FBETA is
zero and consequently may be used to do regular dynamics, actually it is
the same routine which does both dynamics.
In using Langevin dynamics care must be taken to first initialize
the array FBETA by using the scalar commands e.g.,
CHARMM >SCALAR FBETA SET <real> <atom selection>
Failure to do this just means you are doing regular dynamics so no warning is
given by CHARMM.
Contents of a dynamics output
Note: This description of the output of a command is not
normally going to be part of the documentation of commands. The dynamics
output is sufficiently confusing that this description is necessary.
The first part of CHARMM's output after a dynamics command lists
all of the options that apply to that part of the run. Then, any
information about velocity assignments (temperature changes) follows.
Any time the velocities are changed in an anisotropic way, the motion of
and about the center of mass will be stopped. This results in a printout
both before and after this operation of the "DETAILS ABOUT CENTRE OF
MASS". Its position and velocity are output followed by the components
of the angular momentum. The last line gives the translation kinetic
energy of the system, and thus one should expect a drop in the total
energy and temperature of the system afterwards.
Non-bonded interaction and hydrogen bond updates will appear
intermittently and are cleared labeled.
Every NPRINT steps, the total energy and various contributions
will be printed. This output is preceded by a title which gives the
correspondence of numbers to energy names. After IPRFRQ steps will
appear the averages and RMS fluctuations. After the second such printout
of averages and RMS fluctuations, the averages and RMS fluctuations for
the run upto the last turning of the molecule will be given. This gives
you longer range statistics. Such a calculation will not be done if
IPRFRQ equals NTRFRQ. The ratio of total energy to kinetic energy
fluctuations is an excellent measure of the accuracy of the run.
After the averages are printed, a least squares fit of the total
energy against the step number will be made to look for drift in the
energy. Two such values are printed, one for the last IPRFRQ steps, and
one to the previous turn. Next, the initial energy for the statistics,
both short range and long, are printed. Finally, the correlation
coefficient of the energy versus step is given for both ranges. A value
close to zero indicates no systematic drift; a magnitude near 1 means
you have a real problem with the dynamics.
This process of printout continues until the end of the run is
reached. Just before the last energy is printed will appear a message
about the writing of coordinates and velocities to their respective
files.
Reading and writing trajectory frames with direct commands
This facility allows the creation or manipulation of trajectory files
The main uses of this facility are;
1) creating artificial trajectory files from coordinate frames
2) reading an existing trajectory frame by frame for analysis that
requires access to a variety of CHARMM commands
3) modifying an existing trajectory (copy with changes) such as
minimizing each frame or other operations.
[Syntax TRAJectory command]
===================================================================
There are three commands that comprise this facility.
1) Initializing trajectory I/O
TRAJectory {read-spec} {write-spec}
read-spec:== [IREAd unit] [NREAd int] [SKIP int]
[BEGIN INT] [STOP INT]
write-spec:== [IWRIte unit] [NWRIte int] [NFILE int] [EXPAnd]
[NOTHer int] [DELTa real] [SKIP int]
IREAd - first unit to read from (default: do not read)
NREAd - number of units to read from (default:1)
SKIP - skip value for both reading and writing (default:1)
IWRIte - first unit to write to (default: do not write)
NWRIte - number of units to write to (default:1)
NFILe - number of frames on each output file (default: total)
EXPAnd - flag to free fixed atoms in copying (only if reading)
NOTHer - number of frames in previous files (if not reading) (d:0)
DELTa - output delta value (if not reading) (default:0.001)
2) Reading a frame
TRAJectory READ [COMP]
3) Writing a frame
TRAJectory WRITe [COMP]
The reading and writing commands do not have any specifiers
other than whether the comparison or main coordinates will be used.
===================================================================
There are three modes of operation;
1) Create a new trajectory.
The IWRIte and NFILe keywords must be used. The default
values for the others are listed above. If several files
will be made in different CHARMM runs that will be linked together
later, the NOTHer keyword value should be increased by NFILe on
each subsequent run.
EXAMPLE: Create a "movie" trajectory that involves the rotation
of a single sidechain (residue 21).
COOR AXIS SELE ATOM * 21 CA END SELE ATOM * 21 CB
OPEN WRITE UNIT 22 FILE NAME TYR21.ROT
TRAJECTORY IWRITE 22 NWRITE 1 NFILE 360 SKIP 1
* trajectory showing the rotation of sidechain 21
*
SET 1 1
LABEL LOOP
COOR ROTATE AXIS PHI 1.0 SELE ATOM * 21 * .AND. .NOT. ( TYPE C -
.OR. TYPE N .OR. TYPE H ) END
TRAJ WRITE
INCR 1 BY 1
IF 1 LT 360.5 GOTO LOOP
STOP
===================================================================
2) Reading an existing trajectory
The IREAD keyword must be used. The default NFILe value is 1
and the remaining values if not specified will be read from the
trajectory file.
EXAMPLE: find the structure with the lowest energy and save it.
UPDATE ...
OPEN READ UNIT 22 FILE NAME DYN1.TRJ
OPEN READ UNIT 23 FILE NAME DYN2.TRJ
TRAJECTORY IREAD 22 NREAD 2 SKIP 100
SET 1 1
SET 9 9999.0
LABEL LOOP
TRAJ READ
GETE
IF 9 LT ?ENER GOTO NEXT
SET 8 @1
COOR COPY
SET 9 ?ENER
LABEL NEXT
INCR 1 BY 1
IF 1 LT 1000.5 GOTO LOOP
OPEN WRITE CARD UNIT 12 NAME LOWE.CRD
WRITE COOR COMP CARD UNIT 12
* structure with the lowest energy
* frame number @8 with energy @9
*
STOP
===================================================================
3) Copying from one trajectory to another.
The operation of this command works in the same mode as
the MERGE command, except a variety of CHARMM commands can
be executed between reading and writing of frames.
EXAMPLES: Create a new trajectory where every frame is minimized
for 200 steps.
OPEN READ UNIT 22 FILE NAME DYN.TRJ
OPEN WRITE UNIT 32 FILE NAME DYN.MIN
TRAJECTORY IREAD 22 SKIP 100 IWRITE 32
* minimized trajectory
*
SET 1 1
LABEL LOOP
TRAJ READ
MINI ABNR NSTEP 200
TRAJ WRITE
INCR 1 BY 1
IF 1 LT 1000.5 GOTO LOOP
STOP
Merges or breaks up a trajectory into different numbers of files
Frequently, one generates a trajectory into small files to
minimize the CPU time of one job. However, so many files are usually
hard to manage so it is desirable to merge said files into larger units.
This command provides that capacity. In addition, it is possible to
break up the trajectory into smaller pieces and to sample the trajectory
less frequently than originally generated.
Another option is to optionally rotate the structure at each
frame to least squares fix a reference structure.
[Syntax MERGE dynamics trajectories]
MERGE [ COOR ] [FIRSTU unit-number] [NUNIT integer] [SKIP integer]
[ VEL ] [OUTPutu unit-number] [NFILE integer]
[ DRAW ] [BEGIN integer] [STOP integer]
[first-atom-selection]
[ ORIEnt [MASS] [WEIGht] [NOROt] [PRINT] second-atom-selection ]
Keyword table
Option Default Purpose
[COOR] COOR Specification of the type of trajectory file. COOR is
[VEL ] coordinates; VEL is velocities.
[DRAW] Make a CHARMM movie (do not write any files, just display)
FIRSTU 51 The first unit of the trajectory to be read.
NUNIT 1 The number of units to be read starting with FIRSTU
SKIP 1 Only those coordinate whose dynamics step number
modulo SKIP will be reoriented and written out.
OUTPUTU 61 The first unit number of the output trajectory
NFILE The number of coordinate sets written to each output
merged file. If left out, this will be set to the number
of coordinates in the first input file times the number of
input files. WARNING: This default will generate a bad
trajectory file if SKIP is not set to the interval
actually present in the trajectories. Further, if you
set its value to be larger than the number of
coordinates that are actually written in any output
file, you will have problems. The error that is
generated results from the control array in the
beginning specifying that there are more coordinates
than actually exist in the file. EOF errors will result
when the trajectory is read.
BEGIN First step number to start reading from
STOP Last step number to read
first-atom-sel Selection of atoms to include in the output file.
ORIEnt Flag to specify best fit rotation and translations.
MASS Use mass weighting in best fit.
WEIGht Use weighting array for best fit weights.
NOROt Only translate in the best fit.
PRINT Print the details of best fit
second-atom-sel Selection of atoms to use in the best fit.
The title of the output trajectory will be copied from the input
trajectory.
Reorienting a coordinate trajectory
If one is interested in reorienting every set of coordinates
found in a dynamics trajectory with respect to some reference structure,
one can use the ORIEnt option in conjunction with the MERGe command.
The process of reorienting a coordinate trajectory works as
follows: A series of files containing the trajectory are assigned to
successive units prior to a CHARMM run. The coordinates stored therein
are presumed to have been written every NSAVC steps. CHARMM will read
each coordinate, select some periodically, reorient them, and write them
to successive units where each output file will have a user specified
number of coordinates. The following table lists the options involved:
Option Default Purpose
ORIE .false. Specify that a least squares RMS fit will be done.
MASS .false. Use a mass weighting in the fit
WEIGH .false. Use the weighting array (wmain) in the fit
NOROt .false. Just shift the centers to best fit.
PRINt .false. Print what happened to each coordinate set.
atom-selection all Select which atom to use in the fit.
If atoms were fixed during the dynamics, the new trajectory
produced will not have fixed atoms because the rotations applied to
each coordinate set will be different thereby yielding different
coordinates for the fixed atoms. Fixing the coordinates leads to a
large space reductions, so the reorientation process will therefore
result in potentially much larger trajectory files.
See *note fix: (cons.doc)Fixed Atom.
Computes the RMS difference between two trajectory files
and make a matrix of results. Large files should be reduced with the
MERGe command before processing this command.
RMSDynmics [IREAd unit-number] [JREAd unit-number] [IWRIte unit-number]
[BEGIn integer] [STOP integer] [IMAGes]
ORIEnt [MASS] [WEIGht] [NOROt] [RMS] atom-selection
IREAd int - unit number of the first trajectory file.
JREAd int - unit number of the second trajectory file.
IWRIte int - Unit for the output matrix.
BEGIn int - Starting step number (default: first)
STOP int - Ending step number (default: last)
IMAGes - Use image atoms for the analysis
ORIEnt - Do best fit of structures
MASS - Use a mass weighting in best fit.
WEIGht - Use the weighting array in best fit.
NOROt - Best fit without letting the structures rotate.
RMS - Do RMS fit between structures, otherwise,
align structures with the axis.
atom-selection
- Atoms to use in the fitting procedure.
Format or unformat a dynamics trajectory
DYNAmics FORMat FIRStunit <unit> NUNIt <int> BEGIn <int>
SKIP <int> STOP <int> OUTPut <unit>
OFFSet <int> SCALe <int> MODE <FORTRAN-FORMAT>
DYNAmics UNFOrmat INPUt <unit> OUTPut <unit>
These commands allow to convert binary trajectory files into a machine
independent yet compact format and to convert them back into binary files.
The defaults for OFFSet, SCALe and MODE are: OFFSet=600, SCALE=10000, and
MODE=12Z6. The trajectory is converted into positive integers according to
the formula <integer>=INT(<real>+OFFSET)*SCALE). The user has to make
sure that all coordinates of the trajectory are within OFFSET angstroms.
The precision may be increased by choosing a larger SCALE and FORTRAN-format,
e.g. MODE=11Z7 OFFSET=100000. ("Z" is the hexadecimal format and is available
on most machines.)
NIH/DCRT/Laboratory for Structural Biology
FDA/CBER/OVRR Biophysics Laboratory