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

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₁₀ (eq. 10b) and f₁₂ (eq. 11b) in the original ReaxFF paper. The new expression for each term in f₁₀ is \(\left(1 - e^{-2 \lambda_{23} \text{BO}^2} \right)\) and in f₁₂ 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 potential.

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