Other new features and fixes

Repeating the unit cells

A simple way to increase the size of the simulation box is by multiplying it a number of times in the a, b, and c directions. In the ReaxFF program this can be achieved by specifying the reptx, repty, and reptz keys (for the a, b, and c direction, repsectively) in the control file. The specified value must be strictly positive and the default is 1 (one) for each of them.

Molecular charges

In order to freeze the total charge of a group of atoms (for example, a molecule) one can add a MOLCHARGE key to the geo file in the BGF format. The format of the key is as follows:

MOLCHARGE  n1 n2 q

where n1 and n2 specify indices of the first and the last atoms of the molecule and q is the charge. All atoms between n1 and n2 will be included in the molecule. When specifying the MOLCHARGE key all atoms present in the system must be assigned to one of the the molecules. For example, if your system consists of 10 water molecules and the charge of the first of them is constrained to be +1 and the total charge of the system is zero then the following two MOLCHARGE keys must be specified:

MOLCHARGE  1  3  1.0
MOLCHARGE  4 30 -1.0

Note

The MOLCHARGE is implemented as a simple constraint in the charge equilibration procedure and it does not take any bonding information into account. This means that specifying molecular charges using the MOLCHARGE key makes sense only if you are sure that the constitution of the molecule will not change during calculation. This is useful, for example, when fitting EEM-related parameters in a force-field.

External electric fields

External electric field regimes are specified in an external file eregime.in. The presence of this file in the ReaxFF run directory will activate the electric fields.

The format of this file is as follows:

#Electric field regimes
#start #V direction1 Magnitude1(V/A) direction2 Magnitude2(V/A)
 0000  1      x        0.010000
 1000  1      x       -0.010000
 2000  1      y        0.010000
 3000  1      y       -0.010000
 4000  2      x       -0.010000           y       -0.0100
 5000  2      x        0.010000           y        0.0100

Note

The electrostatic potential is sawtooth-like, which requires a vacuum zone (we suggest 10Å) on the box plane perpendicular to the efield vector. Otherwise one risks abnormal polarization for atoms on two sides of the plane. A good strategy to ensure the vacuum zone is to define an elastic wall restraint (see below).

Example:

Adding the following lines to your run script:

cat > eregime.in <<eor
#Electric field regimes
#start #V direction Magnitude(V/Angstrom)
0000 1 x  0.010000
1000 1 x -0.010000
eor

will overlay an external electric field with magnitude 0.01 V/Å in x-direction for the duration of 1000 iterations, then the magnitude is switched to -0.01 V/Å and kept for the rest of the simulation.

Elastic wall restraint

The wall is implemented as a \(cos^2(x)\) function (one period) between Coord-HalfWidth and Coord+HalfWidth with a maximum at Coord. The Height value is given in kcal/mole. At this moment, the elastic wall restraint can only be used with an orthorombic unit cell.

The wall is specified in the geo file as:

EWALL RESTRAINT AxisIndex Coord HalfWidth Height

For example, for a wall that crosses the Z axis at z=7 Angstrom, is 2.0 Angstom thick and 100 kcal/mole high, the input would look like this:

EWALL RESTRAINT  3  7.0  1.0  100.0

Note

We recommended to set the icentr option in the control file to 0 to avoid inconsistencies between the initial atomic coordinates and the wall position.

ACKS2: Atom-condensed Kohn-Sham DFT approximated to second order

The ACKS2 charge equilibration scheme has been implemented following the original paper “ACKS2: Atom-condensed Kohn-Sham DFT approximated to second order” by T. Verstraelen et al. J. Chem. Phys. 138 (2013) 074108.

Using the ACKS2 scheme requires a suitable force-field, which is recognized by “[ acks2 ]” at the start of the first line of force field file (note: the spaces around “acks2” are important!). Besides, the icharg parameter in the control file must be set to 7. In addition to the general EEM parameters the ACKS2 scheme needs the general force-field parameter #35 (“Xamp”) and the atomic cut-off parameter #23 (“softcut”).

Correction for torsion angles asymptotics

There is a discontinuity problem for small bond orders in the expression for torsion angles and conjugation contributions f(BO). These terms correspond to expressions for 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} BO^2} \right)\) and in f12 the new expression is \(sin(\frac{\pi}{3} BO)^4\). The new expressions ensure correct asymptotic behavior for the dE/dBO for BO \(\rightarrow\) 0. Using this makes geometry optimizations more stable and improves conservation of energy during MD.

The correction can be enabled by setting the tors13 flag in the control file to 1.

LG dispersion

The LG dispersion correction was implemented following the paper Liu et al., *ReaxFF-lg: Correction of the ReaxFF Reactive Force Field for London Dispersion, with Applications to the Equations of State for Energetic Materials*, J. Phys. Chem. A, 2011, 115 (40), pp 11016–11022

The LG dispersion correction is turned on when using a suitable forcefield, which is recognized by the “[ lgDispersion=1 ]” key in the file header.

The rxffutil utility

Usage:

rxffutil command [arguments]

The ‘reformat’ command tells rxffutil to read a (possibly mis-formatted) force-field file from standard input and write it in the format required by the ReaxFF program:

rxffutil reformat < input_ffield_file

The result is written as a ‘ffield’ file in the current directory. The input can be almost free-format but the following points need to be taken into account:

  • If you copy-paste force-field parameters from a pdf document (for example, from a publication’s supporting info), then make sure to remove any page numbers that may have been copied together.
  • The atom types section header (with the number of atom types in it) must have three comment lines following the line with the number of atom types.
  • Likewise, the bond types section must have one comment line following the line with the number of bond types. Any comments provided on input are ignored and the standard comments are written instead.
  • Numbers must be separated by at least one space