Molecular dynamics

Molecular dynamics (MD) can be used to simulate the evolution of a system in time.

See also

Examples

To perform a MD simulation, first select the corresponding Task:

Task MolecularDynamics

All aspects of the simulation can then be configured using the MolecularDynamics block.

MolecularDynamics
   Barostat
      BulkModulus float
      Duration integer_list
      Equal [None | XYZ | XY | YZ | XZ]
      Pressure float_list
      Scale [XYZ | Shape | X | Y | Z | XY | YZ | XZ]
      Tau float
      Type [None | Berendsen | MTK]
   End
   CalcPressure [True | False]
   Checkpoint
      Frequency integer
   End
   InitialVelocities
      File string
      Temperature float
      Type [Zero | Random | FromFile | Input]
      Values # Non-standard block. See details.
         ...
      End
   End
   NSteps integer
   Preserve
      AngularMomentum [True | False]
      CenterOfMass [True | False]
      Momentum [True | False]
   End
   Print
      System [True | False]
      Velocities [True | False]
   End
   Restart string
   Thermostat
      BerendsenApply [Local | Global]
      ChainLength integer
      Duration integer_list
      FirstAtom integer
      LastAtom integer
      Tau float
      Temperature float_list
      Type [None | Berendsen | NHC]
   End
   TimeStep float
   Trajectory
      SamplingFreq integer
   End
End

General

The time evolution of the system is simulated by numerically integrating the equations of motion. A velocity Verlet integrator is used with a time step set by the TimeStep key. The MD driver will perform NSteps timesteps in total.

Because the overall computational cost depends on NSteps but not on TimeStep, it is desirable to set the timestep as large as possible to maximize the sampled timescale with a given computational budget. However, numerical integration errors grow rapidly as the timestep increases. These errors will cause a loss of energy conservation, crashes, and other artifacts. It is thus important to set the TimeStep value carefully, as its optimal value strongly depends on the studied system and simulated conditions.

As a rule of thumb, reasonable timesteps for systems not undergoing chemical reactions are 10-20 times lower than the period of the fastest vibration mode. Systems containing hydrogen atoms at room temperature can thus be accurately simulated using a 1 fs timestep. Longer timesteps can be safely used for systems containing only heavy atoms (vibration periods scale with the square root of the atomic mass). Conversely, the timestep needs to be made shorter for high-temperature simulations. The same also applies to simulations of chemical reactions, which are usually accompanied by significant transient local heating. The default timestep of 0.25 fs should work for most of these cases.

MolecularDynamics
NSteps
Type:Integer
Default value:1000
Description:The number of steps to be taken in the MD simulation.
TimeStep
Type:Float
Default value:0.25
Unit:Femtoseconds
Description:The time difference per step.

During a long simulation, numerical integration errors will cause some system-wide quantities to drift from their exact values. For example, the system may develop a nonzero net linear velocity, causing an overall translation or flow. Non-periodic (molecular) systems may also develop nonzero angular momentum (overall rotation) and a Brownian motion of their center of mass through space. These problems are corrected by periodically removing any accumulated drift. This feature can be controlled using the Preserve key.

MolecularDynamics
Preserve
Type:Block
Description:Periodically remove numerical drift accumulated during the simulation to preserve different whole-system parameters.
AngularMomentum
Type:Bool
Default value:True
Description:Remove overall angular momentum of the system. This option is ignored for 3D-periodic systems.
CenterOfMass
Type:Bool
Default value:False
Description:Translate the system to keep its center of mass at the coordinate origin. This option is not very useful for 3D-periodic systems.
Momentum
Type:Bool
Default value:True
Description:Remove overall (linear) momentum of the system.

(Re-)Starting a simulation

The state of a system at the beginning of a simulation is defined by the positions and momenta of all atoms. The positions can be set in the input or loaded from a file as described under System definition. Initial velocities are then supplied using the InitialVelocities block.

Probably the most common way to start up a simulation is to draw the initial velocities from a Maxwell-Boltzmann distribution by setting Type=Random and Temperature to a suitable value. Alternatively, velocities can be loaded from an ams.rkf file produced by an earlier simulation using Type=FromFile and File. This is the recommended way to start a production simulation from the results of a short preparation/equilibration run.

