Merge pull request #2211 from athomps/compute-mliap

Add compute style mliap to MLIAP package

doc/src/Commands_compute.rst

@@ -80,6 +80,7 @@ KOKKOS, o = USER-OMP, t = OPT.
    * :doc:`ke/eff <compute_ke_eff>`
    * :doc:`ke/rigid <compute_ke_rigid>`
    * :doc:`mesont <compute_mesont>`
+   * :doc:`mliap <compute_mliap>`
    * :doc:`momentum <compute_momentum>`
    * :doc:`msd <compute_msd>`
    * :doc:`msd/chunk <compute_msd_chunk>`

doc/src/compute.rst

@@ -225,6 +225,7 @@ The individual style names on the :doc:`Commands compute <Commands_compute>` doc
 * :doc:`ke/atom/eff <compute_ke_atom_eff>` - per-atom translational and radial kinetic energy in the electron force field model
 * :doc:`ke/eff <compute_ke_eff>` - kinetic energy of a group of nuclei and electrons in the electron force field model
 * :doc:`ke/rigid <compute_ke_rigid>` - translational kinetic energy of rigid bodies
+* :doc:`mliap <compute_mliap>` - gradients of energy and forces w.r.t. model parameters and related quantities for training machine learning interatomic potentials
 * :doc:`momentum <compute_momentum>` - translational momentum
 * :doc:`msd <compute_msd>` - mean-squared displacement of group of atoms
 * :doc:`msd/chunk <compute_msd_chunk>` - mean-squared displacement for each chunk
@@ -274,7 +275,7 @@ The individual style names on the :doc:`Commands compute <Commands_compute>` doc
 * :doc:`smd/ulsph/strain/rate <compute_smd_ulsph_strain_rate>` -
 * :doc:`smd/ulsph/stress <compute_smd_ulsph_stress>` - per-particle Cauchy stress tensor and von Mises equivalent stress in Smooth Mach Dynamics
 * :doc:`smd/vol <compute_smd_vol>` - per-particle volumes and their sum in Smooth Mach Dynamics
-* :doc:`snap <compute_sna_atom>` - bispectrum components and related quantities for a group of atoms
+* :doc:`snap <compute_sna_atom>` - gradients of SNAP energy and forces w.r.t. linear coefficients and related quantities for fitting SNAP potentials
 * :doc:`sna/atom <compute_sna_atom>` - bispectrum components for each atom
 * :doc:`snad/atom <compute_sna_atom>` - derivative of bispectrum components for each atom
 * :doc:`snav/atom <compute_sna_atom>` - virial contribution from bispectrum components for each atom

doc/src/compute_mliap.rst

