Molecular dynamics

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

See also

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
   AddMolecules
      AtomTemperature float
      ContactDistance float
      Coords float_list
      CoordsBox float_list
      CoordsSigma float_list
      DeviationAngle float
      Energy float
      EnergySigma float
      FractionalCoords float_list
      FractionalCoordsBox float_list
      FractionalCoordsSigma float_list
      Frequency integer
      MinDistance float
      NumAttempts integer
      Rotate Yes/No
      StartStep integer
      StopStep integer
      System string
      Temperature float
      TemperatureSigma float
      Velocity float
      VelocityDirection float_list
      VelocitySigma float
   End
   Barostat
      BulkModulus float
      ConstantVolume Yes/No
      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
   BondBoost
      Chain
         AtomNames string
         MaxDistances float_list
         MinDistances float_list
      End
      DistanceRestraint string
      NSteps integer
      NumInstances integer
   End
   CRESTMTD
      AddEnergy Yes/No
      GaussianScaling
         ScaleGaussians Yes/No
         ScalingSlope float
      End
      Height float
      NGaussiansMax integer
      NSteps integer
      RestartFile string
      Width float
   End
   CVHD
      Bias
         DampingTemp float
         Delta float
         Height float
      End
      ColVarBB
         at1
            region string
            symbol string
         End
         at2
            region string
            symbol string
         End
         cutoff float
         p integer
         rmax float
         rmin float
      End
      Frequency integer
      StartStep integer
      StopStep integer
      WaitSteps integer
   End
   CalcPressure Yes/No
   Checkpoint
      Frequency integer
      WriteProperties Yes/No
   End
   Deformation
      LatticeVelocity # Non-standard block. See details.
         ...
      End
      LengthRate float_list
      LengthVelocity float_list
      Period float
      ScaleAtoms Yes/No
      StartStep integer
      StopStep integer
      StrainRate # Non-standard block. See details.
         ...
      End
      TargetLattice # Non-standard block. See details.
         ...
      End
      TargetLength float_list
      Type [Linear | Exponential | Sine | Cosine]
   End
   Gravity
      Acceleration float
   End
   HeatExchange
      HeatingRate float
      Method [Simple | HEX | eHEX]
      Sink
         AtomList integer_list
         Box
            Amax float
            Amin float
            Bmax float
            Bmin float
            Cmax float
            Cmin float
         End
         Region string
      End
      Source
         AtomList integer_list
         Box
            Amax float
            Amin float
            Bmax float
            Bmin float
            Cmax float
            Cmin float
         End
         Region string
      End
      StartStep integer
      StopStep integer
   End
   InitialVelocities
      File string
      Temperature float
      Type [Zero | Random | FromFile | Input]
      Values # Non-standard block. See details.
         ...
      End
   End
   NSteps integer
   Plumed
      Input # Non-standard block. See details.
         ...
      End
      Parallel
         nCoresPerGroup integer
         nGroups integer
         nNodesPerGroup integer
      End
   End
   Preserve
      AngularMomentum Yes/No
      CenterOfMass Yes/No
      Momentum Yes/No
   End
   Print
      System Yes/No
      Velocities Yes/No
   End
   Remap
      Type [None | Atoms]
   End
   RemoveMolecules
      Formula string
      Frequency integer
      SafeBox
         Amax float
         Amin float
         Bmax float
         Bmin float
         Cmax float
         Cmin float
         FractionalCoordsBox float_list
      End
      SinkBox
         Amax float
         Amin float
         Bmax float
         Bmin float
         Cmax float
         Cmin float
         FractionalCoordsBox float_list
      End
      StartStep integer
      StopStep integer
   End
   ReplicaExchange
      AllowWrongResults Yes/No
      EWMALength integer
      SwapFrequency integer
      TemperatureFactors float_list
      Temperatures float_list
      nReplicas integer
   End
   Restart string
   Thermostat
      BerendsenApply [Local | Global]
      ChainLength integer
      Duration integer_list
      Region string
      Tau float
      Temperature float_list
      Type [None | Berendsen | NHC]
   End
   TimeStep float
   Trajectory
      PrintFreq integer
      SamplingFreq integer
      TProfileGridPoints integer
      WriteBonds Yes/No
      WriteCharges Yes/No
      WriteGradients Yes/No
      WriteMolecules Yes/No
      WriteVelocities Yes/No
   End
   fbMC
      Frequency integer
      MassRoot float
      NSteps integer
      PrintFreq integer
      StartStep integer
      StepLength float
      StopStep integer
      Temperature float
   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
GUI name:Number of steps
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) and 1D-periodic 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:Yes
GUI name:: Angular momentum
Description:Remove overall angular momentum of the system. This option is ignored for 2D and 3D-periodic systems.
CenterOfMass
Type:Bool
Default value:No
GUI name:: Center of mass
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:Yes
GUI name:Preserve: Total momentum
Description:Remove overall (linear) momentum of the system.

Constrained molecular dynamics

It is possible to keep part of the system frozen in place during MD. This is achieved using the Atom or AtomList keys of the Constraints top-level input block.

Constraints
Atom
Type:Integer
Recurring:True
Description:Fix the position of an atom. Just one integer referring to the index of the atom in the [System%Atoms] block.
AtomList
Type:Integer List
Recurring:True
Description:Fix positions of the specified atoms. A list of integers referring to indices of atoms in the [System%Atoms] block.

Note

During simulations with a changing simulation box (NpT, NpH), the absolute Cartesian coordinates of the frozen atoms cannot be kept fixed. In this case, their fractional cell coordinates are maintained at the original values.

