ReaxFF input
************

This section describes the input keywords to the ReaxFF AMS engine.

See also:

  * AMS driver's tasks and properties


Force field specification
=========================

The only input key required by the engine is "ForceField", used to
select the force field file. Force fields included in the Amsterdam
Modeling Suite can be easily accessed using their file name, such as
"CHO.ff".

"ForceField"
   Type:
      String

   Description:
      Path to the force field parameter file. The path can be absolute
      or relative. Relative paths starting with ./ are considered
      relative to the directory in which the calculation is started,
      otherwise they are considered relative to
      $AMSRESOURCES/ForceFields/ReaxFF.


Recommended lattice convention
==============================

The ReaxFF engine supports molecular (free boundary), 1D-, 2D-, and
3D-periodic systems. Non-orthorhombic lattices are supported in an
arbitrary orientation. However, the engine is slightly more
computationally efficient when the cell is oriented according to the
convention used in standalone ReaxFF, i.e. lattice vector "c" aligned
with the "z" axis and vector "b" in the "yz" plane (zero "x"
component). The "Lattice" block in the system definition then looks
like this:

   System
      Lattice
        xx xy xz
         0 yy yz
         0  0 zz
      End
   End


Smoothened potential energy surface
===================================

The keywords below can be used to enable the **tapered bond orders**
and/or **improved torsion angle potentials**. Although the original
ReaxFF torsion potential is the default to preserve backward
compatibility, the corrected potentials eliminate energy
discontinuities and work well with existing force fields.

Using the tapered bond order (the "TaperBO" key) does not change the
potential form at the chemically relevant distances so it can be used
with any force-field. It may improve the energy conservation during MD
and make geometry optimizations with ReaxFF converge to much tighter
criteria.  The discontinuity and the correction for it are described
in detail in J. Phys. Chem. Lett. 10 (2019) 7215.

"TaperBO"
   Type:
      Bool

   Default value:
      No

   GUI name:
      Taper bond orders

   Description:
      Use tapered bond orders by Furman & Wales (DOI:
      10.1021/acs.jpclett.9b02810).

The discontinuity at small bond orders in the expression for torsion
angles and conjugation contributions can alternatively be corrected
for using the "Torsions 2013" correction. The corresponding terms are
given by f_10 (eq. 10b) and f_12 (eq. 11b) in the original ReaxFF
paper. The new expression for each term in f_10 is \left(1 - e^{-2
\lambda_{23} \text{BO}^2} \right) and in f_12 the new expression is
\sin(\frac{\pi}{3} \text{BO})^4. The new expressions ensure correct
asymptotic behavior for the \frac{dE}{d\text{BO}} for BO \rightarrow
0.

Another discontinuity in the torsion angle term is when one or both
valence angles approach 180 degrees. It is described in detail in J.
Chem. Phys. 153 (2020) 021102, and can be enabled with "FurmanTorsions
Yes".

See also: Geometry optimization issues

"Torsions"
   Type:
      Multiple Choice

   Default value:
      Original

   Options:
      [Original, 2013]

   Description:
      Version of the torsion potential expression.

"FurmanTorsions"
   Type:
      Bool

   Default value:
      No

   Description:
      Use (sin(Theta_ijk)*sin(Theta_jkl))^3 instead of
      sin(Theta_ijk)*sin(Theta_jkl) in the torsion energy term to
      remove discontinuity in the corresponding force.


Bond order and distance cutoffs
===============================

"BondDistanceCutoff"
   Type:
      Float

   Default value:
      5.0

   Unit:
      Angstrom

   Description:
      Maximum distance between two atoms to be considered when
      searching for possible bonds.

"BondOrderCutoff"
   Type:
      Float

   Default value:
      0.001

   Description:
      Minimum bond order required for a bond to be considered during
      the evaluation of the 3- and 4-body potentials.

"StrongBondCutoff"
   Type:
      Float

   Default value:
      0.3

   Description:
      Minimum bond order required for a bond to be returned to the
      driver for bonding analysis and molecule detection. Bonds below
      this threshold are only used to evaluate the potential and not
      written to result files.


Non-reactive mode
=================

The engine can also be switched to a special non-reactive mode useful
mainly for initial preparation of molecular dynamics simulations. This
mode greatly reduces the occurrence of unwanted reactions when
starting from an unrelaxed geometry. In these situations, we recommend
running a short simulation with the "NonReactive" key to relieve the
initial conformational strain and then restarting the MD run without
this key.

Note that if you want to resume or extend an interrupted "NonReactive"
run, it is recommended to also use the "EngineRestart" AMS key to
supply the last ReaxFF ".rkf" file from the previous run. This enables
the engine to load the bonding topology used during the previous run
and ensure that the simulation is seamlessly restarted. If the
"EngineRestart" key is not used, bonds will be re-detected in the
first step and then preserved during the rest of the simulation.

