# ReaxFF input¶

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

## 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 Path to the force field parameter file. Absolute path if starting with / or ./, relative to \$AMSRESOURCES/ForceFields/ReaxFF otherwise.

## 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 No Taper bond orders 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 f10 (eq. 10b) and f12 (eq. 11b) in the original ReaxFF paper. The new expression for each term in f10 is $$\left(1 - e^{-2 \lambda_{23} \text{BO}^2} \right)$$ and in f12 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.

Torsions
Type: Multiple Choice Original [Original, 2013] Version of the torsion potential expression.
FurmanTorsions
Type: Bool No 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 5.0 Angstrom Maximum distance between two atoms to be considered when searching for possible bonds.
BondOrderCutoff
Type: Float 0.001 Minimum bond order required for a bond to be considered during the evaluation of the potential.
StrongBondCutoff
Type: Float 0.3 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 No Non-reactive 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 Settings for the polarizable charge model.
Constraint
Type: Block True Constrain the net charge of a given region.
Charge
Type: Float 0.0 Desired net charge of the region.
Region
Type: String Name of the region to be constrained.
Converge
Type: Block Controls the convergence criteria for charge equilibration.
Charge
Type: Float 1e-06 Requested upper bound on the sum of squared charge residuals.
DisableChecks
Type: Bool No Disable checks for suspicious or unphysical charges.
Predictor
Type: Block Settings for the prediction of new charges before running the solver.
Method
Type: Multiple Choice Simple [None, Simple] Method used to predict the charges.
Solver
Type: Multiple Choice SparseCG [Direct, CG, MINRESQLP, SparseCG, None] 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