(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 Geometry, 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
GUI name:Initial temperature
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. AMSinput will use the first temperature of the first thermostat as default.
Type
Type:Multiple Choice
Default value:Random
Options:[Zero, Random, FromFile, Input]
GUI name:Initial velocities
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 AMSinput.
Values
Type:Non-standard block
Description:This block specifies the velocity of each atom, in Angstrom/fs, 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
GUI name:Restart from
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 non-overlapping regions of the system at the same time. This is done by first defining appropriate Regions in the System block and then specifying the Thermostat block multiple times with the Region key of each thermostat set to an appropriate region expression.

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]
GUI name:Apply Berendsen
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
GUI name:NHC chain length
Description:Number of individual thermostats forming the NHC thermostat
Duration
Type:Integer List
GUI name:Duration(s)
Description:Specifies how many steps should a transition from a particular temperature to the next one in sequence take.
Region
Type:String
Default value:*
Description:The identifier of the region to thermostat. The default ‘*’ applies the thermostat to the entire system. The value can by a plain region name, or a region expression, e.g. ‘*-myregion’ to thermostat all atoms that are not in myregion, or ‘regionA+regionB’ to thermostat the union of the ‘regionA’ and ‘regionB’. Note that if multiple thermostats are used, their regions may not overlap.
Tau
Type:Float
Unit:Femtoseconds
GUI name:Damping constant
Description:The time constant of the thermostat.
Temperature
Type:Float List
Unit:Kelvin
GUI name:Temperature(s)
Description:The target temperature of the thermostat. You can specify multiple temperatures (separated by spaces). In that case the Duration field specifies how many steps to use for the transition from one T to the next T (using a linear ramp). For NHC thermostat, the temperature may not be zero.
Type
Type:Multiple Choice
Default value:None
Options:[None, Berendsen, NHC]
GUI name:Thermostat
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.
ConstantVolume
Type:Bool
Default value:No
Description:Keep the volume constant while allowing the box shape to change. This is currently supported only by the MTK barostat.
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. You can specify multiple pressures (separated by spaces). In that case the Duration field specifies how many steps to use for the transition from one p to the next p (using a linear ramp).
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
GUI name:Damping constant
Description:Specifies the time constant of the barostat.
Type
Type:Multiple Choice
Default value:None
Options:[None, Berendsen, MTK]
GUI name:Barostat
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. By default, this frequency is also used to print basic thermodynamic parameters of the simulation to the output and log file. Set PrintFreq to override this.

Frequently sampling a long trajectory can generate large volumes of data. If the space usage becomes a concern, one can selectively disable storing some parts of the trajectory to save space using the Write* keys. Note however that this will make it impossible to use some analysis methods on the resulting trajectory:

  • WriteBonds is necessary for reaction network analysis (ChemTraYzer). Disabling WriteBonds also makes AMSmovie show only guessed bonds instead of those calculated by the engine.
  • WriteMolecules is required by the Molecule Fractions panel in AMSmovie.
  • WriteVelocities is required to calculate the velocity autocorrelation functions needed for diffusivity and IR spectra.

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.
PrintFreq
Type:Integer
GUI name:Printing frequency
Description:Print current thermodynamic properties to the output every N steps. By default this is done every SamplingFreq steps.
SamplingFreq
Type:Integer
Default value:100
GUI name:Sample frequency
Description:Write the the molecular geometry (and possibly other properties) to the .rkf file once every N steps.
TProfileGridPoints
Type:Integer
Default value:0
Description:Number of points in the temperature profile. If TProfileGridPoints > 0, a temperature profile along each of the three unit cell axes will be written to the .rkf file. By default, no profile is generated.
WriteBonds
Type:Bool
Default value:Yes
Description:Write detected bonds to the .rkf file.
WriteCharges
Type:Bool
Default value:Yes
Description:Write current atomic point charges (if available) to the .rkf file. Disable this to reduce trajectory size if you do not need to analyze charges.
WriteGradients
Type:Bool
Default value:No
Description:Write gradients (negative of the atomic forces) to the .rkf file.
WriteMolecules
Type:Bool
Default value:Yes
Description:Write the results of molecule analysis to the .rkf file.
WriteVelocities
Type:Bool
Default value:Yes
Description:Write velocities to the .rkf file. Disable this to reduce trajectory size if you do not need to analyze the velocities.
Checkpoint
Type:Block
Description:Sets the frequency for storing the entire MD state necessary for restarting the calculation.
Frequency
Type:Integer
Default value:1000
GUI name:Checkpoint frequency
Description:Write the MD state and engine-specific data to the respective .rkf files once every N steps.
WriteProperties
Type:Bool
Default value:No
Description:Write the properties from the properties section to the ChecoPoint file once every N steps.
CalcPressure
Type:Bool
Default value:No
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:No
Description:Print the chemical system before and after the simulation.
Velocities
Type:Bool
Default value:No
Description:Print the atomic velocities before and after the simulation.

Lattice deformations (volume regimes)

The Deformation block can be used to gradually deform the periodic lattice of the system during a MD simulation. This block can be repeated to define multiple deformations, which will be applied on every MD step in the order in which they are listed in the input.

See also

Example demonstrating various lattice deformations

Warning

While a Deformation and a Barostat can be used at the same time, remember to set the Scale parameter of the barostat so that no dimension is simultaneously being deformed and barostatted.

Each Deformation block will be active on MD steps between StartStep and StopStep.

The time dependence of the lattice parameters is defined by the Type key:

Linear
This deformation type adds the same constant amount to the selected lattice parameters on every simulation step. When used with TargetLattice, LatticeVelocity, or StrainRate, the lattice matrix \(H\) evolves as \(H(t) = H_0 + \Delta H \cdot t\), where \(H_0\) is the lattice at StartStep.
Exponential
When used with StrainRate, this type strains the lattice by the given strain matrix \epsilon on every step, so that the lattice matrix \(H\) evolves as \(H(t) = H_0 (1 + \epsilon)^t\). When used with LengthRate, the length of each lattice vector evolves as \(l(t) = l_0 (1 + r)^t\).
Sine
This is a periodic deformation going from the starting value of the selected lattice parameters to a set target, and then with the same amplitude to the opposite direction from the starting lattice.
Cosine
This periodic deformation oscillates between the starting lattice and a defined target.

The period of the oscillation for the Sine and Cosine types must be set using the Period key.

The extent of the deformation is defined by setting one of the six mutually exclusive input keys. These belong to two groups, depending on whether they operate on the lattice matrix as a whole, or just on the lengths of the individual lattice vectors:

