# Molecular dynamics¶

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

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
MoleFraction float
NumAttempts integer
Rotate Yes/No
StartStep integer
StopStep integer
System string
Temperature float
TemperatureSigma float
Velocity float
VelocityDirection float_list
VelocitySigma float
End
ApplyForce
Force float_list
PerAtom Yes/No
Region string
End
ApplyVelocity
Components [X | Y | Z | XY | YZ | XZ | XYZ]
Region string
Velocity float_list
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
BinLog
BiasEnergy Yes/No
BoostFactor Yes/No
ConservedEnergy Yes/No
Density Yes/No
DipoleMoment Yes/No
Hypertime Yes/No
KineticEnergy Yes/No
MaxBiasEnergy Yes/No
MaxBoostFactor Yes/No
PotentialEnergy Yes/No
Pressure Yes/No
PressureTensor Yes/No
Step Yes/No
Temperature Yes/No
Time Yes/No
TotalEnergy Yes/No
Volume Yes/No
End
BondBoost
Chain
AtomNames string
MaxDistances float_list
MinDistances float_list
End
DistanceRestraint string
Moving Yes/No
NSteps integer
NumInstances integer
Units [Default | MD]
End
CRESTMTD
AddEnergy Yes/No
GaussianScaling
ScaleGaussians Yes/No
ScalingSlope float
End
Height float
NGaussiansMax integer
NSteps integer
Region string
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
ColVarBM
at1
region string
symbol string
End
at2
region string
symbol string
End
cutoff float
p integer
rmax float
rmin float
End
Frequency integer
MaxEvents integer
StartStep integer
StopStep integer
WaitSteps integer
End
CalcPressure Yes/No
Checkpoint
Frequency integer
WriteProperties Yes/No
End
CopyRestartTrajectory Yes/No
CosineShear
Acceleration float
Enabled Yes/No
FlowDirection float_list
ProfileAxis [X | Y | Z]
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
RandomVelocitiesMethod [Exact | Boltzmann | Gromacs]
Temperature float
Type [Zero | Random | FromFile | Input]
Values # Non-standard block. See details.
...
End
End
MovingRestraints
Change [Linear | Sine | Cosine]
Distance
Atom1 integer
Atom2 integer
EndValue float
StartValue float
End
Erf
ForceConstant float
MaxForce float
End
GaussianWell
Sigma float
WellDepth float
End
Harmonic
ForceConstant float
End
Hyperbolic
ForceConstant float
MaxForce float
End
Name string
Period float
RestraintType [None | Harmonic | Hyperbolic | Erf | GaussianWell]
StartStep integer
StopStep integer
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
ReactionBoost
BondBreakingRestraints
Erf
ForceConstant float
MaxForce float
End
GaussianWell
Sigma float
WellDepth float
End
Harmonic
ForceConstant float
End
Hyperbolic
ForceConstant float
MaxForce float
End
Taper
Enabled Yes/No
MaxDistance float
MinDistance float
End
Type [None | Harmonic | Hyperbolic | Erf | GaussianWell]
End
BondMakingRestraints
Erf
ForceConstant float
MaxForce float
End
GaussianWell
Sigma float
WellDepth float
End
Harmonic
ForceConstant float
End
Hyperbolic
ForceConstant float
MaxForce float
End
Type [None | Harmonic | Hyperbolic | Erf | GaussianWell]
End
BondedRestraints
Harmonic
ForceConstant float
End
Type [None | Harmonic]
End
Change [TargetCoordinate | Force | LogForce]
InitialFraction float
InterEquilibrationSteps integer
MinBondChange float
MinBondStrength float
NSteps integer
NonBondedRestraints
Exponential
Epsilon float
Sigma float
End
Type [None | Exponential]
End
PreEquilibrationSteps integer
RMSDRestraint
Erf
ForceConstant float
MaxForce float
End
GaussianWell
Sigma float
WellDepth float
End
Harmonic
ForceConstant float
End
Hyperbolic
ForceConstant float
MaxForce float
End
Type [None | Harmonic | Hyperbolic | Erf | GaussianWell]
End
Region string
TargetSystem string
Type [None | Pair | RMSD]
End
Reactor
ForceConstant float
MassScaled Yes/No
NSteps integer
Radius float
End
ReflectiveWall
Axis float_list
Region string
Threshold float
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
Shake
All string
ConvergeR2 float
ConvergeRV float
Iterations integer
ShakeInitialCoordinates Yes/No
End
Thermostat
BerendsenApply [Local | Global]
ChainLength integer
Duration integer_list
Region string
Tau float
Temperature float_list
Type [None | Berendsen | NHC]
End
TimeStep float
Trajectory
ExitConditionFreq integer
PrintFreq integer
SamplingFreq integer
TProfileGridPoints integer
WriteBonds Yes/No
WriteCharges Yes/No
WriteCoordinates Yes/No
WriteEngineGradients Yes/No
WriteMolecules Yes/No
WriteVelocities Yes/No
End
fbMC
Frequency integer
MassRoot float
MolecularMoves
Enabled Yes/No
RotationStepAngle float
TranslationStepLength float
End
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, and disabled by default for systems which are not translationally invariant (for example when frozen atoms are present).