@@ -0,0 +1,170 @@
+.. index:: compute mliap
+
+compute mliap command
+=====================
+
+Syntax
+""""""
+
+.. code-block:: LAMMPS
+
+   compute ID group-ID mliap ... keyword values ...
+
+* ID, group-ID are documented in :doc:`compute <compute>` command
+* mliap = style name of this compute command
+* two keyword/value pairs must be appended
+* keyword = *model* or *descriptor*
+
+  .. parsed-literal::
+
+       *model* values = style Nelems Nparams
+         style = *linear* or *quadratic*
+         Nelems = number of elements
+         Nparams = number of parameters per element 
+       *descriptor* values = style filename
+         style = *sna*
+         filename = name of file containing descriptor definitions
+
+Examples
+""""""""
+
+.. code-block:: LAMMPS
+
+   compute mliap model linear 2 31 descriptor sna Ta06A.mliap.descriptor
+
+Description
+"""""""""""
+
+Compute style *mliap* provides a general interface to the gradient
+of machine-learning interatomic potentials w.r.t. model parameters. 
+It is used primarily for calculating the gradient of energy, force, and
+stress components w.r.t. model parameters, which is useful when training
+:doc:`mliap pair_style <pair_mliap>` models to match target data.
+It provides separate 
+definitions of the interatomic potential functional form (*model*)
+and the geometric quantities that characterize the atomic positions
+(*descriptor*). By defining *model* and *descriptor* separately, 
+it is possible to use many different models with a given descriptor,
+or many different descriptors with a given model. Currently, the 
+compute supports just two models, *linear* and *quadratic*,
+and one descriptor, *sna*, the SNAP descriptor used by 
+:doc:`pair_style snap <pair_snap>`, including the linear, quadratic,
+and chem variants. Work is currently underway to extend
+the interface to handle neural network energy models,
+and it is also straightforward to add new descriptor styles.
+
+The compute *mliap* command must be followed by two keywords
+*model* and *descriptor* in either order. 
+
+The *model* keyword is followed by a model style, currently limited to
+either *linear* or *quadratic*. In both cases,
+this is followed by two arguments. *nelems* is the number of elements.
+It must be equal to the number of LAMMPS atom types. *nparams*
+is the number of parameters per element for this model i.e.
+the number of parameter gradients for each element. Note these definitions
+are identical to those of *nelems* and *nparams* in the 
+:doc:`pair_style mliap <pair_mliap>` model file.
+
+The *descriptor* keyword is followed by a descriptor style, and additional arguments.
+Currently the only descriptor style is *sna*, indicating the bispectrum component 
+descriptors used by the Spectral Neighbor Analysis Potential (SNAP) potentials of 
+:doc:`pair_style snap <pair_snap>`.
+The \'p\' in SNAP is dropped, because keywords that match pair_styles are silently stripped 
+out by the LAMMPS command parser. A single additional argument specifies the descriptor filename 
+containing the parameters and setting used by the SNAP descriptor. 
+The descriptor filename usually ends in the *.mliap.descriptor* extension.
+The format of this file is identical to the descriptor file in the 
+:doc:`pair_style mliap <pair_mliap>`, and is described in detail
+there. 
+
+.. note::
+
+   The number of LAMMPS atom types (and the value of *nelems* in the model)
+   must match the value of *nelems* in the descriptor file. 
+
+Compute *mliap* calculates a global array containing gradient information.
+The number of columns in the array is :math:`nelems \times nparams + 1`.
+The first row of the array contain the derivative of potential energy w.r.t. to
+each parameter and each element. The last six rows
+of the array contain the corresponding derivatives of the
+virial stress tensor, listed in Voigt notation: *pxx*, *pyy*, *pzz*,
+*pyz*, *pxz*, *pxy*. In between the energy and stress rows are
+the 3\*\ *N* rows containing the derivatives of the force components. 
+See section below on output for a detailed description of how 
+rows and columns are ordered. 
+
+The element in the last column of each row contains
+the potential energy, force, or stress, according to the row.
+These quantities correspond to the user-specified reference potential
+that must be subtracted from the target data when training a model.
+The potential energy calculation uses the built in compute *thermo_pe*.
+The stress calculation uses a compute called *mliap_press* that is
+automatically created behind the scenes, according to the following
+command:
+
+.. code-block:: LAMMPS
+
+   compute mliap_press all pressure NULL virial
+
+See section below on output for a detailed explanation of the data
+layout in the global array.
+
+Atoms not in the group do not contribute to this compute. 
+Neighbor atoms not in the group do not contribute to this compute.
+The neighbor list needed to compute this quantity is constructed each
+time the calculation is performed (i.e. each time a snapshot of atoms
+is dumped).  Thus it can be inefficient to compute/dump this quantity
+too frequently.
+
+.. note::
+
+   If the user-specified reference potentials includes bonded and
+   non-bonded pairwise interactions, then the settings of
+   :doc:`special_bonds <special_bonds>` command can remove pairwise 
+   interactions between atoms in the same bond, angle, or dihedral.  This
+   is the default setting for the :doc:`special_bonds <special_bonds>`
+   command, and means those pairwise interactions do not appear in the
+   neighbor list.  Because this fix uses the neighbor list, it also means
+   those pairs will not be included in the calculation. The :doc:`rerun <rerun>`
+   command is not an option here, since the reference potential is required
+   for the last column of the global array. A work-around is to prevent
+   pairwise interactions from being removed by explicitly adding a 
+   *tiny* positive value for every pairwise interaction that would otherwise be
+   set to zero in the :doc:`special_bonds <special_bonds>` command.
+
+----------
+
+**Output info:**
+
+Compute *mliap* evaluates a global array.
+The columns are arranged into
+*nelems* blocks, listed in order of element *I*\ . Each block
+contains one column for each of the *nparams* model parameters. 
+A final column contains the corresponding energy, force component
+on an atom, or virial stress component. The rows of the array appear
+in the following order:
+
+* 1 row: Derivatives of potential energy w.r.t. each parameter of each element.
+* 3\*\ *N* rows: Derivatives of force components. x, y, and z components of force on atom *i* appearing in consecutive rows. The atoms are sorted based on atom ID.
+* 6 rows: Derivatives of virial stress tensor  w.r.t. each parameter of each element. The ordering of the rows follows Voigt notation: *pxx*, *pyy*, *pzz*, *pyz*, *pxz*, *pxy*.
+
+These values can be accessed by any command that uses a global array
+from a compute as input.  See the :doc:`Howto output <Howto_output>` doc
+page for an overview of LAMMPS output options. To see how this command
+can be used within a Python workflow to train machine-learning interatomic
+potentials, see the examples in `FitSNAP <https://github.com/FitSNAP/FitSNAP>`_.
+
+Restrictions
+""""""""""""
+
+This compute is part of the MLIAP package.  It is only enabled if
+LAMMPS was built with that package.  In addition, building LAMMPS with the MLIAP package
+requires building LAMMPS with the SNAP package.
+See the :doc:`Build package <Build_package>` doc page for more info.
+
+Related commands
+""""""""""""""""
+
+:doc:`pair_style mliap <pair_mliap>`
+
+**Default:** none

