Go to the previous, next section.

Minimization and Dynamics

An important function of CONGEN is the evaluation and manipulation of the potential energy of a macromolecular system. Two types of manipulation are supported. First, one can minimize the energy by adjusting the coordinates of all the atoms in order to reduce its value. Several minimization algorithms are provided. Second, one can solve Newton's equations of motion using the potential energy in order to obtain a dynamics trajectory.

The commands which invoke these three types of manipulation are all parsed by the same routine, and any common actions are done therein.

Syntax for Energy Manipulation Commands

{ MINImize } 
{ DYNAmics } repeat(keyword) repeat(keyword value)
{ ENERgy   } 

value ::= {  real   }
          { integer }

Syntactic ordering: All keywords and keyword value pairs can appear in any order.

Table of Keywords with references

Keyword
Reference (where to read about it)

ABNR
See section Minimization Options.
AKMASTP
See section Running Molecular Dynamics.
CONJ
See section Minimization Options.
CRASHU
See section Running Molecular Dynamics.
DYNAmics
See section Running Molecular Dynamics.
EIGRNG
See section Minimization Options.
FINALT
See section Running Molecular Dynamics.
EKETEST
See section Running Molecular Dynamics.
ETETEST
See section Running Molecular Dynamics.
FIRSTT
See section Running Molecular Dynamics.
hbond-spec
See section Generation of Hydrogen Bonds.
IASORS
See section Running Molecular Dynamics.
IASVEL
See section Running Molecular Dynamics.
ICHECW
See section Running Molecular Dynamics.
IEQFRQ
See section Running Molecular Dynamics.
IHBFRQ
See section Options Common to Minimization and Dynamics, and see section Symmetry and Molecular Images
IHTFRQ
See section Running Molecular Dynamics.
IMAXP
See section Options Common to Minimization and Dynamics.
INBFRQ
See section Options Common to Minimization and Dynamics.
IPRFRQ
See section Running Molecular Dynamics.
IPRINT
See section Options Common to Minimization and Dynamics.
ISCALE
See section Running Molecular Dynamics.
ISCVEL
See section Running Molecular Dynamics.
ISEED
See section Running Molecular Dynamics.
IUNCRD
See section Running Molecular Dynamics.
IUNREA
See section Running Molecular Dynamics.
IUNVEL
See section Running Molecular Dynamics.
IUNWRI
See section Running Molecular Dynamics.
KUNIT
See section Running Molecular Dynamics.
MASS
See section Minimization Options.
MINDIM
See section Minimization Options.
MINImize
See section Minimization Options.
NCGCYC
See section Minimization Options.
non-bond-spec
See section Generation of Non-bonded Interactions.
NPRINT
See section Running Molecular Dynamics.
NRAP
See section Minimization Options.
NSAVC
See section Running Molecular Dynamics.
NSAVV
See section Running Molecular Dynamics.
NSTEP
See section Options Common to Minimization and Dynamics.
NTRFRQ
See section Running Molecular Dynamics.
PCUT
See section Minimization Options.
PRTMIN
See section Minimization Options, and section Vibrational Analysis.
SCALE
See section Running Molecular Dynamics.
SD
See section Minimization Options.
STEP
See section Minimization Options.
STPLIM
See section Minimization Options.
STRICT
See section Minimization Options.
TEMINC
See section Running Molecular Dynamics.
TFREQ
See section Vibrational Analysis.
TIMESTP
See section Running Molecular Dynamics.
TLIMIT
See section Running Molecular Dynamics.
TOL
See section Running Molecular Dynamics.
TOLENR
See section Minimization Options.
TOLGRD
See section Minimization Options.
TOLITR
See section Minimization Options.
TOLSTP
See section Minimization Options.
TSTRUC
See section Running Molecular Dynamics.
TWINDH
See section Running Molecular Dynamics.
TWINDL
See section Running Molecular Dynamics.
VIBRation
See section Vibrational Analysis.

Requirements for Energy Manipulations

Before the energy of a system can be evaluated and manipulated, a number of data structures must be present.

First, a PSF must be present.

Second, a parameter set must be present. It must contain all parameters which are required by the PSF being used.