`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. This is disabled by default for systems which are not translationally invariant (for example when frozen atoms are present).

## Constrained molecular dynamics¶

There are two types of constraint available for MD: frozen atoms and Shake/Rattle.
Atoms can be kept frozen using the `Atom`

, `AtomList`

and `FixedRegion`

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.

`FixedRegion`

- Type
String

- Recurring
True

- Description
Fix positions of all atoms in a region.

The Shake/Rattle constraints are implemented following H.C. Andersen, J Comp. Phys., 52 (1983) 24. The procedure consists of two steps: update of the coordinates to satisfy the constraints (the Shake part) and orthogonalization of atomic velocities to them (the Rattle part). The Shake part is implemented as an additional correcting force that brings each bond length exactly to its constrained value upon the subsequent velocity Verlet integration step. Application of the constraints is an iterative procedure and the convergence greatly depends on how intertwined the constraints are. For diatomic molecules, usually no more than two iterations are necessary, while tens of iterations may be required for rings (a triangle constraint is treated as a 3-membered ring).

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. Also the Shake constraints cannot be satisfied exactly during such simulations but the difference is usually small enough.

The Shake/Rattle constraints are specified using the Shake input block:

`MolecularDynamics`

`Shake`

- Type
Block

- Description
Parameters of the Shake/Rattle algorithm.

`All`

- Type
String

- Recurring
True

- GUI name
Constrain all

- Description
Constraint description in one the following formats: All [bondOrder] bonds at1 at2 [to distance] All triangles at1 at2 at3 The first option constrains all bonds between atoms at1 at2 to a certain length, while the second - bonds at1-at2 and at2-at3 and the angle between them. The [bondOrder] can be a number or a string such as single, double, triple or aromatic. If it’s omitted then all bonds between specified atoms will be constrained. Atom names are case-sensitive and they must be as they are in the Atoms block, or an asterisk ‘*’ denoting any atom. The distance, if present, must be in Angstrom. If it is omitted then the bond length from the initial geometry is used. Important: only the bonds present in the system at certain points of the simulation (at the start or right after adding/removing atoms) can be constrained, which means that the bonds may need to be specified in the System block. Warning: the triangles constraint should be used with care because each constrained bond or angle means removing one degree of freedom from the dynamics. When there are too many constraints (for example, “All triangles H C H” in methane) some of them may be linearly dependent, which will lead to an error in the temperature computation. Valid examples: All single bonds C C to 1.4 All bonds O H to 0.98 All bonds O H All bonds H * All triangles H * H

`ConvergeR2`

- Type
Float

- Default value
1e-08

- Description
Convergence criterion on the max squared difference, in atomic units.

`ConvergeRV`

- Type
Float

- Default value
1e-08

- Description
Convergence criterion on the orthogonality of the constraint and the relative atomic velocity, in atomic units.

`Iterations`

- Type
Integer

- Default value
100

- Description
Number of iterations.

`ShakeInitialCoordinates`

- Type
Bool

- Default value
Yes

- Description
Apply constraints before computing the first energy and gradients.

Note

The list of constrained bonds is built by matching the input against the bonds as specified in the System block. This means that bonds not specified in the System block will *not* be constrained.

## (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.

`RandomVelocitiesMethod`

- Type
Multiple Choice

- Default value
Exact

- Options
[Exact, Boltzmann, Gromacs]

- GUI name
Velocity randomization method

- Description
Specifies how are random velocities generated. Three methods are available. Exact: Velocities are scaled to exactly match set random velocities temperature. Boltzmann: Velocities are not scaled and sample Maxwell-Boltzmann distribution. However, the distribution is not corrected for constraints. Gromacs: Velocities are scaled to match set random velocities temperature, but removal of net momentum is performed only after the scaling. Resulting kinetic energy is lower based on how much net momentum the system had.

`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:

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

`ExitConditionFreq`

- Type
Integer

- GUI name
Exit condition frequency

- Description
Check the exit conditions every N steps. By default this is done every SamplingFreq steps.

`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 lattice axes will be written to the .rkf file. The temperature at a given profile point is calculated as the total temperature of all atoms inside the corresponding slice of the simulation box, time-averaged over all MD steps since the previous snapshot. 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.

`WriteCoordinates`

- Type
Bool

- Default value
Yes

- Description
Write atomic coordinates to the .rkf file.

`WriteEngineGradients`

- Type
Bool

- Default value
No

- Description
Write atomic gradients (negative of the atomic forces, as calculated by the engine) to the History section of ams.rkf.

`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

- GUI name
Calculate pressure

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

In rare cases, one may require some properties from every MD step, for example, the
pressure tensor for the Green-Kubo viscosity calculation. In this case one can
use the `BinLog`

(binary log) feature. When enabled, the selected MD variables will be stored
in the `BinLog`

section of `ams.rkf`

in the format similar to that of the `MDHistory`

section
with a few differences: the blocking factor is larger (1000 instead of 100)
and the pressure tensor is stored per component, for example, P_{xx} is stored in
`PressureTensor_xx`

, and so on. Note that none of the atomic properties (coordinates,
velocities) are included in the binlog.

`MolecularDynamics`

`BinLog`

- Type
Block

- Description
This block controls writing the BinLog section in ams.rkf, which contains the selected MD state scalars and tensors from every MD step. No per-atom data is written. If you need the per-atom data then you can set the sampling frequency to 1 and get the required data from the MDHistory section. The tensors are written per component, that is, the pressure tensor is written as six variables: PressureTensor_xx, PressureTensor_yy, etc.. To reduce the file size, all data is written in blocks.

`BiasEnergy`

- Type
Bool

- Default value
No

- Description
Write the CVDH bias energy.

`BoostFactor`

- Type
Bool

- Default value
No

- Description
Write the CVDH boost factor.

`ConservedEnergy`

- Type
Bool

- Default value
No

- Description
Write the conserved energy value.

`Density`

- Type
Bool

- Default value
No

- Description
Write the density.

`DipoleMoment`

- Type
Bool

- Default value
No

- Description
Write the dipole moment. Each component of the tensor is written in its own variable.

`Hypertime`

- Type
Bool

- Default value
No

- Description
Write the CVDH hypertime.

`KineticEnergy`

- Type
Bool

- Default value
No

- Description
Write the kinetic energy value.

`MaxBiasEnergy`

- Type
Bool

- Default value
No

- Description
Write the max CVDH bias energy.

`MaxBoostFactor`

- Type
Bool

- Default value
No

- Description
Write the max CVDH boost factor.

`PotentialEnergy`

- Type
Bool

- Default value
No

- Description
Write the potential energy value.

`Pressure`

- Type
Bool

- Default value
No

- Description
Write the pressure.

`PressureTensor`

- Type
Bool

- Default value
No

- Description
Write the pressure tensor in Voigt notation. Each component of the tensor is written in its own variable.

`Step`

- Type
Bool

- Default value
No

