Molecular dynamics
******************

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

See also:

  * Examples

  * AMS GUI Tutorial

  * MD trajectory analysis tool

To perform an 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
      CheckTemperature
         Enabled Yes/No
         MaxTemperature float
      End
      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
         TargetDensity float
         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
         Thickness float
         Threshold float
      End
      Remap
         Range [Auto | ZeroToOne | PlusMinusHalf | AroundCenter]
         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
         ExcludeRegion string
         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
         EngineResultsFreq integer
         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.

      "ExcludeRegion"
         Type:
            String

         Default value:
         Description:
            Exclude all atoms covered by this region expression from
            all constraints. If both atoms of a bond belong to
            ExcludeRegion, the bond will be kept unconstrained,
            regardless of Shake%All or engine-specific constraints.

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

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" (or a System block
containing the most recent geometry to start from) 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's behavior.

      "BerendsenApply"
         Type:
            Multiple Choice

         Default value:
            Global

         Options:
            [Local, Global]

         GUI name:
            Apply Berendsen

         Description:
            Select how to apply the scaling correction for the
            Berendsen thermostat:                         - per-atom-
            velocity (Local)                         - on the
            molecular system as a whole (Global).

      "ChainLength"
         Type:
            Integer

         Default value:
            10

         GUI name:
            NHC chain length

         Description:
            Number of individual thermostats forming the NHC
            thermostat

      "Duration"
         Type:
            Integer List

         GUI name:
            Duration(s)

         Description:
            Specifies how many steps should a transition from a
            particular temperature to the next one in sequence take.

      "Region"
         Type:
            String

         Default value:
            *

         Description:
            The identifier of the region to thermostat. The default
            *** applies the thermostat to the entire system. The value
            can by a plain region name, or a region expression, e.g.
            **-myregion* to thermostat all atoms that are not in
            myregion, or *regionA+regionB* to thermostat the union of
            the 'regionA' and 'regionB'. Note that if multiple
            thermostats are used, their regions may not overlap.

      "Tau"
         Type:
            Float

         Unit:
            Femtoseconds

         GUI name:
            Damping constant

         Description:
            The time constant of the thermostat.

      "Temperature"
         Type:
            Float List

         Unit:
            Kelvin

         GUI name:
            Temperature(s)

         Description:
            The target temperature of the thermostat.
            You can specify multiple temperatures (separated by
            spaces). In that case the Duration field specifies how
            many steps to use for the transition from one T to the
            next T (using a linear ramp). For NHC thermostat, the
            temperature may not be zero.

      "Type"
         Type:
            Multiple Choice

         Default value:
            None

         Options:
            [None, Berendsen, NHC]

         GUI name:
            Thermostat

         Description:
            Selects the type of the thermostat.

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

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

Berendsen
   The Berendsen friction-like isobaric ensemble method rescales the
   system in each step to drive the pressure towards a target value.
   Similarly to the "Berendsen" thermostat, the relaxation is
   exponential with a time constant "Tau". Similar considerations for
   the choice of "Tau" apply as in the case of the thermostat, but the
   value of "Tau" for the barostat is usually at least several times
   higher than the corresponding "Tau" used for the thermostat. This
   barostat does not have any conserved quantity.

MTK
   This enables the Martyna-Tobias-Klein extended Lagrangian barostat,
   which generates a true isobaric ensemble by integrating the cell
   parameters as additional degrees of freedom. This barostat is
   derived from the Andersen-Hoover isotropic barostat and the
   Parrinello-Rahman-Hoover anisotropic barostat. Like the "NHC"
   thermostat, it exhibits oscillatory relaxation unsuitable for
   systems far from equilibrium. This barostat must always be combined
   with a "NHC" thermostat. One instance of such thermostat coupled to
   the atoms as usual, while a second instance is created internally
   and coupled to the cell degrees of freedom.

