$MD group                       (relevant if RUNTYP=MD)
 
This group controls the molecular dynamics trajectory for a
collection of quantum mechanical atoms and/or Effective
Fragment Potential particles.
 
A typical MD simulation starts with an equilibration phase,
running long enough to produce a randomized structure and
velocity distribution.  Typically equilibration is done
with an NVT ensemble, allowing the system to equilibrate to
a desired temperature.  A production run restarts with the
positions and the velocity and quaternion data from the
equilibration run, might use either a NVE or NVT ensemble,
and collects radial distribution functions and other
properties.
 
Only a few properties are computed from the MD trajectory,
apart from correct radial distribution functions.  In
particular, the pressures, diffusion constants, and heats
of vaporization that appear on the printout (presently only
for pure EFP runs) are from a preliminary code, which has
not yet been verified.
 
If the system contains only EFP particles, it may be placed
in a periodic box, according to the minimum image
convention.  The optional periodic boundary conditions,
along with cut-offs, are given in the $EFRAG input.  See
also the $EWALD input group for long-range electrostatic
treatment if PBC is used.
 
            The first keywords relate to the steps:
 
MDINT  = MD integrator selection.
       = FROG (leapfrog).  This is less accurate, and lacks
         the special ensemble stepper option NVTNH.
       = VVERLET (velocity Verlet) - default.
 
DT     = MD time step size, in seconds, default=1.0d-15,
         which is a femtosecond.
 
NVTNH    selects a integrator step appropriate to the
         desired ensemble.  This is only implemented for
         velocity Verlet.
       = 0 means use NVE Verlet stepping
       = 1 means use NVT Verlet stepping
       = 2 means use Nose/Hoover chain NVT Verlet stepping
         The default is 2 if either NVT option RSTEMP or
         RSRAND is chosen, but is 0 otherwise.
 
NSTEPS = number of MD time steps to be found in this run,
         default=10000.
 
TTOTAL = total time elapsed in the previous part of a MD
         trajectory which is being restarted (READ=.TRUE.).
         The default means this trajectory is a new one, or
         perhaps the start of a production phase of the MD.
         (default=0.0 seconds)
 
                      * * *
 
BATHT  = bath temperature, in Kelvin (default=300.0)
         This value is used during NVT runs, or if the
         MD is initialized to a Maxwell-Boltzmann velocity
         distribution. In REMD, BATHT is an array with the
         size given by the number of replicas.
 
                      * * *
 
        Two options exist to create NVT runs, to bring
        the system to a desired bath temperature.
        If neither is chosen, the ensemble is NVE:
 
RSTEMP = flag to rescale the temperature.  default=.FALSE.
 
DTEMP  = temperature range for the RSTEMP option.  The
         velocities are rescaled to the bath temperature
         if T < (BATHT-DTEMP) or T > (BATHT+DTEMP).
         The default is DTEMP=100.0 degrees.
 
RSRAND = flag to reset to Maxwell-Boltzmann distribution,
         using random numbers (same algorithm as MBT and
         MBR) to choose individual velocity magnitudes and
         directions.  default=.FALSE.
 
NRAND  = number of steps for the RSRAND option.  Reassign
         velocities (translational and rotational) every
         NRAND time steps.  Default=1000.
 
NVTOFF = step number at which to turn off either NVT
         thermostat, and switch to NVE.  At this point, the
         NVTNH parameter will be reset to 0, and the PROD
         flag will be turned on, so that the production
         run will start (gathering and printing the RDF
         information to .log file).  This keyword is also
         useful in NVE runs to postpone the accumulation of
         production information.  The default means no
         switch to NVE (default=0).
 
JEVERY = report simulation quantities (write info such as
         energies, temps, etc. to .log file) and collect
         RDF info each JEVERY time step.  Default=10
 
KEVERY = write coordinates (to log and TRAJECT files),
         velocity/quaternion restart info (to the TRAJECT
         file and RDFs (to log file) at each KEVERY step.
         default=100
 
PROD   = production run, at present this means only that
         information for radial distribution functions is
         collected, and printed.  default=.FALSE.
 
DELR   = spacing for radial bins in RDF calculations,
         default=0.02 Angstroms.
 
NPROP  = step number at which to begin collecting data for
         the other properties, such as pressure and
         diffusion constants.  This should be a value
         between 1 and NSTEPS, as it counts off the current
         run's steps.  Default=0.
 
PBCOUT = print PBC coordinates in the end of simulation
         (i.e. all molecules will be contained in one box)
         Default=.FALSE.
 
                         * * *
 
SSBP   = flag to add spherical solvent boundary condition
         using the harmonic restraint potential(V)
              V=0.5*SFORCE*(R-DROFF)**2.
         (Default=.FALSE)
 
SFORCE = the force constant for SSBP in kcal/mol-A**2
         (default: 0.0).
 
CCMS   = flag to add a harmonic potential to constrain
         the center of mass of the QM subsystem. This will
         keep the QM subsystem in QM/MM-MD with spherical
         boundary conditions near the center of sphere.
         (Default=.FALSE)
 
CFORCE = the force constant for CCMS in kcal/mol-A**2
         (default: 0.0).
 
DROFF  = is an array of distances such that V=0 if R