- Description
Write the step index during the simulation.

`Temperature`

- Type
Bool

- Default value
No

- Description
Write the temperature.

`Time`

- Type
Bool

- Default value
No

- Description
Write the simulation time (fs).

`TotalEnergy`

- Type
Bool

- Default value
No

- Description
Write the total energy value.

`Volume`

- Type
Bool

- Default value
No

- Description
Write the simulation cell volume, area or length, depending on the system periodicity.

## 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

The ReactionsDiscovery workflow uses lattice deformations to predict chemical reactions

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.

`MoleFraction`

- Type
Float

- Description
Defines a mixture to be deposited using one AddMolecules block per component. AMS will randomly alternate between any guns that have MoleFraction set. These need to all have the same settings for StartStep, StopStep and Frequency. Any additional AddMolecules blocks without MoleFraction will remain completely independent.

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

Since the 2022 release of the AMS driver, it is possible to generate molecule mixtures by using multiple guns that each have the `MoleFraction`

keyword set: AMS will then randomly alternate between any of these guns. The cooperating guns need to all have the same settings for `StartStep`

, `StopStep`

and `Frequency`

. (Any additional AddMolecules blocks without MoleFraction will remain completely independent.) The following example shows the generation of a 80:20 mixture of `molA`

and `molB`

:

```
MolecularDynamics
AddMolecules
System molA
MoleFraction 0.8
Frequency 2000
StartStep 1
End
AddMolecules
System molB
MoleFraction 0.2
Frequency 2000
StartStep 1
End
End
System anthracene
...
End
System benzene
...
End
```

Note that the exact number of inserted molecules of each type is not deterministic: running the same job twice may give you slightly different ratios. However, when running the job many times, the relative ratios averaged over the different runs should approach the fractions as requested in the input.

## 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 exponentially from 0 to 1
with the simulation time. The keyword `ScalingSlope`

in the input block `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.

By default, all atoms are included in the RMSD calculation. However, the selection can be restricted to specific regions
using the `Region`

keyword. This can be a single region or a combination thereof using region expressions (for example, `reg1+reg2`

for a union).

`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

`Region`

- Type
String

- Default value
*

- Description
Restrict the range of atoms for RMSD calculation to the specified region.

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

`ColVarBM`

- Type
Block

- Recurring
True

- GUI name
Collective Variable

- Description
Description of a bond-making collective variable (CV). A collective variable may consist of multiple ColVar blocks.

`at1`

- Type
Block

- Description
Specifies selection criteria for the first atom of a pair in the collective variable.

`region`

- Type
String

- Default value
*

- Description
Restrict the selection 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 pair. 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 selection criteria for the second atom of a pair in the collective variable.

`region`

- Type
String

- Default value
*

- Description
Restrict the selection 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 pair. The value is allowed to be the same as [at1], in which case pairs of 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 above this value are ignored when creating the initial atom-pair list for the CV. The list does not change during lifetime of the variable even if some bond orders rise above 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 much larger than the corresponding Rmin.

`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 the transition-state 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].

`MaxEvents`

- Type
Integer

- Default value
0

- Description
Max number of events to allow during dynamics. When this number is reached no new bias will be added for this input block.

`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 optimum distance in Angstrom, the first parameter and, optionally, the profile type and the second parameter. 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 atom index indicates a position of the corresponding atom in the AtomNames key. Currently recognized restraint profile types: Harmonic (default), Hyperbolic, Erf, GaussianWell. The first parameter is the force constant for the Harmonic, Hyperbolic, and Erf profiles, or well depth for GaussianWell. The second parameter is the asymptotic force value F(Inf) for Hyperbolic and Erf profiles, or the factor before x^2 in the exponent for GaussianWell. Note: the GaussianWell restraint should be used with the Moving flag.

`Moving`

- Type
Bool

- Default value
No

- GUI name
Move restraint

- Description
Move the restraints created with this BondBoost. The restraint value will start at the current coordinate’s value and will move towards the optimum during the restraint’s lifetime. The increment is calculated from the initial deviation and the [NSteps] parameter. This feature should be used with the GaussianWell restraint types.

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

`Units`

- Type
Multiple Choice

- Default value
Default

- Options
[Default, MD]

- GUI name
Restr. parameter units

- Description
Change energy, force and force constant units in the DistanceRestraint key from the default (atomic units) to those often used in the MD community (based on kcal/mol and Angstrom). Units for the optimum distances are not affected.

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.

### Reaction boost (Targeted MD)¶

The principle of targeted (or steered) MD is simple: during molecular dynamics, the system is pushed towards a certain known configuration in small steps. The steering can be in the form of holonomic constraints, then it is the standard Targeted MD (TMD). When restraints instead of constraints are used, then some authors call it steered MD 6.

The ReactionBoost feature uses restraints to push the system towards the target configuration. The restraints can be one of the following two types: atom-pair or (mass-weighted) RMSD. In both cases, the user can select a region to which the steering will be limited. AMS will then apply the restraints only to the atoms of the region(s). The RMSD type of restraints is suitable for very large systems or regions while the pair restraints can be used for small to moderately-sized ones. The atom-pair restraints can be further sub-divided into the following types, depending on the mutual bonding situation of the two atoms in the initial and the final state: bonded, non-bonded, bond-breaking and bond-making. The type is determined automatically at the time the steering is started based on the interatomic distances and the covalent radii of the atoms. A non-bonded restraint is a static repulsive potential, it is intended to keep atoms from forming a bond. The other three are moving: the restraint’s target distance changes between its initial and final values as the MD progresses. Setting the profile to `None`

disables the corresponding pair restraint type.

The following keys define parameters of the method:

`MolecularDynamics`

`ReactionBoost`

- Type
Block

- GUI name
Reaction Boost

- Description
Define a series of transitions between different states of the system. Each transition is defined by a TargetSystem and by a set of restraints that force the transition.

`BondBreakingRestraints`

- Type
Block

- Description
Define parameters for moving restraints that are added for pairs of atoms that become disconnected during the transition. It is intended to make sure the corresponding bonds get broken, although this may not always be required because forming other bonds will likely get these bonds broken.