Velocities of all atoms in units of Å/fs can also be explicitly defined in the Values block after setting Type=Input. This is mainly useful to repeat or extend simulations done by other programs. For example, velocities can be extracted from the vels or moldyn.vel files used by the standalone ReaxFF program. A simple AWK script is supplied in scripting/standalone/reaxff-ams/vels2ams.awk to help with the conversion.

MolecularDynamics
InitialVelocities
Type:Block
Description:Sets the frequency for printing to stdout and storing the molecular configuration on the .rkf file.
File
Type:String
Description:AMS RKF file containing the initial velocities.
Temperature
Type:Float
Unit:Kelvin
Description:Sets the temperature for the Maxwell-Boltzmann distribution when the type of the initial velocities is set to random, in which case specifying this key is mandatory. ADFinput will use the thermostat temperature as default.
Type
Type:Multiple Choice
Default value:Random
Options:[Zero, Random, FromFile, Input]
Description:Specifies the initial velocities to assign to the atoms. Three methods to assign velocities are available. Zero: All atom are at rest at the beginning of the calculation. Random: Initial atom velocities follow a Maxwell-Boltzmann distribution for the temperature given by the [MolecularDynamics%InitialVelocities%Temperature] keyword. FromFile: Load the velocities from a previous ams result file. Input: Atom’s velocities are set to the values specified in the [MolecularDynamics%InitialVelocities%Values] block, which can be accessed via the Expert AMS panel in ADFinput.
Values
Type:Non-standard block
Description:This block specifies the velocity of each atom when [MolecularDynamics%InitialVelocities%Type] is set to Input. Each row must contain three floating point values (corresponding to the x,y,z component of the velocity vector) and a number of rows equal to the number of atoms must be present, given in the same order as the [System%Atoms] block.

The MD module also supports exact restarts of interrupted simulations by pointing the Restart key to an ams.rkf file. This will restore the entire state of the MD module from the last available checkpoint (if the previous simulation was interrupted) or from the final state (if the previous simulation ended after NSteps). An earlier trajectory can thus be seamlessly extended by increasing NSteps and using Restart.

Note

Restart should be combined with LoadSystem from the same ams.rkf to restore the atomic positions.

Warning

The Restart feature is only intended for exact restarts, so the rest of the MolecularDynamics settings should be the same as in the original run. Only NSteps and engine settings (contents of the Engine block) can always be changed safely across restarts.

Although some MD settings (such as the trajectory sampling options) can in practice be changed without problems, changing others (such as thermostat or barostat settings) will cause the restart to fail or produce physically incorrect results. It is thus strongly recommended to only use Restart for exact continuation and InitialVelocities Type=FromFile together with LoadSystem otherwise.

MolecularDynamics
Restart
Type:String
Description:The path to the ams.rkf file from which to restart the simulation.

Thermostats and barostats

By default, the MD simulation samples the microcanonical (NVE) ensemble. Although this is useful to check energy conservation and other basic physical properties, it does not directly map to common experimental conditions. The canonical (NVT) ensemble can be sampled instead by applying a Thermostat, which serves as a simulated heat bath around the system, keeping its average temperature at a set value.

AMS offers two thermostats with drastically different properties, mode of operation, and applicability, selected using the Type key:

Berendsen

The Berendsen friction thermostat drives the system to a particular target temperature by rescaling the velocities of all atoms in each step. This ensures rapid (exponential) convergence of the temperature with a time constant Tau. However, this thermostat produces an incorrect velocity distribution and should thus be avoided in all situations where correct energy fluctuations are important. Additionally, using a too short time constant Tau tends to cause incorrect equipartition of energy between different degrees of freedom in the system, leading to the “flying ice cube” phenomenon. The time constant Tau should thus be set as large as possible to limit these artifacts while still providing sufficient temperature control. Common values of Tau for condensed-phase systems lie between 100 fs (strong damping, rapid convergence) and 10 ps (weak coupling with minimal artifacts).

