SeqQuest input - Molecular dynamics
Table of Contents
- Input options
- Details and defaults
- MD methods and thermostats
- Diagnostic tools
This page gives a description of the input section that controls
and modifies how a molecular dynamics run is done.
The MD input section is optional; all the user
need do is add a "do dynamics" instruction in the
command options and the code will
automatically launch a molecular dynamics run, with an NVE
simulation at room temperature.
The default units for the dynamics section are Ry for energies,
bohr for distances, Kelvin for temperature, and femtoseconds for time.
Current constraints, important version notes:
- The atom files must have atomic mass data in them OR
the input file must specify masses for each atom type.
- Default integrator is leapfrog Verlet
The molecular dynamics input section can appear anywhere in the
"run phase" or "geometry relaxation" section of the input file,
begins with the "dynamics data" keyword and
ends with the "end dynamics" keyword.
The following are the common input keywords in this section.
All the keywords are optional, and can appear in any order within this section.
The keywords must be left-justified, but all data input is free-format.
- dynamics input - begin MD input section
- timestep - time step (femtoseconds) for dynamics
- md_steps - number of (thermalized stage) MD steps
- md_method - molecular dynamics method
- mdtype ( NVE | BERENDSEN | HOOVER | TSCALE )
- md_temperature - molecular dynamics temperature (Kelvin)
- md_parameter - parameter defining MD method (see below) (2.61)
- equilibration steps - number of equilibration steps
- eq_method - equilibration method (2.61)
- eqtype ( NVE | BERENDSEN | HOOVER | TSCALE )
- eq_temperature - equilibration target temperature (Kelvin)
- eq_parameter - parameter defining eq method (see below) (2.61)
- end dynamics input - end of MD input section
Details and defaults
The dynamics section enables the SeqQuest program to run molecular
dynamics on the Born-Oppenheimer potential energy surface.
The units of temperature are Kelvin, the unit of time is fs, and
the length of simulation is in number of timesteps.
The default temperature for the dynamics is 300 K, room temperature,
unless the user specifies an input md_temperature.
For initialization (if there is no MD restart file),
velocities are randomly assigned from a
Boltzmann distribution oriented in random directions.
The initial velocities are consistent with a temperature twice
the simulation temperature.
Note: this implicitly assumes the atomic positions are at, or near,
equilibrium positions, such that half the energy will go into
potential modes and thereby leading to an equilibrated temperature
half that of the temperature defined by the initial velocities.
If the atoms are away from initial positions, the internal energy
they possess will result in a higher equilibrated temperature than
the target simulation temperature.
The equilibration phase takes the initialized position and
velocities and runs a number of equilibration steps specified
by the user with the eq_steps keyword.
The default number of equilibration steps is 0,
but 100-1000 steps are typical values.
The default equilibration scheme is TSCALE.
The user can select a any valid thermostat with the
eq_method keyword, and input a parameter to tune that
thermostat with eq_parameter keyword.
After the system has been equilibrated, the MD phase begins.
The user may specify the number of time steps to be run
with the md_steps keyword.
The default number of time steps is 1000000, to effectively run
until the user kills the job, or the job times out.
The user may specify the type of molecular dynamics
to be run (see below) with the "md_method" keyword.
The default dynamics method is NVE.
To continue a dynamics run from a previous set of
positions and velocities, keep the last "lcao.vxyz" file from the
previous run and place it in the running directory of the SeqQuest job.
This MD job will then continue from the last point in the previous run
(with the latest coordinates, velocities, and thermostat conditions).
The vxyz file begins with a line specifying the number of atoms in
the molecule, followed by a title line. These lines are followed by
one line for each atom in the unit cell containing a symbol for the
atom (ignored by SeqQuest), and six floating point fields
(input in free format) containing the x, y, and z coordinates, and
the vx, vy, and vz velocities.
At the end of the vxyz file is a set of additional data necessary to
consistently continue the thermostat.
The vxyz file is a variant of
the XMol XYZ
Format, except that instead of records for only the coordinates,
the velocity fields are added.
If the vxyz file is present, the positions of the atoms in the
SeqQuest input file will be overwritten by the values in the vxyz
MD methods and thermostats
The update scheme used in the molecular dynamics is leapfrog verlet.
The default thermostat for dynamics is NVE, and the default thermostat
for the equilibration is TSCALE.
The available thermostats are:
- NVE - microcanonical ensemble
The default dynamics is NVE, the microcanonical ensemble.
In principle, this method conserves total (electronic + kinetic energy
of ions) energy.
- BERENDSEN - NVT
The BERENDSEN keyword selects the method of Berendsen et al
[J. Chem. Phys. 81, 3684 (1984)] for running molecular
dynamics coupled to an external heat bath to conserve the
The Berendsen friction is expressed as s=(1/(2*tau))(T_c-T)/T_c,
where T is the target temperature, T_c is the current temperature,
and tau is parameter corresponding to a relaxation time, fs time units.
This parameter is set to 400fs as a default when BERENDSEN is invoked.
- HOOVER - NVT
The HOOVER keyword selects the Hoover thermostat,
[Phys. Rev. A 31, 1695 (1985)], which, in principle,
produces states in the canonical ensemble.
The time derivative of the friction coefficient is implemented
as ds=(1/tau^2)(T_c-T)/T_c, where T is the target temperature,
T_c is the current temperature, and tau is a parameter that scales
the connection to the heat bath, again a relaxation time in fs.
This parameter is set to 100fs as a default when HOOVER is invoked.
However, the value of this parameter should be tailored to the
specific application rather than depending upon the default.
This is a poor means to equilibrate a system, and one should
use some temperature-scaling method (TSCALE/BERENDSEN/LANGEVIN)
to equilibrate the system to near the target temperature before
starting a MD simulation with a Hoover thermostat.
- LANGEVIN - NVT
The Langevin thermostat can be very useful for equilibration.
It is not yet available.
Temperature scaled dynamics scales the velocities to match a
temperature at every time step.
Temperature scaled dynamics do not reproduce good physics, and should
only be used for equilibrating velocities at the beginning of an MD run.
The scaling is implemented with a parameter f, such that the new
target temperature T_i at each step is T_i = f*T_c - (1-f)*T,
where T_c is the current simulation temperature and T is the target.
To scale the temperature to the target at every step, set f=0.
The limit f=1 recovers the NVE microcanonical ensemble.
An intermediate value allows ramping the temperature more gently.
The default f is 0.98, to allow for a modest quenching of the system
toward the desired temperature.
There are a few built-in diagnostic tools in the MD output:
- grep output for "MDtemperature" to monitor simulation temperature
- Total energy
- grep output for "MDenergy" to monitor total (KE+electronic) energy
- grep "TEMP" to view step-wise atom-type-resolved temperatures.
- grep "temps" to view current job-average of atom-type-resolved temperatures.
- grep "last100T" for atom-type-resolved temperature averages over 100 step intervals.
- final entry is overall temperature, preceeding entries are temps over atom-types.
- Temperatures computed four separate ways (TEMP#,temps#,#=1-4):
- temperature using v(t) (i.e., in-synch with positions/forces)
- at time t+1/2 (out-of-synch with positions, in-synch with velocities)
- time-averaged KE around time t, using v(t-1/2) to v(t+1/2)
- approximately time-averaged KE around time t+1/2
The dynamics module will print some warning messages if parameters
are input that appear suspect:
Other problems that may occur:
- WARNING: Quest MD -- Temperature high
- indicates that the input temperature is outside of the range
that SeqQuest believes will work.
- WARNING: Quest MD -- Large time step entered
- indicates that the time step entered is too large for SeqQuest
to accurately integrate Newton's equations
- WARNING: Quest MD -- Small time step entered
- indicates that the time step entered is very small, and that
too many time steps will be required to obtain meaningful
- WARNING: likely uninitialized mass
- suggests that an atom file may have been input without field
describing the atomic mass (note: the default mass is 1.0,
and Hydrogen in your system may be incorrectly flagged as
not being initialized).
Return to Top
or Input Manual
or User Guides
or SeqQuest Home
- With NVE, a larger energy drift is observed
- Reduce the size of the time step in the MD.
- Slow equilibration for NVE (or Hoover, Berendsen):
- Use an aggressive temperature scaling method to pre-equilibrate
(e.g., a 100fs Berendsen, 20fs Hoover, TSCALE, etc.).
- Poor equipartition (uneven temperature distribution among atoms)
- Check that size of time step is small enough (5-10 steps, or
more, per fastest vibration period in problem is supposedly needed).
- Check that thermostat relaxation times are (at least)
a factor of 2-10 larger than that fastest vibrational period.
- A too long relaxation time in the thermostat can lead
to poor equilibration that looks like poor equipartition.
Send questions and comments to:
Last updated: January 17, 2014