`Erf`

- Type
Block

- Description
Define parameters for the Int(erf) potential V = alpha*(beta*x*erf(beta*x) + (exp(-(beta*x)**2) - 1)/sqrt(PI)). The alpha and beta parameters are computed from the user-defined ForceConstant and MaxForce.

`ForceConstant`

- Type
Float

- Default value
0.5

- Unit
Hartree/Bohr^2

- Description
The force constant (second derivative of the potential) at the optimum point.

`MaxForce`

- Type
Float

- Default value
0.05

- Unit
Hartree/Bohr

- Description
Asymptotic value of the force at the infinity.

`GaussianWell`

- Type
Block

- Description
Define parameters in the Gaussian well potential V=-WellDepth*exp(-Sigma*(r-r0)^2).

`Sigma`

- Type
Float

- Default value
1.0

- Unit
1/Bohr^2

- Description
Sigma parameter in the potential expression.

`WellDepth`

- Type
Float

- Default value
1.0

- Unit
Hartree

- Description
WellDepth parameter in the potential expression.

`Harmonic`

- Type
Block

- Description
Define parameters for the harmonic potential V=0.5*FC*(r-r0)^2.

`ForceConstant`

- Type
Float

- Default value
0.5

- Unit
Hartree/Bohr^2

- Description
The FC parameter of the harmonic potential.

`Hyperbolic`

- Type
Block

- Description
Define parameters for the hyperbolic potential V=alpha*(sqrt(1 + beta*x^2) - 1). The alpha and beta parameters are computed from the user-defined ForceConstant and MaxForce: beta=ForceConstant/MaxForce, alpha=MaxForce/beta

`ForceConstant`

- Type
Float

- Default value
0.5

- Unit
Hartree/Bohr^2

- Description
The force constant (second derivative of the potential) at the optimum point.

`MaxForce`

- Type
Float

- Default value
0.05

- Unit
Hartree/Bohr

- Description
Asymptotic value of the force at the infinity.

`Taper`

- Type
Block

- Description

`Enabled`

- Type
Bool

- Default value
No

- GUI name
Tapering

- Description
Enable tapering of the restraint potential and force between the given range of bond distances. A 7-th order tapering function on the actual (not target!) distance will be used. The MaxDistance must be greater than MinDistance.

`MaxDistance`

- Type
Float

- Default value
0.0

- Unit
Angstrom

- GUI name
End tapering at

- Description
Bond length at which the restraint potential and force decays to zero.

`MinDistance`

- Type
Float

- Default value
0.0

- Unit
Angstrom

- GUI name
Start tapering at

- Description
Bond length at which the restraint potential and force will start decaying to zero.

`Type`

- Type
Multiple Choice

- Default value
Erf

- Options
[None, Harmonic, Hyperbolic, Erf, GaussianWell]

- GUI name
Bond breaking restraints

- Description
Select type of the moving restraint profile. Harmonic: V=0.5*FC*(r-r0)^2 Hyperbolic: V=alpha*(sqrt(1 + beta*x^2) - 1) Erf: V = alpha*(beta*x*erf(beta*x) + (exp(-(beta*x)**2) - 1)/sqrt(PI)) GaussianWell: V=-WellDepth*exp(-Sigma*(r-r0)^2) Here beta=ForceConstant/MaxForce, alpha=MaxForce/beta. The force for Hyperbolic and Erf is bounded by a user-defined value, the latter converging to it faster than the former. The GaussianWell has a finite depth so it is suitable for cases when crossing a high reaction barrier is not desirable. Moving restraints are added for pairs of atoms that become disconnected during the transition. It is intended to make sure the corresponding bonds get broken, although this may not always be required because forming other bonds will likely get these bonds broken.

`BondMakingRestraints`

- Type
Block

- Description
Define parameters for moving restraints that are added for pairs of atoms that become connected during the transition. It is intended to make sure the bonds are created as required.

`Erf`

- Type
Block

- Description
Define parameters for the Int(erf) potential V = alpha*(beta*x*erf(beta*x) + (exp(-(beta*x)**2) - 1)/sqrt(PI)). The alpha and beta parameters are computed from the user-defined ForceConstant and MaxForce.

`ForceConstant`

- Type
Float

- Default value
0.5

- Unit
Hartree/Bohr^2

- Description
The force constant (second derivative of the potential) at the optimum point.

`MaxForce`

- Type
Float

- Default value
0.05

- Unit
Hartree/Bohr

- Description
Asymptotic value of the force at the infinity.

`GaussianWell`

- Type
Block

- Description
Define parameters in the Gaussian well potential V=-WellDepth*exp(-Sigma*(r-r0)^2).

`Sigma`

- Type
Float

- Default value
1.0

- Unit
1/Bohr^2

- Description
Sigma parameter in the potential expression.

`WellDepth`

- Type
Float

- Default value
1.0

- Unit
Hartree

- Description
WellDepth parameter in the potential expression.

`Harmonic`

- Type
Block

- Description
Define parameters for the harmonic potential V=0.5*FC*(r-r0)^2.

`ForceConstant`

- Type
Float

- Default value
0.5

- Unit
Hartree/Bohr^2

- Description
The FC parameter of the harmonic potential.

`Hyperbolic`

- Type
Block

- Description
Define parameters for the hyperbolic potential V=alpha*(sqrt(1 + beta*x^2) - 1). The alpha and beta parameters are computed from the user-defined ForceConstant and MaxForce: beta=ForceConstant/MaxForce, alpha=MaxForce/beta

`ForceConstant`

- Type
Float

- Default value
0.5

- Unit
Hartree/Bohr^2

- Description
The force constant (second derivative of the potential) at the optimum point.

`MaxForce`

- Type
Float

- Default value
0.05

- Unit
Hartree/Bohr

- Description
Asymptotic value of the force at the infinity.

`Type`

- Type
Multiple Choice

- Default value
Erf

- Options
[None, Harmonic, Hyperbolic, Erf, GaussianWell]

- GUI name
Bond making restraints

