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
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₁₀ (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.

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