Molecular dynamics

Overview

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

Input options

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
     dt_md
md_steps - number of (thermalized stage) MD steps
     n_step_md
md_method - molecular dynamics method
     mdtype ( NVE | BERENDSEN | HOOVER | TSCALE )
md_temperature - molecular dynamics temperature (Kelvin)
     temp_md
md_parameter - parameter defining MD method (see below) (2.61)
     param_md
equilibration steps - number of equilibration steps
     n_step_eq
eq_method - equilibration method (2.61)
     eqtype ( NVE | BERENDSEN | HOOVER | TSCALE )
eq_temperature - equilibration target temperature (Kelvin)
     temp_eq
eq_parameter - parameter defining eq method (see below) (2.61)
     param_eq
...
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.

Temperature. The default temperature for the dynamics is 300 K, room temperature, unless the user specifies an input md_temperature.

Initialization. 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.

Equilibration. 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.

Dynamics. 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.

Restart. 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 file.

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 temperature. 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.
  • TSCALE
    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.

Diagnostic tools

There are a few built-in diagnostic tools in the MD output:

  • Temperature – grep output for “MDtemperature” to monitor simulation temperature
  • Total energy – grep output for “MDenergy” to monitor total (KE+electronic) energy
  • Equipartition
    • 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):
      1. temperature using v(t) (i.e., in-synch with positions/forces)
      2. at time t+1/2 (out-of-synch with positions, in-synch with velocities)
      3. time-averaged KE around time t, using v(t-1/2) to v(t+1/2)
      4. approximately time-averaged KE around time t+1/2

Troubleshooting

The dynamics module will print some warning messages if parameters are input that appear suspect:WARNING: Quest MD — Temperature highindicates that the input temperature is outside of the range that SeqQuest believes will work.WARNING: Quest MD — Large time step enteredindicates that the time step entered is too large for SeqQuest to accurately integrate Newton’s equationsWARNING: Quest MD — Small time step enteredindicates that the time step entered is very small, and that too many time steps will be required to obtain meaningful statistics.WARNING: likely uninitialized masssuggests 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).Other problems that may occur:With NVE, a larger energy drift is observedReduce 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.