Third, coordinates must be defined for every atom in the system. An undefined coordinate has a particular value, and if two coordinates have the same value, division by zero will occur in the evaluation of the energy. If the positions of hydrogens are required, the hydrogen bond generation routine, see section HBUILD Command), must be called before the energy is evaluated. In the case of energy evaluations, CONGEN will ignore interactions involving undefined atom positions. However, minimization and dynamics will fail if any positions are unknown.

Fourth, provisions must be made for having a hydrogen bond list and a non-bonded interaction list. Having non-zero frequencies for updating this lists is one way, one can also read these lists in, see section READ -- Reads Data from External Sources, or generate them with separate commands, see section Generation of Hydrogen Bonds, or section Generation of Non-bonded Interactions.

All of these requirements can be met through the use of a restart file, see section Running Molecular Dynamics.

Energy Manipulation Options

There exist several commands which can modify the way the potential energy is calculated or can affect the way energy manipulations are performed.

The CONSTRAINT command, see section Constraints, can be used to set constraints of various kinds. First, it can be used to set flags for particular atoms which will prevent them from being moved during minimization or dynamics. Second, it can be used to add positional constraint term to the potential energy. This term will be harmonic about some reference position. The user is free to set the force constant. Third, the user can place a harmonic constraint on the value of particular torsion angles in an attempt to force the geometry of a molecule. Fourth, atoms can be fixed in place. When atoms are fixed, they are not allowed to move at all. Any energy terms involving nothing but fixed atom are ignored, because they will have no effect. In addition the trajectories written by the dynamics command only have one copy of the fixed atom coordinates written. Thus, such trajectories take less disk space.

The SHAKE command, see section SHAKE -- Fixing Bond Lengths Or Angles in Dynamics, is used to set constraints on bond lengths and also bond angles during dynamics. It is very valuable in that it permits a larger step size to be used during dynamics. This is vital for dynamics where hydrogens are explicitly represented as the low mass and high force constant of bonds involving hydrogen require a ridiculously small step size.

The WEIGHT command, see section Set Weight Command Variables -- WEIGHT, is used to change the weights on each of the terms in the potential energy function. This is useful during dynamics with NMR constraints when one is interested in getting to the final solution as quick as possible.

The user interface commands, see section Interfacing to CONGEN, can be used to modify the calculation of the potential and to add another term to the potential energy.

Options Common to Minimization and Dynamics

The following table describes the keywords which apply to both minimization and dynamics. See section Generation of Non-bonded Interactions, for a more detailed description of the non-bonded options. See section Generation of Hydrogen Bonds, for more information on hydrogen bonds.

Keyword  Default  Purpose


NSTEP    100      The number of steps to be taken. For minimization, this
                  is the number of cycles of minimization, not the
                  number of energy evaluations. For dynamics, 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.


IHBFRQ   50       The frequency of regenerating the hydrogen bond list.
                  Analogous to INBFRQ.

non-bond-         The specifications for generating the non-bonded list.
-spec             

hbond-            The specifications for generating the hydrogen bond list.
-spec             


IPRINT   0        If given a non-zero value, causes the printing of
                  all contributions to the energy whose atoms are less
                  than IMAXP, more or less, every IPRINT cycles. Useful
                  in debugging although somewhat obsolete since the
                  analysis facility can provide the same information in
                  a more readable form.  There may be bugs in the
                  frequency of when the printing occurs.


IMAXP    0        In conjunction with IPRINT, causes the print out of 
                  contributions to the energy. This keyword determines
                  the largest atom number which has its contributions
                  printed. In some cases the check is made on all atoms
                  in the interaction; in others, the check is made on
                  just the first.


NPRINT   1        The frequency with which energies are printed during
                  course of dynamics or ABNR minimization.

Minimization Options

There are several different algorithms for minimizing the energy of the system. They all involve calculating the derivative of the potential energy, and possibly the second derivative, and using that information to adjust the coordinates in order to find a lower energy. In the descriptions of the algorithms below, a capitalized keyword at the beginning of each paragraph is the key word used to invoke that method. After the descriptions is a listing of all keywords associated with minimization.

