In AMS2020, ReaxFF integrated with the AMS driver is the default. This gives consistent input for different modules and enables new properties through the central driver. It also helps users more easily switch between modules and use ReaxFF in multi-layer hybrid simulations.

Some ReaxFF features have not yet been ported to the integrated ReaxFF, and can be used from the standalone ReaxFF.

Click on i next to Force field input in ReaxFF Input to see which force fields are included with ReaxFF. Alternatively, open $AMSHOME/atomicdata/ForceFields/ReaxFF/Readme, which also contains the original references and sometimes a short description. The list of currently available ReaxFF force fields are available in the documentation.

You are strongly advised to check the original publications to see how the force field has been trained and to what type of systems it has been applied.

If no force field is available for your elements, check the literature or parametrize it yourself, see the other FAQs.

This is hard work, which takes quite a bit of training and experience.

The ReaxFF parameters are typically trained against DFT data and parameters should be optimized for all atom types in your system. In principle, parameters are transferable across systems, but including additional atoms requires also to include new interactions with all new element combinations. This is not trivial, but we try to help you with our ParAMS ReaxFF parametrization tools.

RxFF consulting offers consulting services for optimizing force field parameters for new materials.

There are two elements in your simulation box for which no bonded parameters are present in the force field file. Please check the description and the paper on the ReaxFF force field page.

You can find more info on which bond, angle, and dihedral terms exist, and which ones are missing by using the checkreaxffforcefield.py utility:


$AMSBIN/startpython $AMSHOME/scripting/standalone/reaxff/checkreaxffforcefield.py /path/to/yourreaxff.ff

You’ll probably find that this force field was not designed for the elements you want to study. Perhaps there is another force field more suitable for your application.

You may try to our reparametrize the force field.

In the GUI: build a ReaxFF input with these elements X, Y, and Z, and click on the folder icon next to Force field:. A list of force fields included in the AMS/ReaxFF distribution that contain these atom types will be shown. If force fields are available, consult the Readme (click ‘i’) to check the publication to make sure this is the force field for your application.

If no force field shows up in the list, no suitable force field is available in your distribution. Use google or web of knowledge to find out if a ReaxFF force field exists for your application. If it does, and the force field is printed in the supporting information, you can add the ReaxFF force field file (see FAQ) to the existing ones in the distribution.

Note that not all terms may be available for all element combinations. To get more insight in which force field terms are available, use the python ReaxFF force field checking tool:


$AMSBIN/startpython $AMSHOME/scripting/standalone/reaxff/checkreaxffforcefield.py /path/to/yourreaxff.ff

You can find a description of the general structure of the force field file and the parameters from the original ReaxFF documentation.

We try to include as much information on the origin of parameters the ReaxFF manual as well as in the GUI help files.

Many published ReaxFF force fields are building upon a previous set of parameters, and so only the elements studied in that paper are being refined, and the other parameters are just kept in place. You could try to trace back the origin of some of these parameters by looking at the existing force fields and finding out where the same parameters for a particular element also appear. Any force field that has matching atomic parameters might have some information about the origins of the parameters in the corresponding publication or in the information in the force field description file.

However, you need to be careful by just combining and using elements / parameters from different force fields. See also the FAQ about ReaxFF parameter transferability. You should also be aware that certain bonded terms may not be available for a particular unpublished element and/or combination of elements, so you are strongly advised to check this with our ReaxFF force field checking tool: $AMSHOME/scripting/checkreaxffforcefield.py. The original ReaxFF documentation describes the general structure and correspondence of the columns with the parameters.

In general the parameters are quite transferable, but some forcefields were developed for specialized cases such as explosives research (high pressure). It is good practice to scan through the publication of the forcefield (you can find the information by clicking the “i” button behind the forcefield ) before using it extensively. It is also wise to do some testing simulations to see if your structures are not falling apart or folding into weird geometries.

This is most easily done from the GUI using the File → Export Trajectory As: option in AMSJobs. You can export as Gromacs, PDB trajectory, or xmol file. For the .gro export you will get an additional dialogue window where you can include the timestep and the velocities.

Alternatively, you can use the scripting tool amsreport, e.g. amsreport yourreaxff.rxkf -r gro-t0.00025 > yourreaxff.gro
The options for amsreport in your AMS distribution are displayed if you type amsreport -h

During a non-reactive ReaxFF simulation the connection table is determined only once at the beginning of a simulation and is then kept fixed. Other than this the same ReaxFF procedure with the same forces and energy expressions is used. Therefore, if the system is far from equilibrium, or if the temperature is too high, bonds may still be broken or formed.

You can also think about simulating the system at zero Kelvin with a strong thermostat to try and get your system to the nearest local minimum without changing the bonds too much.

This depends on what parameters are available in the force field file. A ReaxFF force field consists of some generic parameters, a number of atomic parameters (1 set for every atom type), bond parameters (combination of 2 atom types), angle parameters (combination of 3 atom types), dihedral parameters (combination of 4 atom types) and possibly hydrogen bond parameters (see original ReaxFF documentation for more details).

If a force field has elements X, Y an Z it means that it contains atomic parameters for X, Y and Z. However, it does not necessarily contain all combinations of the bond/angle/dihedral parameters. For example, X-X, Y-Y, Z-Z, X-Y and X-Z might have bond parameters, but Y-Z could be left out. This means that Y and Z should not be used together in the simulation, as their interaction is not defined.