"MolecularDynamics"
   "Barostat"
      Type:
         Block

      Description:
         This block allows to specify the use of a barostat during the
         simulation.

      "BulkModulus"
         Type:
            Float

         Default value:
            2200000000.0

         Unit:
            Pascal

         Description:
            An estimate of the bulk modulus (inverse compressibility)
            of the system for the Berendsen barostat.
            This is only used to make Tau correspond to the true
            observed relaxation time constant. Values are commonly on
            the order of 10-100 GPa (1e10 to 1e11) for solids and 1
            GPa (1e9) for liquids (2.2e9 for water). Use 1e9 to match
            the behavior of standalone ReaxFF.

      "ConstantVolume"
         Type:
            Bool

         Default value:
            No

         Description:
            Keep the volume constant while allowing the box shape to
            change.                         This is currently
            supported only by the MTK barostat.

      "Duration"
         Type:
            Integer List

         Description:
            Specifies how many steps should a transition from a
            particular pressure to the next one in sequence take.

      "Equal"
         Type:
            Multiple Choice

         Default value:
            None

         Options:
            [None, XYZ, XY, YZ, XZ]

         Description:
            Enforce equal scaling of the selected set of dimensions.
            They will be barostatted as one dimension according to the
            average pressure over the components.

      "Pressure"
         Type:
            Float List

         Unit:
            Pascal

         Description:
            Specifies the target pressure.
            You can specify multiple pressures (separated by spaces).
            In that case the Duration field specifies how many steps
            to use for the transition from one p to the next p (using
            a linear ramp).

      "Scale"
         Type:
            Multiple Choice

         Default value:
            XYZ

         Options:
            [XYZ, Shape, X, Y, Z, XY, YZ, XZ]

         Description:
            Dimensions that should be scaled by the barostat to
            maintain pressure. Selecting Shape means that all three
            dimensions and also all the cell angles are allowed to
            change.

      "Tau"
         Type:
            Float

         Unit:
            Femtoseconds

         GUI name:
            Damping constant

         Description:
            Specifies the time constant of the barostat.

      "Type"
         Type:
            Multiple Choice

         Default value:
            None

         Options:
            [None, Berendsen, MTK]

         GUI name:
            Barostat

         Description:
            Selects the type of the barostat.


Temperature and pressure regimes
--------------------------------

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

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

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

1. Increase linearly from 0 to 300 K over 100 steps.

2. Stay constant at 300 K for 200 steps.

3. Increase linearly from 300 to 500 K over 100 steps.

4. Stay constant at 500 K for 200 steps.

5. Decrease linearly from 500 to 300 K over 100 steps.

6. Stay constant at 300 K for the rest of the simulation.


Trajectory sampling and output
==============================

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

AMS writes the trajectory to the "History" and "MDHistory" sections of
"ams.rkf", according to the settings in the "Trajectory" block. A
snapshot of the system and various MD variables is stored every
"SamplingFreq" timesteps. By default, this frequency is also used to
print basic thermodynamic parameters of the simulation to the output
and log file. Set "PrintFreq" to override this.

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

* "WriteBonds" is necessary for reaction network analysis
  (ChemTraYzer). Disabling "WriteBonds" also makes AMSmovie show only
  guessed bonds instead of those calculated by the engine.

* "WriteMolecules" is required by the Molecule Fractions panel in
  AMSmovie.

* "WriteVelocities" is required to calculate the velocity
  autocorrelation functions needed for diffusivity and IR spectra.

The trajectory itself contains only the data needed for subsequent
analysis of the dynamics of the system. However, much more data is
usually generated on every integration step. This includes, for
example, the internal data used by an engine when evaluating the
energies and forces. This information is normally discarded after most
steps, because it is often very large. However, every
"EngineResultsFreq" steps this engine-specific data is written to a
new engine result file called "MDStep*.rkf". These files can be
inspected afterwards to extract various engine-specific details about
the system, such as the orbitals for QM engines.

For some engines such as ForceField or ReaxFF, the data on engine
results files is rarely needed for subsequent analysis. In that case,
setting "EngineResultsFreq" to zero can be used to disable writing
engine results files and save disk space.

To enable restarting interrupted MD simulations or extending them to a
larger number of steps, a "Checkpoint" containing the complete
internal state of the MD driver is stored every "Frequency" steps. An
interrupted simulation can then be restarted from this checkpoint
using the "Restart" keyword.