This thermostat is mainly useful for systems far from equilibrium, for example during the initial preparation and equilibration phase of a simulation. The NHC thermostat should be preferred where possible.

NHC
This enables a chain of coupled Nosé-Hoover thermostats. This method introduces artificial degrees of freedom representing the heat bath and ensures correct sampling of the canonical ensemble. The combined total energy of the system and the heat bath is conserved and shown in the GUI as Conserved Energy. Checking this quantity for drift and artifacts thus offers a valuable test of the correctness of the simulation. This thermostat exhibits oscillatory relaxation with a period of Tau. It is thus not well suited for systems starting far from equilibrium, because the oscillations may take long to settle. The time constant Tau should be at least comparable to the period of some natural oscillation of the system to ensure efficient energy transfer. It is commonly on the order of hundreds of femtoseconds, although higher values may be used if weak coupling is desired.

Multiple independent thermostats can be used to separately control different regions of the system at the same time. This is done by specifying the Thermostat block multiple times and setting the FirstAtom and/or LastAtom keys to the desired range of atoms. Care needs to be taken to avoid defining thermostats with overlapping atom ranges.

MolecularDynamics
Thermostat
Type:Block
Recurring:True
Description:This block allows to specify the use of a thermostat during the simulation. Depending on the selected thermostat type, different additional options may be needed to characterize the specific thermostat’ behavior.
BerendsenApply
Type:Multiple Choice
Default value:Global
Options:[Local, Global]
Description:Select how to apply the scaling correction for the Berendsen thermostat: - per-atom-velocity (Local) - on the molecular system as a whole (Global).
ChainLength
Type:Integer
Default value:10
Description:Number of individual thermostats forming the NHC thermostat
Duration
Type:Integer List
Description:Specifies how many steps should a transition from a particular temperature to the next one in sequence take.
FirstAtom
Type:Integer
Default value:1
Description:Index of the first atom to be thermostatted
LastAtom
Type:Integer
Default value:0
Description:Index of the last atom to be thermostatted. A value of zero means the last atom in the system.
Tau
Type:Float
Unit:Femtoseconds
Description:The time constant of the thermostat.
Temperature
Type:Float List
Unit:Kelvin
Description:The target temperature of the thermostat.
Type
Type:Multiple Choice
Default value:None
Options:[None, Berendsen, NHC]
Description:Selects the type of the thermostat.

Just like using a Thermostat to control the temperature of the system, a Barostat can be applied to keep the pressure constant by adjusting the volume. This enables sampling the isenthalpic-isobaric (NpH) ensemble by using only a barostat or the isothermal-isobaric (NpT) ensemble by combining a barostat and a thermostat. Unlike thermostats, a barostat always applies to the entire system and there can thus be at most one barostat defined.

AMS offers two barostats with similar properties to the related thermostats:

Berendsen
The Berendsen friction-like isobaric ensemble method rescales the system in each step to drive the pressure towards a target value. Similarly to the Berendsen thermostat, the relaxation is exponential with a time constant Tau. Similar considerations for the choice of Tau apply as in the case of the thermostat, but the value of Tau for the barostat is usually at least several times higher than the corresponding Tau used for the thermostat. This barostat does not have any conserved quantity.
MTK
This enables the Martyna-Tobias-Klein extended Lagrangian barostat, which generates a true isobaric ensemble by integrating the cell parameters as additional degrees of freedom. This barostat is derived from the Andersen-Hoover isotropic barostat and the Parrinello-Rahman-Hoover anisotropic barostat. Like the NHC thermostat, it exhibits oscillatory relaxation unsuitable for systems far from equilibrium. This barostat must always be combined with a NHC thermostat. One instance of such thermostat coupled to the atoms as usual, while a second instance is created internally and coupled to the cell degrees of freedom.
MolecularDynamics
Barostat
Type:Block
Description:This block allows to specify the use of a barostat during the simulation.
BulkModulus
Type:Float
Default value:2200000000.0
Unit:Pascal
Description:An estimate of the bulk modulus (inverse compressibility) of the system for the Berendsen barostat. This is only used to make Tau correspond to the true observed relaxation time constant. Values are commonly on the order of 10-100 GPa (1e10 to 1e11) for solids and 1 GPa (1e9) for liquids (2.2e9 for water). Use 1e9 to match the behavior of standalone ReaxFF.
Duration
Type:Integer List
Description:Specifies how many steps should a transition from a particular pressure to the next one in sequence take.
Equal
Type:Multiple Choice
Default value:None
Options:[None, XYZ, XY, YZ, XZ]
Description:Enforce equal scaling of the selected set of dimensions. They will be barostatted as one dimension according to the average pressure over the components.
Pressure
Type:Float List
Unit:Pascal
Description:Specifies the target pressure.
Scale
Type:Multiple Choice
Default value:XYZ
Options:[XYZ, Shape, X, Y, Z, XY, YZ, XZ]
Description:Dimensions that should be scaled by the barostat to maintain pressure. Selecting Shape means that all three dimensions and also all the cell angles are allowed to change.
Tau
Type:Float
Unit:Femtoseconds
Description:Specifies the time constant of the barostat.
Type
Type:Multiple Choice
Default value:None
Options:[None, Berendsen, MTK]
Description:Selects the type of the barostat.