TargetLattice, Velocity, and Rate
These input keys expect a “lattice-like” matrix of numbers, consisting of up to three rows containing up to three values each. Each row contains the components of a single lattice vector and corresponds to a row of the Lattice block. For systems with 1D or 2D periodicity the matrix may be padded to 3x3 with zeros.
TargetLength, LengthVelocity, and LengthRate
These input keys expect a list of up to three values, defining the desired length of each lattice vector or the (absolute or relative) rate of its change.
MolecularDynamics
Deformation
Type:Block
Recurring:True
Description:Deform the periodic lattice of the system during the simulation.
LatticeVelocity
Type:Non-standard block
Description:Velocity of individual lattice vector components in Angstrom/fs. The format is identical to the System%Lattice block. For Type Sine and Cosine, this defines the maximum velocity (at the inflection point).
LengthRate
Type:Float List
Default value:[0.0, 0.0, 0.0]
Description:Relative rate of change of each lattice vector per step.
LengthVelocity
Type:Float List
Default value:[0.0, 0.0, 0.0]
Unit:Angstrom/fs
Description:Change the length of each lattice vector with this velocity. With Type=Exponential, LengthVelocity is divided by the current lattice vector lengths on StartStep to determine a LengthRate, which is then applied on all subsequent steps. For Type Sine and Cosine, this defines the maximum velocity (at the inflection point).
Period
Type:Float
Unit:Femtoseconds
Description:Period of oscillation for Type Sine and Cosine.
ScaleAtoms
Type:Bool
Default value:Yes
Description:Scale the atomic positions together with the lattice vectors. Disable this to deform only the lattice, keeping the coordinates of atoms unchanged.
StartStep
Type:Integer
Default value:1
Description:First step at which the deformation will be applied.
StopStep
Type:Integer
Default value:0
Description:Last step at which the deformation will be applied. If unset or zero, nSteps will be used instead.
StrainRate
Type:Non-standard block
Description:Strain rate matrix to be applied on every step. The format is identical to the System%Lattice block.
TargetLattice
Type:Non-standard block
Description:Target lattice vectors to be achieved by StopStep. The format is identical to the System%Lattice block.
TargetLength
Type:Float List
Default value:[0.0, 0.0, 0.0]
Unit:Angstrom
Description:Target lengths of each lattice vector to be achieved by StopStep. The number of values should equal the periodicity of the system. If a value is zero, the corresponding lattice vector will not be modified.
Type
Type:Multiple Choice
Default value:Linear
Options:[Linear, Exponential, Sine, Cosine]
Description:Function defining the time dependence of the deformed lattice parameters. Linear increments the lattice parameters by the same absolute amount every timestep. Exponential multiplies the lattice parameters by the same factor every timestep. Only StrainRate, LengthRate, and LengthVelocity are supported for Type=Exponential. Sine deforms the system from the starting lattice to TargetLattice/TargetLength and then by the same amount to the opposite direction, while Cosine deforms the system from the starting lattice to the target and back.

Example: Transition from the initial lattice to a 10 Å cube over 1000 steps:

MolecularDynamics
   Deformation
      StopStep 1000
      TargetLattice
         10.0  0.0  0.0
          0.0 10.0  0.0
          0.0  0.0 10.0
      End
   End
End

Example: Oscillate the length of the c lattice vector between the initial value and 20 Å, leaving the a and b vectors unchanged:

MolecularDynamics
   Deformation
      Type Cosine
      Period 100
      TargetLength 0 0 20
   End
End

Example: Stretch the box in the “z” direction by a true exponential strain of 10 ppm per timestep while barostatting the remaining dimensions:

MolecularDynamics
   Barostat Type=MTK Pressure=1e5 Tau=1000 Scale=XY
   Deformation
      Type Exponential
      StrainRate
          0.0  0.0  0.0
          0.0  0.0  0.0
          0.0  0.0 1e-5
      End
   End
End

Molecule Gun: adding molecules during simulation

The molecule gun allows you to “shoot” (add with velocity) a molecule into the simulation box.

See also

The GUI tutorial on the molecule gun.

Molecules can be continuously added to the simulation or only once. The initial position can be pre-set or be random within the simulation box or a part thereof. It can be defined either in the Cartesian or fractional coordinates. The initial velocity can be specified either directly (in Angstrom per femtosecond) or as translational temperature or kinetic energy. Possible applications of the molecule gun include e.g. the simulation of enforced collisions or deposition processes on surfaces.

MolecularDynamics
AddMolecules
Type:Block
Recurring:True
GUI name:Add molecules
Description:This block controls adding molecules to the system (a.k.a. the Molecule Gun). Multiple occurrences of this block are possible. By default, molecules are added at random positions in the simulation box with velocity matching the current system temperature. The initial position can be modified using one of the following keywords: Coords, CoordsBox, FractionalCoords, FractionalCoordsBox. The Coords and FractionalCoords keys can optionally be accompanied by CoordsSigma or FractionalCoordsSigma, respectively.
AtomTemperature
Type:Float
Default value:0.0
Unit:Kelvin
Description:Add random velocity corresponding to the specified temperature to individual atoms of the molecule. This only affects rotational and internal degrees of freedom, not the net translational velocity of the inserted molecule as set by the other options.
ContactDistance
Type:Float
Default value:0.0
Unit:Angstrom
Description:Translate the bullet along the velocity vector until it comes within ContactDistance of any other atom.
Coords
Type:Float List
Unit:Angstrom
Description:Place molecules at or around the specified Cartesian coordinates. This setting takes precedence over other ways to specify initial coordinates of the molecule: [CoordsBox], [FractionalCoords], and [FractionalCoordsBox].
CoordsBox
Type:Float List
Unit:Angstrom
Description:Place molecules at random locations inside the specified box in Cartesian coordinates. Coordinates of the box corners are specified as: Xmin, Xmax, Ymin, Ymax, Zmin, Zmax. This setting is ignored if Coords is used. In AMSinput, if this field is not empty it will be used instead of the default Coords.
CoordsSigma
Type:Float List
Unit:Angstrom
Description:Sigma values (one per Cartesian axis) for a Gauss distribution of the initial coordinates. Can only be used together with Coords.
DeviationAngle
Type:Float
Default value:0.0
Unit:Degree
Description:Randomly tilt the shooting direction up to this angle away from the VelocityDirection vector.
Energy
Type:Float
Unit:Hartree
Description:Initial kinetic energy of the molecule in the shooting direction.
EnergySigma
Type:Float
Default value:0.0
Unit:Hartree
Description:Sigma value for the Gauss distribution of the initial kinetic energy around the specified value. Should only be used together with Energy.
FractionalCoords
Type:Float List
Description:Place molecules at or around the specified fractional coordinates in the main system’s lattice. For non-periodic dimensions a Cartesian value in Angstrom is expected. This setting is ignored if [Coords] or [CoordsBox] is used.
FractionalCoordsBox
Type:Float List
Description:Place molecules at random locations inside the box specified as fractional coordinates in the main system’s lattice. Coordinates of the box corners are specified as: Xmin, Xmax, Ymin, Ymax, Zmin, Zmax. For non-periodic dimensions the Cartesian value in Angstrom is expected. This setting is ignored if [Coords], [CoordsBox], or [FractionalCoords] is used.
FractionalCoordsSigma
Type:Float List
Description:Sigma values (one per axis) for a Gauss distribution of the initial coordinates. For non-periodic dimensions the Cartesian value in Angstrom is expected. Can only be used together with FractionalCoords.
Frequency
Type:Integer
Default value:0
Description:A molecule is added every [Frequency] steps after the StartStep. There is never a molecule added at step 0.
MinDistance
Type:Float
Default value:0.0
Unit:Angstrom
Description:Keep the minimal distance to other atoms of the system when adding the molecule.
NumAttempts
Type:Integer
Default value:10
Description:Try adding the molecule up to the specified number of times or until the MinDistance constraint is satisfied. If all attempts fail a message will be printed and the simulation will continue normally.
Rotate
Type:Bool
Default value:No
Description:Rotate the molecule randomly before adding it to the system.
StartStep
Type:Integer
Default value:0
Description:Step number when the first molecule should be added. After that, molecules are added every Frequency steps. For example, ff StartStep=99 and Frequency=100 then a molecule will be added at steps 99, 199, 299, etc… No molecule will be added at step 0, so if StartStep=0 the first molecule is added at the step number equal to [Frequency].
StopStep
Type:Integer
Description:Do not add this molecule after the specified step.
System
Type:String
Description:String ID of the [System] that will be added with this ‘gun’. The lattice specified with this System is ignored and the main system’s lattice is used instead. AMSinput adds the system at the coordinates of the System (thus setting Coords to the center of the System).
Temperature
Type:Float
Unit:Kelvin
Description:Initial energy of the molecule in the shooting direction will correspond to the given temperature.
TemperatureSigma
Type:Float
Default value:0.0
Unit:Kelvin
Description:Sigma value for the Gauss distribution of the initial temperature the specified value. Should only be used together with Temperature.
Velocity
Type:Float
Unit:Angstrom/fs
Description:Initial velocity of the molecule in the shooting direction.
VelocityDirection
Type:Float List
Description:Velocity direction vector for aimed shooting. It will be random if not specified. In AMSinput add one or two atoms (which may be dummies). One atom: use vector from center of the system to add to that atom. Two atoms: use vector from the first to the second atom.
VelocitySigma
Type:Float
Default value:0.0
Unit:Angstrom/fs
Description:Sigma value for the Gauss distribution of the initial velocity around the specified value. Should only be used together with Velocity.

