8.3.1.2. ReaxFF
***************

The "ReaxFFParameters" class provides a parameterizable interface to
the ReaxFF reactive force field. It can be initialized by passing a
force field path to the constructor:

   ff = ReaxFFParameters('path/to/ffield.ff')

By default, all relevant parameters will be marked as
*is_active==True* (*i.e.* marked for parameterization). See the API
for parameter naming conventions.

See also:

  * ReaxFF force field format specification

  * "Detailed SCM developer notes" about all ReaxFF parameters


Contents of this Page
^^^^^^^^^^^^^^^^^^^^^

* Reading Legacy params files

* Addition of Parameters

* Trim (remove elements) from a ReaxFF force field

* Apply recommended parameter ranges (constraints)

* Parameter Databases and Value Suggestions

* ReaxFF Conversion Functions

* ReaxFF Interface API

* Full List of Parameter Names

* Parameter Categories

======================================================================


8.3.1.2.1. Reading Legacy params files
======================================

Outside of the ParAMS formalism, a params file is common to mark a set
of parameters for optimization, and set their respective ranges.
Support for these files is implemented in the
"ReaxFFParameters.read_paramsfile()" method:

   ff = ReaxFFParameters('path/to/ffield.ff')
   ff.read_paramsfile('path/to/params')

The method will mark all matching parameters as *is_active=True* (and
deactivate all remaining ones). It will raise an error if blocks that
are not present in the force field are referenced in the *params*
file. Note that the 4th column in the *params* file (delta value) will
be ignored.


8.3.1.2.2. Addition of Parameters
=================================

The ReaxFF interface allows for the addition of new parameters through
the "add_blocks()" method.


8.3.1.2.3. Trim (remove elements) from a ReaxFF force field
===========================================================

Use the "trim()" method.


8.3.1.2.4. Apply recommended parameter ranges (constraints)
===========================================================

Note: This is only about setting technically meaningful upper/lower
bounds on parameter values. Even when you apply these recommendations
you are likely to have to modify the ranges yourself for your
particular problem.

* GUI: Select the parameters in the parameters table for which you
  would like to apply the recommended ranges. Then select **Parameters
  → Apply Recommended Constraints**. The recommended constraints are
  documented in the "apply_default_constraints()" method.

* Python: Use the "apply_default_constraints()" method.


8.3.1.2.5. Parameter Databases and Value Suggestions
====================================================

The class can suggest suitable parameter values for a number of
queried parameters. See the "ReaxFFParameters.params_from_db()" method
for more information. By default, all force fields in
*$AMSHOME/atomicdata/ForceFields/ReaxFF* will be used to return
suggestions. Users can create their own databases from fewer (or more)
force fields through the "ReaxFFParameters.generate_paramsdb()"
method.


8.3.1.2.6. ReaxFF Conversion Functions
======================================

Outside of ParAMS, ReaxFF is parameterized with a training set defined
in *trainset.in* and *geo* files. There are two functions available to
help with the conversion to the ParAMS format:

trainset_to_params(trainsetfile, use_weights=True, default_sigmas=None) -> DataSet

   Converts a ReaxFF *trainset.in* file to a "DataSet" object.

   Parameters:
   trainsetfile : str
      Path to *trainset.in*.

   use_weights : bool
      If False, all weights will be set to 1.0, and Sigma will be set
      to the value from trainset.in

      If True, then new_weight/(new_sigma**2) ==
      1/(value_from_trainset_in**2).

      new_sigma is taken from the standard Sigma value for the
      pertinent extractor, or from the default_sigmas dictionary.

   default_sigmas : dict
      A dictionary of the type

      { 'energy': sigma_energy, 'forces': sigma_forces, 'charges':
      sigma_charges, 'angles': sigma_angles,
         'distance': sigma_distance, 'dihedral': sigma_dihedral }

         The sigma values should be given in the classic ReaxFF units,
         so kcal/mol for energy, angstrom for distance, etc.

         The sigma values are taken from this dictionary only if
         use_weights == True

   Returns:
   ds : "DataSet" instance

