ReaxFF input¶
This section describes the input keywords to the ReaxFF AMS engine.
See also
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.
ForceFieldType: String Description: Path to the force field parameter file. Absolute path if starting with / or ./, relative to $AMSRESOURCES/ForceFields/ReaxFF otherwise.
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.
TaperBOType: 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 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.
See also
TorsionsType: Multiple Choice Default value: Original Options: [Original, 2013] Description: Version of the torsion potential expression.
FurmanTorsionsType: 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¶
BondDistanceCutoffType: Float Default value: 5.0 Unit: Angstrom Description: Maximum distance between two atoms to be considered when searching for possible bonds.
BondOrderCutoffType: Float Default value: 0.001 Description: Minimum bond order required for a bond to be considered during the evaluation of the potential.
StrongBondCutoffType: 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.
NonReactiveType: 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
ChargesType: Block Description: Settings for the polarizable charge model. ConstraintType: Block Recurring: True Description: Constrain the net charge of a given region. ChargeType: Float Default value: 0.0 Description: Desired net charge of the region. RegionType: String Description: Name of the region to be constrained.
ConvergeType: Block Description: Controls the convergence criteria for charge equilibration. ChargeType: Float Default value: 1e-06 Description: Requested upper bound on the sum of squared charge residuals.
DisableChecksType: Bool Default value: No Description: Disable checks for suspicious or unphysical charges. PredictorType: Block Description: Settings for the prediction of new charges before running the solver. MethodType: Multiple Choice Default value: Simple Options: [None, Simple] Description: Method used to predict the charges.
SolverType: 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