Removing molecules during simulation

Molecules detected by the AMS driver can also be removed from the system. This feature can, for example, be used to remove reaction products.

MolecularDynamics
RemoveMolecules
Type:Block
Recurring:True
GUI name:Remove molecules
Description:This block controls removal of molecules from the system. Multiple occurrences of this block are possible.
Formula
Type:String
Description:Molecular formula of the molecules that should be removed from the system. The order of elements in the formula is very important and the correct order is: C, H, all other elements in the strictly alphabetic order. Element names are case-sensitive, spaces in the formula are not allowed. Digit ‘1’ must be omitted. Valid formula examples: C2H6O, H2O, O2S. Invalid formula examples: C2H5OH, H2O1, OH, SO2. Invalid formulas are silently ignored. Use * to remove any molecule, which must be combined with SinkBox or SafeBox.
Frequency
Type:Integer
Default value:0
Description:The specified molecules are removed every so many steps after the StartStep. There is never a molecule removed at step 0.
SafeBox
Type:Block
Description:Part of the simulation box where molecules may not be removed. Only one of the SinkBox or SafeBox blocks may be present. If this block is present the molecule will not be removed if any of its atoms is within the box. For a periodic dimension it is given as a fraction of the simulation box (the full 0 to 1 range by default). For a non-periodic dimension it represents absolute Cartesian coordinates in Angstrom.
Amax
Type:Float
Description:Coordinate of the upper bound along the first axis.
Amin
Type:Float
Description:Coordinate of the lower bound along the first axis.
Bmax
Type:Float
Description:Coordinate of the upper bound along the second axis.
Bmin
Type:Float
Description:Coordinate of the lower bound along the second axis.
Cmax
Type:Float
Description:Coordinate of the upper bound along the third axis.
Cmin
Type:Float
Description:Coordinate of the lower bound along the third axis.
FractionalCoordsBox
Type:Float List
GUI name:Safe box
Description:Do not remove molecules that are (partly) inside the safe box. Borders of the safe box specified as: Amin, Amax, Bmin, Bmax, Cmin, Cmax. For periodic dimensions fractional coordinates between 0 and 1 and for non-periodic dimensions Cartesian values in Angstrom are expected.
SinkBox
Type:Block
Description:Part of the simulation box where matching molecules will be removed. By default, molecules matching the formula will be removed regardless of their location. If this block is present then such a molecule will only be removed if any of its atoms is within the box. For a periodic dimension it is given as a fraction of the simulation box (the full 0 to 1 range by default). For a non-periodic dimension it represents absolute Cartesian coordinates in Angstrom.
Amax
Type:Float
Description:Coordinate of the upper bound along the first axis.
Amin
Type:Float
Description:Coordinate of the lower bound along the first axis.
Bmax
Type:Float
Description:Coordinate of the upper bound along the second axis.
Bmin
Type:Float
Description:Coordinate of the lower bound along the second axis.
Cmax
Type:Float
Description:Coordinate of the upper bound along the third axis.
Cmin
Type:Float
Description:Coordinate of the lower bound along the third axis.
FractionalCoordsBox
Type:Float List
GUI name:Sink box
Description:Remove molecules that are (partly) inside the sink box. Borders of the sink box specified as: Amin, Amax, Bmin, Bmax, Cmin, Cmax. For periodic dimensions fractional coordinates between 0 and 1 and for non-periodic dimensions Cartesian values in Angstrom are expected.
StartStep
Type:Integer
Default value:0
Description:Step number when molecules are removed for the first time. After that, molecules are removed every [Frequency] steps. For example, if StartStep=99 and Frequency=100 then molecules will be removed at steps 99, 199, 299, etc… No molecule will be removed at step 0, so if StartStep=0 the first molecules are removed at the step number equal to [Frequency].
StopStep
Type:Integer
Description:Do not remove the specified molecules after this step.

Warning

When there is a Molecules%AdsorptionSupportRegion defined, the molecule formulas depend on whether the molecule is adsorbed or not.

Accelerated dynamics

The PLUMED library support in AMS

PLUMED is a plugin that works with various MD programs and is also available in AMS. It can be used for on-the-fly analysis of the dynamics, or to perform a wide variety of free energy methods. The interface with the plugin is really simple: you just need to specify the PLUMED input in the MolecularDynamics%Plumed%Input block and it will be passed to the library “as is”. At each MD step, the current state of the system will be passed to the plugin to be updated according to the PLUMED input.

MolecularDynamics
Plumed
Type:Block
Description:Input for PLUMED. The parallel option is still experimental.
Input
Type:Non-standard block
Description:Input for PLUMED. Contents of this block is passed to PLUMED as is.
Parallel
Type:Block
Description:Options for double parallelization, which allows to split the available processor cores into groups working through all the available tasks in parallel, resulting in a better parallel performance. The keys in this block determine how to split the available processor cores into groups working in parallel.
nCoresPerGroup
Type:Integer
GUI name:Cores per group
Description:Number of cores in each working group.
nGroups
Type:Integer
GUI name:Number of groups
Description:Total number of processor groups. This is the number of tasks that will be executed in parallel.
nNodesPerGroup
Type:Integer
GUI name:Nodes per group
Description:Number of nodes in each group. This option should only be used on homogeneous compute clusters, where all used compute nodes have the same number of processor cores.

