pair_style reax command
pair_style reax hbcut hbnewflag tripflag precision
- hbcut = hydrogen-bond cutoff (optional) (distance units)
- hbnewflag = use old or new hbond function style (0 or 1) (optional)
- tripflag = apply stabilization to all triple bonds (0 or 1) (optional)
- precision = precision for charge equilibration (optional)
pair_style reax pair_style reax 10.0 0 1 1.0e-5 pair_coeff * * ffield.reax 3 1 2 2 pair_coeff * * ffield.reax 3 NULL NULL 3
Style reax computes the ReaxFF potential of van Duin, Goddard and co-workers. ReaxFF uses distance-dependent bond-order functions to represent the contributions of chemical bonding to the potential energy. There is more than one version of ReaxFF. The version implemented in LAMMPS uses the functional forms documented in the supplemental information of the following paper: (Chenoweth). The version integrated into LAMMPS matches the most up-to-date version of ReaxFF as of summer 2010.
WARNING: pair style reax is now deprecated and will soon be retired. Users should switch to pair_style reax/c. The reax style differs from the reax/c style in the lo-level implementation details. The reax style is a Fortran library, linked to LAMMPS. The reax/c style was initially implemented as stand-alone C code and is now integrated into LAMMPS as a package.
LAMMPS requires that a file called ffield.reax be provided, containing the ReaxFF parameters for each atom type, bond type, etc. The format is identical to the ffield file used by van Duin and co-workers. The filename is required as an argument in the pair_coeff command. Any value other than “ffield.reax” will be rejected (see below).
LAMMPS provides several different versions of ffield.reax in its potentials dir, each called potentials/ffield.reax.label. These are documented in potentials/README.reax. The default ffield.reax contains parameterizations for the following elements: C, H, O, N.
We do not distribute a wide variety of ReaxFF force field files with LAMMPS. Adri van Duin’s group at PSU is the central repository for this kind of data as they are continuously deriving and updating parameterizations for different classes of materials. You can submit a contact request at the Materials Computation Center (MCC) website https://www.mri.psu.edu/materials-computation-center/connect-mcc, describing the material(s) you are interested in modeling with ReaxFF. They can tell you what is currently available or what it would take to create a suitable ReaxFF parameterization.
The format of these files is identical to that used originally by van Duin. We have tested the accuracy of pair_style reax potential against the original ReaxFF code for the systems mentioned above. You can use other ffield files for specific chemical systems that may be available elsewhere (but note that their accuracy may not have been tested).
The hbcut, hbnewflag, tripflag, and precision settings are optional arguments. If none are provided, default settings are used: hbcut = 6 (which is Angstroms in real units), hbnewflag = 1 (use new hbond function style), tripflag = 1 (apply stabilization to all triple bonds), and precision = 1.0e-6 (one part in 10^6). If you wish to override any of these defaults, then all of the settings must be specified.
Two examples using pair_style reax are provided in the examples/reax sub-directory, along with corresponding examples for pair_style reax/c. Note that while the energy and force calculated by both of these pair styles match very closely, the contributions due to the valence angles differ slightly due to the fact that with pair_style reax/c the default value of thb_cutoff_sq is 0.00001, while for pair_style reax it is hard-coded to be 0.001.
Use of this pair style requires that a charge be defined for every atom since the reax pair style performs a charge equilibration (QEq) calculation. See the atom_style and read_data commands for details on how to specify charges.
The thermo variable evdwl stores the sum of all the ReaxFF potential energy contributions, with the exception of the Coulombic and charge equilibration contributions which are stored in the thermo variable ecoul. The output of these quantities is controlled by the thermo command.
This pair style tallies a breakdown of the total ReaxFF potential energy into sub-categories, which can be accessed via the compute pair command as a vector of values of length 14. The 14 values correspond to the following sub-categories (the variable names in italics match those used in the ReaxFF FORTRAN library):
- eb = bond energy
- ea = atom energy
- elp = lone-pair energy
- emol = molecule energy (always 0.0)
- ev = valence angle energy
- epen = double-bond valence angle penalty
- ecoa = valence angle conjugation energy
- ehb = hydrogen bond energy
- et = torsion energy
- eco = conjugation energy
- ew = van der Waals energy
- ep = Coulomb energy
- efi = electric field energy (always 0.0)
- eqeq = charge equilibration energy
To print these quantities to the log file (with descriptive column headings) the following commands could be included in an input script:
compute reax all pair reax variable eb equal c_reax variable ea equal c_reax ... variable eqeq equal c_reax thermo_style custom step temp epair v_eb v_ea ... v_eqeq
Only a single pair_coeff command is used with the reax style which specifies a ReaxFF potential file with parameters for all needed elements. These are mapped to LAMMPS atom types by specifying N additional arguments after the filename in the pair_coeff command, where N is the number of LAMMPS atom types:
- N indices = mapping of ReaxFF elements to atom types
The specification of the filename and the mapping of LAMMPS atom types recognized by the ReaxFF is done differently than for other LAMMPS potentials, due to the non-portable difficulty of passing character strings (e.g. filename, element names) between C++ and Fortran.
The filename has to be “ffield.reax” and it has to exist in the directory you are running LAMMPS in. This means you cannot prepend a path to the file in the potentials dir. Rather, you should copy that file into the directory you are running from. If you wish to use another ReaxFF potential file, then name it “ffield.reax” and put it in the directory you run from.
In the ReaxFF potential file, near the top, after the general parameters, is the atomic parameters section that contains element names, each with a couple dozen numeric parameters. If there are M elements specified in the ffield file, think of these as numbered 1 to M. Each of the N indices you specify for the N atom types of LAMMPS atoms must be an integer from 1 to M. Atoms with LAMMPS type 1 will be mapped to whatever element you specify as the first index value, etc. If a mapping value is specified as NULL, the mapping is not performed. This can be used when a ReaxFF potential is used as part of the hybrid pair style. The NULL values are placeholders for atom types that will be used with other potentials.
Currently the reax pair style cannot be used as part of the hybrid pair style. Some additional changes still need to be made to enable this.
As an example, say your LAMMPS simulation has 4 atom types and the elements are ordered as C, H, O, N in the ffield file. If you want the LAMMPS atom type 1 and 2 to be C, type 3 to be N, and type 4 to be H, you would use the following pair_coeff command:
pair_coeff * * ffield.reax 1 1 4 2
Mixing, shift, table, tail correction, restart, rRESPA info:
This pair style does not support the pair_modify mix, shift, table, and tail options.
This pair style does not write its information to binary restart files, since it is stored in potential files. Thus, you need to re-specify the pair_style and pair_coeff commands in an input script that reads a restart file.
This pair style can only be used via the pair keyword of the run_style respa command. It does not support the inner, middle, outer keywords.
The ReaxFF potential files provided with LAMMPS in the potentials directory are parameterized for real units. You can use the ReaxFF potential with any LAMMPS units, but you would need to create your own potential file with coefficients listed in the appropriate units if your simulation doesn’t use “real” units.
The keyword defaults are hbcut = 6, hbnewflag = 1, tripflag = 1, precision = 1.0e-6.
(Chenoweth_2008) Chenoweth, van Duin and Goddard, Journal of Physical Chemistry A, 112, 1040-1053 (2008).