Temperature and pressure regimes

Arbitrary temperature and pressure regimes can be generated by setting Temperature or Pressure to a list of values, corresponding to the successive set points. This needs to be accompanied by a Duration key specifying the length of each regime segment in steps:

Thermostat
   Temperature 0     300     300     500     500     300
   Duration      100     200     100     200     100
End

Note that there is always N-1 Duration values for N Temperature values. The target temperature of the thermostat in this example will evolve as follows:

  1. Increase linearly from 0 to 300 K over 100 steps.
  2. Stay constant at 300 K for 200 steps.
  3. Increase linearly from 300 to 500 K over 100 steps.
  4. Stay constant at 500 K for 200 steps.
  5. Decrease linearly from 500 to 300 K over 100 steps.
  6. Stay constant at 300 K for the rest of the simulation.

Trajectory sampling and output

A basic principle of the numerical integration of motion in MD is that the changes in the state of the system between successive time steps are small. This means that storing the results of every step is not useful, because all the data is strongly correlated. Instead, a snapshot of the system is taken every N steps, where N is set low enough to still capture the fastest motion of interest but high enough to avoid wasting space due to correlations. The resulting sequence of snapshots is then commonly called the trajectory.

AMS writes the trajectory to the History and MDHistory sections of ams.rkf, according to the settings in the Trajectory block. A snapshot of the system and various MD variables is stored every SamplingFreq timesteps.

The trajectory itself contains only the data needed for subsequent analysis of the dynamics of the system. However, much more data is usually generated on every integration step. This includes, for example, the internal data used by an engine when evaluating the energies and forces. This information is normally discarded after each step, because it is often very large. However, a Checkpoint containing the complete internal state of the MD driver together with a result file generated by the engine is stored every Frequency steps. An interrupted simulation can then be restarted from this checkpoint using the Restart keyword. Additionally, the engine result files called MDStep*.rkf can also be used to extract various engine-specific details about the system, such as the orbitals for QM engines.

MolecularDynamics
Trajectory
Type:Block
Description:Sets the frequency for printing to stdout and storing the molecular configuration on the .rkf file.
SamplingFreq
Type:Integer
Default value:100
Description:Write the the molecular geometry (and possibly other properties) to the .rkf file once every N steps.
Checkpoint
Type:Block
Description:Sets the frequency for storing the entire MD state necessary for restarting the calculation.
Frequency
Type:Integer
Default value:1000
Description:Write the MD state and engine-specific data to the respective .rkf files once every N steps.
CalcPressure
Type:Bool
Default value:False
Description:Calculate the pressure in periodic systems. This may be computationally expensive for some engines that require numerical differentiation. Some other engines can calculate the pressure for negligible additional cost and will always do so, even if this option is disabled.
Print
Type:Block
Description:This block controls the printing of additional information to stdout.
System
Type:Bool
Default value:False
Description:Print the chemical system before and after the simulation.
Velocities
Type:Bool
Default value:False
Description:Print the atomic velocities before and after the simulation.