Metadynamics for Conformer-Rotamer Ensemble Sampling (CREST-MTD)

This is a very specific implementation of metadynamics that is only meant for exploration of conformer space, as used in the Conformer-Rotamer Ensemble Sampling (CREST) approach. It is an RMSD-based metadynamics that places a new 1-dimensional Gaussian on the potential energy surface every NSteps steps. The Gaussian is a function of the RMSD from the current geometry, which will always be zero at the moment of placement. All Gaussians have the same maximum height (Height) and the same width (Width). There is an upper limit of NGaussiansMax to the number of Gaussians present in the system, and when this is exceeded the oldest Gaussian is removed. By default the Gaussians are gradually introduced, using a scaling function that increases from 0 to 1 with the simulation time. The keyword ScalingSlope in the inputblock GaussianScaling determines the slope of the scaling function with respect to time. The default value of 0.03 \(steps^{-1}\) yields a scaling factor of nearly 1 after 100 steps. The keyword ScaleGaussians in the inputblock GaussianScaling determines whether the Gaussians are scaled at all.

MolecularDynamics
CRESTMTD
Type:Block
GUI name:CREST_MTD
Description:Input for CREST metadynamics simulation.
AddEnergy
Type:Bool
Default value:No
Description:Add the bias energy to the potential energy (to match the gradients)
GaussianScaling
Type:Block
Description:Options for gradual introduction of the Gaussians
ScaleGaussians
Type:Bool
Default value:Yes
Description:Introduce the Gaussians gradually, using a scaling function
ScalingSlope
Type:Float
Default value:0.03
Description:Slope of the scaling function for the Gaussians with respect to time
Height
Type:Float
Unit:Hartree
Description:The height of the Gaussians added
NGaussiansMax
Type:Integer
Description:Maximum number of Gaussians stored
NSteps
Type:Integer
Description:Interval of Gaussian placement
RestartFile
Type:String
Description:Filename for file from which to read data on Gaussians placed previously.
Width
Type:Float
Unit:Bohr
Description:The width of the Gaussians added in terms of the RMSD

By default the gradients of the Gaussians with respect to the atom coordinates are added to the state gradients, but the value of the Gaussians is not added to the energy, which is common in metadynamics. For testing purposes it can be useful to add the Gaussian value to the energy, and this can be done with the keyword AddEnergy. Irrespective of this choice, the energy value of the Gaussian at each geometry is printed in the ams.rkf file in the section CrestMTDHistory. It is also possible to use Gaussians from an earlier CREST-MTD simulation, using the keyword RestartFile.

Collective Variable-driven HyperDynamics (CVHD)

The Collective Variable-driven HyperDynamics is a molecular dynamics acceleration method that allows observation of rare events by filling energy minima with a bias potential. In this sense it is similar to metadynamics. The difference of the hyperdynamics is that it ensures that the bias disappears in the transition state region. This difference allows hyperdynamics to calculate the rate of slow processes, for example the ignition phase of combustion.

See also

The GUI tutorial on CVHD.

The CVHD implementation in AMS follows the algorithm described in K.M. Bal, E.C. Neyts, JCTC, 11 (2015)

The StartStep, Frequency, StopStep, and WaitSteps keys define when and how often the bias potential is added, and when it is removed. The Bias block defines parameters of the bias potential peaks and the ColVarBB block describes parameters of the bond-breaking collective variable.

MolecularDynamics
CVHD
Type:Block
Recurring:True
GUI name:CVHD
Description:Input for the Collective Variable-driven HyperDynamics (CVHD).
Bias
Type:Block
Description:The bias is built from a series of Gaussian peaks deposited on the collective variable axis every [Frequency] steps during MD. Each peak is characterized by its (possibly damped) height and the RMS width (standard deviation).
DampingTemp
Type:Float
Default value:0.0
Unit:Kelvin
GUI name:Bias damping T
Description:During well-tempered hyperdynamics the height of the added bias is scaled down with an exp(-E/kT) factor [PhysRevLett 100, 020603 (2008)], where E is the current value of the bias at the given CV value and T is the damping temperature DampingTemp. If DampingTemp is zero then no damping is applied.
Delta
Type:Float
Description:Standard deviation parameter of the Gaussian bias peak.
Height
Type:Float
Unit:Hartree
Description:Height of the Gaussian bias peak.
ColVarBB
Type:Block
Recurring:True
GUI name:Collective Variable
Description:Description of a bond-breaking collective variable (CV) as described in [Bal & Neyts, JCTC, 11 (2015)]. A collective variable may consist of multiple ColVar blocks.
at1
Type:Block
Description:Specifies the first bonded atom in the collective variable.
region
Type:String
Default value:*
Description:Restrict the selection of bonded atoms to a specific region. If this is not set, atoms anywhere in the system will be selected.
symbol
Type:String
Description:Atom type name of the first atom of the bond. The name must be as it appears in the System block. That is, if the atom name contains an extension (e.g C.1) then the full name including the extension must be used here.
at2
Type:Block
Description:Specifies the second bonded atom in the collective variable.
region
Type:String
Default value:*
Description:Restrict the selection of bonded atoms to a specific region. If this is not set, atoms anywhere in the system will be selected.
symbol
Type:String
Description:Atom type name of the second atom of the bond. The value is allowed to be the same as [at1], in which case bonds between atoms of the same type will be included.
cutoff
Type:Float
Default value:0.3
GUI name:Bond order cutoff
Description:Bond order cutoff. Bonds with BO below this value are ignored when creating the initial bond list for the CV. The bond list does not change during lifetime of the variable even if some bond orders drop below the cutoff.
p
Type:Integer
Default value:6
GUI name:Exponent p
Description:Exponent value p used to calculate the p-norm for this CV.
rmax
Type:Float
Unit:Angstrom
GUI name:R max
Description:Max bond distance parameter Rmax used for calculating the CV. It should be close to the transition-state distance for the corresponding bond.
rmin
Type:Float
Unit:Angstrom
GUI name:R min
Description:Min bond distance parameter Rmin used for calculating the CV. It should be close to equilibrium distance for the corresponding bond.
Frequency
Type:Integer
Description:Frequency of adding a new bias peak, in steps. New bias is deposited every [Frequency] steps after [StartStep] if the following conditions are satisfied: the current CV value is less than 0.9 (to avoid creating barriers at the transition state), the step number is greater than or equal to [StartStep], and the step number is less than or equal to [StopStep].
StartStep
Type:Integer
Description:If this key is specified, the first bias will be deposited at this step. Otherwise, the first bias peak is added at the step number equal to the Frequency parameter. The bias is never deposited at step 0.
StopStep
Type:Integer
Description:No bias will be deposited after the specified step. The already deposited bias will continue to be applied until the reaction event occurs. After that no new CVHD will be started. By default, the CVHD runs for the whole duration of the MD calculation.
WaitSteps
Type:Integer
Description:If the CV value becomes equal to 1 and remains at this value for this many steps then the reaction event is considered having taken place. After this, the collective variable will be reset and the bias will be removed.