The simplest minimization algorithm is steepest descents (SD). In each step of this iterative procedure, we adjust the coordinates in the negative direction of the gradient. It has one adjustable parameter, the step size, which determines how far to shift the coordinates at each step. The step size is adjusted depending on whether a step results in a lower energy. I.e., if the energy drops, we increase the step size by 20% to accelerate the convergence. If the energy rises, we overshot a minimum, so the step size is halved. Steepest descents does not converge in general, but it will rapidly improve a very poor conformation.

A second method is the conjugate gradient technique (CONJ) which has better convergence characteristics (R. Flectcher & C.M. Reeves, The Computer Journal 7:149 (1964)). The method is iterative and makes use of the previous history of minimization steps as well as the current gradient to determine the next step. It can be shown that the method converges to the minimum energy in N steps for a quadratic energy surface where N is the number of degrees of freedom in the energy. Since several terms in the potential are quadratic, it requires less evaluations of the energy and gradient to achieve the same reduction in energy in comparison to steepest descents. Its main drawback is that with very poor conformations, it is more likely to generate numerical overflows than steepest descents. The algorithm used in CONGEN has a slightly better interpolation scheme and automatic step size selection.

A third method involves solving the Newton-Raphson minimization equations iteratively (NRAP). This procedure requires the computation of the derivative of the gradient which is a matrix of size O(n**2). The procedure here is to find a point where the gradient will be zero (hopefully a minimum in energy) assuming that the potential is quadratic. The Newton-Raphson equations can be solved by a number of means, but the method adopted for CONGEN involves diagonalizing the second derivative matrix and then finding the optimum step size along each eigenvector. When there are one or more negative eigenvalues, a blind application of the equations will find a saddle point in the potential. To overcome this problem, a single additional energy and gradient determination is performed along the eigenvector displacement for each small or negative eigenvalue. From this additional data, the energy function is approximated by a cubic potential and the step size that minimizes this function is adopted. The advantages of this method are that the geometry cannot remain at a saddle point, as sometimes occurs with the previous procedures, and that the procedure converges rapidly when the potential is nearly quadratic (or cubic). The major disadvantage is that this procedure can require excessive storage requirements, o(n**2), and computation time, o(n**3), for large molecules. Thus we are currently restricted to systems with about 200 atoms or less for this method. This method may not be used when Images are used (see section Symmetry and Molecular Images).

The fourth method available is an adopted basis Newton-Raphson method (ABNR) (D. J. States). This routine performs energy minimization using a Newton-Raphson algorithm applied to a subspace of the coordinate vector spanned by the displacement coordinates of the last (MINDIM) positions. The second derivative matrix is constructed numerically from the change in the gradient vectors, and is inverted by an eigen vector analysis allowing the routine to recognize and avoid saddle points in the energy surface. At each step the residual gradient vector is calculated and used to add a steepest descent step onto the Newton-Raphson step, incorporating the new direction into the basis set.

A fifth method involving an amalgam of conjugate gradient and steepest descents (GCSD) is currently unavailable due to a bug in the command parser.

In the table which follows, keywords enclosed in square brackets means that one can choose one in the set. Such enclosed keywords do not expect a value after them. All other keywords are used for specifying values, See section Syntax for Energy Manipulation Commands. The method column shows which method the keyword affects. See section Options Common to Minimization and Dynamics, for common variables.

Keyword  Default  Method  Purpose






[ CONJ ] CONJ             Do conjugate gradient minimization.
[ SD   ]                  Do steepest descent minimization.
[ NRAP ]                  Do Newton-Raphson minimization.
[ ABNR ] [MASS]           Do Adopted Basis Newton-Raphson minimization,
                             with mass weighted forces if specified.
[ CGSD ]                  This is to combine CONJ and SD (Not available)


STEP     .02      ALL     Initial step size for the minimization algorithms.
                          Reasonable values for the various methods are best
                          determined by trial and error.

                  
PRTMIN   1        ALL     A flag indicating how much to print during
                          minimization.
                  SD      No effect
                  CONJ    If less than 2, the energy is printed only once
                          each cycle. A setting of 2 shows the energy for
                          each evaluation plus variables used in the method.
                  NRAP    if greater than 1, a brief overview of
                          the least squares cubic fitting procedure
                          given for eigenvalues less than TFREQ.
                  ABNR    If less than 2, the energy is printed out only for
                          successful steps (improvements in total
                          energy). Some description of how the step was
                          chosen is printed if it is set to 2, and a
                          very verbose description is given for values
                          of 3 or higher.