- Description
Select type of the moving restraint profile. Harmonic: V=0.5*FC*(r-r0)^2 Hyperbolic: V=alpha*(sqrt(1 + beta*x^2) - 1) Erf: V = alpha*(beta*x*erf(beta*x) + (exp(-(beta*x)**2) - 1)/sqrt(PI)) GaussianWell: V=-WellDepth*exp(-Sigma*(r-r0)^2) Here beta=ForceConstant/MaxForce, alpha=MaxForce/beta. The force for Hyperbolic and Erf is bounded by a user-defined value. The GaussianWell has a finite depth so it is suitable for cases when crossing a high reaction barrier is not desirable. Moving restraints are added for pairs of atoms that become connected during the transition. It is intended to make sure the bonds are created as required.

`BondedRestraints`

- Type
Block

- Description
Define parameters for bonded restraints. A bonded restraint is added for each pair of atoms that are bonded both in the current and in the final state. It is intended to make sure they remain bonded during simulation.

`Harmonic`

- Type
Block

- Description
Define parameters for the harmonic potential V=0.5*FC*(r-r0)^2.

`ForceConstant`

- Type
Float

- Default value
0.5

- Unit
Hartree/Bohr^2

- Description
The FC parameter of the harmonic potential.

`Type`

- Type
Multiple Choice

- Default value
None

- Options
[None, Harmonic]

- GUI name
Bonded restraints

- Description
Select type of the bonded restraints: Harmonic: V=0.5*FC*(r-r0)^2 A bonded restraint is added for each pair of atoms that are bonded both in the current and in the final state. It is intended to make sure they remain bonded during simulation.

`Change`

- Type
Multiple Choice

- Default value
TargetCoordinate

- Options
[TargetCoordinate, Force, LogForce]

- GUI name
Move type

- Description
Select what to change during dynamics. By default, once the restraints are switched on, AMS will change the restraint’s target coordinate towards its final value. If the [Force] or [LogForce] option is selected then the target coordinate is set to its final value immediately and instead the restraint force is gradually scaled from 0 to 1. The scaling is either linear (Force) or logarithmic (LogForce).

`InitialFraction`

- Type
Float

- Default value
0.0

- Description
Initial fraction of the boost variable. At the first boosting step, the restraint’s target value (or force or log(force)) is equal to InitialFraction + 1/NSteps.

`InterEquilibrationSteps`

- Type
Integer

- Default value
0

- Description
Number of equilibration steps after reaching a target before setting up restraints for the next one.

`MinBondChange`

- Type
Float

- Default value
1.0

- Unit
Bohr

- Description
Minimal change in the distance for an individual restraint to be considered bond-breaking/making vs bonded.

`MinBondStrength`

- Type
Float

- Default value
0.5

- Description
Minimum strength (usually ranges from 0 to 1) for a bond to be considered.

`NSteps`

- Type
Integer

- Default value
500

- GUI name
Steps per target

- Description
Number of steps per target the restraints should be active for.

`NonBondedRestraints`

- Type
Block

- Description
Define parameters for non-bonded restraints. A non-bonded restraint is added for each pair of atoms that are bonded neither in the current nor in the final state. It is intended to keep them from forming a bond unintentionally. They are represented by a repulsive potential

`Exponential`

- Type
Block

- Description
Define parameters for the repulsive potential V=Epsilon*exp(-Sigma*r).

`Epsilon`

- Type
Float

- Default value
1.0

- Unit
Hartree

- Description
Epsilon parameter in the repulsive potential expression.

`Sigma`

- Type
Float

- Default value
1.0

- Unit
1/Bohr

- Description
Sigma parameter in the repulsive potential expression.

`Type`

- Type
Multiple Choice

- Default value
None

- Options
[None, Exponential]

- GUI name
Non-bonded restraints

- Description
Select type of the non-bonded restraints: Exponential: V=Epsilon*exp(-Sigma*r) A non-bonded restraint is added for each pair of atoms that are bonded neither in the current nor in the final state. It is intended to keep them from forming a bond unintentionally. They are represented by a repulsive potential.

`PreEquilibrationSteps`

- Type
Integer

- Default value
0

- Description
Number of steps before enabling the first set of restraints.

`RMSDRestraint`

- Type
Block

- GUI name
RMSD restraint

- Description
Define a static restraint that pulls each atom to its position in the target system, but in contrast to the individual restraints, the force for this one depends on the total mass-weighted root-mean-squared distance (RMSD) between the two structures.

`Erf`

- Type
Block

- Description
Define parameters for the Int(erf) potential V = alpha*(beta*x*erf(beta*x) + (exp(-(beta*x)**2) - 1)/sqrt(PI)). The alpha and beta parameters are computed from the user-defined ForceConstant and MaxForce.

`ForceConstant`

- Type
Float

- Default value
0.5

- Unit
Hartree/Bohr^2

- Description
The force constant (second derivative of the potential) at the optimum point.

`MaxForce`

- Type
Float

- Default value
0.05

- Unit
Hartree/Bohr

- Description
Asymptotic value of the force at the infinity.

`GaussianWell`

- Type
Block

- Description
Define parameters in the Gaussian well potential V=-WellDepth*exp(-Sigma*(r-r0)^2).

`Sigma`

- Type
Float

- Default value
1.0

- Unit
1/Bohr^2

- Description
Sigma parameter in the potential expression.

`WellDepth`

- Type
Float

- Default value
1.0

- Unit
Hartree

- Description
WellDepth parameter in the potential expression.

`Harmonic`

- Type
Block

- Description
Define parameters for the harmonic potential V=0.5*FC*(r-r0)^2.

`ForceConstant`

- Type
Float

- Default value
0.5

- Unit
Hartree/Bohr^2

- Description
The FC parameter of the harmonic potential.

`Hyperbolic`

- Type
Block

- Description
Define parameters for the hyperbolic potential V=alpha*(sqrt(1 + beta*x^2) - 1). The alpha and beta parameters are computed from the user-defined ForceConstant and MaxForce: beta=ForceConstant/MaxForce, alpha=MaxForce/beta