During a CVHD calculation, the following variables are saved to the MDHistory section of the RKF file, in addition to other MD properties:

  • BiasEnergy - value the bias energy at the current MD step, in Hartree.
  • MaxBiasEnergy - max BiasEnergy since the last sampling step.
  • BoostFactor - the boost factor at the given MD step. The boost factor is calculated at each MD step as \(boost = e^{E_{bias}/kT}\), where T is the MD ensemble temperature.
  • MaxBoostFactor - max BoostFactor value since the last sampling step.
  • HyperTime - boosted MD time, in femtoseconds, which is a sum of the hyper-time steps calculated from the current boost factor and the MD time step as \(\Delta t_{boost} = boost * \Delta t\). In hyperdynamics, the hyper-time value is directly related to the rate of the process boosted by the corresponding collective variable.

Temperature Replica Exchange

Sampling of rare events can be accelerated using the Replica Exchange Molecular Dynamics (REMD) method, also known as Parallel Tempering. This method runs multiple replicas (copies) of the simulated system in parallel, each in a different ensemble. In the case of Temperature REMD, these ensembles are all NVT or NpT, each at a different temperature. Periodically, Monte Carlo swaps are attempted between neighboring ensembles. If the current configuration of replica A has a sufficient Boltzmann probability in the ensemble of replica B (and vice versa), the two configurations will be swapped. This causes high-energy configuration to migrate into the high-temperature replicas while low-energy configurations eventually end up in the coldest ensemble. This facilitates the crossing of energy barriers in the high-temperature ensembles while keeping the coldest replica at a given temperature of interest. Because each replica always samples an unbiased ensemble, any property can be calculated using standard MD analysis methods without special preparation.

The method is controlled using the ReplicaExchange block:

MolecularDynamics
ReplicaExchange
Type:Block
Description:This block is used for (temperature) Replica Exchange MD (Parallel Tempering) simulations.
AllowWrongResults
Type:Bool
Default value:No
Description:Allow combining Replica Exchange with other features when the combination is known to produce physically incorrect results.
EWMALength
Type:Integer
Default value:10
Description:Length of the exponentially weighted moving average used to smooth swap probabilities for monitoring. This value is equal to the inverse of the EWMA mixing factor.
SwapFrequency
Type:Integer
Default value:100
Description:Attempt an exchange every N steps.
TemperatureFactors
Type:Float List
Description:This is the ratio of the temperatures of two successive replicas. The first value sets the temperature of the second replica with respect to the first replica, the second value sets the temperature of the third replica with respect to the second one, and so on. If there are fewer values than nReplicas, the last value of TemperatureFactor is used for all the remaining replicas.
Temperatures
Type:Float List
Description:List of temperatures for all replicas except for the first one. This is mutually exclusive with TemperatureFactors. Exactly nReplicas-1 temperature values need to be specified, in increasing order. The temperature of the first replica is given by [Thermostat%Temperature].
nReplicas
Type:Integer
Default value:1
GUI name:Number of replicas
Description:Number of replicas to run in parallel.

The number of replicas set by nReplicas must never exceed the total number of processors used for the simulation. If possible, the total number of processors should be an integer multiple of nReplicas to ensure good load balancing.

The temperature of the base (coldest) replica is determined by the Thermostat input block, just like in an ordinary MD simulation. There are two ways to set the temperatures of the remaining replicas, either using Temperatures or TemperatureFactors. The latter is typically more convenient, as it makes it easy to set up the optimal geometric progression of temperatures. In the simplest case, it is enough to supply just a single value in TemperatureFactors, setting the common ratio of temperatures of any two adjacent replicas.

SwapFrequency should be set as low as practical for maximum efficiency. The value of this parameter isn’t critical because it doesn’t affect the validity of the results. However, setting it too high will decrease overall acceleration by missing some opportunities to exchange. Conversely, using a value that is too low will increase the communication overhead and lead to useless back-and-forth swaps between adjacent replicas. Ideally, SwapFrequency should be comparable to the correlation time of the system to ensure that individual exchange attempts are uncorrelated.

The trajectory of each replica is written to a separate RKF file: ams.rkf for the base replica and replicaX.rkf for the other replicas. One can easily switch between these files in the GUI using File → Related Files. In addition to data present in any MD trajectory, these files also contain an extra section ReplicaExchangeHistory with the following data items written every SwapFrequency steps:

  • AvgSwapProbability – Average swap acceptance for each pair of replicas, smoothed using an exponentially weighted moving average with a mixing factor equal to the inverse of EWMALength.
  • {Min,Max,Mean,StdDev}PotentialEnergy – Statistics of the potential energy for each ensemble over the last SwapFrequency steps.
  • SystemInEnsemble i – Identifies the system (continuous trajectory) currently running in ensemble (replica) i.
  • EnsembleOfSystem i – Inverse mapping of SystemInEnsemble, giving the current replica number in which the system number i runs.
  • TemperatureOfSystem i – Equivalent to EnsembleOfSystem using temperatures instead of integer numbers to identify ensembles.

These data items can be plotted using the MD Replica Exchange menu in AMSmovie. For example, plotting TemperatureOfSystem or EnsembleOfSystem is useful to visualize the migration of each system through the space of ensembles, where each curve represents one continuous trajectory. Plotting potential energy statistics or average acceptances facilitates tuning the number of replicas and their temperatures to achieve efficient acceleration. The replica exchange method can only work when the potential energy distributions of adjacent ensembles have a sufficient overlap. This can be easily seen by comparing MaxPotentialEnergy of ensemble i with MinPotentialEnergy of ensemble i+1. The optimal degree of overlap is such that leads to approximately 20 % of swap attempts getting accepted. The acceptance of swaps can be monitored by plotting AvgSwapProbability and the corresponding TemperatureFactors can then be adjusted to keep it near the optimal value.

Bond Boost Method

The bond boost method implemented here is described in [3]. In this method, the distances between atoms that are relevant to the reaction of interest are calculated to determine the orientation of the reactant molecules. If a suitable initial configuration is recognized, an additional restraint energy (possibly consisting of more than one term) is added to the system that is intended to stretch or compress bonds at a pre-defined rate such that this additional energy can help achieve the energy to cross the reaction barrier. A single term of the restraint energy is depends on the restraints type and its parameters, see restraint definitions for details. If more than one suitable configuration is found then the one with the smallest sum of distances is used to create the restraints.

See also

The GUI tutorial on the Bond Boost Method.