Note:

  In previous AMS versions, users sometimes increased
  "Checkpoint%Frequency" to suppress unwanted engine results files.
  This is no longer necessary because "EngineResultsFreq" now defaults
  to 0. The "Checkpoint%Frequency" should be set based only on how
  often the MD driver should save its state for restarts. The default
  should work well for fast engines, but it might need to be adjusted
  downwards for expensive (DFT) engines to ensure frequent
  checkpointing to avoid losing many steps since the last checkpoint
  in case the simulation gets interrupted.

"MolecularDynamics"
   "Trajectory"
      Type:
         Block

      Description:
         Sets the frequency for printing to stdout and storing the
         molecular configuration on the .rkf file.

      "EngineResultsFreq"
         Type:
            Integer

         Default value:
            0

         GUI name:
            Engine results frequency

         Description:
            Write MDStep*.rkf files with engine-specific results once
            every N steps. Setting this to zero disables writing
            engine results files except for one file at the end of the
            simulation. If this is set to a negative value, the value
            of Checkpoint%Frequency will be used instead.

      "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 molecular geometry (and possibly other
            properties) to the ams.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 to the Molecule and MDResults sections
            on ams.rkf once every N steps. Setting this to zero
            disables checkpointing during the simulation, but one
            checkpoint at the very end of the simulation will always
            be written.                         NOTE: This no longer
            controls the frequency of writing MDStep*.rkf files, set
            Trajectory%EngineResultsFreq instead.

      "WriteProperties"
         Type:
            Bool

         Default value:
            No

         Description:
            Honor top-level Properties input block when writing engine
            results files.                         DEPRECATED: This
            input option will be removed in AMS2026 and will be
            considered always enabled. If you rely on it being
            disabled, please contact SCM support.

   "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 CVHD bias energy at every time step.

      "BoostFactor"
         Type:
            Bool

         Default value:
            No

         Description:
            Write the CVHD boost factor at every time step.

      "ConservedEnergy"
         Type:
            Bool

         Default value:
            No

         Description:
            Write the conserved energy value at every time step.

      "Density"
         Type:
            Bool

         Default value:
            No

         Description:
            Write the density at every time step.

      "DipoleMoment"
         Type:
            Bool

         Default value:
            No

         GUI name:
            Dipole moment binlog

         Description:
            Write the dipole moment at every time step. Each component
            of the tensor is written in its own variable.

      "Hypertime"
         Type:
            Bool

         Default value:
            No

         Description:
            Write the CVHD hypertime at every time step.

      "KineticEnergy"
         Type:
            Bool

         Default value:
            No

         Description:
            Write the kinetic energy value at every time step.

      "MaxBiasEnergy"
         Type:
            Bool

         Default value:
            No

         Description:
            Write the max CVHD bias energy at every time step.

      "MaxBoostFactor"
         Type:
            Bool

         Default value:
            No

         Description:
            Write the max CVHD boost factor at every time step.

      "PotentialEnergy"
         Type:
            Bool

         Default value:
            No

         Description:
            Write the potential energy value at every time step.

      "Pressure"
         Type:
            Bool

         Default value:
            No

         Description:
            Write the pressure at every time step.

      "PressureTensor"
         Type:
            Bool

         Default value:
            No

         GUI name:
            Pressure tensor binlog

         Description:
            Write the pressure tensor in Voigt notation at every time
            step. 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 at every time
            step.

      "Temperature"
         Type:
            Bool

         Default value:
            No

         Description:
            Write the temperature at every time step.

      "Time"
         Type:
            Bool

         Default value:
            No

         GUI name:
            Time (fs) binlog

         Description:
            Write the simulation time (fs) at every time step.

      "TotalEnergy"
         Type:
            Bool

         Default value:
            No

         Description:
            Write the total energy value at every time step.

      "Volume"
         Type:
            Bool

         Default value:
            No

         Description:
            Write the simulation cell volume, area or length,
            depending on the system periodicity, at every time step.

Every time a snapshot of the system is stored in the trajectory file,
atomic coordinates are remapped (wrapped) across periodic boundaries
into the unit cell. This can be disabled by setting "Remap Type=None".
By default, the MD driver tries to detect which of the two common PBC
conventions the system uses, based on the atomic positions at the
start of the simulation. This can be customized using the "Remap
Range" setting.