geo_to_params(geofile: str, normal_run_settings: Settings | str) -> JobCollection

   Convert a ReaxFF *geo* file to a "JobCollection" instance.

   * Each entry's ID will be based on the *DESCRP* line. Duplicates
     will not be added.

   * Will guess the settings for each job based on the *RUTYPE* line.
     Can either be *SinglePoint* or *GeometryOptimization* (optionally
     with *MaxIterations*).

   Parameters:
   geofile : str
      Path to *geo* file.

   normal_run_settings : plams.Settings, str
      A "plams.Settings" instance representing the *NORMAL RUN* or the
      path to a *control* file from which the settings will be
      extracted

   Returns:
   jc : "JobCollection"


8.3.1.2.7. ReaxFF Interface API
===============================

class ReaxFFParameters(ffieldfile, bounds_scale=1.2, settings=None)

   Interface to the ReaxFF force field format.

   Attributes:
   In addition to the ones already availabe in "BaseParameters"

   header : dict
      A dictionary with metadata. *header['head']* contains the first
      line of the force field

   path : str
      The path to the last written ffield file.

   blocks : List[str]
      Strings that mark valid parameter blocks. See ReaxFF manual for
      more information.

      Note:

        This is a class attribute, meaning that an instance can
        contain fewer blocks.

   allnames : Dict[str : List[str]]
      A dictionary with *BLOCK_NAME:List[PARAMETER_NAMES]*, where
      *BLOCK_NAME* is one of the elements from "blocks". See also
      "name_to_index()" and "names_in_block()".

      Note:

        This is a class attribute, meaning that an instance can
        contain fewer parameters.

   shape : Dict[str : List[int]]
      A dictionary of *BLOCK_NAME:[NBLOCKS, N_PARAMS_IN_BLOCK]*, where
      *BLOCK_NAME* is one of the elements from "blocks", *NBLOCKS* is
      the number of *BLOCK_NAME* blocks and *N_PARAMS_IN_BLOCK* is the
      number of parameters per block. See ReaxFF manual for more
      information.

   flags : dict[str, str]
      This attribute stores flags for special force field flavours.
      The data is read from the first line of the force field file, if
      it starts with *[ X ]*, where *X* is a series of flag strings.
      Currently, the following flags are understood: *ereaxff*,
      *ereaxff2*, *acks2*, *lgDispersion*.

   **Parameter Naming Conventions:**

   Naming follows the tables as documented in the ReaxFF manual.  The
   name string is generated from the table's *Name in Eq* column.  All
   atom-specific parameter names will start with the atom prefix,
   followed by the name: *At1.At2.At3:NAME*. For example, the name of
   the first parameter in the atoms block for a Hydrogen will be
   "H:r_0^sigma".  See below for a full list of names.  Each parameter
   in the ReaxFF interface has the following additional attributes:
   "block", "block_index", "equation", "description",
   "atoms_reversible". For easy filitering of parameters, see "get()".
   For conversions between block indices and name strings, see the
   "blocks", "name_to_index()" and "names_in_block()".

   class Parameter(block, block_index, atoms, bounds_scale=0, *a, **kw)

      __init__(block, block_index, atoms, bounds_scale=0, *a, **kw)

   __init__(ffieldfile, bounds_scale=1.2, settings=None)

      Initialize the class instance, reading a ffied from
      *ffieldfile*.

      Parameters:
      ffieldfile : str
         Path to the force field file.

      bounds_scale : float >= 0
         Apply arbitrary bounds, by providing a factor which will
         determine lower/upper bounds by "b = x0 ±
         (bounds_scale-1)*abs(x0)".

         Note:

           Bounds for any "x0 == 0" will be set to "+-1".

      settings : optional, Settings
         A *plams.Settings* instance that will be added to the output
         of "get_engine()". Must have the *input.ReaxFF* keys defined.

   __str__()

      String representation of self

   get_engine(parameters=None, path=None, write=True)

      Note:

        When called without *path*, will store a ffield file in a
        temporary directory for as long as the returned "Engine"
        instance is referenced.

      Returns:
         A ReaxFF "Engine" class matching the passed *params*, or the
         current set of parameters if *None*.

   write(file=None, parameters=None)

      Writes the current parameters to *file*.

      Note:

        When called without *file*, will store the ffield file in a
        temporary directory.  Use the "path" attribute to get the file
        path.

   get_constraints() -> List[Constraint]

      Returns a list of physically reasonable Constraints for this
      instance to be passed to the Optimization.

      Will constrain the covalent radii and bond lengths r_0 for the
      \sigma, \pi, \pi\pi bonds to r_0^\sigma >= r_0^\pi >= r_0^\pi\pi
      and the Van der Waals radius r_\mathrm{VdW} to r_\mathrm{VdW} >=
      r_0^sigma`.

   read_paramsfile(paramsfile)

      Reads a legacy params file, sets all matching parameters to
      *is_active=True* (*i.e.* flags them as paramters to be fitted),
      and adjusts each parameter's ranges to the ones from
      *paramsfile*.

   add_blocks(params: List[str], allow_duplicates=False)

      Adds new parameters to the current instance.  The addition
      happens block-wise by providing a list of strings, where each
      string represents the addition of one block. Strings must have
      the form *At1.At2.At3.At4*, where *At* represents the chemical
      symbol of the desired element (separated by dots) and the number
      of elements determines the block to be added (1 for the atomic
      block, 2 for bonds, 3 for angles, 4 for torsions). For example,

      >>> ffield.add_blocks(['C.O', 'C.O.H'])

      adds the C-O bonds block and the C-O-H angles block to the
      *ffield* instance. Off-diagonals and hydrogen-bonds can be added
      by prefixing a string with *OFD:* and *HBD:* respectively:

      >>> ffield.add_blocks(['OFD:C.O', 'HBD:O.H.C'])

      The above will add the C-O off-diagonal and O-H-C hydrogen-bond
      block.

      Note:

        The following is of note when adding new blocks:

        * The addition of multiatomic blocks to a force field that
          does not yet contain all respective elements will
          automatically add the atomic block for that element. *E.g.*,
          the addition of C-O to a force field that does not yet
          contain parameters for oxygen also implies the addition of O

        * Parameters added by this method have a value of *0.* and a
          range of *(-1, 1)*

        * Already present blocks will not be overwritten, however, a
          duplicate block can be added by setting
          *allow_duplicates=True*

   classmethod generate_paramsdb(paths=['$AMSHOME'], dbfile='reaxffdb.gz')

      Generates a parameter database based on all force field files
      found in *paths* and stores it in *dbfile*.

      Parameters:
      paths : List[str]
         When generating the database, each entry in *paths* can be
         either a directory, or a full path to a force field file. If
         a directory, will scan recursively and consider all files
         with the *.ff* ending.

      dbfile : optional, str
         Stores the generated database in *dbfile*. Defaults to
         *reaxffdb.gz*.

   classmethod params_from_db(params: List[str], dbfile=None, **subst_atoms) -> List[namedtuple]

      Given a list of *n* parameter queries in *params*, scans a force
      field database stored in *dbfile* and returns *n* named tuples,
      which hold all parameters matching the query.

      Properties:
      params : List[str]
         Scans the data base for all entries in *params*. Each element
         is a string with either of the following syntaxes:

         * *'BLOCK:At1.At2.At3.At4:INDEX'* **or**

         * *'At1.At2.At3.At4:NAME'*

         where *At1.At2.At3.At4* are chemical symbols (one to four).
         For the general block, this is replaced by *None*.

         In the former case, *BLOCK* has to be any of *{'GEN', 'ATM',
         'BND', 'OFD', 'ANG', 'TOR', 'HBD'}* referring to the General,
         Atoms, Bonds, Off-Diagonals, Angles, Torsions, Hydrogen-Bonds
         blocks respectively. And *INDEX* is an integer specifying the
         parameter number for the block (starting with *0*).

         In the latter case, *NAME* refers to the string name of the
         desired parameter, following the same naming as in the
         "names" attribute (see Full List of Parameter Names).

      dbfile: optional, str
         Path to the database file to be used for the search. Defaults
         to a database generated from
         *$AMSHOME/atomicdata/ForceFields/ReaxFF*. See
         "generate_paramsdb()" for more information.

      subst_atoms: optional keywords
         Additional keywords will be interpreted as chemical elements
         (*keys*) that can be substitued by a number of other elements
         (*values*). This extends each query by matching all valid
         combinations of chemical elements for a specific parameter,
         rather than limiting it to one specific set of elements.
         Useful when queries should also include chemically similar
         elements.

         Example:

         >>> cls.params_from_db(['BND:C.O:0'], C=('Si', 'Ge'), O=('S',))

         Will search all first entries of the bonds block considering
         the follwing combinaitons of elements: C.O, C.S, Ge.O, G.S,
         Si.O, Si.S.

      Returns:
      Returns a list of named 4-tuples. Each element in the list
      corresponds to a query. Each named tuple unpacks to *(query,
      parametername, values, atoms, sources)*. Elements within the
      tuple are stored in the respective attribute (*e.g.*
      "tup.values").

      query : str
         String used for the query, Same as in *params*

      parametername : str
         String with the requested parameter's name (excluding the
         chemical symbols prefix)

      values : array
         Array of all parameter values matching the query

      atoms : array of str
         Array of chemical elements corresponding to each returned
         parameter value

      sources : array of str
         Array of force field paths corresponding to the origin of
         each parameter value

   classmethod block_from_db(block: str, dbfile=None)

      block : str
         Should include the atoms.  'GEN', 'BND:C.Cl', 'ANG:C.C.H',
         etc.

      dbfile: optional, str
         Path to the database file to be used for the search. Defaults
         to a database generated from
         *$AMSHOME/atomicdata/ForceFields/ReaxFF*. See
         "generate_paramsdb()" for more information.

      Returns a 2-tuple containing a dictionary and a list of
      parameter names. The keys of the dictionary are the sources
      (paths to .ff files), and the values are lists of parameter
      values.

   names_in_block(block: str | int, atoms: str = '') -> List[str]

      Returns the names of the all parameters in *block*, optionally
      prepended by an *{atoms}:* string. *block* can either be a an
      int in the range of *[0, 6]* or any of *{'GEN', 'ATM', 'BND',
      'OFD', 'ANG', 'TOR', 'HBD'}* referring to the General, Atoms,
      Bonds, Off-Diagonals, Angles, Torsions, Hydrogen-Bonds blocks
      respectively.

      Note:

        Index numbering starts from zero

      See also: "block_to_names()"

   classmethod name_to_index(name: str, block_as='str') -> Tuple[int | str, int]

      Given a parameter name, returns the indices at which this name
      is located within the force field blocks. The *name* argument
      should be a string of *At1.At2:NAME*, where *At* are the atoms
      (one to four) the parameter affects.  If *block_as=='str'*, the
      returned value is a *(str, int)*, otherwise *(int, int)*, where
      the first element specifies the block and the second one the
      index at which this is parameter located within the block.
      Blocks are represented by an int in the range of *[0, 6]* or
      strings of *{'GEN', 'ATM', 'BND', 'OFD', 'ANG', 'TOR', 'HBD'}*
      referring to the General, Atoms, Bonds, Off-Diagonals, Angles,
      Torsions, Hydrogen-Bonds blocks respectively.

      Note:

        Index numbering starts from zero.

   copy_block(other_interface: ReaxFFParameters, from_block: str, to_block: str | None = None, copy_range: bool = True, allow_duplicates=False)

      Copies the values (and optionally ranges) of a whole parameter
      block from *other_interface* to this instance.

      Parameters:
      other_ff: ReaxFFParameters
         Copy the parameters in the block from this other instance

      from_block: str
         'BND:C.C' or similar. It is also possible to copy the
         complete block, _i.e._ all values in 'BND' or 'OFD' by
         passing only the three-letter block name.

      to_block : str or None
         'BND:C.O' or similar. If None, set to the same value as
         from_block

      copy_range: bool
         Whether to copy the .range attribute from the other interface

   get(name: str | None = None, atoms: str | None = None, block: str | int | None = None, index: int | None = None) -> List[Parameter]

      Returns a list of parameters in self based on a set of filters.
      Filtering can be based on:

      Parameters:
      name : str
         Match a all parameters with *name*. Can be exact in the form
         of *At1.At2:NAME*, where *At* are the atoms for a parameter,
         or fuzzy in the form of *NAME*, meaning that all paraemters
         with *NAME* (regardless of atoms) will be matched

      atoms : str
         A string of the form 'C.H', 'C.C.H', etc

      block : str or int
         Parameter block name (str) or index (int), matching
         *:attr:blocks*

      index : int
         The parameter at *index* within a block.

   trim(keepatoms: List[str] | None = None)

      Returns a new ReaxFFParameters instance, only keeping parameters
      describing interactions between atoms in the "keepatoms" list.

      keepatoms: list of str
         Elements to keep, e.g. ['C', 'H'] or ['C', 'H', '0']

      NOTE: to save wildcard parameters, you need to explicitly
      include '0' as one of the keepatoms.

      Example:

      >>> old = ReaxFFParameters('CHON2017_weak_bb.ff')
      >>> new_interf = old.trim(keepatoms=['C', 'H', 'O'])

      will keep the GEN block, ATM:C, ..., BND:C.H, ...., ANG:C.C.H,
      ...., etc.

   apply_default_constraints(only_active: bool = False, verbose: bool = False, parameter_names: List[str] | None = None)

      Changes the parameter values and ranges to agree with the
      default constraints.

      Example: Below the p_bo4 parameter is set to a value greater
      than or equal to 1.

      >>> ff = ReaxFFParameters('CHO.ff')
      >>> p = ff['H.H:p_bo4']
      >>> p.is_active = True
      >>> print(p.value, p.range, p.is_active)
      0.0 [-1.0, 1.0] True
      >>> ff.apply_default_constraints(only_active=True)
      >>> print(p.value, p.range, p.is_active)
      1.0 [1.0, 1.0] True

      The following values/ranges will be altered:

         ('ATM', 'alpha_ij') >= 0
         ('OFD', 'alpha_ij') >= 0
         ('ATM', 'D_ij') >= 0
         ('OFD', 'D_ij', >= 0
         ('ATM', 'p_boc3') >= 0
         ('ATM', 'p_boc4') >= 0
         ('ATM', 'p_boc5') >= 0
         ('ATM', 'r_0^sigma') >= 0
         ('OFD', 'r_0^sigma') >= 0
         ('ATM', 'r_vdW') >= 0
         ('ATM', 'gamma_i') >= 0
         ('ATM', 'Val_i^e') >= 0
         ('BND', 'D_e^pi') >= 0
         ('BND', 'D_e^pipi') >= 0
         ('BND', 'D_e^sigma') >= 0
         ('BND', "Val'_i^boc") >= 0
         ('BND', 'p_be2') >= 0
         ('BND', 'p_ovun1') >= 0
         ('ANG', 'Theta_0,0') >= 0
         ('ANG', 'p_val2') >= 0
         ('ANG', 'p_val7') >= 0
         ('HBD', 'r_hb^0') >= 0

         ('ATM', '1/gamma_w') > 0.5

         ('BND', 'p_bo6') >= 1
         ('BND', 'p_bo4') >= 1
         ('BND', 'p_bo2') >= 1
         ('ANG', 'p_val4') >= 1

         ('ANG', 'Theta_0,0') < 180

      only_active: bool
         If True, make changes only to the active parameter subset

      verbose: bool
         Prints all changed parameters and the reason for changing
         them to stdout

      parameter_names: list of str
         Names of the parameters to apply default constraints to. If
         not set, all parameters will be considered. Note: the value
         of "only_active" is still considered.

   get_default_constraints(only_active=True)

      Returns a list of constraints that can be used in the input
      Constraints block in ParAMS (if each entry is separated by a
      newline).

      only_active : bool
         Whether to only give constraints for the active subset.

   classmethod yaml_load(yamlfile, *a, **kw)

      Loads a parameter interface from a (compressed) human-readable
      *yamlfile*. Child classes should call super() in this method
      before implementing their own routines.


8.3.1.2.8. Full List of Parameter Names
=======================================

*General* block, no prefix
   1. "p_boc1"

   2. "p_boc2"

   3. "-p_coa2"

   4. "p_trip4"

   5. "p_trip3"

   6. "k_c2"

   7. "p_ovun6"

   8. "p_trip2"

   9. "p_ovun7"

   10. "p_ovun8"

   11. "p_trip1"

   12. "nonb_low,swa"

   13. "R_cut"

   14. "p_fe1"

   15. "p_val6"

   16. "p_lp1"

   17. "p_val9"

   18. "p_val10"

   19. "p_fe2"

   20. "p_pen2"

   21. "p_pen3"

   22. "p_pen4"

   23. "p_fe3"

   24. "p_tor2"

   25. "p_tor3"

   26. "p_tor4"

   27. "p_elho" **only if eReaxFF is enabled**

   28. "p_cot2"

   29. "p_vdW1"

   30. "cutoff*100"

   31. "p_coa4"

   32. "p_ovun4"

   33. "p_ovun3"

   34. "p_val8"

   35. "X_soft"

   36. "n/a 1"

   37. "p_val" **only if eReaxFF is enabled**

   38. "n/a 2"

   39. "p_coa3"

   40. "n/a 3"

   41. "n/a 4" **only if eReaxFF is enabled**

*Atoms* block, prefix: "At1:"
   1. "r_0^sigma"

   2. "Val_i"

   3. "n/a 1"

   4. "r_vdW"

   5. "D_ij"

   6. "gamma_i"

   7. "r_0^pi"

   8. "Val_i^e"

   9. "alpha_ij"

   10. "1/gamma_w"

   11. "Val_j^angle"

   12. "p_ovun5"

   13. "p_i^xel2" **only if eReaxFF is enabled**

   14. "chi_i"

   15. "eta_i"

   16. "n/a 2"

   17. "r_0^pipi"

   18. "p_lp2"

   19. "n/a 3"

   20. "p_boc4"

   21. "p_boc3"

   22. "p_boc5"

   23. "C_i"

   24. "alpha_i" **only if eReaxFF is enabled**

   25. "p_ovun2"

   26. "p_val3"

   27. "beta_i" **only if eReaxFF is enabled**

   28. "Val_i^'boc"

   29. "p_val5"

   30. "p_c1"

   31. "p_c2"

   32. "p_c3"

   33. "C_i" **only if Lg dispersion is enabled**

   34. "R_eij" **only if Lg dispersion is enabled**

*Bonds* block, prefix: "At1.At2:" (equivalent to "At2.At1:")
   1. "D_e^sigma"

   2. "D_e^pi"

   3. "D_e^pipi"

   4. "p_be1"

   5. "p_bo5"

   6. "Val'_i^boc"

   7. "p_bo6"

   8. "p_ovun1"

   9. "p_be2"

   10. "p_bo3"

   11. "p_bo4"

   12. "n/a"

   13. "p_bo1"

   14. "p_bo2"

   15. "delta'_i"

   16. "p_ij^xel1" **only if eReaxFF is enabled**

*Off-diagonals* block, prefix: "At1.At2:" (equivalent to "At2.At1:")
   1. "D_ij"

   2. "r_vdW"

   3. "alpha_ij"

   4. "r_0^sigma"

   5. "r_0^pi"

   6. "r_0^pipi"

   7. "C_i,C_lg,ij" **only if Lg dispersion is enabled**

*Angles* block, prefix "At1.At2.At3:" (equivalent to "At3.At2.At1",
At2 is the central atom)
   1. "Theta_0,0"

   2. "p_val1"

   3. "p_val2"

   4. "p_coa1"

   5. "p_val7"

   6. "p_pen1"

   7. "p_val4"

*Torsions* block, prefix "At1.At2.At3.At4:" (equivalent to
"At4.At3.At2.A1:")
   1. "V_1"

   2. "V_2"

   3. "V_3"

   4. "p_tor1"

   5. "p_cot1"

   6. "n/a 1"

   7. "n/a 2"

*Hydrogen Bonds* block, prefix "At1.H.At2:" (NOT equivalent to
"At2.H.At1:")
   1. "r_hb^0"

   2. "p_hb1"

   3. "-p_hb2"

   4. "-p_hb3"


8.3.1.2.9. Parameter Categories
===============================

To guide your selection we have split the ReaxFF parameters into three
categories:

1. **Standard**: Parameters which are generally safe to optimize.

2. **Expert**: Parameters which require some expert knowledge/insight.
   Generally, these are parameters which should be handled with more
   care and should only be tuned at the end of the parametrization
   process, if at all.

3. **DoNotOptimize**: Parameters which should never be activated and
   changed through optimization. This final category generally
   represents parameters which act as boolean flags or switches, are
   not used by the force field, or contain physical constants like
   atomic mass.

Please note that these categorizations are entirely empirical and
based only on in-house experience and expert knowledge. They are not
definitive groupings, and serve only to better inform and guide your
own decisions.

The category for each ReaxFF parameter is shown in a separate column
of the GUI. When scripting it can be accessed via the "metadata"
attribute of the "Parameter" class e.g. "p.metadata.category" or
"p.metadata['category']".