MolecularDynamics
BondBoost
Type:Block
Recurring:True
Description:Forced reaction (bond boost) definitions. Multiple BondBoost blocks may be specified, which will be treated independently.
Chain
Type:Block
Description:Specifications of a chain of atoms. When a chain is detected the distance restraints will be activated. No other chain of this type will be detected while any restraints for this chain is active.
AtomNames
Type:String
Description:Atom names specifying the chain. An atom name can optionally be followed by ‘@’ and a region name, in this case only atoms of this type from the given region will be matched. A leading ‘@’ followed by a number indicates that this position in the chain must be occupied by the atom found earlier at the specified position in the chain. For example “O H N C @1” indicates that the last atom in the chain of the five atoms must be the first oxygen, thus defining a 4-membered ring. This is the only way to define a ring because implicit rings will not be detected. For example, “O H N C O” does not include rings.
MaxDistances
Type:Float List
Unit:Angstrom
Description:Maximum distances for each pair of atoms in the chain. The number of distances must be one less than the number of AtomNames.
MinDistances
Type:Float List
Unit:Angstrom
Description:Minimum distances for each pair of atoms in the chain. The number of distances must be one less than the number of AtomNames.
DistanceRestraint
Type:String
Recurring:True
Description:Specify two atom indices followed by the distance in Angstrom, the ForceConstant (in a.u.) and, optionally, the profile type and F(Inf) (in a.u.). This restraint will try to keep the distance between the two specified atoms at the given value. For periodic systems this restraint follows the minimum image convention. Each index indicates position of the corresponding atom in the AtomNames key. Currently recognized restraint profile types: Harmonic (default), Hyperbolic, Erf.
NSteps
Type:Integer
GUI name:Boost lifetime
Description:Number of steps the restraints will remain active until removed. Atoms participating in one reaction are not available for the given number of steps.
NumInstances
Type:Integer
Default value:1
GUI name:Number of instances
Description:Number of reactions of this type taking place simultaneously.

For example:

MolecularDynamics
   BondBoost
      NSteps 10000
      Chain
         AtomNames     N.1   C.t    O     H     @1
         MinDistances     3.0   1.2   1.5   0.8
         MaxDistances     3.8   3.0   2.5   1.5
      End
#     Restraints         i  j   R_0    FC
      DistanceRestraint  1  2   1.5   0.05   Hyperbolic   0.7
      DistanceRestraint  2  3   3.5   0.03   Hyperbolic   0.5
      DistanceRestraint  3  4   1.0   0.03   Hyperbolic   0.5
   End
End

The AtomNames, MinDistances, and MaxDistances keys constitute the atom chain definitions for the initial configuration. Thus, the example above defines a chain of atoms N.1-C.t-O-H-N.1 with R(N.1-C) in the (3.0,3.8) range, R(C.t-O) in the (1.2,3.0) range, etc.. In this example, the last atoms in the chain is required to be the same as the fist one, thus defining a ring. The specified restraints will push pairs of atoms C-N and O-H close together, which will hopefully let them form a bond, and pull atoms C.t and O away from each other, thus breaking the C-O bond. The restraint type is set to Hyperbolic to avoid very large forces that would otherwise result from a harmonic potential at a large deviation.

Note

When detecting coordinates, the program uses the full atom name and not just the element name. An atom name consists of the element name optionally followed by a period and a suffix, just like N.1 and C.t in the example above. Using extended names for some atoms one may allow only a subset of bonds to be boosted.

Force bias Monte Carlo (fbMC)

Configurational sampling and structural relaxation can be accelerated by the Force bias Monte Carlo method. Unlike molecular dynamics, this approach does not use the equations of motion or atomic velocities. Instead, all atoms are randomly displaced in each fbMC steps. AMS uses an uniform-acceptance force bias Monte Carlo flavour of fbMC called tfMC, as published by Mees et al. [4] Unlike other well-known Monte Carlo methods such as the Metropolis algorithm, fbMC does not generate completely random moves and then test whether these should be accepted or rejected. Instead, the atomic displacements are generated so that atoms preferably move in the direction of the forces acting on them. This special biasing enables fbMC to accept every move, sampling the correct canonical ensemble with high efficiency.

In AMS the fbMC method is implemented as an add-on to molecular dynamics, enabling mixed MC-MD simulations. After every Frequency MD steps, the MD procedure is stopped and fbMC takes over, executing NSteps Monte Carlo steps. Once the fbMC procedure completes, the MD simulation continues from the new geometry generated by fbMC while reusing atomic velocities from before the start of fbMC.

The fbMC module requires three main input settings. Apart from the NSteps and Frequency defining the mix of molecular dynamics and Monte Carlo, Temperature needs to be set as desired to sample a particular canonical (NVT) ensemble. Although the fbMC procedure always runs in constant temperature mode, the rest of the MD simulation can in principle use any ensemble (with or without a thermostat or barostat).

Note

To run a pure fbMC simulation with no molecular dynamics component, leave Frequency set to 1 (the default), set the MD TimeStep and InitialVelocities to zero, and do not use other MD features such as a thermostat or a barostat.

Note

AMS currently cannot write sampled geometries to the trajectory during the fbMC procedure. To see the evolution of the system during a long fbMC run, split it into multiple shorter fbMC passes to let the MD driver sample the current state between them. To do this, set the NSteps for fbMC to the desired trajectory sampling interval, keep Frequency set to 1 and set Trajectory SamplingFreq to 1.

The amount of acceleration provided by fbMC depends on the StepLength parameter, which defines the maximum displacement of the lightest atom in the system. The default value of 0.1 Å is relatively conservative. Increasing StepLength leads to faster evolution of the system, but setting it too high can lead to inaccurate sampling or cause crashes by generating geometries that are too distorted.

An estimate of the timescale associated with fbMC (tfMC) steps will be printed to the output.

MolecularDynamics
fbMC
Type:Block
Recurring:True
GUI name:fbMC
Description:This block sets up force bias Monte Carlo interleaved with the molecular dynamics simulation.
Frequency
Type:Integer
Default value:1
Description:Run the fbMC procedure every Frequency MD steps.
MassRoot
Type:Float
Default value:2.0
Description:Inverse of the exponent used to mass-weight fbMC steps.
NSteps
Type:Integer
GUI name:Number of steps
Description:Number of fbMC steps to perform on every invocation of the procedure.
PrintFreq
Type:Integer
GUI name:Printing frequency
Description:Print current thermodynamic properties to the output every N fbMC steps. This defaults to the PrintFreq set in the Trajectory block. Setting this to zero disables printing fbMC steps.
StartStep
Type:Integer
Default value:1
Description:First step at which the fbMC procedure may run.
StepLength
Type:Float
Default value:0.1
Unit:Angstrom
Description:Maximum allowed displacement of the lightest atom in the system in each Cartesian coordinate in one fbMC step.
StopStep
Type:Integer
Default value:0
Description:Last step at which the fbMC procedure may run. If unset or zero, there is no limit.
Temperature
Type:Float
Unit:Kelvin
Description:Temperature used for fbMC.