To get more insight in which terms are available in a ReaxFF force field, use the python ReaxFF force field checking tool: $AMSHOME/scripting/checkreaxffforcefield.py

ReaxFF consists of two parts: the program and the force field files:

– The ReaxFF program is an implementation of the ReaxFF forcefield equations (see the Required citations ReaxFF page for the original papers).
– The ReaxFF force field files contain the parameters needed by the force field equations when running your simulation, they consist of 5 or 6 parts (see original ReaxFF documentation for more details):
* generic parameters
* atom parameters (per element)
* atom pairs / bond parameters (combination of two elements)
* angle parameters (combination of three elements)
* dihedrals (combination of four elements)
* and sometimes hydrogen bonds (per element)

A ReaxFF force field is build using a training set (see other FAQs) and is usually targeted for a specific simulation and typically includes about 4 to 12 elements, as the number of parameters scales to the 4th power with the number of elements included. Because ReaxFF uses pair/angle/dihedral parameters, it is not possible to simply combine two forcefield files.

The AMS GUI can tell you if a force field is available for the elements in your simulation (see other FAQs), but it will not be able to tell if the force field is conditioned for your specific simulation type. For this you will have to check the publication of the force field (doi link is usually included in the Readme file or the first line of the force field file)

The AMS implementation of ReaxFF still takes the original format. If a publication has the parameters printed in the supporting information in native format, just copy that data and save it as a file, e.g. YourNewReaxFF.ff, in $AMSHOME/atomicdata/ForceFields/ReaxFF. Then you can start using that force field by selecting it in the GUI, by clicking other after you’ve clicked the folder icon next to Force field: in the input panel.

If the parameters are printed in tables in the manuscript, but not as a comprehensive SI file, you can still compile your own ReaxFF force field file from this data in the same original ReaxFF format.

It is advisable to use the checkreaxffforcefield.py script in $AMSHOME/scripting/ to check your force field. If you have modified or downloaded a force field you may want to clean it up:

Open a shell in a terminal with AMS settings loaded. Then run the rxffutil tool

rxffutil reformat < your_ff.ff

The resulting cleaned up ReaxFF force field is stored in a file ‘ffield’.

For the NHC themostat/MTK barostat, the simulation samples the correct ensemble no matter what the damping time is set to, but it is recommended to set the damping time to roughly correspond to the correlation time of the system (how long it takes for the system to “forget” a disturbance). This is 2-3 picoseconds in liquid water, up to a couple hundred fs in solids. Setting it any longer doesn’t improve the accuracy of the system anymore, just makes the thermostatting/barostatting less efficient (taking longer to settle).

Berendsen thermostats/barostats only approximate the correct ensemble and this approximation is the worse the shorter the damping time. It is thus crucial to set the damping time as high as practical to still reach the target value in acceptable time. Especially for the Berendsen thermostat, setting the damping time too low is known to result in severe equipartition artifacts (all kinetic energy is transferred from vibrations into rotational degrees of freedom, leading to the so-called “flying ice cube problem“).

Setting the damping constant/relaxation time to a lower value means more aggressive thermostating/barostating. The system will be driven to reach the target temperature/pressure faster, which may be necessary if it starts far from equilibrium or if there’s something in the system releasing or consuming a lot of energy (a chemical reaction). However, using a short relaxation time may also lead to artifacts (e.g. rapid changes of the volume by an aggressive barostat may disturb the structure of a condensed-phase system). Conversely, setting the relaxation time too high will make the thermostat/barostat very weak, affecting the system only minimally but also taking a lot of time to reach the target temperature/pressure.

The computational cost of any thermostat/barostat in AMS doesn’t depend on the damping at all.

See also the webinar by Tomas Trnka on MD basics (slides)

It is possible to constrain the net charge of a particular atom or molecule to a particular value, this can be done in the GUI under Model→Charge Constraints.

For the old standalone ReaxFF, charge constraints are defined by putting lines like the following into the “geo” file:

MOLCHARGE firstAtom lastAtom netCharge

where firstAtom and lastAtom define the range of atoms to be constrained (inclusive, numbering starts from 1). An unfortunate limitation of this approach is that all atoms in the system need to be covered by one of the MOLCHARGE lines because the ranges need to be continuous, this often requires some reordering of the atoms to be possible. 

As an example, let’s say there are two ions that we want to constrain to +2 charge:

MOLCHARGE 1 1 2.0 
 MOLCHARGE 2 2 2.0
 MOLCHARGE 3 1000 -4.0

We have moved the ions to the beginning of the system to be able to constrain the rest of the system and make sure the net charge of the whole simulation box is zero.

From on AMS2020, the AMS driver supports charge constraints for ReaxFF.

They can be defined in the ReaxFF Engine block as follows:

Engine ReaxFF
Charges
Constraint Region=ion1 Charge=2.0 
Constraint Region=ion2 Charge=2.0
End
EndEngine

One thus only needs to put the ions into separate regions (for example using the GUI) and use the appropriate region names on the Constraint lines. There’s no need to reorder atoms here or to constrain the rest of the system (this happens naturally to keep the total charge at the value set in the System block).