NCGCYC     100    CONJ    The number of conjugate gradient cycles executed 
                          before the algorithm restarts. This number
                          will be automatically lowered to the shortest
                          hydrogen bond or non-bonded list update
                          frequency.  The algorithm will fail if the
                          potential is changed will it is running.


PCUT     .9999    CONJ    If the cosine of the angle between the old and new 
                          P vector is greater than PCUT, the algorithm will be
                          restarted. This prevents the algorithm from plodding
                          down the same path repeatedly. If PRTMIN is less
                          than 2, one effect of the restart is that the step
                          size will go its initial value. If this happens many
                          times, you've converged.

      
EIGRNG   .0005    ABNR    The smallest eigenvalue (relative to the largest)
                          that will be considered nonsingular.

                 
MINDIM   5        ABNR    The dimension of the basis set stored.

                
STPLIM   1.0      ABNR    The maximum Newton Raphson step that will
                          be allowed.


STRICT   0.1      ABNR    The strictness of descent.  The energy of a new step
                          must be less than the previous best energy + STRICT
                          for the new step to be accepted.

                          
MASS    --       ABNR    Use unweighted forces by default or mass weighted
                          if specified.  Mass weights converge more slowly but
                          allow association with normal mode frequencies.


TFREQ    1.0      NRAP    The smallest eigenvalue that is considered to be
                          non-negative (i.e. do cubic fitting on all
                          eigenvalues smaller than this).


NFREQ   NATOM*3   NRAP    The number of degrees of freedom to be
                          optimized (the number of lowest eigenvalues).
                          Use the default whenever practical.


TOLENR  0.0       ABNR    A tolerance applied to the change in total energy
                          change during a cycle of minimization (NCYCLE steps)
                          If the energy change is less than or equal to
                          TOLENR, the minimization routine will exit.


TOLGRD  0.0       ABNR    A tolerance applied to the average gradient during
                          a cycle of minimization.  If the average gradient
                          is less than or equal to TOLGRD, the routine
                          will exit.


TOLITR  100       ABNR    The maximum number of energy evaluations allowed
                  CONJ    for a single step of minimization.


TOLSTP  0.0       ABNR    A tolerance applied to the average step size during
                          a cycle of minimization.  If the average step size
                          is less than or equal to TOLSTP, the routine
                          will exit.

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. 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 section SHAKE -- Fixing Bond Lengths Or Angles in Dynamics. 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 section 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 section Set Debugging Variables -- DEBUG.

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 section Dynamics Output, 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 section Fixing Atoms in Place, 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 section The Analysis Facility of CONGEN). In addition, trajectory files can be merged, broken in smaller pieces, and sampled at different intervals. See section Manipulating Trajectories, for details. Likewise, said operations can be performed on coordinate trajectories while rotating the coordinates to match a given coordinate set. See section 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 section Syntax for Energy Manipulation Commands). See section Options Common to Minimization and Dynamics, for other variables which apply. See section The CONGEN System of Units: 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.

Dynamics Output

The first part of CONGEN'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 up to 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.

Manipulating Trajectories

The trajectories produced from a dynamics calculation generally require some sort of additional processing to make them easier to analyze. Currently, two such commands are provided, MERGE and ROTATE. MERGE is described below, and ROTATE is described under the COMPARE command, see section Reorienting a Coordinate Trajectory.

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.

Merge Trajectory Syntax

     MERGE DYN [ COOR ] [FIRSTU unit-number] [NUNIT integer] [SKIP integer]
               [ VEL  ]       [OUTPUTU unit-number] [NFILE integer]

Merge Trajectory Options

Option  Default  Purpose



[COOR]    COOR   Specification of the type of trajectory file. COOR is 
[VEL ]           coordinates; VEL is velocities.






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 coordinates written to each output 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.
The title of the output trajectory will be copied from the input trajectory.k

Go to the previous, next section.