Non-equilibrium MD (NEMD)

T-NEMD for thermoconductivity: heat exchange

There are different methods to study thermal conductivity using non-equilibrium molecular dynamics (NEMD). A common feature of these methods is that they require the system to be divided into three or more zones, each with its own thermostat and other properties. One method maintains the temperature of the heat source and the heat sink zones at the given temperature using two different thermostats and measures the amount of heat transferred. These method does not require any special features besides a standard thermostat and a possibility to calculate the amount of heat accumulated by the thermostat per unit of time. The accumulated thermostat energies are available in the MDHistory section of ams.rkf file, in variables called ‘XXXXstat#Energy’, where XXXX is a 4-letter abbreviation of the thermo-/barostat (‘BerT’ for a Berendset thermostat, ‘NHCT’ for an NHC thermostat, ‘NHTB’ for an MTK barostat, etc.) and ‘#’ is a 1-digit index of the thermo-/barostat.

In the other method, the heat flow is constant and the induced temperature gradient is measured. This method is implemented in AMS in three variants: a simple heat exchange, HEX [1] and eHEX [2]. In the simple heat exchange method the atom velocities are scaled up (or down) by a factor corresponding to the amount of heat deposited to the heat source (or withdrawn from the heat sink) without paying attention to the conservation of the total momentum of the heat source (or sink). In the HEX method the velocities are scaled in such a way that the total momentum is conserved. This, however, introduces a small but measurable drift in the total energy. The eHEX algorithm improves upon the HEX by adding a third-order time-integration correction to remove the drift.

This method is controlled by the HeatExchange sub-block of the MolecularDynamics block:

MolecularDynamics
HeatExchange
Type:Block
Recurring:True
GUI name:Heat exchange
Description:Input for the heat-exchange non-equilibrium MD (T-NEMD).
HeatingRate
Type:Float
Unit:Hartree/fs
Description:Rate at which the energy is added to the Source and removed from the Sink. A heating rate of 1 Hartree/fs equals to about 0.00436 Watt of power being transferred through the system.
Method
Type:Multiple Choice
Default value:Simple
Options:[Simple, HEX, eHEX]
Description:Heat exchange method used. Simple: kinetic energy of the atoms of the source and sink regions is modified irrespective of that of the center of mass (CoM) of the region (recommended for solids). HEX: kinetic energy of the atoms of these regions is modified keeping that of the corresponding CoM constant. eHEX: an enhanced version of HEX that conserves the total energy better (recommended for gases and liquids).
Sink
Type:Block
Description:Defines the heat sink region (where the heat will be removed).
AtomList
Type:Integer List
GUI name:Sink region
Description:The atoms that are part of the sink. This key is ignored if the [Box] block or [Region] key is present.
Box
Type:Block
Description:Part of the simulation box (in fractional cell coordinates) defining the heat sink. If this block is specified, then by default, the whole box in each of the three dimensions is used, which usually does not make much sense. Normally, you will want to set the bounds along one of the axes.
Amax
Type:Float
Default value:1.0
Description:Coordinate of the upper bound along the first axis.
Amin
Type:Float
Default value:0.0
Description:Coordinate of the lower bound along the first axis.
Bmax
Type:Float
Default value:1.0
Description:Coordinate of the upper bound along the second axis.
Bmin
Type:Float
Default value:0.0
Description:Coordinate of the lower bound along the second axis.
Cmax
Type:Float
Default value:1.0
Description:Coordinate of the upper bound along the third axis.
Cmin
Type:Float
Default value:0.0
Description:Coordinate of the lower bound along the third axis.
Region
Type:String
GUI name:Sink region
Description:The region that is the sink. This key is ignored if the [Box] block is present.
Source
Type:Block
Description:Defines the heat source region (where the heat will be added).
AtomList
Type:Integer List
GUI name:Source region
Description:The atoms that are part of the source. This key is ignored if the [Box] block or [Region] key is present.
Box
Type:Block
Description:Part of the simulation box (in fractional cell coordinates) defining the heat source. If this block is specified, then by default, the whole box in each of the three dimensions is used, which usually does not make much sense. Normally, you will want to set the bounds along one of the axes. This block is mutually exclusive with the FirstAtom/LastAtom setting.
Amax
Type:Float
Default value:1.0
Description:Coordinate of the upper bound along the first axis.
Amin
Type:Float
Default value:0.0
Description:Coordinate of the lower bound along the first axis.
Bmax
Type:Float
Default value:1.0
Description:Coordinate of the upper bound along the second axis.
Bmin
Type:Float
Default value:0.0
Description:Coordinate of the lower bound along the second axis.
Cmax
Type:Float
Default value:1.0
Description:Coordinate of the upper bound along the third axis.
Cmin
Type:Float
Default value:0.0
Description:Coordinate of the lower bound along the third axis.
Region
Type:String
GUI name:Source region
Description:The region that is the source. This key is ignored if the [Box] block is present.
StartStep
Type:Integer
Default value:0
Description:Index of the MD step at which the heat exchange will start.
StopStep
Type:Integer
Description:Index of the MD step at which the heat exchange will stop.

One should be careful when choosing a value for the HeatingRate because a too large value may lead to pyrolysis of the heat source or to an abnormal termination when all the kinetic energy of the heat sink has been drained. The optimal value depends on the size of the system, its heat conductivity and the desired temperature gradient value. The thermal conductivity k can be calculated by dividing the heat flow rate W by the temperature gradient grad(T) and by the flow cross-section area S: \(k = W/(S \cdot grad(T))\). See the trajectory sampling section above on how to obtain the temperature profile from which the grad(T) can be computed.

References

[1]T. Ikeshoji and B. Hafskjold, Non-equilibrium molecular dynamics calculation of heat conduction in liquid and through liquid-gas interface Molecular Physics 81, 251-261 (1994)
[2]P. Wirnsberger, D. Frenkel, and C. Dellago, An enhanced version of the heat exchange algorithm with excellent energy conservation properties Journal of Chemical Physics 143, 124104 (2015)
[3]A. Vashisth, C. Ashraf, W. Zhang, C.E. Bakis, and A.C.T. van Duin, Accelerated ReaxFF simulations for describing the reactive cross-linking of polymers, J. Phys. Chem. A, 122, 6633-6642 (2018)
[4]M.J. Mees, G. Pourtois, E. C. Neyts, B. J. Thijsse, and André Stesmans, Uniform-acceptance force-bias Monte Carlo method with time scale to study solid-state diffusion, Phys. Rev. B 85, 134301 (2012)