"MolecularDynamics"
   "Remap"
      Type:
         Block

      Description:
         Control periodic remapping (backtranslation) of atoms into
         the PBC box.

      "Range"
         Type:
            Multiple Choice

         Default value:
            Auto

         Options:
            [Auto, ZeroToOne, PlusMinusHalf, AroundCenter]

         Description:
            Select the range of fractional coordinates to which atoms
            are remapped.                         Auto: Use
            PlusMinusHalf if the geometrical center of the starting
            positions of all atoms lies between -0.25 and 0.25 in
            fractional coordinates in all periodic directions,
            ZeroToOne otherwise.
            PlusMinusHalf: Remap atoms to lie between -0.5 and 0.5 in
            fractional coordinates.                         ZeroToOne:
            Remap atoms to lie between 0 and 1 in fractional
            coordinates.                         AroundCenter: Remap
            atoms to lie within 0.5 in fractional coordinates around
            the geometrical center of the starting positions of all
            atoms.

      "Type"
         Type:
            Multiple Choice

         Default value:
            Atoms

         Options:
            [None, Atoms]

         Description:
            Select the method used to remap atoms into the unit cell.
            None: Disable remapping completely.
            Atoms: Remap any atoms that leave the unit cell.


Lattice deformations (volume regimes)
=====================================

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

See also:

  * Example demonstrating various lattice deformations

  * 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 seven
mutually exclusive input keys. These belong to two main 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.

TargetDensity
   This input key serves as a shorthand for setting TargetLattice to
   the current lattice at "StartStep" isotropically scaled by a factor
   needed to reach a particular target density at "StopStep".

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

      "TargetDensity"
         Type:
            Float

         Unit:
            g/cm3

         Description:
            Target density of the system to be achieved by StopStep by
            isotropically rescaling the lattice.
            This can only be used for 3D-periodic systems. Note that
            the selected Type determines how the lengths of the
            lattice vectors will evolve (linearly or as a
            sine/cosine), the evolution of the volume or density is
            thus necessarily non-linear.

      "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
====================

AMS provides several methods for accelerating the simulation of rare
events, such as chemical reactions or conformational changes that
would be impractical to observe with conventional molecular dynamics.

See also:

  For path sampling simulations using the Replica Exchange Transition
  Interface Sampling (RETIS) method, see the InfRETIS manual.


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 (version 2.5.0). 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 target distance
            in Angstrom, the first parameter and, optionally, the
            profile type and the second parameter.
            For periodic systems, restraints follow 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 profile name can optionally be prefixed by
            "OneSided-", thus making this particular restraint one-
            sided. A one-sided restraint will only push or pull the
            distance towards its target value depending on the sign of
            the difference between the initial and the target
            distance. For example, if the initial distance is larger
            than the target one then the restraint will only push the
            value down. This is especially useful with moving
            restraints because then the restraint will effectively
            disappear as soon as the system crosses a barrier along
            the reaction path. By default, the restraints are two-
            sided, which means they will try to keep the distance
            between the two specified atoms near the (possibly moving)
            target value. For moving restraints, this means that after
            crossing a reaction barrier the system will slide slowly
            towards products held back by the restraints.
            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*(1-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.

See also: Example demonstrating a NEMD friction simulation

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

See also: Example demonstrating a NEMD viscosity simulation

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.

The optional "Thickness" key defines how far beyond the reflective
face particles are still detected and reflected. This enables the
placement of two walls back-to-back (separated by two times
"Thickness" so that their zones of operation do not overlap) to form a
wall that reflects in both directions. The default value of
"Thickness" should be appropriate for most systems, but may need to be
increased for very high-speed particles (such as light atoms at high
temperatures) to prevent them from tunneling through the wall in a
single timestep.

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

      "Thickness"
         Type:
            Float

         Default value:
            0.1

         Unit:
            Angstrom

         Description:
            Any particle not deeper than Thickness beyond the
            reflective face of the wall will be reflected. Particles
            deeper than this will be ignored by the wall. This setting
            allows the placement of two walls back-to-back to create a
            two-sided wall that reflects in both directions (offset
            the position of each individual wall by Thickness so that
            their areas of influence do not overlap). Systems with
            very high speed particles (light atoms at high
            temperatures) need a higher Thickness to stop atoms from
            tunneling through the wall in a single timestep.

      "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

[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
