Molecular dynamics (MD) can be used to simulate the evolution of a system in time.
To perform a MD simulation, first select the corresponding
All aspects of the simulation can then be configured using the
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
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
TimeStep key. The MD driver will perform
NSteps timesteps in
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.
Type: Integer Default value: 1000 Description: The number of steps to be taken in the MD simulation.
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
Type: Block Description: Periodically remove numerical drift accumulated during the simulation to preserve different whole-system parameters.
Type: Bool Default value: True Description: Remove overall angular momentum of the system. This option is ignored for 3D-periodic systems.
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.
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
Probably the most common way to start up a simulation is to draw the initial
velocities from a Maxwell-Boltzmann distribution by setting
Temperature to a suitable value. Alternatively, velocities can be loaded
ams.rkf file produced by an earlier simulation using
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
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.
Type: Block Description: Sets the frequency for printing to stdout and storing the molecular configuration on the .rkf file.
Type: String Description: AMS RKF file containing the initial velocities.
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: 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.
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
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 should be combined with
LoadSystem from the same
to restore the atomic positions.
Restart feature is only intended for exact restarts, so the rest of
MolecularDynamics settings should be the same as in the original run.
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
exact continuation and
InitialVelocities Type=FromFile together with
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
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
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
Tautends to cause incorrect equipartition of energy between different degrees of freedom in the system, leading to the “flying ice cube” phenomenon. The time constant
Taushould thus be set as large as possible to limit these artifacts while still providing sufficient temperature control. Common values of
Taufor 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
NHCthermostat should be preferred where possible.
- 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
Taushould 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
LastAtom keys to the desired range of atoms. Care needs to be taken to avoid
defining thermostats with overlapping atom ranges.
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.
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).
Type: Integer Default value: 10 Description: Number of individual thermostats forming the NHC thermostat
Type: Integer List Description: Specifies how many steps should a transition from a particular temperature to the next one in sequence take.
Type: Integer Default value: 1 Description: Index of the first atom to be thermostatted
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.
Type: Float Unit: Femtoseconds Description: The time constant of the thermostat.
Type: Float List Unit: Kelvin Description: The target temperature of the thermostat.
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:
- The Berendsen friction-like isobaric ensemble method rescales the system in
each step to drive the pressure towards a target value. Similarly to the
Berendsenthermostat, the relaxation is exponential with a time constant
Tau. Similar considerations for the choice of
Tauapply as in the case of the thermostat, but the value of
Taufor the barostat is usually at least several times higher than the corresponding
Tauused for the thermostat. This barostat does not have any conserved quantity.
- 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
NHCthermostat, it exhibits oscillatory relaxation unsuitable for systems far from equilibrium. This barostat must always be combined with a
NHCthermostat. 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.
Type: Block Description: This block allows to specify the use of a barostat during the simulation.
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.
Type: Integer List Description: Specifies how many steps should a transition from a particular pressure to the next one in sequence take.
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.
Type: Float List Unit: Pascal Description: Specifies the target pressure.
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.
Type: Float Unit: Femtoseconds Description: Specifies the time constant of the barostat.
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
Pressure to a list of values, corresponding to the
successive set points. This needs to be accompanied by a
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
The target temperature of the thermostat in this example will evolve as follows:
- Increase linearly from 0 to 300 K over 100 steps.
- Stay constant at 300 K for 200 steps.
- Increase linearly from 300 to 500 K over 100 steps.
- Stay constant at 500 K for 200 steps.
- Decrease linearly from 500 to 300 K over 100 steps.
- 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
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
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
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.
Type: Block Description: Sets the frequency for printing to stdout and storing the molecular configuration on the .rkf file.
Type: Integer Default value: 100 Description: Write the the molecular geometry (and possibly other properties) to the .rkf file once every N steps.
Type: Block Description: Sets the frequency for storing the entire MD state necessary for restarting the calculation.
Type: Integer Default value: 1000 Description: Write the MD state and engine-specific data to the respective .rkf files once every N steps.
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. Type: Block Description: This block controls the printing of additional information to stdout.
Type: Bool Default value: False Description: Print the chemical system before and after the simulation.
Type: Bool Default value: False Description: Print the atomic velocities before and after the simulation.