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.

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.

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

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