doc/src/compute_sna_atom.rst

@@ -391,7 +391,9 @@ of :math:`K N_{elem}^3` columns.
 
 These values can be accessed by any command that uses per-atom values
 from a compute as input.  See the :doc:`Howto output <Howto_output>` doc
-page for an overview of LAMMPS output options.
+page for an overview of LAMMPS output options. To see how this command
+can be used within a Python workflow to train SNAP potentials, 
+see the examples in `FitSNAP <https://github.com/FitSNAP/FitSNAP>`_.
 
 Restrictions
 """"""""""""

doc/src/pair_mliap.rst

@@ -8,7 +8,19 @@ Syntax
 
 .. code-block:: LAMMPS
 
-   pair_style mliap
+   pair_style mliap ... keyword values ...
+
+* two keyword/value pairs must be appended
+* keyword = *model* or *descriptor*
+
+  .. parsed-literal::
+
+       *model* values = style filename
+         style = *linear* or *quadratic*
+         filename = name of file containing model definitions
+       *descriptor* values = style filename
+         style = *sna*
+         filename = name of file containing descriptor definitions
 
 Examples
 """"""""
@@ -22,8 +34,8 @@ Examples
 Description
 """""""""""
 
-Pair style *mliap* provides a general interface to families of
-machine-learning interatomic potentials. It provides separate
+Pair style *mliap* provides a general interface to families of 
+machine-learning interatomic potentials. It allows separate 
 definitions of the interatomic potential functional form (*model*)
 and the geometric quantities that characterize the atomic positions
 (*descriptor*). By defining *model* and *descriptor* separately,
@@ -34,6 +46,9 @@ and one descriptor, *sna*, the SNAP descriptor used by :doc:`pair_style snap <pa
 and chem variants. Work is currently underway to extend
 the interface to handle neural network energy models,
 and it is also straightforward to add new descriptor styles.
+In order to train a model, it is useful to know the gradient or derivative
+of energy, force, and stress w.r.t. model parameters. This information
+can be accessed using the related :doc:`compute mliap <compute_mliap>` command.
 
 The pair_style *mliap* command must be followed by two keywords
 *model* and *descriptor* in either order. A single
@@ -46,10 +61,10 @@ where N is the number of LAMMPS atom types.
 
 The *model* keyword is followed by a model style, currently limited to
 either *linear* or *quadratic*. In both cases,
-this is followed by a single argument specifying the model filename containing the
-linear or quadratic coefficients for a set of elements.
+this is followed by a single argument specifying the model filename containing the 
+parameters for a set of elements. 
 The model filename usually ends in the *.mliap.model* extension.
-It may contain coefficients for many elements. The only requirement is that it
+It may contain parameters for many elements. The only requirement is that it
 contain at least those element names appearing in the
 *pair_coeff* command.
 
@@ -58,10 +73,10 @@ but follows a strict format after that. The first non-blank non-comment
 line must contain two integers:
 
 * nelems  = Number of elements
-* ncoeff = Number of coefficients
+* nparams = Number of parameters
 
 This is followed by one block for each of the *nelem* elements.
-Each block consists of *ncoeff* coefficients, one per line.
+Each block consists of *nparams* parameters, one per line.
 Note that this format is similar, but not identical to that used
 for the :doc:`pair_style snap <pair_snap>` coefficient file.
 Specifically, the line containing the element weight and radius is omitted,
@@ -131,6 +146,6 @@ See the :doc:`Build package <Build_package>` doc page for more info.
 Related commands
 """"""""""""""""
 
-:doc:`pair_style snap  <pair_snap>`,
+:doc:`pair_style snap  <pair_snap>`, :doc:`compute mliap <compute_mliap>`
 
 **Default:** none

doc/utils/sphinx-config/false_positives.txt

@@ -2083,6 +2083,7 @@ Novint
 np
 Npair
 Npairs
+nparams
 nparticle
 npernode
 nph

examples/mliap/compute.mliap.descriptor

@@ -0,0 +1,19 @@
+# LAMMPS SNAP parameters for MLIAP compute example
+
+# required
+rcutfac 1.0
+twojmax 2
+
+# elements
+
+nelems 2
+elems A B 
+radelems 2.3 2.0
+welems 1.0 0.96
+
+# optional
+
+rfac0 0.99363
+rmin0 0
+bzeroflag 0
+switchflag 0