"NonReactive"
   Type:
      Bool

   Default value:
      No

   GUI name:
      Non-reactive

   Description:
      Enable the non-reactive mode. Bonds are determined only once at
      the beginning and subsequent steps only update their bond
      orders. Thus, no new bonds can form during the simulation, but
      existing bonds can still stretch and dissociate.


Charge equilibration
====================

Details of the charge equilibration (electronegativity equalization
method, EEM) procedure can be adjusted using the "Charges" block.

   Charges
      Constraint
         Charge float
         Region string
      End
      Converge
         Charge float
      End
      DisableChecks Yes/No
      Predictor
         Method [None | Simple]
      End
      Solver [Direct | CG | MINRESQLP | SparseCG | None]
   End

"Charges"
   Type:
      Block

   Description:
      Settings for the polarizable charge model.

   "Constraint"
      Type:
         Block

      Recurring:
         True

      Description:
         Constrain the net charge of a given region.

      "Charge"
         Type:
            Float

         Default value:
            0.0

         Description:
            Desired net charge of the region.

      "Region"
         Type:
            String

         Description:
            Name of the region to be constrained.

   "Converge"
      Type:
         Block

      Description:
         Controls the convergence criteria for charge equilibration.

      "Charge"
         Type:
            Float

         Default value:
            1e-06

         Description:
            Requested upper bound on the sum of squared charge
            residuals.

   "DisableChecks"
      Type:
         Bool

      Default value:
         No

      Description:
         Disable checks for suspicious or unphysical charges.

   "Predictor"
      Type:
         Block

      Description:
         Settings for the prediction of new charges before running the
         solver.

      "Method"
         Type:
            Multiple Choice

         Default value:
            Simple

         Options:
            [None, Simple]

         Description:
            Method used to predict the charges.

   "Solver"
      Type:
         Multiple Choice

      Default value:
         SparseCG

      Options:
         [Direct, CG, MINRESQLP, SparseCG, None]

      Description:
         Algorithm used to solve the charge equilibration equations.


Charge constraints
------------------

The net charge of an arbitrary group of atoms can be constrained to a
particular value using the "Constraint" block. This block can be
repeated as needed to constrain multiple non-overlapping parts of the
system. To define charge constraints, first define appropriate
regions, in the "System" block and then set the "Region" key inside
each "Constraint" block accordingly.

Note:

  Unlike the similar MOLCHARGE constraints in standalone ReaxFF, it is
  not necessary for the constrained regions to span a consecutive
  range of atoms. It is also not necessary to define constraints for
  all atoms in the system. The necessary sum of charges of any
  unconstrained atoms will be determined from the overall net charge
  of the entire system, as set by the "Charge" key in the "System"
  block.

In the following example, we constrain the net charge of one water
molecule in a dimer while the other molecule automatically assumes the
opposite charge to keep the whole system neutral:

   System
       Charge 0.0
       Atoms
           O -0.0509 -0.2754  0.6371 region=donor
           H  0.0157  0.5063  0.0531 region=donor
           H -0.0055 -1.0411  0.0658 region=donor
           O  0.0981  1.7960 -1.2550 region=acceptor
           H -0.6686  2.2908 -1.5343 region=acceptor
           H  0.8128  2.3488 -1.5619 region=acceptor
       End
   End

   Engine ReaxFF
       ForceField Water2017.ff
       Charges
           Constraint Region=donor Charge=0.1
           # The following constraint is implied and need not be specified explicitly.
           # It is only shown here as an example of multiple constraints in a single system.
           Constraint Region=acceptor Charge=-0.1
       End
   EndEngine

Sometimes it may be useful to disable the charge equilibration
altogether and set the charges from input. This can be done using the
"reaxff.charge" atom property, as shown in the example below:

   System
       Atoms
           O -0.0509 -0.2754  0.6371 reaxff.charge=-0.8
           H  0.0157  0.5063  0.0531 reaxff.charge=0.4
           H -0.0055 -1.0411  0.0658 reaxff.charge=0.4
           O  0.0981  1.7960 -1.2550 reaxff.charge=-0.8
           H -0.6686  2.2908 -1.5343 reaxff.charge=0.4
           H  0.8128  2.3488 -1.5619 reaxff.charge=0.4
       End
   End


Atomic stress (per-atom stress tensor)
======================================

The per-atom stress tensor is calculated according to Thompson,
Plimpton, Mattson, J Chem Phys, 131, 154107 (2009) . It does not
include any kinetic contribution (i.e., it is atomic stress, not
atomic pressure).

The calculated stress tensor is stored in the engine .rkf file in the
"RxfAtomData%AtomicStressIso" variable in MPa. It is calculated as
Sαβ/V, where V is atomic volume calculated using the Voronoi
partitioning scheme.

"ComputeAtomicPressure"
   Type:
      Bool

   Default value:
      No

   Description:
      Compute the virial part of the atomic pressure (the kinetic part
      cannot be computed by the engine).
