Next: Dynamics Outputs, Previous: Minimization, Up: Energy Manipulations
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. CONGEN can solve this equation numerically for all atoms in the molecule. It has two different methods available, a simple second order predictor two step method due to Verlet and a fifth order predictor - corrector multivalue method described by Gear.
In either algorithm, the choice of 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 either picoseconds (using the TIMESTP keyword) or the internal AKMA units (using the AKMASTP keyword).
CONGEN 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 Shake Command. 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 zero velocity, a uniform distribution of kinetic energy along each coordinate with random sign of the motion along each axis (IASORS 0) or a Gaussian distribution of velocities (IASORS 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.
The heating of a 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 based 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, CONGEN 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.
When dynamics is used for simulated annealing, it is useful to use the TLIMIT option. This option applies a velocity damping to any atom whose temperature exceeds the limit. It is especially useful for solving structures using NMR constraints, see NMR Constraints. Consider the case of two atoms which are supposed to be close in space, but are not at the beginning of the simulation. When they approach, they will both acquire substantial kinetic energy as they travel down the potential well. With the TLIMIT option, their velocity will be reduced so that they will not accelerate to such high speeds that the numerical integration of their motion will fail because the step size is too large. In addition, the velocity damping will help to hold them in position. You can see details of the damping process by turning on the TLIMIT debugging variable, see Debug Command.
The use of a restart file is essential for running dynamics. Since running dynamics requires storing various derivatives of the position with respect to time, this information must be stored for the life of the dynamics run. This capability is provided by the restart file. When the run is initiated, a restart file must be written using the IUNWRI keyword. As the dynamics routine complete NCYCLE 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 CONGEN 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. See Dynamics Outputs, for a description of these variables.
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 Fixed Atoms, for more information. 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.
The trajectory produced by the dynamics procedure can be analyzed in great detail using the analysis facility (see Analysis). In addition, trajectory files can be merged, broken in smaller pieces, and sampled at different intervals. See Manipulating Trajectories, for details. Likewise, said operations can be performed on coordinate trajectories while rotating the coordinates to match a given coordinate set. See Reorienting a Coordinate Trajectory, for details.
In the table of keywords which follow, one can select on keyword from a set of keywords enclosed in square brackets and such keywords take no values after them. All other keywords must be specified with a value (see Energy Manipulation Syntax). See General Energy Operands, for other variables which apply. See AKMA, for a description of the AKMA units system.
Keyword Default Purpose
[ VERL ] VERL Verlet algorithm is used for integration in dynamics.
[ GEAR ] Gear ( 6 vector ) algorithm is used for integration.
If SHAKE is used, GEAR option is overridden and Verlet
algorithm is used.
[ 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.
AKMASTP .02 Time step for Dynamics in AKMA units. The AKMASTP
TIMESTP keyword is used to enter a step size in AKMA units.
TIMESTP is used for picoseconds. The default value is
0.02 AKMA units (0.000977 picoseconds).
TOL 1.0E-6 Shake tolerance, i.e. the maximum relative error allowed
in the constraining of a SHAKEn bond length or bond angle.
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.
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 5 The step frequency for writing coordinates.
NSAVV 5 The step frequency for writing velocities.
NPRINT 1 The step frequency for storing on KUNIT as well as printing
on unit 6, the total energy data of the dynamics run.
IPRFRQ 50 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.
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 298.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 5.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 -5.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.
This number must specified as a negative number to
be meaningful.
TLIMIT 0.0 The temperature limit for atoms. When this option is
positive, CONGEN will compute the velocity that
corresponds to this temperature for all atoms. After
every dynamics time step, the new velocity for each
atom will be checked against this limit. If it is
exceeded, the atom's velocity will be lowered to the
limit, and a new corresponding position will be calculated.
The use of this option can result in a loss of energy
in the system. The limit should be set above the desired
temperature of the system, but CONGEN makes no check to
see that the limit is reasonable. Also, this option
only works with Verlet integration. It will also work
with the SHAKE algorithm.
IASORS 0 The option for scaling or assigning of velocities during
heating ( every IHTFRQ steps) or equilibration
(every IEQFRQ steps).
.eq. 0 - scale velocities.
.ne. 0 - assign velocities.
IASVEL 1 The option for different assignments of velocities.
.eq. 0 - zero velocity assignment
.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.0 Scale factor for the previous option.
ETETEST 20.0 This variable is used for the total energy conservation
test in dynamics. If the total energy varies by more
than this amount and EKETEST multiplied by the kinetic
energy, the run will be terminated. This check is turned
off if TLIMIT is set.
EKETEST 0.1 See ETETEST above.