`ForceConstant`

- Type
Float

- Default value
0.5

- Unit
Hartree/Bohr^2

- Description
The force constant (second derivative of the potential) at the optimum point.

`MaxForce`

- Type
Float

- Default value
0.05

- Unit
Hartree/Bohr

- Description
Asymptotic value of the force at the infinity.

`Type`

- Type
Multiple Choice

- Default value
None

- Options
[None, Harmonic, Hyperbolic, Erf, GaussianWell]

- GUI name
Type

- Description
Select type of the RMSD restraint profile: Harmonic: V=0.5*FC*(r-r0)^2 Hyperbolic: V=alpha*(sqrt(1 + beta*x^2) - 1) Erf: V = alpha*(beta*x*erf(beta*x) + (exp(-(beta*x)**2) - 1)/sqrt(PI), GaussianWell: V=-WellDepth*exp(-Sigma*(r-r0)^2) Here beta=ForceConstant/MaxForce, alpha=MaxForce/beta. The Harmonic profile can be problematic at large deviations as it may result in large forces. The force for Hyperbolic and Erf is bounded by a user-defined value. The GaussianWell has a finite depth so it is suitable for cases when crossing a high reaction barrier is not desirable.

`Region`

- Type
String

- Default value
*

- GUI name
Region

- Description
Region to which the restraints should be limited.

`TargetSystem`

- Type
String

- Recurring
True

- GUI name
Target system

- Description
The target system’s name for this transition. Multiple targets can be specified to request multiple transitions in one simulation. Note that only the lattice and the atomic coordinates of the target system are used and other properties (bonds, charge, etc.) are ignored. The target system’s lattice is used only to determine connections and it cannot be restrained.

`Type`

- Type
Multiple Choice

- Default value
None

- Options
[None, Pair, RMSD]

- GUI name
Restraint set type

- Description
Reaction Boost uses a series of transitions between different states of the system. Each transition is defined by a TargetSystem and by a set of restraints that force the transition. Select the type of the restraint set: -None: no Reaction Boost - Pair: use pair restraints - RMSD: use RMSD restraints. Pair restraints are defined per atom pair while the RMSD defines one collective restraint for all atoms and is thus suitable for very large systems. The pair restraints are further divided into four sub-types: bonding, non-bonding, bond-breaking and bond-making. The sub-type of restraints for each pair is determined automatically depending on whether the two atoms are bonded in the initial/final state. Parameters of the pair restraints are defined by [NonBondedRestraints], [BondedRestraints], [BondBreakingRestraints] and [BondMakingRestraints] blocks, while those of the RMSD restraint by the [RMSDRestraint] block.

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

`MolecularMoves`

- Type
Block

- Description
Move molecules as rigid bodies in addition to normal atomic moves.

`Enabled`

- Type
Bool

- Default value
No

- GUI name
Enable molecular moves

- Description
Enable moving molecules as rigid bodies based on net forces and torques. Ordinary per-atom displacements will then be based on residual atomic forces.

`RotationStepAngle`

- Type
Float

- Default value
0.1

- Unit
Radian

- Description
Maximum allowed angle of rotation of each molecule in one fbMC step.

`TranslationStepLength`

- Type
Float

- Default value
0.1

- Unit
Angstrom

- Description
Maximum allowed displacement of each molecule in each Cartesian coordinate in one fbMC step.

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

### Nanoreactor¶

See also

The ReactionsDiscovery workflow uses the nanoreactor to predict chemical reactions.

Nanoreactor is implemented following the article by L.-P. Wang et al. 5. The nanoreactor is a spherical region of space centered
at the origin and surrounded by an elastic wall (represented by a harmonic potential), its size changing over time. The idea of the
nanoreactor is to have two phases: reaction and relaxation. In the following, each nanoreactor phase is referred to as a *reactor*.
During the relaxation phase, the reactor is large and if the temperature is high enough the reactants are distributed more or less
evenly throughout its volume. In the reaction phase, the reactor shrinks and the atoms that now fall outside the reactor’s
sphere accelerate towards its center. During this phase, high-energy collisions may lead to reactions. After the reaction phase is over,
the relaxation phase begins again and the reactor is expanded. This cycle is repeated indefinitely.

Each phase is specified by its own `Reactor`

input block. There should be at least two of them but there can be more. The `NSteps`

key specifies the duration of the phase in MD steps. The radius of the reactor and the force constant of its wall are set by the
`Radius`

and `ForceConstant`

keys, respectively. To minimize additional strain on internal degrees of freedom of the affected molecules, the force
applied by the wall on an atom is proportional to the atom’s mass by default. This means that atoms situated at the same distance
from the center of the reactor will receive the same acceleration irrespective of their mass. This should be taken into account when
choosing a value for the ForceConstant. Alternatively, this feature can be disabled using the `MassScaled`

key.

`MolecularDynamics`

`Reactor`

- Type
Block

- Recurring
True

- Description
Define one phase of the nanoreactor. A reactor is a region of space surrounded by an elastic wall. Atoms inside the region are not affected. Atoms outside it will be pushed back with force depending on the [ForceConstant] and the [MassScaled] flag.

`ForceConstant`

- Type
Float

- GUI name
Reactor force constant

- Description
Force constant of the reactor wall in Hartree/Bohr^2 (or Hartree/Bohr^2/Dalton if [MassScaled] is true).

`MassScaled`

- Type
Bool

- Default value
Yes

- GUI name
Scale force by mass

- Description
If this flag is disabled the force on an atom outside of the reactor depends only on the atomic coordinates and the force constant. Otherwise, the force is also multiplied by the mass of the atom. This means that atoms at the same distance from the wall will receive the same accelerate due to the wall potential.

`NSteps`

- Type
Integer

- GUI name
Reactor lifetime

- Description
Number of steps for which the reactor will remain active until disabled. The next reactor will be activated immediately after this. After the last reactor is disabled the cycle will repeat.

`Radius`

- Type
Float

- Unit
Angstrom

- GUI name
Reactor radius

- Description
Radius of the reactor sphere.

## 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 Berendsen 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.

### NEMD for tribology: Applying velocities / forces¶

Non-equilibrium Molecular Dynamics (NEMD) is a commonly used tool to study phenomena like friction & wear in tribochemical systems.
In these MD simulations, an constant external force or velocity is applied to a specific subgroup of atoms.
For example, when studying sliding solid surfaces, the `ApplyForce`

feature can be used to exert a compressive normal force in the -Z direction on one of the solid slabs while the `ApplyVelocity`

block simultaneously slides the same group of atoms along the X or Y axis.

In the case of `ApplyVelocity`

, the specified group of atoms (defined as a Region) will move with the specified constant net `Velocity`

.
The imposed velocity is automatically excluded from thermostatting and system-wide momentum removal.
The `Components`

input key can be used to select just a subset of Cartesian components of the net velocity of the group to be controlled, so that other components can be left to evolve freely.
In the sliding example in the previous paragraph, one would use `Components=XY`

to let the Z motion of the sliding slab to freely follow the applied normal force.

In the case of `ApplyForce`

, the selected Region (group of atoms) will experience an additional specified constant `Force`

on each timestep of the simulation.
By default, the `Force`

parameter defines the total (net) force acting on the whole Region.
This can be changed by setting the `PerAtom`

option to apply the `Force`

to each atom in the Region separately.

`MolecularDynamics`

`ApplyVelocity`

- Type
Block

- Recurring
True

- Description
The ApplyVelocity keyword can be used to move an arbitrary group of atoms in the system with a constant net velocity

`Components`

- Type
Multiple Choice

- Default value
XY

- Options
[X, Y, Z, XY, YZ, XZ, XYZ]

- Description
Select which components of the Velocity vector are used to set the corresponding components of the net velocity of the specified set of atoms. Any other components of Velocity are ignored and the motion of the selected atoms in those directions is unaffected by ApplyVelocity.

`Region`

- Type
String

- Recurring
True

- Description
Applies the defined velocity to all atoms in this region.

`Velocity`

- Type
Float List

- Default value
[0.0, 0.0, 0.0]

- Unit
Angstrom/fs

- Recurring
False

- Description
The constant velocity that will be applied to the specified atoms.

`MolecularDynamics`

`ApplyForce`

- Type
Block

- Recurring
True

- Description
The ApplyForce keyword can be used to apply an arbitrary constant force to a certain (subgroups of) atoms in the system

`Force`

- Type
Float List

- Default value
[0.0, 0.0, 0.0]

- Unit
Hartree/Bohr

- Description
Defines the constant force vector

`PerAtom`

- Type
Bool

- Default value
No

- Description
If enabled, the Force vector is applied separately to every atom in the selected Region, so that the total net force on the Region equals the value of Force times the number of atoms in Region. This was the behavior of ApplyForce in AMS2022. By default, with PerAtom disabled, the Force vector defines the total net force on the Region, so the force applied to each atom equals the value of Force divided by the number of atoms in Region.

`Region`

- Type
String

- Recurring
True

- Description
Apply the constant force to all atoms in this region.

When the `ApplyVelocity`

block is used to slide a region, the measured net force experienced by that region is written to the MDHistory section on the trajectory.
This force is calculated directly from the gradients returned by an Engine and thus excludes the effect of any `ApplyForce`

blocks.
The ratio of the tangential (sliding) and normal (compression) component of the mean force is thus equal to the friction coefficient.
Three variables will be written by any `ApplyVelocity`

block:

*EngXFrcRgnA*: instantaneous engine force on region*A*, Cartesian component*X**MeanEngXFrcRgnA*: average of*EngXFrcRgnA*over all timesteps since the previous trajectory frame*StdevEngXFrcRgnA*: standard deviation of*EngXFrcRgnA*over all timesteps since the previous trajectory frame

Note

Only the first eight characters of the region name are used to construct the *EngXFrcRgn* variables. When using multiple `ApplyVelocity`

instances with different regions, make sure their names are not ambiguous after shortening to eight characters.

### NEMD viscosity via cosine shear perturbation¶

In order to calculate the viscosity using NEMD, apply a periodic perturbation to the system using the `CosineShear`

input block.
This will add an amount of acceleration to the atoms in the system according to a cosine function running along the specified Cartesian `ProfileAxis`

.
The acceleration applied to each atom will have the direction of `FlowDirection`

and magnitude equal to \(a_x(z) = A \cos(\frac{2 \pi z}{L_z})\), where *A* is the amplitude defined by the `Acceleration`

input key and *L_z* is the length of the periodic box in the Z axis (assuming `ProfileAxis=Z`

and `FlowDirection="1 0 0"`

.
It is typically the most efficient to use a box roughly three times taller in the `ProfileAxis`

direction than in the other two axes.
AMS currently requires the `ProfileAxis`

to be aligned with one of the lattice vectors and perpendicular to the remaining lattice vectors.

Warning

`CosineShear`

simulations are typically performed in the NVT ensemble.
While using a `Barostat`

is theoretically possible, such a combination has not been physically validated by SCM.

`MolecularDynamics`

`CosineShear`

- Type
Block

- Description
Apply an external acceleration to all atoms of a fluid using a periodic (cosine) function along a selected coordinate axis. This induces a periodic shear flow profile which can be used to determine the viscosity.

`Acceleration`

- Type
Float

- Default value
5e-06

- Unit
Angstrom/fs^2

- Description
Amplitude of the applied cosine shear acceleration profile. The default value should be a rough first guess for water and it needs to be adjusted by experimentation for other systems.

`Enabled`

- Type
Bool

- Default value
No

- GUI name
Enable cosine shear

- Description
Apply a cosine shear acceleration profile for a NEMD calculation of viscosity.

`FlowDirection`

- Type
Float List

- Default value
[1.0, 0.0, 0.0]

- Description
The direction in which to apply the shear acceleration, in Cartesian coordinates. The magnitude of this vector is ignored (AMS will normalize it internally). FlowDirection has to be perpendicular to ProfileAxis.

`ProfileAxis`

- Type
Multiple Choice

- Default value
Z

- Options
[X, Y, Z]

- Description
The Cartesian coordinate axis along which the cosine wave runs

This generates a velocity profile in the steady state that looks like: \(v_x(z) = V \cos(\frac{2 \pi z}{L_z})\).
After the system reaches a steady state, the velocity amplitude *V* will be approximately constant, although typically relatively noisy.
Averaging *V* over a long enough trajectory yields a reasonable estimate of the steady-state flow velocity, which can then be used to calculate the viscosity as \(\eta = \frac{A \rho L_z^2}{4 \pi^2}\).
The measured value of *V* averaged over every `Trajectory SamplingFreq`

timesteps is written to the trajectory as *MeanCosineShearVelocity* in the units of Bohr/fs.

The calculated viscosity value will be different from the equilibrium viscosity of the studied liquid due to shear thinning and heating caused by the applied acceleration *A*.
It is thus necessary to measure the viscosity at several different values of *A* and extrapolate the obtained values to the zero-shear limit \(A = 0\).
To do this, it is useful to calculate the “shear parameter” \(S = \frac{A^2 m \rho L_z^2}{8 \pi^2 N_f k_B}\), where *m* is the total mass of the system and *N_f* is the number of degrees of freedom.
The measured viscosity \(\eta\) is then approximately a linear function of *S*, enabling a straightforward extrapolation to zero *S*.
Typical values of the shear parameter *S* in SI units are on the order of 10⁹ - 10¹⁰ kg m s⁻² K.

## Applying reflective walls¶

Note

This feature is still experimental. It may contain bugs and some input keys are still subject to change in future releases of AMS. If you are interested in trying it out, monitor the bug fix changelog for updates and feel free to get in touch with us at support@scm.com if you run into any issues.

The ReflectiveWall keyword can be used to define impenetrable planes in space, causing molecules to bounce off by an elastic collision. Therefore, the main usage of reflective walls is to make sure that molecules will stay within a certain area in space.

The reflective wall is defined by a normal vector (`Axis`

) and a `Threshold`

.
The plane of the wall is perpendicular to `Axis`

and passes through a point given by multiplying `Axis`

by `Threshold`

.
Any atom moving in the direction of `Axis`

will be reflected back, but the opposite does not apply.
Particles moving in the direction opposite to `Axis`

will not be affected (the wall works like a one-way mirror).
The `Axis`

vector does not need to be normalized to unity, AMS will normalize it internally and adjust `Threshold`

accordingly.

`MolecularDynamics`

`ReflectiveWall`

- Type
Block

- Recurring
True

- Description
Apply a reflective wall in space

`Axis`

- Type
Float List

- Unit
Angstrom

- Description
Defines the normal vector perpendicular to the plane of the reflective wall. Any particle moving in this direction will be reflected back.

`Region`

- Type
String

- Recurring
True

- Description
Apply the reflective wall to all atoms in this region.

`Threshold`

- Type
Float

- Unit
Angstrom

- Description
Defines the threshold value determining the position of the reflective wall. If the dot product of a position of a particle with Axis exceeds Threshold, the particle will be reflected. This means that the plane of the wall passes through a point given by Axis times Threshold.

## Exit Conditions¶

It is possible to define exit conditions to end a molecular dynamics calculations before all timesteps have been performed.
If an exit condition is met, the variable `History%ExitConditionMsg`

is populated with a message describing why the molecular dynamics task was stopped.
When an exit condition is met, its status is considered to be `NORMAL TERMINATION`

.

```
ExitCondition
AtomsTooClose
MinimumDistance float
PairCalculation [NeighborList | DistanceMatrix]
End
EngineEnergyUncertainty
MaxUncertainty float
Normalization float
End
EngineGradientsUncertainty
MaxUncertainty float
End
Type [AtomsTooClose | EngineEnergyUncertainty | EngineGradientsUncertainty]
End
```

`ExitCondition`

- Type
Block

- Recurring
True

- Description
If any of the specified exitconditions are met, the AMS driver will exit cleanly.

`AtomsTooClose`

- Type
Block

- Description
If any pair of atoms is closer than the specified minimum value, the program will exit cleanly.

`MinimumDistance`

- Type
Float

- Default value
0.7

- Unit
Angstrom

- Description
Two atoms closer than this threshold value are considered too close.

`PairCalculation`

- Type
Multiple Choice

- Default value
NeighborList

- Options
[NeighborList, DistanceMatrix]

- Description
Two atoms closer than this threshold value are considered too close.

`EngineEnergyUncertainty`

- Type
Block

- Description
If the engine reports an uncertainty that is too high, the program will exit cleanly.

`MaxUncertainty`

- Type
Float

- Default value
0.001

- Unit
Hartree

- Description
Threshold for Engine Energy Uncertainty divided by Normalization (by default the number of atoms)

`Normalization`

- Type
Float

- Value Range
value >= 0.0

- Description
Divide the reported Engine Energy Uncertainty by this normalization. Will divide by the number of atoms if unset.

`EngineGradientsUncertainty`

- Type
Block

- Description
If the engine reports an uncertainty in the magnitude of the nuclear gradient of any atom that is too high, the program will exit cleanly.

`MaxUncertainty`

- Type
Float

- Default value
0.01580221

- Unit
Hartree/Angstrom

- Description
Threshold for Engine Gradients Uncertainty.

`Type`

- Type
Multiple Choice

- Default value
AtomsTooClose

- Options
[AtomsTooClose, EngineEnergyUncertainty, EngineGradientsUncertainty]

- Description
The type of exitcondition specified

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)- 5
L.-P. Wang, A. Titov, R. McGibbon, F. Liu, V.S. Pande and T.J. Martínez,

*Discovering chemistry with an ab initio nanoreactor*, Nature Chemistry 6, 1044 (2014)- 6
He Huang, E. Ozkirimli, C.B. Post,

*A Comparison of Three Perturbation Molecular Dynamics Methods for Modeling Conformational Transitions*, J Chem Theory Comput. 5 (2009) 1301–1314