Merge pull request #3079 from akohlmey/next_patch_release

Step version strings for the next patch release

cmake/CMakeLists.txt

@@ -339,7 +339,6 @@ pkg_depends(ML-IAP ML-SNAP)
 pkg_depends(MPIIO MPI)
 pkg_depends(ATC MANYBODY)
 pkg_depends(LATBOLTZ MPI)
-pkg_depends(PHONON KSPACE)
 pkg_depends(SCAFACOS MPI)
 pkg_depends(DIELECTRIC KSPACE)
 pkg_depends(DIELECTRIC EXTRA-PAIR)
@@ -609,7 +608,7 @@ endif()
 # packages which selectively include variants based on enabled styles
 # e.g. accelerator packages
 ######################################################################
-foreach(PKG_WITH_INCL CORESHELL QEQ OPENMP DPD-SMOOTH KOKKOS OPT INTEL GPU)
+foreach(PKG_WITH_INCL CORESHELL DPD-SMOOTH PHONON QEQ OPENMP KOKKOS OPT INTEL GPU)
   if(PKG_${PKG_WITH_INCL})
     include(Packages/${PKG_WITH_INCL})
   endif()

cmake/Modules/Packages/GPU.cmake

@@ -306,12 +306,12 @@ elseif(GPU_API STREQUAL "HIP")
 
         if(HIP_COMPILER STREQUAL "clang")
           add_custom_command(OUTPUT ${CUBIN_FILE}
-            VERBATIM COMMAND ${HIP_HIPCC_EXECUTABLE} --genco --offload-arch=${HIP_ARCH} -O3 -ffast-math -DUSE_HIP -D_${GPU_PREC_SETTING} -DLAMMPS_${LAMMPS_SIZES} -I${LAMMPS_LIB_SOURCE_DIR}/gpu -o ${CUBIN_FILE} ${CU_CPP_FILE}
+            VERBATIM COMMAND ${HIP_HIPCC_EXECUTABLE} --genco --offload-arch=${HIP_ARCH} -O3 -DUSE_HIP -D_${GPU_PREC_SETTING} -DLAMMPS_${LAMMPS_SIZES} -I${LAMMPS_LIB_SOURCE_DIR}/gpu -o ${CUBIN_FILE} ${CU_CPP_FILE}
             DEPENDS ${CU_CPP_FILE}
             COMMENT "Generating ${CU_NAME}.cubin")
         else()
           add_custom_command(OUTPUT ${CUBIN_FILE}
-            VERBATIM COMMAND ${HIP_HIPCC_EXECUTABLE} --genco -t="${HIP_ARCH}" -f=\"-O3 -ffast-math -DUSE_HIP -D_${GPU_PREC_SETTING} -DLAMMPS_${LAMMPS_SIZES} -I${LAMMPS_LIB_SOURCE_DIR}/gpu\" -o ${CUBIN_FILE} ${CU_CPP_FILE}
+            VERBATIM COMMAND ${HIP_HIPCC_EXECUTABLE} --genco -t="${HIP_ARCH}" -f=\"-O3 -DUSE_HIP -D_${GPU_PREC_SETTING} -DLAMMPS_${LAMMPS_SIZES} -I${LAMMPS_LIB_SOURCE_DIR}/gpu\" -o ${CUBIN_FILE} ${CU_CPP_FILE}
             DEPENDS ${CU_CPP_FILE}
             COMMENT "Generating ${CU_NAME}.cubin")
         endif()

cmake/Modules/Packages/PHONON.cmake

@@ -0,0 +1,9 @@
+# fix phonon may only be installed if also the FFT wrappers from KSPACE are installed
+if(NOT PKG_KSPACE)
+  get_property(LAMMPS_FIX_HEADERS GLOBAL PROPERTY FIX)
+  list(REMOVE_ITEM LAMMPS_FIX_HEADERS ${LAMMPS_SOURCE_DIR}/PHONON/fix_phonon.h)
+  set_property(GLOBAL PROPERTY FIX "${LAMMPS_FIX_HEADERS}")
+  get_target_property(LAMMPS_SOURCES lammps SOURCES)
+  list(REMOVE_ITEM LAMMPS_SOURCES ${LAMMPS_SOURCE_DIR}/PHONON/fix_phonon.cpp)
+  set_property(TARGET lammps PROPERTY SOURCES "${LAMMPS_SOURCES}")
+endif()

doc/lammps.1

@@ -1,4 +1,4 @@
-.TH LAMMPS "1" "14 December 2021" "2021-12-14"
+.TH LAMMPS "1" "7 January 2022" "2022-1-7"
 .SH NAME
 .B LAMMPS
 \- Molecular Dynamics Simulator.

doc/src/Bibliography.rst

@@ -1123,9 +1123,12 @@ Bibliography
 **(Sun)**
    Sun, J. Phys. Chem. B, 102, 7338-7364 (1998).
 
-**(Surblys)**
+**(Surblys2019)**
    Surblys, Matsubara, Kikugawa, Ohara, Phys Rev E, 99, 051301(R) (2019).
 
+**(Surblys2021)**
+   Surblys, Matsubara, Kikugawa, Ohara, J Appl Phys 130, 215104 (2021).
+
 **(Sutmann)**
    Sutmann, Arnold, Fahrenberger, et. al., Physical review / E 88(6), 063308 (2013)
 

doc/src/Commands_compute.rst

@@ -28,6 +28,7 @@ KOKKOS, o = OPENMP, t = OPT.
    * :doc:`angle <compute_angle>`
    * :doc:`angle/local <compute_angle_local>`
    * :doc:`angmom/chunk <compute_angmom_chunk>`
+   * :doc:`ave/sphere/atom (k) <compute_ave_sphere_atom>`
    * :doc:`basal/atom <compute_basal_atom>`
    * :doc:`body/local <compute_body_local>`
    * :doc:`bond <compute_bond>`

doc/src/Developer_utils.rst

@@ -56,11 +56,11 @@ String to number conversions with validity check
 
 These functions should be used to convert strings to numbers. They are
 are strongly preferred over C library calls like ``atoi()`` or
-``atof()`` since they check if the **entire** provided string is a valid
+``atof()`` since they check if the **entire** string is a valid
 (floating-point or integer) number, and will error out instead of
 silently returning the result of a partial conversion or zero in cases
-where the string is not a valid number.  This behavior allows to more
-easily detect typos or issues when processing input files.
+where the string is not a valid number.  This behavior improves
+detecting typos or issues when processing input files.
 
 Similarly the :cpp:func:`logical() <LAMMPS_NS::utils::logical>` function
 will convert a string into a boolean and will only accept certain words.
@@ -76,19 +76,34 @@ strings for compliance without conversion.
 
 ----------
 
-.. doxygenfunction:: numeric
+.. doxygenfunction:: numeric(const char *file, int line, const std::string &str, bool do_abort, LAMMPS *lmp)
    :project: progguide
 
-.. doxygenfunction:: inumeric
+.. doxygenfunction:: numeric(const char *file, int line, const char *str, bool do_abort, LAMMPS *lmp)
    :project: progguide
 
-.. doxygenfunction:: bnumeric
+.. doxygenfunction:: inumeric(const char *file, int line, const std::string &str, bool do_abort, LAMMPS *lmp)
    :project: progguide
 
-.. doxygenfunction:: tnumeric
+.. doxygenfunction:: inumeric(const char *file, int line, const char *str, bool do_abort, LAMMPS *lmp)
    :project: progguide
 
-.. doxygenfunction:: logical
+.. doxygenfunction:: bnumeric(const char *file, int line, const std::string &str, bool do_abort, LAMMPS *lmp)
+   :project: progguide
+
+.. doxygenfunction:: bnumeric(const char *file, int line, const char *str, bool do_abort, LAMMPS *lmp)
+   :project: progguide
+
+.. doxygenfunction:: tnumeric(const char *file, int line, const std::string &str, bool do_abort, LAMMPS *lmp)
+   :project: progguide
+
+.. doxygenfunction:: tnumeric(const char *file, int line, const char *str, bool do_abort, LAMMPS *lmp)
+   :project: progguide
+
+.. doxygenfunction:: logical(const char *file, int line, const std::string &str, bool do_abort, LAMMPS *lmp)
+   :project: progguide
+
+.. doxygenfunction:: logical(const char *file, int line, const char *str, bool do_abort, LAMMPS *lmp)
    :project: progguide
 
 

doc/src/Developer_write.rst

@@ -55,7 +55,7 @@ of each timestep. First of all, implement a constructor:
      if (narg < 4)
        error->all(FLERR,"Illegal fix print/vel command");
 
-     nevery = force->inumeric(FLERR,arg[3]);
+     nevery = utils::inumeric(FLERR,arg[3],false,lmp);
      if (nevery <= 0)
        error->all(FLERR,"Illegal fix print/vel command");
    }

doc/src/Errors_messages.rst

@@ -7772,9 +7772,6 @@ keyword to allow for additional bonds to be formed
    The system size must fit in a 32-bit integer to use this dump
    style.
 
-*Too many atoms to dump sort*
-   Cannot sort when running with more than 2\^31 atoms.
-
 *Too many elements extracted from MEAM library.*
    Increase 'maxelt' in meam.h and recompile.
 

doc/src/Modify_pair.rst

@@ -12,24 +12,24 @@ includes some optional methods to enable its use with rRESPA.
 
 Here is a brief description of the class methods in pair.h:
 
-+---------------------------------+-------------------------------------------------------------------+
-| compute                         | workhorse routine that computes pairwise interactions             |
-+---------------------------------+-------------------------------------------------------------------+
-| settings                        | reads the input script line with arguments you define             |
-+---------------------------------+-------------------------------------------------------------------+
-| coeff                           | set coefficients for one i,j type pair                            |
-+---------------------------------+-------------------------------------------------------------------+
-| init_one                        | perform initialization for one i,j type pair                      |
-+---------------------------------+-------------------------------------------------------------------+
-| init_style                      | initialization specific to this pair style                        |
-+---------------------------------+-------------------------------------------------------------------+
-| write & read_restart            | write/read i,j pair coeffs to restart files                       |
-+---------------------------------+-------------------------------------------------------------------+
-| write & read_restart_settings   | write/read global settings to restart files                       |
-+---------------------------------+-------------------------------------------------------------------+
-| single                          | force and energy of a single pairwise interaction between 2 atoms |
-+---------------------------------+-------------------------------------------------------------------+
-| compute_inner/middle/outer      | versions of compute used by rRESPA                                |
-+---------------------------------+-------------------------------------------------------------------+
++---------------------------------+---------------------------------------------------------------------+
+| compute                         | workhorse routine that computes pairwise interactions               |
++---------------------------------+---------------------------------------------------------------------+
+| settings                        | reads the input script line with arguments you define               |
++---------------------------------+---------------------------------------------------------------------+
+| coeff                           | set coefficients for one i,j type pair                              |
++---------------------------------+---------------------------------------------------------------------+
+| init_one                        | perform initialization for one i,j type pair                        |
++---------------------------------+---------------------------------------------------------------------+
+| init_style                      | initialization specific to this pair style                          |
++---------------------------------+---------------------------------------------------------------------+
+| write & read_restart            | write/read i,j pair coeffs to restart files                         |
++---------------------------------+---------------------------------------------------------------------+
+| write & read_restart_settings   | write/read global settings to restart files                         |
++---------------------------------+---------------------------------------------------------------------+
+| single                          | force/r and energy of a single pairwise interaction between 2 atoms |
++---------------------------------+---------------------------------------------------------------------+
+| compute_inner/middle/outer      | versions of compute used by rRESPA                                  |
++---------------------------------+---------------------------------------------------------------------+
 
 The inner/middle/outer routines are optional.

doc/src/Packages_details.rst

@@ -1880,6 +1880,12 @@ MPIIO library.  It adds :doc:`dump styles <dump>` with a "mpiio" in
 their style name.  Restart files with an ".mpiio" suffix are also
 written and read in parallel.
 
+.. warning::
+
+   The MPIIO package is currently unmaintained and has become
+   unreliable. Use with caution.
+
+
 **Install:**
 
 The MPIIO package requires that LAMMPS is build in :ref:`MPI parallel mode <serial>`.

doc/src/angle_class2.rst

@@ -64,34 +64,44 @@ These are the 4 coefficients for the :math:`E_a` formula:
 radians internally; hence the various :math:`K` are effectively energy
 per radian\^2 or radian\^3 or radian\^4.
 
-For the :math:`E_{bb}` formula, each line in a :doc:`angle_coeff <angle_coeff>`
-command in the input script lists 4 coefficients, the first of which
-is "bb" to indicate they are BondBond coefficients.  In a data file,
-these coefficients should be listed under a "BondBond Coeffs" heading
-and you must leave out the "bb", i.e. only list 3 coefficients after
-the angle type.
+For the :math:`E_{bb}` formula, each line in a :doc:`angle_coeff
+<angle_coeff>` command in the input script lists 4 coefficients, the
+first of which is "bb" to indicate they are BondBond coefficients.  In
+a data file, these coefficients should be listed under a "BondBond
+Coeffs" heading and you must leave out the "bb", i.e. only list 3
+coefficients after the angle type.
 
 * bb
 * :math:`M` (energy/distance\^2)
 * :math:`r_1` (distance)
 * :math:`r_2` (distance)
 
-For the :math:`E_{ba}` formula, each line in a :doc:`angle_coeff <angle_coeff>`
-command in the input script lists 5 coefficients, the first of which
-is "ba" to indicate they are BondAngle coefficients.  In a data file,
-these coefficients should be listed under a "BondAngle Coeffs" heading
-and you must leave out the "ba", i.e. only list 4 coefficients after
-the angle type.
+For the :math:`E_{ba}` formula, each line in a :doc:`angle_coeff
+<angle_coeff>` command in the input script lists 5 coefficients, the
+first of which is "ba" to indicate they are BondAngle coefficients.
+In a data file, these coefficients should be listed under a "BondAngle
+Coeffs" heading and you must leave out the "ba", i.e. only list 4
+coefficients after the angle type.
 
 * ba
-* :math:`N_1` (energy/distance\^2)
-* :math:`N_2` (energy/distance\^2)
+* :math:`N_1` (energy/distance)
+* :math:`N_2` (energy/distance)
 * :math:`r_1` (distance)
 * :math:`r_2` (distance)
 
 The :math:`\theta_0` value in the :math:`E_{ba}` formula is not specified,
 since it is the same value from the :math:`E_a` formula.
 
+.. note::
+
+   It is important that the order of the I,J,K atoms in each angle
+   listed in the Angles section of the data file read by the
+   :doc:`read_data <read_data>` command be consistent with the order
+   of the :math:`r_1` and :math:`r_2` BondBond and BondAngle
+   coefficients.  This is because the terms in the formulas for
+   :math:`E_{bb}` and :math:`E_{ba}` will use the I,J atoms to compute
+   :math:`r_{ij}` and the J,K atoms to compute :math:`r_{jk}`.
+
 ----------
 
 .. include:: accel_styles.rst

doc/src/compute.rst

@@ -174,6 +174,7 @@ The individual style names on the :doc:`Commands compute <Commands_compute>` pag
 * :doc:`angle <compute_angle>` - energy of each angle sub-style
 * :doc:`angle/local <compute_angle_local>` - theta and energy of each angle
 * :doc:`angmom/chunk <compute_angmom_chunk>` - angular momentum for each chunk
+* :doc:`ave/sphere/atom <compute_ave_sphere_atom>` - compute local density and temperature around each atom
 * :doc:`basal/atom <compute_basal_atom>` - calculates the hexagonal close-packed "c" lattice vector of each atom
 * :doc:`body/local <compute_body_local>` - attributes of body sub-particles
 * :doc:`bond <compute_bond>` - energy of each bond sub-style

doc/src/compute_ave_sphere_atom.rst

@@ -0,0 +1,101 @@
+.. index:: compute ave/sphere/atom
+.. index:: compute ave/sphere/atom/kk
+
+compute ave/sphere/atom command
+================================
+
+Accelerator Variants: *ave/sphere/atom/kk*
+
+Syntax
+""""""
+
+.. parsed-literal::
+
+   compute ID group-ID ave/sphere/atom keyword values ...
+
+* ID, group-ID are documented in :doc:`compute <compute>` command
+* ave/sphere/atom = style name of this compute command
+* one or more keyword/value pairs may be appended
+
+  .. parsed-literal::
+
+     keyword = *cutoff*
+       *cutoff* value = distance cutoff
+
+Examples
+""""""""
+
+.. code-block:: LAMMPS
+
+   compute 1 all ave/sphere/atom
+
+   compute 1 all ave/sphere/atom cutoff 5.0
+   comm_modify cutoff 5.0
+
+Description
+"""""""""""
+
+Define a computation that calculates the local density and temperature
+for each atom and neighbors inside a spherical cutoff.
+
+The optional keyword *cutoff* defines the distance cutoff
+used when searching for neighbors. The default value is the cutoff
+specified by the pair style. If no pair style is defined, then a cutoff
+must be defined using this keyword. If the specified cutoff is larger than
+that of the pair_style plus neighbor skin (or no pair style is defined),
+the *comm_modify cutoff* option must also be set to match that of the
+*cutoff* keyword.
+
+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 you have a bonded system, 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 order parameter.  This
+   difficulty can be circumvented by writing a dump file, and using the
+   :doc:`rerun <rerun>` command to compute the order parameter for
+   snapshots in the dump file.  The rerun script can use a
+   :doc:`special_bonds <special_bonds>` command that includes all pairs in
+   the neighbor list.
+
+----------
+
+
+.. include:: accel_styles.rst
+
+
+----------
+
+Output info
+"""""""""""
+
+This compute calculates a per-atom array with two columns: density and temperature.
+
+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.
+
+Restrictions
+""""""""""""
+
+This compute is part of the EXTRA-COMPUTE package.  It is only enabled if
+LAMMPS was built with that package.  See the :doc:`Build package <Build_package>` page for more info.
+
+Related commands
+""""""""""""""""
+
+:doc:`comm_modify <comm_modify>`
+
+Default
+"""""""
+
+The option defaults are *cutoff* = pair style cutoff
+

doc/src/compute_heat_flux.rst

@@ -89,13 +89,20 @@ included in the calculation.
 .. warning::
 
    The compute *heat/flux* has been reported to produce unphysical
-   values for angle, dihedral and improper contributions
+   values for angle, dihedral, improper and constraint force contributions
    when used with :doc:`compute stress/atom <compute_stress_atom>`,
-   as discussed in :ref:`(Surblys) <Surblys2>` and :ref:`(Boone) <Boone>`.
-   You are strongly advised to
+   as discussed in :ref:`(Surblys2019) <Surblys3>`, :ref:`(Boone) <Boone>`
+   and :ref:`(Surblys2021) <Surblys4>`. You are strongly advised to
    use :doc:`compute centroid/stress/atom <compute_stress_atom>`,
    which has been implemented specifically for such cases.
 
+.. warning::
+
+   Due to an implementation detail, the :math:`y` and :math:`z`
+   components of heat flux from :doc:`fix rigid <fix_rigid>`
+   contribution when computed via :doc:`compute stress/atom <compute_stress_atom>`
+   are highly unphysical and should not be used.
+
 The Green-Kubo formulas relate the ensemble average of the
 auto-correlation of the heat flux :math:`\mathbf{J}`
 to the thermal conductivity :math:`\kappa`:
@@ -232,10 +239,14 @@ none
 
 ----------
 
-.. _Surblys2:
+.. _Surblys3:
 
-**(Surblys)** Surblys, Matsubara, Kikugawa, Ohara, Phys Rev E, 99, 051301(R) (2019).
+**(Surblys2019)** Surblys, Matsubara, Kikugawa, Ohara, Phys Rev E, 99, 051301(R) (2019).
 
 .. _Boone:
 
 **(Boone)** Boone, Babaei, Wilmer, J Chem Theory Comput, 15, 5579--5587 (2019).
+
+.. _Surblys4:
+
+**(Surblys2021)** Surblys, Matsubara, Kikugawa, Ohara, J Appl Phys 130, 215104 (2021).

doc/src/compute_stress_atom.rst

@@ -87,6 +87,10 @@ Tersoff 3-body interaction) is assigned in equal portions to each atom
 in the set.  E.g. 1/4 of the dihedral virial to each of the 4 atoms,
 or 1/3 of the fix virial due to SHAKE constraints applied to atoms in
 a water molecule via the :doc:`fix shake <fix_shake>` command.
+As an exception, the virial contribution from
+constraint forces in :doc:`fix rigid <fix_rigid>` on each atom
+is computed from the constraint force acting on the corresponding atom
+and its position, i.e. the total virial is not equally distributed.
 
 In case of compute *centroid/stress/atom*, the virial contribution is:
 
@@ -103,13 +107,25 @@ atom :math:`I` due to the interaction and the relative position
 :math:`\mathbf{r}_{I0}` of the atom :math:`I` to the geometric center
 of the interacting atoms, i.e. centroid, is used.  As the geometric
 center is different for each interaction, the :math:`\mathbf{r}_{I0}`
-also differs.  The sixth and seventh terms, Kspace and :doc:`fix
-<fix>` contribution respectively, are computed identical to compute
-*stress/atom*.  Although the total system virial is the same as
+also differs. The sixth term, Kspace contribution,
+is computed identically to compute *stress/atom*.
+The seventh term is handed differently depending on
+if the constraint forces are due to :doc:`fix shake <fix_shake>`
+or :doc:`fix rigid <fix_rigid>`.
+In case of SHAKE constraints, each distance constraint is
+handed as a pairwise interaction.
+E.g. in case of a water molecule, two OH and one HH distance
+constraints are treated as three pairwise interactions.
+In case of :doc:`fix rigid <fix_rigid>`,
+all constraint forces in the molecule are treated
+as a single many-body interaction with a single centroid position.
+In case of water molecule, the formula expression would become
+identical to that of the three-body angle interaction.
+Although the total system virial is the same as
 compute *stress/atom*, compute *centroid/stress/atom* is know to
-result in more consistent heat flux values for angle, dihedrals and
-improper contributions when computed via :doc:`compute heat/flux
-<compute_heat_flux>`.
+result in more consistent heat flux values for angle, dihedrals,
+improper and constraint force contributions
+when computed via :doc:`compute heat/flux <compute_heat_flux>`.
 
 If no extra keywords are listed, the kinetic contribution all of the
 virial contribution terms are included in the per-atom stress tensor.
@@ -134,7 +150,8 @@ contribution for the cluster interaction is divided evenly among those
 atoms.
 
 Details of how compute *centroid/stress/atom* obtains the virial for
-individual atoms is given in :ref:`(Surblys) <Surblys1>`, where the
+individual atoms are given in :ref:`(Surblys2019) <Surblys1>` and
+:ref:`(Surblys2021) <Surblys2>`, where the
 idea is that the virial of the atom :math:`I` is the result of only
 the force :math:`\mathbf{F}_I` on the atom due to the interaction and
 its positional vector :math:`\mathbf{r}_{I0}`, relative to the
@@ -235,10 +252,10 @@ between the pair of particles.  All bond styles are supported.  All
 angle, dihedral, improper styles are supported with the exception of
 INTEL and KOKKOS variants of specific styles.  It also does not
 support models with long-range Coulombic or dispersion forces,
-i.e. the kspace_style command in LAMMPS.  It also does not support the
-following fixes which add rigid-body constraints: :doc:`fix shake
-<fix_shake>`, :doc:`fix rattle <fix_shake>`, :doc:`fix rigid
-<fix_rigid>`, :doc:`fix rigid/small <fix_rigid>`.
+i.e. the kspace_style command in LAMMPS.  It also does not implement the
+following fixes which add rigid-body constraints:
+:doc:`fix rigid/* <fix_rigid>` and the OpenMP accelerated version of :doc:`fix rigid/small <fix_rigid>`,
+while all other :doc:`fix rigid/*/small <fix_rigid>` are implemented.
 
 LAMMPS will generate an error if one of these options is included in
 your model.  Extension of centroid stress calculations to these force
@@ -270,4 +287,8 @@ none
 
 .. _Surblys1:
 
-**(Surblys)** Surblys, Matsubara, Kikugawa, Ohara, Phys Rev E, 99, 051301(R) (2019).
+**(Surblys2019)** Surblys, Matsubara, Kikugawa, Ohara, Phys Rev E, 99, 051301(R) (2019).
+
+.. _Surblys2:
+
+**(Surblys2021)** Surblys, Matsubara, Kikugawa, Ohara, J Appl Phys 130, 215104 (2021).

doc/src/dump.rst

@@ -137,7 +137,7 @@ Examples
    dump myDump all atom/gz 100 dump.atom.gz
    dump myDump all atom/zstd 100 dump.atom.zst
    dump 2 subgroup atom 50 dump.run.bin
-   dump 2 subgroup atom 50 dump.run.mpiio.bin
+   dump 2 subgroup atom/mpiio 50 dump.run.mpiio.bin
    dump 4a all custom 100 dump.myforce.* id type x y vx fx
    dump 4b flow custom 100 dump.%.myforce id type c_myF[3] v_ke
    dump 4b flow custom 100 dump.%.myforce id type c_myF[*] v_ke
@@ -169,11 +169,12 @@ or multiple smaller files).
 
 .. note::
 
-   Because periodic boundary conditions are enforced only on
-   timesteps when neighbor lists are rebuilt, the coordinates of an atom
-   written to a dump file may be slightly outside the simulation box.
-   Re-neighbor timesteps will not typically coincide with the timesteps
-   dump snapshots are written.  See the :doc:`dump_modify pbc <dump_modify>` command if you with to force coordinates to be
+   Because periodic boundary conditions are enforced only on timesteps
+   when neighbor lists are rebuilt, the coordinates of an atom written
+   to a dump file may be slightly outside the simulation box.
+   Re-neighbor timesteps will not typically coincide with the
+   timesteps dump snapshots are written.  See the :doc:`dump_modify
+   pbc <dump_modify>` command if you with to force coordinates to be
    strictly inside the simulation box.
 
 .. note::
@@ -189,20 +190,21 @@ or multiple smaller files).
    multiple processors, each of which owns a subset of the atoms.
 
 For the *atom*, *custom*, *cfg*, and *local* styles, sorting is off by
-default.  For the *dcd*, *xtc*, *xyz*, and *molfile* styles, sorting by
-atom ID is on by default. See the :doc:`dump_modify <dump_modify>` doc
-page for details.
+default.  For the *dcd*, *xtc*, *xyz*, and *molfile* styles, sorting
+by atom ID is on by default. See the :doc:`dump_modify <dump_modify>`
+doc page for details.
 
-The *atom/gz*, *cfg/gz*, *custom/gz*, *local/gz*, and *xyz/gz* styles are identical
-in command syntax to the corresponding styles without "gz", however,
-they generate compressed files using the zlib library. Thus the filename
-suffix ".gz" is mandatory. This is an alternative approach to writing
-compressed files via a pipe, as done by the regular dump styles, which
-may be required on clusters where the interface to the high-speed network
-disallows using the fork() library call (which is needed for a pipe).
-For the remainder of this doc page, you should thus consider the *atom*
-and *atom/gz* styles (etc) to be inter-changeable, with the exception
-of the required filename suffix.
+The *atom/gz*, *cfg/gz*, *custom/gz*, *local/gz*, and *xyz/gz* styles
+are identical in command syntax to the corresponding styles without
+"gz", however, they generate compressed files using the zlib
+library. Thus the filename suffix ".gz" is mandatory. This is an
+alternative approach to writing compressed files via a pipe, as done
+by the regular dump styles, which may be required on clusters where
+the interface to the high-speed network disallows using the fork()
+library call (which is needed for a pipe).  For the remainder of this
+doc page, you should thus consider the *atom* and *atom/gz* styles
+(etc) to be inter-changeable, with the exception of the required
+filename suffix.
 
 Similarly, the *atom/zstd*, *cfg/zstd*, *custom/zstd*, *local/zstd*,
 and *xyz/zstd* styles are identical to the gz styles, but use the Zstd
@@ -219,6 +221,11 @@ you should thus consider the *atom* and *atom/mpiio* styles (etc) to
 be inter-changeable.  The one exception is how the filename is
 specified for the MPI-IO styles, as explained below.
 
+.. warning::
+
+   The MPIIO package is currently unmaintained and has become
+   unreliable. Use with caution.
+
 The precision of values output to text-based dump files can be
 controlled by the :doc:`dump_modify format <dump_modify>` command and
 its options.
@@ -275,10 +282,11 @@ This bounding box is convenient for many visualization programs.  The
 meaning of the 6 character flags for "xx yy zz" is the same as above.
 
 Note that the first two numbers on each line are now xlo_bound instead
-of xlo, etc, since they represent a bounding box.  See the :doc:`Howto triclinic <Howto_triclinic>` page for a geometric description
-of triclinic boxes, as defined by LAMMPS, simple formulas for how the
-6 bounding box extents (xlo_bound,xhi_bound,etc) are calculated from
-the triclinic parameters, and how to transform those parameters to and
+of xlo, etc, since they represent a bounding box.  See the :doc:`Howto
+triclinic <Howto_triclinic>` page for a geometric description of
+triclinic boxes, as defined by LAMMPS, simple formulas for how the 6
+bounding box extents (xlo_bound,xhi_bound,etc) are calculated from the
+triclinic parameters, and how to transform those parameters to and
 from other commonly used triclinic representations.
 
 The "ITEM: ATOMS" line in each snapshot lists column descriptors for
@@ -310,23 +318,24 @@ written to the dump file.  This local data is typically calculated by
 each processor based on the atoms it owns, but there may be zero or
 more entities per atom, e.g. a list of bond distances.  An explanation
 of the possible dump local attributes is given below.  Note that by
-using input from the :doc:`compute property/local <compute_property_local>` command with dump local,
-it is possible to generate information on bonds, angles, etc that can
-be cut and pasted directly into a data file read by the
-:doc:`read_data <read_data>` command.
+using input from the :doc:`compute property/local
+<compute_property_local>` command with dump local, it is possible to
+generate information on bonds, angles, etc that can be cut and pasted
+directly into a data file read by the :doc:`read_data <read_data>`
+command.
 
 Style *cfg* has the same command syntax as style *custom* and writes
-extended CFG format files, as used by the
-`AtomEye <http://li.mit.edu/Archive/Graphics/A/>`_ visualization
-package.  Since the extended CFG format uses a single snapshot of the
-system per file, a wildcard "\*" must be included in the filename, as
-discussed below.  The list of atom attributes for style *cfg* must
-begin with either "mass type xs ys zs" or "mass type xsu ysu zsu"
-since these quantities are needed to write the CFG files in the
-appropriate format (though the "mass" and "type" fields do not appear
-explicitly in the file).  Any remaining attributes will be stored as
-"auxiliary properties" in the CFG files.  Note that you will typically
-want to use the :doc:`dump_modify element <dump_modify>` command with
+extended CFG format files, as used by the `AtomEye
+<http://li.mit.edu/Archive/Graphics/A/>`_ visualization package.
+Since the extended CFG format uses a single snapshot of the system per
+file, a wildcard "\*" must be included in the filename, as discussed
+below.  The list of atom attributes for style *cfg* must begin with
+either "mass type xs ys zs" or "mass type xsu ysu zsu" since these
+quantities are needed to write the CFG files in the appropriate format
+(though the "mass" and "type" fields do not appear explicitly in the
+file).  Any remaining attributes will be stored as "auxiliary
+properties" in the CFG files.  Note that you will typically want to
+use the :doc:`dump_modify element <dump_modify>` command with
 CFG-formatted files, to associate element names with atom types, so
 that AtomEye can render atoms appropriately. When unwrapped
 coordinates *xsu*, *ysu*, and *zsu* are requested, the nominal AtomEye
@@ -452,6 +461,11 @@ use the :doc:`read_dump <read_dump>` command or perform other
 post-processing, just as if the dump file was not written using
 MPI-IO.
 
+.. warning::
+
+   The MPIIO package is currently unmaintained and has become
+   unreliable. Use with caution.
+
 Note that MPI-IO dump files are one large file which all processors
 write to.  You thus cannot use the "%" wildcard character described
 above in the filename since that specifies generation of multiple

doc/src/dump_modify.rst

@@ -17,7 +17,7 @@ Syntax
 * one or more keyword/value pairs may be appended
 
 * these keywords apply to various dump styles
-* keyword = *append* or *at* or *buffer* or *delay* or *element* or *every* or *fileper* or *first* or *flush* or *format* or *header* or *image* or *label* or *maxfiles* or *nfile* or *pad* or *pbc* or *precision* or *region* or *refresh* or *scale* or *sfactor* or *sort* or *tfactor* or *thermo* or *thresh* or *time* or *units* or *unwrap*
+* keyword = *append* or *at* or *buffer* or *delay* or *element* or *every* or *every/time* or *fileper* or *first* or *flush* or *format* or *header* or *image* or *label* or *maxfiles* or *nfile* or *pad* or *pbc* or *precision* or *region* or *refresh* or *scale* or *sfactor* or *sort* or *tfactor* or *thermo* or *thresh* or *time* or *units* or *unwrap*
 
   .. parsed-literal::
 
@@ -32,6 +32,9 @@ Syntax
        *every* arg = N
          N = dump every this many timesteps
          N can be a variable (see below)
+       *every/time* arg = Delta
+         Delta = dump every this interval in simulation time (time units)
+         Delta can be a variable (see below)
        *fileper* arg = Np
          Np = write one file for every this many processors
        *first* arg = *yes* or *no*
@@ -197,11 +200,19 @@ will be accepted.
 
 ----------
 
-The *every* keyword changes the dump frequency originally specified by
-the :doc:`dump <dump>` command to a new value.  The every keyword can be
-specified in one of two ways.  It can be a numeric value in which case
-it must be > 0.  Or it can be an :doc:`equal-style variable <variable>`,
-which should be specified as v_name, where name is the variable name.
+The *every* keyword can be used with any dump style except the *dcd*
+and *xtc* styles.  It does two things.  It specifies that the interval
+between dump snapshots will be set in timesteps, which is the default
+if the *every* or *every/time* keywords are not used.  See the
+*every/time* keyword for how to specify the interval in simulation
+time, i.e. in time units of the :doc:`units <units>` command.  The
+*every* keyword also sets the interval value, which overrides the dump
+frequency originally specified by the :doc:`dump <dump>` command.
+
+The *every* keyword can be specified in one of two ways.  It can be a
+numeric value in which case it must be > 0.  Or it can be an
+:doc:`equal-style variable <variable>`, which should be specified as
+v_name, where name is the variable name.
 
 In this case, the variable is evaluated at the beginning of a run to
 determine the next timestep at which a dump snapshot will be written
@@ -210,11 +221,12 @@ determine the next timestep, etc.  Thus the variable should return
 timestep values.  See the stagger() and logfreq() and stride() math
 functions for :doc:`equal-style variables <variable>`, as examples of
 useful functions to use in this context.  Other similar math functions
-could easily be added as options for :doc:`equal-style variables <variable>`.  Also see the next() function, which allows
-use of a file-style variable which reads successive values from a
-file, each time the variable is evaluated.  Used with the *every*
-keyword, if the file contains a list of ascending timesteps, you can
-output snapshots whenever you wish.
+could easily be added as options for :doc:`equal-style variables
+<variable>`.  Also see the next() function, which allows use of a
+file-style variable which reads successive values from a file, each
+time the variable is evaluated.  Used with the *every* keyword, if the
+file contains a list of ascending timesteps, you can output snapshots
+whenever you wish.
 
 Note that when using the variable option with the *every* keyword, you
 need to use the *first* option if you want an initial snapshot written
@@ -255,14 +267,103 @@ in file tmp.times:
 
 ----------
 
+The *every/time* keyword can be used with any dump style except the
+*dcd* and *xtc* styles.  It does two things.  It specifies that the
+interval between dump snapshots will be set in simulation time,
+i.e. in time units of the :doc:`units <units>` command.  This can be
+useful when the timestep size varies during a simulation run, e.g. by
+use of the :doc:`fix dt/reset <fix_dt_reset>` command.  The default is
+to specify the interval in timesteps; see the *every* keyword.  The
+*every/time* command also sets the interval value.
+
+.. note::
+
+   If you wish dump styles *atom*, *custom*, *local*, or *xyz* to
+   include the simulation time as a field in the header portion of
+   each snapshot, you also need to use the dump_modify *time* keyword
+   with a setting of *yes*.  See its documentation below.
+
+Note that since snapshots are output on simulation steps, each
+snapshot will be written on the first timestep whose associated
+simulation time is >= the exact snapshot time value.
+
+As with the *every* option, the *Delta* value can be specified in one
+of two ways.  It can be a numeric value in which case it must be >
+0.0.  Or it can be an :doc:`equal-style variable <variable>`, which
+should be specified as v_name, where name is the variable name.
+
+In this case, the variable is evaluated at the beginning of a run to
+determine the next simulation time at which a dump snapshot will be
+written out.  On that timestep the variable will be evaluated again to
+determine the next simulation time, etc.  Thus the variable should
+return values in time units.  Note the current timestep or simulation
+time can be used in an :doc:`equal-style variables <variable>` since
+they are both thermodynamic keywords.  Also see the next() function,
+which allows use of a file-style variable which reads successive
+values from a file, each time the variable is evaluated.  Used with
+the *every/time* keyword, if the file contains a list of ascending
+simulation times, you can output snapshots whenever you wish.
+
+Note that when using the variable option with the *every/time*
+keyword, you need to use the *first* option if you want an initial
+snapshot written to the dump file.  The *every/time* keyword cannot be
+used with the dump *dcd* style.
+
+For example, the following commands will write snapshots at successive
+simulation times which grow by a factor of 1.5 with each interval.
+The dt value used in the variable is to avoid a zero result when the
+initial simulation time is 0.0.
+
+.. code-block:: LAMMPS
+
+   variable        increase equal 1.5*(time+dt)
+   dump            1 all atom 100 tmp.dump
+   dump_modify     1 every/time v_increase first yes
+
+The following commands would write snapshots at the times listed in
+file tmp.times:
+
+.. code-block:: LAMMPS
+
+   variable        f file tmp.times
+   variable        s equal next(f)
+   dump            1 all atom 100 tmp.dump
+   dump_modify     1 every/time v_s
+
+.. note::
+
+   When using a file-style variable with the *every/time* keyword, the
+   file of timesteps must list a first time that is beyond the time
+   associated with the current timestep (e.g. it cannot be 0.0).  And
+   it must list one or more times beyond the length of the run you
+   perform.  This is because the dump command will generate an error
+   if the next time it reads from the file is not a value greater than
+   the current time.  Thus if you wanted output at times 0,15,100 of a
+   run of length 100 in simulation time, the file should contain the
+   values 15,100,101 and you should also use the dump_modify first
+   command.  Any final value > 100 could be used in place of 101.
+
+----------
+
 The *first* keyword determines whether a dump snapshot is written on
 the very first timestep after the dump command is invoked.  This will
-always occur if the current timestep is a multiple of N, the frequency
-specified in the :doc:`dump <dump>` command, including timestep 0.  But
-if this is not the case, a dump snapshot will only be written if the
-setting of this keyword is *yes*\ .  If it is *no*, which is the
+always occur if the current timestep is a multiple of $N$, the
+frequency specified in the :doc:`dump <dump>` command or
+:doc:`dump_modify every <dump_modify>` command, including timestep 0.
+It will also always occur if the current simulation time is a multiple
+of *Delta*, the time interval specified in the doc:`dump_modify
+every/time <dump_modify>` command.
+
+But if this is not the case, a dump snapshot will only be written if
+the setting of this keyword is *yes*\ .  If it is *no*, which is the
 default, then it will not be written.
 
+Note that if the argument to the :doc:`dump_modify every
+<dump_modify>` doc:`dump_modify every/time <dump_modify>` commands is
+a variable and not a numeric value, then specifying *first yes* is the
+only way to write a dump snapshot on the first timestep after the dump
+command is invoked.
+
 ----------
 
 The *flush* keyword determines whether a flush operation is invoked
@@ -342,10 +443,10 @@ The *fileper* keyword is documented below with the *nfile* keyword.
 
 ----------
 
-The *header* keyword toggles whether the dump file will include a header.
-Excluding a header will reduce the size of the dump file for fixes such as
-:doc:`fix pair/tracker <fix_pair_tracker>` which do not require the information
-typically written to the header.
+The *header* keyword toggles whether the dump file will include a
+header.  Excluding a header will reduce the size of the dump file for
+fixes such as :doc:`fix pair/tracker <fix_pair_tracker>` which do not
+require the information typically written to the header.
 
 ----------
 
@@ -561,7 +662,9 @@ The dump *local* style cannot be sorted by atom ID, since there are
 typically multiple lines of output per atom.  Some dump styles, such
 as *dcd* and *xtc*, require sorting by atom ID to format the output
 file correctly.  If multiple processors are writing the dump file, via
-the "%" wildcard in the dump filename, then sorting cannot be
+the "%" wildcard in the dump filename and the *nfile* or *fileper*
+keywords are set to non-default values (i.e. the number of dump file
+pieces is not equal to the number of procs), then sorting cannot be
 performed.
 
 .. note::
@@ -639,16 +742,20 @@ threshold criterion is met.  Otherwise it is not met.
 
 ----------
 
-The *time* keyword only applies to the dump *atom*, *custom*, and
-*local* styles (and their COMPRESS package versions *atom/gz*,
-*custom/gz* and *local/gz*\ ). If set to *yes*, each frame will will
-contain two extra lines before the "ITEM: TIMESTEP" entry:
+The *time* keyword only applies to the dump *atom*, *custom*, *local*,
+and *xyz* styles (and their COMPRESS package versions *atom/gz*,
+*custom/gz* and *local/gz*\ ).  For the first 3 styles, if set to
+*yes*, each frame will will contain two extra lines before the "ITEM:
+TIMESTEP" entry:
 
 .. parsed-literal::
 
    ITEM: TIME
    \<elapsed time\>
 
+For the *xyz* style, the simulation time is included on the same line
+as the timestep value.
+
 This will output the current elapsed simulation time in current
 time units equivalent to the :doc:`thermo keyword <thermo_style>` *time*\ .
 This is to simplify post-processing of trajectories using a variable time

doc/src/fix_dt_reset.rst

@@ -78,13 +78,20 @@ outer loop (largest) timestep, which is the same timestep that the
 
 Note that the cumulative simulation time (in time units), which
 accounts for changes in the timestep size as a simulation proceeds,
-can be accessed by the :doc:`thermo_style time <thermo_style>` keyword.
+can be accessed by the :doc:`thermo_style time <thermo_style>`
+keyword.
+
+Also note that the :doc:`dump_modify every/time <dump_modify>` option
+allows dump files to be written at intervals specified by simulation
+time, rather than by timesteps.  Simulation time is in time units;
+see the :doc:`units <units>` doc page for details.
 
 Restart, fix_modify, output, run start/stop, minimize info
 """""""""""""""""""""""""""""""""""""""""""""""""""""""""""
 
-No information about this fix is written to :doc:`binary restart files <restart>`.  None of the :doc:`fix_modify <fix_modify>` options
-are relevant to this fix.
+No information about this fix is written to :doc:`binary restart files
+<restart>`.  None of the :doc:`fix_modify <fix_modify>` options are
+relevant to this fix.
 
 This fix computes a global scalar which can be accessed by various
 :doc:`output commands <Howto_output>`.  The scalar stores the last
@@ -93,7 +100,8 @@ timestep on which the timestep was reset to a new value.
 The scalar value calculated by this fix is "intensive".
 
 No parameter of this fix can be used with the *start/stop* keywords of
-the :doc:`run <run>` command.  This fix is not invoked during :doc:`energy minimization <minimize>`.
+the :doc:`run <run>` command.  This fix is not invoked during
+:doc:`energy minimization <minimize>`.
 
 Restrictions
 """"""""""""
@@ -102,7 +110,7 @@ Restrictions
 Related commands
 """"""""""""""""
 
-:doc:`timestep <timestep>`
+:doc:`timestep <timestep>`, :doc:`dump_modify every/time <dump_modify>`
 
 Default
 """""""

doc/src/pair_nm.rst

@@ -195,5 +195,4 @@ none
 
 .. _Dietz:
 
-**(Dietz)** J.D. Dietz, R.S. Hoy, "Facile equilibration of well-entangled
-semi-flexible bead-spring polymer melts" arXiv:2109.11001
+**(Dietz)** Dietz and Hoy, J. Chem Phys, 156, 014103 (2022).

doc/src/pair_python.rst

@@ -126,11 +126,11 @@ and *compute_energy*, which both take 3 numerical arguments:
 * itype = the (numerical) type of the first atom
 * jtype = the (numerical) type of the second atom
 
-This functions need to compute the force and the energy, respectively,
-and use the result as return value. The functions need to use the
-*pmap* dictionary to convert the LAMMPS atom type number to the symbolic
-value of the internal potential parameter data structure. Following
-the *LJCutMelt* example, here are the two functions:
+This functions need to compute the (scaled) force and the energy,
+respectively, and use the result as return value. The functions need
+to use the *pmap* dictionary to convert the LAMMPS atom type number
+to the symbolic value of the internal potential parameter data structure.
+Following the *LJCutMelt* example, here are the two functions:
 
 .. code-block:: python
 
@@ -154,10 +154,10 @@ the *LJCutMelt* example, here are the two functions:
 
    for consistency with the C++ pair styles in LAMMPS, the
    *compute_force* function follows the conventions of the Pair::single()
-   methods and does not return the full force, but the force scaled by
-   the distance between the two atoms, so this value only needs to be
-   multiplied by delta x, delta y, and delta z to conveniently obtain the
-   three components of the force vector between these two atoms.
+   methods and does not return the pairwise force directly, but the force
+   divided by the distance between the two atoms, so this value only needs
+   to be  multiplied by delta x, delta y, and delta z to conveniently obtain
+   the three components of the force vector between these two atoms.
 
 ----------
 

lib/gpu/Makefile.cuda_mps

@@ -1,5 +1,5 @@
 # /* ----------------------------------------------------------------------
-#  Generic Linux Makefile for CUDA
+#  Generic Linux Makefile for CUDA with the Multi-Process Service (MPS)
 #     - change CUDA_ARCH for your GPU
 # ------------------------------------------------------------------------- */
 

lib/gpu/Makefile.hip

@@ -39,11 +39,9 @@ HIP_PLATFORM=$(shell $(HIP_PATH)/bin/hipconfig --platform)
 HIP_COMPILER=$(shell $(HIP_PATH)/bin/hipconfig --compiler)
 
 ifeq (hcc,$(HIP_PLATFORM))
-	HIP_OPTS  += -ffast-math
 	# possible values: gfx803,gfx900,gfx906
 	HIP_ARCH = gfx906
 else ifeq (amd,$(HIP_PLATFORM))
-	HIP_OPTS  += -ffast-math
 	# possible values: gfx803,gfx900,gfx906
 	HIP_ARCH = gfx906
 else ifeq (nvcc,$(HIP_PLATFORM))

lib/gpu/Makefile.linux

@@ -1,5 +1,5 @@
 # /* ----------------------------------------------------------------------   
-#  Generic Linux Makefile for CUDA 
+#  Generic Linux Makefile for CUDA
 #     - Change CUDA_ARCH for your GPU
 # ------------------------------------------------------------------------- */
 
@@ -13,7 +13,7 @@ endif
 
 NVCC = nvcc
 
-# obsolete hardware. not supported by current drivers anymore.
+# obsolete hardware. not supported by current drivers and toolkits anymore.
 #CUDA_ARCH = -arch=sm_13
 #CUDA_ARCH = -arch=sm_10 -DCUDA_PRE_THREE
 
@@ -28,11 +28,11 @@ NVCC = nvcc
 #CUDA_ARCH = -arch=sm_37
 
 # Maxwell hardware
-CUDA_ARCH = -arch=sm_50
+#CUDA_ARCH = -arch=sm_50
 #CUDA_ARCH = -arch=sm_52
 
 # Pascal hardware
-#CUDA_ARCH = -arch=sm_60
+CUDA_ARCH = -arch=sm_60
 #CUDA_ARCH = -arch=sm_61
 
 # Volta hardware
@@ -70,7 +70,7 @@ LIB_DIR = ./
 AR = ar
 BSH = /bin/sh
 
-# GPU binning not recommended with modern GPUs
+# GPU binning not recommended for most modern GPUs
 CUDPP_OPT = #-DUSE_CUDPP -Icudpp_mini
 
 include Nvidia.makefile

lib/gpu/Makefile.linux_multi

@@ -1,6 +1,6 @@
 # /* ----------------------------------------------------------------------   
-#  Generic Linux Makefile for CUDA 
-#     - Change CUDA_ARCH for your GPU
+#  Generic Linux Makefile for CUDA complied for multiple compute capabilities
+#     - Add your GPU to CUDA_CODE
 # ------------------------------------------------------------------------- */
 
 # which file will be copied to Makefile.lammps

lib/gpu/Makefile.mpi

@@ -0,0 +1 @@
+Makefile.linux

lib/gpu/Makefile.serial

@@ -1,5 +1,5 @@
 # /* ----------------------------------------------------------------------   
-#  Generic Linux Makefile for CUDA 
+#  Generic Linux Makefile for CUDA without MPI libraries
 #     - Change CUDA_ARCH for your GPU
 # ------------------------------------------------------------------------- */
 
@@ -28,11 +28,11 @@ NVCC = nvcc
 #CUDA_ARCH = -arch=sm_37
 
 # Maxwell hardware
-CUDA_ARCH = -arch=sm_50
+#CUDA_ARCH = -arch=sm_50
 #CUDA_ARCH = -arch=sm_52
 
 # Pascal hardware
-#CUDA_ARCH = -arch=sm_60
+CUDA_ARCH = -arch=sm_60
 #CUDA_ARCH = -arch=sm_61
 
 # Volta hardware
@@ -41,6 +41,10 @@ CUDA_ARCH = -arch=sm_50
 # Turing hardware
 #CUDA_ARCH = -arch=sm_75
 
+# Ampere hardware
+#CUDA_ARCH = -arch=sm_80
+#CUDA_ARCH = -arch=sm_86
+
 # this setting should match LAMMPS Makefile
 # one of LAMMPS_SMALLBIG (default), LAMMPS_BIGBIG and LAMMPS_SMALLSMALL
 

lib/gpu/Makefile.turing

@@ -1,23 +0,0 @@
-NVCC  = $(CUDA_HOME)/bin/nvcc
-EXTRAMAKE = Makefile.lammps.standard
-
-CUDA_ARCH = -arch=sm_75
-CUDA_PRECISION = -D_SINGLE_DOUBLE
-CUDA_INCLUDE = -I$(CUDA_HOME)/include
-CUDA_LIB = -L$(CUDA_HOME)/lib64 -Xlinker -rpath -Xlinker $(CUDA_HOME)/lib64 -lcudart
-CUDA_OPTS = -DUNIX -O3 --use_fast_math --ftz=true
-
-CUDR_CPP = mpic++ -DMPI_GERYON -DUCL_NO_EXIT -I$(CUDA_HOME)/include
-CUDR_OPTS = -O3 -ffast-math -funroll-loops -DMPI_GERYON -DLAMMPS_SMALLBIG
-
-BIN_DIR = .
-OBJ_DIR = obj
-LIB_DIR = .
-AR = ar
-BSH = /bin/sh
-
-# GPU binning not recommended with most modern GPUs
-CUDPP_OPT = #-DUSE_CUDPP -Icudpp_mini
-
-include Nvidia.makefile
-

src/.gitignore

@@ -429,6 +429,8 @@
 /commgrid.h
 /compute_ackland_atom.cpp
 /compute_ackland_atom.h
+/compute_ave_sphere_atom.cpp
+/compute_ave_sphere_atom.h
 /compute_basal_atom.cpp
 /compute_basal_atom.h
 /compute_body_local.cpp

src/CLASS2/angle_class2.cpp

@@ -169,6 +169,8 @@ void AngleClass2::compute(int eflag, int vflag)
 
     // force & energy for bond-angle term
 
+    dr1 = r1 - ba_r1[type];
+    dr2 = r2 - ba_r2[type];
     aa1 = s * dr1 * ba_k1[type];
     aa2 = s * dr2 * ba_k2[type];
 
@@ -459,6 +461,9 @@ double AngleClass2::single(int type, int i1, int i2, int i3)
   double dr2 = r2 - bb_r2[type];
   energy += bb_k[type]*dr1*dr2;
 
+  dr1 = r1 - ba_r1[type];
+  dr2 = r2 - ba_r2[type];
   energy += ba_k1[type]*dr1*dtheta + ba_k2[type]*dr2*dtheta;
+
   return energy;
 }

src/COMPRESS/dump_xyz_gz.cpp

@@ -88,11 +88,16 @@ void DumpXYZGZ::openfile()
   if (multifile) delete[] filecurrent;
 }
 
+/* ---------------------------------------------------------------------- */
+
 void DumpXYZGZ::write_header(bigint ndump)
 {
   if (me == 0) {
-    std::string header = fmt::format("{}\n", ndump);
-    header += fmt::format("Atoms. Timestep: {}\n", update->ntimestep);
+    auto header = fmt::format("{}\n", ndump);
+    if (time_flag) {
+      double tcurrent = update->atime + (update->ntimestep-update->atimestep) + update->dt;
+      header += fmt::format(" Atoms. Timestep: {} Time: {:.6f}\n", update->ntimestep, tcurrent);
+    } else header += fmt::format(" Atoms. Timestep: {}\n", update->ntimestep);
     writer.write(header.c_str(), header.length());
   }
 }

src/COMPRESS/dump_xyz_zstd.cpp

@@ -96,11 +96,16 @@ void DumpXYZZstd::openfile()
   if (multifile) delete[] filecurrent;
 }
 
+/* ---------------------------------------------------------------------- */
+
 void DumpXYZZstd::write_header(bigint ndump)
 {
   if (me == 0) {
-    std::string header = fmt::format("{}\n", ndump);
-    header += fmt::format("Atoms. Timestep: {}\n", update->ntimestep);
+    auto header = fmt::format("{}\n", ndump);
+    if (time_flag) {
+      double tcurrent = update->atime + (update->ntimestep-update->atimestep) + update->dt;
+      header += fmt::format(" Atoms. Timestep: {} Time: {:.6f}\n", update->ntimestep, tcurrent);
+    } else header += fmt::format(" Atoms. Timestep: {}\n", update->ntimestep);
     writer.write(header.c_str(), header.length());
   }
 }

src/Depend.sh

@@ -77,6 +77,10 @@ if (test $1 = "DPD-BASIC") then
   depend INTEL
 fi
 
+if (test $1 = "EXTRA-COMPUTE") then
+  depend KOKKOS
+fi
+
 if (test $1 = "EXTRA-MOLECULE") then
   depend GPU
   depend OPENMP

src/EXTRA-COMPUTE/compute_ave_sphere_atom.cpp

@@ -0,0 +1,278 @@
+/* ----------------------------------------------------------------------
+   LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
+   https://www.lammps.org/, Sandia National Laboratories
+   Steve Plimpton, sjplimp@sandia.gov
+
+   Copyright (2003) Sandia Corporation.  Under the terms of Contract
+   DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains
+   certain rights in this software.  This software is distributed under
+   the GNU General Public License.
+
+   See the README file in the top-level LAMMPS directory.
+------------------------------------------------------------------------- */
+
+#include "compute_ave_sphere_atom.h"
+
+#include "atom.h"
+#include "comm.h"
+#include "error.h"
+#include "force.h"
+#include "group.h"
+#include "memory.h"
+#include "modify.h"
+#include "neigh_list.h"
+#include "neigh_request.h"
+#include "neighbor.h"
+#include "pair.h"
+#include "update.h"
+#include "math_const.h"
+
+#include <cmath>
+#include <cstring>
+
+using namespace LAMMPS_NS;
+using namespace MathConst;
+
+
+/* ---------------------------------------------------------------------- */
+
+ComputeAveSphereAtom::ComputeAveSphereAtom(LAMMPS *lmp, int narg, char **arg) :
+  Compute(lmp, narg, arg),
+  result(nullptr)
+{
+  if (narg < 3 || narg > 5) error->all(FLERR,"Illegal compute ave/sphere/atom command");
+
+  // process optional args
+
+  cutoff = 0.0;
+
+  int iarg = 3;
+  while (iarg < narg) {
+    if (strcmp(arg[iarg],"cutoff") == 0) {
+      if (iarg+2 > narg) error->all(FLERR,"Illegal compute ave/sphere/atom command");
+      cutoff = utils::numeric(FLERR,arg[iarg+1],false,lmp);
+      if (cutoff <= 0.0) error->all(FLERR,"Illegal compute ave/sphere/atom command");
+      iarg += 2;
+    } else error->all(FLERR,"Illegal compute ave/sphere/atom command");
+  }
+
+  peratom_flag = 1;
+  size_peratom_cols = 2;
+  comm_forward = 3;
+
+  nmax = 0;
+}
+
+/* ---------------------------------------------------------------------- */
+
+ComputeAveSphereAtom::~ComputeAveSphereAtom()
+{
+  if (copymode) return;
+
+  memory->destroy(result);
+}
+
+/* ---------------------------------------------------------------------- */
+
+void ComputeAveSphereAtom::init()
+{
+  if (!force->pair && cutoff == 0.0)
+    error->all(FLERR,"Compute ave/sphere/atom requires a cutoff be specified "
+               "or a pair style be defined");
+
+  double skin = neighbor->skin;
+  if (cutoff != 0.0) {
+    double cutghost;            // as computed by Neighbor and Comm
+    if (force->pair)
+      cutghost = MAX(force->pair->cutforce+skin,comm->cutghostuser);
+    else
+      cutghost = comm->cutghostuser;
+
+    if (cutoff > cutghost)
+      error->all(FLERR,"Compute ave/sphere/atom cutoff exceeds ghost atom range - "
+                 "use comm_modify cutoff command");
+  }
+
+  int cutflag = 1;
+  if (force->pair) {
+    if (cutoff == 0.0) {
+      cutoff = force->pair->cutforce;
+    }
+    if (cutoff <= force->pair->cutforce+skin) cutflag = 0;
+  }
+
+  cutsq = cutoff*cutoff;
+  sphere_vol = 4.0/3.0*MY_PI*cutsq*cutoff;
+
+  // need an occasional full neighbor list
+
+  int irequest = neighbor->request(this,instance_me);
+  neighbor->requests[irequest]->pair = 0;
+  neighbor->requests[irequest]->compute = 1;
+  neighbor->requests[irequest]->half = 0;
+  neighbor->requests[irequest]->full = 1;
+  neighbor->requests[irequest]->occasional = 1;
+  if (cutflag) {
+    neighbor->requests[irequest]->cut = 1;
+    neighbor->requests[irequest]->cutoff = cutoff;
+  }
+}
+
+/* ---------------------------------------------------------------------- */
+
+void ComputeAveSphereAtom::init_list(int /*id*/, NeighList *ptr)
+{
+  list = ptr;
+}
+
+/* ---------------------------------------------------------------------- */
+
+void ComputeAveSphereAtom::compute_peratom()
+{
+  int i,j,ii,jj,inum,jnum;
+  double xtmp,ytmp,ztmp,delx,dely,delz,rsq;
+  int *ilist,*jlist,*numneigh,**firstneigh;
+  int count;
+  double vsum[3],vavg[3],vnet[3];
+
+  invoked_peratom = update->ntimestep;
+
+  // grow result array if necessary
+
+  if (atom->nmax > nmax) {
+    memory->destroy(result);
+    nmax = atom->nmax;
+    memory->create(result,nmax,2,"ave/sphere/atom:result");
+    array_atom = result;
+  }
+
+  // need velocities of ghost atoms
+
+  comm->forward_comm_compute(this);
+
+  // invoke full neighbor list (will copy or build if necessary)
+
+  neighbor->build_one(list);
+
+  inum = list->inum;
+  ilist = list->ilist;
+  numneigh = list->numneigh;
+  firstneigh = list->firstneigh;
+
+  // compute properties for each atom in group
+  // use full neighbor list to count atoms less than cutoff
+
+  double **x = atom->x;
+  double **v = atom->v;
+  int *mask = atom->mask;
+
+  for (ii = 0; ii < inum; ii++) {
+    i = ilist[ii];
+
+    if (mask[i] & groupbit) {
+      xtmp = x[i][0];
+      ytmp = x[i][1];
+      ztmp = x[i][2];
+      jlist = firstneigh[i];
+      jnum = numneigh[i];
+
+      // i atom contribution
+
+      count = 1;
+      vsum[0] = v[i][0];
+      vsum[1] = v[i][1];
+      vsum[2] = v[i][2];
+
+      for (jj = 0; jj < jnum; jj++) {
+        j = jlist[jj];
+        j &= NEIGHMASK;
+
+        delx = xtmp - x[j][0];
+        dely = ytmp - x[j][1];
+        delz = ztmp - x[j][2];
+        rsq = delx*delx + dely*dely + delz*delz;
+        if (rsq < cutsq) {
+          count++;
+          vsum[0] += v[j][0];
+          vsum[1] += v[j][1];
+          vsum[2] += v[j][2];
+        }
+      }
+
+      vavg[0] = vsum[0]/count;
+      vavg[1] = vsum[1]/count;
+      vavg[2] = vsum[2]/count;
+
+      // i atom contribution
+
+      count = 1;
+      vnet[0] = v[i][0] - vavg[0];
+      vnet[1] = v[i][1] - vavg[1];
+      vnet[2] = v[i][2] - vavg[2];
+      double ke_sum = vnet[0]*vnet[0] + vnet[1]*vnet[1] + vnet[2]*vnet[2];
+
+      for (jj = 0; jj < jnum; jj++) {
+        j = jlist[jj];
+        j &= NEIGHMASK;
+
+        delx = xtmp - x[j][0];
+        dely = ytmp - x[j][1];
+        delz = ztmp - x[j][2];
+        rsq = delx*delx + dely*dely + delz*delz;
+        if (rsq < cutsq) {
+          count++;
+          vnet[0] = v[j][0] - vavg[0];
+          vnet[1] = v[j][1] - vavg[1];
+          vnet[2] = v[j][2] - vavg[2];
+          ke_sum += vnet[0]*vnet[0] + vnet[1]*vnet[1] + vnet[2]*vnet[2];
+        }
+      }
+      double density = count/sphere_vol;
+      double temp = ke_sum/3.0/count;
+      result[i][0] = density;
+      result[i][1] = temp;
+    }
+  }
+}
+
+/* ---------------------------------------------------------------------- */
+
+int ComputeAveSphereAtom::pack_forward_comm(int n, int *list, double *buf,
+                                        int /*pbc_flag*/, int * /*pbc*/)
+{
+  double **v = atom->v;
+
+  int i,m=0;
+  for (i = 0; i < n; ++i) {
+    buf[m++] = v[list[i]][0];
+    buf[m++] = v[list[i]][1];
+    buf[m++] = v[list[i]][2];
+  }
+
+  return m;
+}
+
+/* ---------------------------------------------------------------------- */
+
+void ComputeAveSphereAtom::unpack_forward_comm(int n, int first, double *buf)
+{
+  double **v = atom->v;
+
+  int i,last,m=0;
+  last = first + n;
+  for (i = first; i < last; ++i) {
+    v[i][0] = buf[m++];
+    v[i][1] = buf[m++];
+    v[i][2] = buf[m++];
+  }
+}
+
+/* ----------------------------------------------------------------------
+   memory usage of local atom-based array
+------------------------------------------------------------------------- */
+
+double ComputeAveSphereAtom::memory_usage()
+{
+  double bytes = (double)2*nmax * sizeof(double);
+  return bytes;
+}

src/EXTRA-COMPUTE/compute_ave_sphere_atom.h

@@ -0,0 +1,67 @@
+/* -*- c++ -*- ----------------------------------------------------------
+   LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
+   https://www.lammps.org/, Sandia National Laboratories
+   Steve Plimpton, sjplimp@sandia.gov
+
+   Copyright (2003) Sandia Corporation.  Under the terms of Contract
+   DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains
+   certain rights in this software.  This software is distributed under
+   the GNU General Public License.
+
+   See the README file in the top-level LAMMPS directory.
+------------------------------------------------------------------------- */
+
+#ifdef COMPUTE_CLASS
+
+ComputeStyle(ave/sphere/atom,ComputeAveSphereAtom)
+
+#else
+
+#ifndef LMP_COMPUTE_AVE_SPHERE_ATOM_H
+#define LMP_COMPUTE_AVE_SPHERE_ATOM_H
+
+#include "compute.h"
+
+namespace LAMMPS_NS {
+
+class ComputeAveSphereAtom : public Compute {
+ public:
+  ComputeAveSphereAtom(class LAMMPS *, int, char **);
+  virtual ~ComputeAveSphereAtom();
+  virtual void init();
+  void init_list(int, class NeighList *);
+  virtual void compute_peratom();
+  int pack_forward_comm(int, int *, double *, int, int *);
+  void unpack_forward_comm(int, int, double *);
+  double memory_usage();
+
+ protected:
+  int nmax;
+  double cutoff,cutsq,sphere_vol;
+  class NeighList *list;
+
+  double **result;
+};
+
+}
+
+#endif
+#endif
+
+/* ERROR/WARNING messages:
+
+E: Illegal ... command
+
+Self-explanatory.  Check the input script syntax and compare to the
+documentation for the command.  You can use -echo screen as a
+command-line option when running LAMMPS to see the offending line.
+
+E: Compute ave/sphere/atom requires a cutoff be specified or a pair style be defined
+
+Self-explanatory.
+
+E: Compute ave/sphere/atom cutoff exceeds ghost atom range - use comm_modify cutoff command
+
+Self-explanatory.
+
+*/

src/EXTRA-DUMP/dump_dcd.cpp

@@ -100,13 +100,17 @@ void DumpDCD::init_style()
   if (sort_flag == 0 || sortcol != 0)
     error->all(FLERR,"Dump dcd requires sorting by atom ID");
 
-  // check that dump frequency has not changed and is not a variable
-  // but only when not being called from the "write_dump" command.
+  // check that dump modify settings are compatible with dcd
+  // but only when not being called from the "write_dump" command
 
   if (strcmp(id,"WRITE_DUMP") != 0) {
     int idump;
     for (idump = 0; idump < output->ndump; idump++)
       if (strcmp(id,output->dump[idump]->id) == 0) break;
+
+    if (output->mode_dump[idump] == 1)
+      error->all(FLERR,"Cannot use every/time setting for dump dcd");
+
     if (output->every_dump[idump] == 0)
       error->all(FLERR,"Cannot use variable every setting for dump dcd");
 

src/EXTRA-DUMP/dump_xtc.cpp

@@ -121,17 +121,24 @@ void DumpXTC::init_style()
 
   if (flush_flag) error->all(FLERR,"Cannot set dump_modify flush for dump xtc");
 
-  // check that dump frequency has not changed and is not a variable
+  // check that dump modify settings are compatible with xtc
+  // but only when not being called from the "write_dump" command
 
-  int idump;
-  for (idump = 0; idump < output->ndump; idump++)
-    if (strcmp(id,output->dump[idump]->id) == 0) break;
-  if (output->every_dump[idump] == 0)
-    error->all(FLERR,"Cannot use variable every setting for dump xtc");
+  if (strcmp(id,"WRITE_DUMP") != 0) {
+    int idump;
+    for (idump = 0; idump < output->ndump; idump++)
+      if (strcmp(id,output->dump[idump]->id) == 0) break;
 
-  if (nevery_save == 0) nevery_save = output->every_dump[idump];
-  else if (nevery_save != output->every_dump[idump])
-    error->all(FLERR,"Cannot change dump_modify every for dump xtc");
+    if (output->mode_dump[idump] == 1)
+      error->all(FLERR,"Cannot use every/time setting for dump xtc");
+
+    if (output->every_dump[idump] == 0)
+      error->all(FLERR,"Cannot use every variable setting for dump xtc");
+
+    if (nevery_save == 0) nevery_save = output->every_dump[idump];
+    else if (nevery_save != output->every_dump[idump])
+      error->all(FLERR,"Cannot change dump_modify every for dump xtc");
+  }
 }
 
 /* ---------------------------------------------------------------------- */

src/EXTRA-MOLECULE/bond_fene_nm.cpp

@@ -47,7 +47,7 @@ void BondFENENM::compute(int eflag, int vflag)
 {
   int i1, i2, n, type;
   double delx, dely, delz, ebond, fbond;
-  double rsq, r0sq, rlogarg, sr2, sr6;
+  double rsq, r0sq, rlogarg, sr6;
   double r;
 
   ebond = sr6 = 0.0;

src/EXTRA-MOLECULE/dihedral_table_cut.cpp

@@ -82,20 +82,14 @@ enum { //GSL status return codes.
   GSL_EBADLEN  = 19
 };
 
-
-
 // cyc_splintD(): Evaluate the deriviative of a cyclic spline at position x,
 //           with n control points at xa[], ya[], with parameters y2a[].
 //           The xa[] must be monotonically increasing and their
 //           range should not exceed period (ie xa[n-1] < xa[0] + period).
 //           x must lie in the range:  [(xa[n-1]-period), (xa[0]+period)]
 //           "period" is typically 2*PI.
-static double cyc_splintD(double const *xa,
-                          double const *ya,
-                          double const *y2a,
-                          int n,
-                          double period,
-                          double x)
+static double cyc_splintD(double const *xa, double const *ya, double const *y2a,
+                          int n, double period, double x)
 {
   int klo = -1;
   int khi = n; // (not n-1)
@@ -490,8 +484,7 @@ void DihedralTableCut::coeff(int narg, char **arg)
   if (tb->ninput < 2)
     error->all(FLERR,"Invalid dihedral table length: {}",arg[5]);
   else if ((tb->ninput == 2) && (tabstyle == SPLINE))
-    error->all(FLERR,"Invalid dihedral spline table length: {} "
-                                 "(Try linear)",arg[5]);
+    error->all(FLERR,"Invalid dihedral spline table length: {} (Try linear)",arg[5]);
 
   // check for monotonicity
   for (int i=0; i < tb->ninput-1; i++) {
@@ -509,12 +502,10 @@ void DihedralTableCut::coeff(int narg, char **arg)
   double phihi = tb->phifile[tb->ninput-1];
   if (tb->use_degrees) {
     if ((phihi - philo) >= 360)
-      error->all(FLERR,"Dihedral table angle range must be < 360 "
-                                   "degrees ({})",arg[5]);
+      error->all(FLERR,"Dihedral table angle range must be < 360 degrees ({})",arg[5]);
   } else {
     if ((phihi - philo) >= MY_2PI)
-      error->all(FLERR,"Dihedral table angle range must be < 2*PI "
-                                   "radians ({})",arg[5]);
+      error->all(FLERR,"Dihedral table angle range must be < 2*PI radians ({})",arg[5]);
   }
 
   // convert phi from degrees to radians
@@ -532,9 +523,9 @@ void DihedralTableCut::coeff(int narg, char **arg)
   // We also want the angles to be sorted in increasing order.
   // This messy code fixes these problems with the user's data:
   {
-    double *phifile_tmp = new double [tb->ninput];  //temporary arrays
-    double *ffile_tmp = new double [tb->ninput];  //used for sorting
-    double *efile_tmp = new double [tb->ninput];
+    double *phifile_tmp = new double[tb->ninput];  //temporary arrays
+    double *ffile_tmp = new double[tb->ninput];  //used for sorting
+    double *efile_tmp = new double[tb->ninput];
 
     // After re-imaging, does the range of angles cross the 0 or 2*PI boundary?
     // If so, find the discontinuity:

src/GPU/fix_nve_gpu.cpp

@@ -37,7 +37,7 @@ using namespace FixConst;
 FixNVEGPU::FixNVEGPU(LAMMPS *lmp, int narg, char **arg) :
   FixNVE(lmp, narg, arg)
 {
-  _dtfm = 0;
+  _dtfm = nullptr;
   _nlocal_max = 0;
 }
 
@@ -57,7 +57,11 @@ void FixNVEGPU::setup(int vflag)
     _respa_on = 1;
   else
     _respa_on = 0;
-  if (atom->ntypes > 1) reset_dt();
+
+  // ensure that _dtfm array is initialized if the group is not "all"
+  // or there is more than one atom type as that re-ordeted array is used for
+  // per-type/per-atom masses and group membership detection.
+  if ((igroup != 0) || (atom->ntypes > 1)) reset_dt();
 }
 
 /* ----------------------------------------------------------------------

src/H5MD/dump_h5md.cpp

@@ -181,6 +181,7 @@ DumpH5MD::DumpH5MD(LAMMPS *lmp, int narg, char **arg) : Dump(lmp, narg, arg)
   // allocate global array for atom coords
 
   bigint n = group->count(igroup);
+  if ((bigint) domain->dimension*n > MAXSMALLINT) error->all(FLERR,"Too many atoms for dump h5md");
   natoms = static_cast<int> (n);
 
   if (every_position>=0)

src/H5MD/dump_h5md.h

@@ -90,6 +90,11 @@ E: Dump h5md requires sorting by atom ID
 
 Use the dump_modify sort command to enable this.
 
+E: Too many atoms for dump h5md
+
+The system size must fit in a 32-bit integer to use this dump
+style.
+
 E: Cannot use variable every setting for dump xtc
 
 The format of this file requires snapshots at regular intervals.

src/KOKKOS/Install.sh

@@ -88,6 +88,8 @@ action comm_kokkos.cpp
 action comm_kokkos.h
 action comm_tiled_kokkos.cpp
 action comm_tiled_kokkos.h
+action compute_ave_sphere_atom_kokkos.cpp compute_ave_sphere_atom.cpp
+action compute_ave_sphere_atom_kokkos.h compute_ave_sphere_atom.h
 action compute_coord_atom_kokkos.cpp
 action compute_coord_atom_kokkos.h
 action compute_orientorder_atom_kokkos.cpp

src/KOKKOS/angle_class2_kokkos.cpp

@@ -224,8 +224,8 @@ void AngleClass2Kokkos<DeviceType>::operator()(TagAngleClass2Compute<NEWTON_BOND
 
   // force & energy for bond-bond term
 
-  const F_FLOAT dr1 = r1 - d_bb_r1[type];
-  const F_FLOAT dr2 = r2 - d_bb_r2[type];
+  F_FLOAT dr1 = r1 - d_bb_r1[type];
+  F_FLOAT dr2 = r2 - d_bb_r2[type];
   const F_FLOAT tk1 = d_bb_k[type] * dr1;
   const F_FLOAT tk2 = d_bb_k[type] * dr2;
 
@@ -241,6 +241,8 @@ void AngleClass2Kokkos<DeviceType>::operator()(TagAngleClass2Compute<NEWTON_BOND
 
   // force & energy for bond-angle term
 
+  dr1 = r1 - d_ba_r1[type];
+  dr2 = r2 - d_ba_r2[type];
   const F_FLOAT aa1 = s * dr1 * d_ba_k1[type];
   const F_FLOAT aa2 = s * dr2 * d_ba_k2[type];
 

src/KOKKOS/atom_vec_angle_kokkos.cpp

@@ -1630,7 +1630,7 @@ void AtomVecAngleKokkos::create_atom(int itype, double *coord)
 ------------------------------------------------------------------------- */
 
 void AtomVecAngleKokkos::data_atom(double *coord, imageint imagetmp,
-                                  char **values)
+                                  const std::vector<std::string> &values)
 {
   int nlocal = atom->nlocal;
   if (nlocal == nmax) grow(0);
@@ -1663,9 +1663,10 @@ void AtomVecAngleKokkos::data_atom(double *coord, imageint imagetmp,
    initialize other atom quantities for this sub-style
 ------------------------------------------------------------------------- */
 
-int AtomVecAngleKokkos::data_atom_hybrid(int nlocal, char **values)
+int AtomVecAngleKokkos::data_atom_hybrid(int nlocal, const std::vector<std::string> &values,
+                                         int offset)
 {
-  h_molecule(nlocal) = utils::inumeric(FLERR,values[0],true,lmp);
+  h_molecule(nlocal) = utils::inumeric(FLERR,values[offset],true,lmp);
   h_num_bond(nlocal) = 0;
   h_num_angle(nlocal) = 0;
   return 1;

src/KOKKOS/atom_vec_angle_kokkos.h

@@ -51,8 +51,8 @@ class AtomVecAngleKokkos : public AtomVecKokkos {
   int pack_restart(int, double *);
   int unpack_restart(double *);
   void create_atom(int, double *);
-  void data_atom(double *, tagint, char **);
-  int data_atom_hybrid(int, char **);
+  void data_atom(double *, imageint, const std::vector<std::string> &);
+  int data_atom_hybrid(int, const std::vector<std::string> &, int);
   void pack_data(double **);
   int pack_data_hybrid(int, double *);
   void write_data(FILE *, int, double **);

src/KOKKOS/atom_vec_atomic_kokkos.cpp

@@ -821,8 +821,8 @@ void AtomVecAtomicKokkos::create_atom(int itype, double *coord)
    initialize other atom quantities
 ------------------------------------------------------------------------- */
 
-void AtomVecAtomicKokkos::data_atom(double *coord, tagint imagetmp,
-                                    char **values)
+void AtomVecAtomicKokkos::data_atom(double *coord, imageint imagetmp,
+                                    const std::vector<std::string> &values)
 {
   int nlocal = atom->nlocal;
   if (nlocal == nmax) grow(0);

src/KOKKOS/atom_vec_atomic_kokkos.h

@@ -44,7 +44,7 @@ class AtomVecAtomicKokkos : public AtomVecKokkos {
   int pack_restart(int, double *);
   int unpack_restart(double *);
   void create_atom(int, double *);
-  void data_atom(double *, tagint, char **);
+  void data_atom(double *, imageint, const std::vector<std::string> &);
   void pack_data(double **);
   void write_data(FILE *, int, double **);
   double memory_usage();

src/KOKKOS/atom_vec_bond_kokkos.cpp

@@ -1056,7 +1056,7 @@ void AtomVecBondKokkos::create_atom(int itype, double *coord)
 ------------------------------------------------------------------------- */
 
 void AtomVecBondKokkos::data_atom(double *coord, imageint imagetmp,
-                                  char **values)
+                                  const std::vector<std::string> &values)
 {
   int nlocal = atomKK->nlocal;
   if (nlocal == nmax) grow(0);
@@ -1088,9 +1088,10 @@ void AtomVecBondKokkos::data_atom(double *coord, imageint imagetmp,
    initialize other atom quantities for this sub-style
 ------------------------------------------------------------------------- */
 
-int AtomVecBondKokkos::data_atom_hybrid(int nlocal, char **values)
+int AtomVecBondKokkos::data_atom_hybrid(int nlocal, const std::vector<std::string> &values,
+                                        int offset)
 {
-  h_molecule(nlocal) = utils::inumeric(FLERR,values[0],true,lmp);
+  h_molecule(nlocal) = utils::inumeric(FLERR,values[offset],true,lmp);
   h_num_bond(nlocal) = 0;
   return 1;
 }

src/KOKKOS/atom_vec_bond_kokkos.h

@@ -45,8 +45,8 @@ class AtomVecBondKokkos : public AtomVecKokkos {
   int pack_restart(int, double *);
   int unpack_restart(double *);
   void create_atom(int, double *);
-  void data_atom(double *, tagint, char **);
-  int data_atom_hybrid(int, char **);
+  void data_atom(double *, imageint, const std::vector<std::string> &);
+  int data_atom_hybrid(int, const std::vector<std::string> &, int);
   void pack_data(double **);
   int pack_data_hybrid(int, double *);
   void write_data(FILE *, int, double **);

src/KOKKOS/atom_vec_charge_kokkos.cpp

@@ -955,7 +955,7 @@ void AtomVecChargeKokkos::create_atom(int itype, double *coord)
 ------------------------------------------------------------------------- */
 
 void AtomVecChargeKokkos::data_atom(double *coord, imageint imagetmp,
-                                    char **values)
+                                    const std::vector<std::string> &values)
 {
   int nlocal = atom->nlocal;
   if (nlocal == nmax) grow(0);
@@ -987,9 +987,10 @@ void AtomVecChargeKokkos::data_atom(double *coord, imageint imagetmp,
    initialize other atom quantities for this sub-style
 ------------------------------------------------------------------------- */
 
-int AtomVecChargeKokkos::data_atom_hybrid(int nlocal, char **values)
+int AtomVecChargeKokkos::data_atom_hybrid(int nlocal, const std::vector<std::string> &values,
+                                          int offset)
 {
-  h_q[nlocal] = utils::numeric(FLERR,values[0],true,lmp);
+  h_q[nlocal] = utils::numeric(FLERR,values[offset],true,lmp);
 
   return 1;
 }

src/KOKKOS/atom_vec_charge_kokkos.h

@@ -46,8 +46,8 @@ class AtomVecChargeKokkos : public AtomVecKokkos {
   int pack_restart(int, double *);
   int unpack_restart(double *);
   void create_atom(int, double *);
-  void data_atom(double *, tagint, char **);
-  int data_atom_hybrid(int , char **);
+  void data_atom(double *, imageint, const std::vector<std::string> &);
+  int data_atom_hybrid(int , const std::vector<std::string> &, int);
   void pack_data(double **);
   int pack_data_hybrid(int, double *);
   void write_data(FILE *, int, double **);

src/KOKKOS/atom_vec_dpd_kokkos.cpp

@@ -1716,8 +1716,8 @@ void AtomVecDPDKokkos::create_atom(int itype, double *coord)
    initialize other atom quantities
 ------------------------------------------------------------------------- */
 
-void AtomVecDPDKokkos::data_atom(double *coord, tagint imagetmp,
-                                    char **values)
+void AtomVecDPDKokkos::data_atom(double *coord, imageint imagetmp,
+                                 const std::vector<std::string> &values)
 {
   int nlocal = atom->nlocal;
   if (nlocal == nmax) grow(0);
@@ -1759,9 +1759,10 @@ void AtomVecDPDKokkos::data_atom(double *coord, tagint imagetmp,
    initialize other atom quantities for this sub-style
 ------------------------------------------------------------------------- */
 
-int AtomVecDPDKokkos::data_atom_hybrid(int nlocal, char **values)
+int AtomVecDPDKokkos::data_atom_hybrid(int nlocal, const std::vector<std::string> &values,
+                                       int offset)
 {
-  h_dpdTheta(nlocal) = utils::numeric(FLERR,values[0],true,lmp);
+  h_dpdTheta(nlocal) = utils::numeric(FLERR,values[offset],true,lmp);
 
   atomKK->modified(Host,DPDTHETA_MASK);
 

src/KOKKOS/atom_vec_dpd_kokkos.h

@@ -54,8 +54,8 @@ class AtomVecDPDKokkos : public AtomVecKokkos {
   int pack_restart(int, double *);
   int unpack_restart(double *);
   void create_atom(int, double *);
-  void data_atom(double *, tagint, char **);
-  int data_atom_hybrid(int, char **);
+  void data_atom(double *, imageint, const std::vector<std::string> &);
+  int data_atom_hybrid(int, const std::vector<std::string> &, int);
   void pack_data(double **);
   int pack_data_hybrid(int, double *);
   void write_data(FILE *, int, double **);

src/KOKKOS/atom_vec_full_kokkos.cpp

@@ -1488,7 +1488,7 @@ void AtomVecFullKokkos::create_atom(int itype, double *coord)
 ------------------------------------------------------------------------- */
 
 void AtomVecFullKokkos::data_atom(double *coord, imageint imagetmp,
-                                       char **values)
+                                  const std::vector<std::string> &values)
 {
   int nlocal = atom->nlocal;
   if (nlocal == nmax) grow(0);
@@ -1525,10 +1525,11 @@ void AtomVecFullKokkos::data_atom(double *coord, imageint imagetmp,
    initialize other atom quantities for this sub-style
 ------------------------------------------------------------------------- */
 
-int AtomVecFullKokkos::data_atom_hybrid(int nlocal, char **values)
+int AtomVecFullKokkos::data_atom_hybrid(int nlocal, const std::vector<std::string> &values,
+                                        int offset)
 {
-  h_molecule(nlocal) = utils::inumeric(FLERR,values[0],true,lmp);
-  h_q(nlocal) = utils::numeric(FLERR,values[1],true,lmp);
+  h_molecule(nlocal) = utils::inumeric(FLERR,values[offset],true,lmp);
+  h_q(nlocal) = utils::numeric(FLERR,values[offset+1],true,lmp);
   h_num_bond(nlocal) = 0;
   h_num_angle(nlocal) = 0;
   h_num_dihedral(nlocal) = 0;

src/KOKKOS/atom_vec_full_kokkos.h

@@ -45,8 +45,8 @@ class AtomVecFullKokkos : public AtomVecKokkos {
   int pack_restart(int, double *);
   int unpack_restart(double *);
   void create_atom(int, double *);
-  void data_atom(double *, tagint, char **);
-  int data_atom_hybrid(int, char **);
+  void data_atom(double *, imageint, const std::vector<std::string> &);
+  int data_atom_hybrid(int, const std::vector<std::string> &, int);
   void pack_data(double **);
   int pack_data_hybrid(int, double *);
   void write_data(FILE *, int, double **);

src/KOKKOS/atom_vec_hybrid_kokkos.cpp

@@ -970,7 +970,8 @@ void AtomVecHybridKokkos::create_atom(int itype, double *coord)
    grow() occurs here so arrays for all sub-styles are grown
 ------------------------------------------------------------------------- */
 
-void AtomVecHybridKokkos::data_atom(double *coord, imageint imagetmp, char **values)
+void AtomVecHybridKokkos::data_atom(double *coord, imageint imagetmp,
+                                    const std::vector<std::string> &values)
 {
   atomKK->sync(Host,X_MASK|TAG_MASK|TYPE_MASK|IMAGE_MASK|MASK_MASK|V_MASK|OMEGA_MASK/*|ANGMOM_MASK*/);
 
@@ -1009,7 +1010,7 @@ void AtomVecHybridKokkos::data_atom(double *coord, imageint imagetmp, char **val
 
   int m = 5;
   for (int k = 0; k < nstyles; k++)
-    m += styles[k]->data_atom_hybrid(nlocal,&values[m]);
+    m += styles[k]->data_atom_hybrid(nlocal,values,m);
 
   atom->nlocal++;
 }
@@ -1018,21 +1019,21 @@ void AtomVecHybridKokkos::data_atom(double *coord, imageint imagetmp, char **val
    unpack one line from Velocities section of data file
 ------------------------------------------------------------------------- */
 
-void AtomVecHybridKokkos::data_vel(int m, char **values)
+void AtomVecHybridKokkos::data_vel(int m, const std::vector<std::string> &values)
 {
   atomKK->sync(Host,V_MASK);
 
-  h_v(m,0) = utils::numeric(FLERR,values[0],true,lmp);
-  h_v(m,1) = utils::numeric(FLERR,values[1],true,lmp);
-  h_v(m,2) = utils::numeric(FLERR,values[2],true,lmp);
+  int ivalue = 1;
+  h_v(m,0) = utils::numeric(FLERR,values[ivalue++],true,lmp);
+  h_v(m,1) = utils::numeric(FLERR,values[ivalue++],true,lmp);
+  h_v(m,2) = utils::numeric(FLERR,values[ivalue++],true,lmp);
 
   atomKK->modified(Host,V_MASK);
 
   // each sub-style parses sub-style specific values
 
-  int n = 3;
   for (int k = 0; k < nstyles; k++)
-    n += styles[k]->data_vel_hybrid(m,&values[n]);
+    ivalue += styles[k]->data_vel_hybrid(m,values,ivalue);
 }
 
 /* ----------------------------------------------------------------------

src/KOKKOS/atom_vec_hybrid_kokkos.h

@@ -57,9 +57,9 @@ class AtomVecHybridKokkos : public AtomVecKokkos {
   int pack_restart(int, double *);
   int unpack_restart(double *);
   void create_atom(int, double *);
-  void data_atom(double *, imageint, char **);
-  int data_atom_hybrid(int, char **) {return 0;}
-  void data_vel(int, char **);
+  void data_atom(double *, imageint, const std::vector<std::string> &);
+  int data_atom_hybrid(int, const std::vector<std::string> &, int) {return 0;}
+  void data_vel(int, const std::vector<std::string> &);
   void pack_data(double **);
   void write_data(FILE *, int, double **);
   void pack_vel(double **);

src/KOKKOS/atom_vec_kokkos.cpp

@@ -1016,12 +1016,13 @@ void AtomVecKokkos::unpack_reverse(int n, int *list, double *buf)
  *    unpack one line from Velocities section of data file
  *    ------------------------------------------------------------------------- */
 
-void AtomVecKokkos::data_vel(int m, char **values)
+void AtomVecKokkos::data_vel(int m, const std::vector<std::string> &values)
 {
   double **v = atom->v;
-  v[m][0] = utils::numeric(FLERR,values[0],true,lmp);
-  v[m][1] = utils::numeric(FLERR,values[1],true,lmp);
-  v[m][2] = utils::numeric(FLERR,values[2],true,lmp);
+  int ivalue = 1;
+  v[m][0] = utils::numeric(FLERR,values[ivalue++],true,lmp);
+  v[m][1] = utils::numeric(FLERR,values[ivalue++],true,lmp);
+  v[m][2] = utils::numeric(FLERR,values[ivalue++],true,lmp);
 
   atomKK->modified(Host,V_MASK);
 }

src/KOKKOS/atom_vec_kokkos.h

@@ -44,7 +44,7 @@ class AtomVecKokkos : public AtomVec {
   virtual void unpack_comm_vel(int, int, double *);
   virtual int pack_reverse(int, int, double *);
   virtual void unpack_reverse(int, int *, double *);
-  virtual void data_vel(int, char **);
+  virtual void data_vel(int, const std::vector<std::string> &);
   virtual void pack_vel(double **);
   virtual void write_vel(FILE *, int, double **);
 

src/KOKKOS/atom_vec_molecular_kokkos.cpp

@@ -1889,7 +1889,7 @@ void AtomVecMolecularKokkos::create_atom(int itype, double *coord)
 ------------------------------------------------------------------------- */
 
 void AtomVecMolecularKokkos::data_atom(double *coord, imageint imagetmp,
-                                       char **values)
+                                       const std::vector<std::string> &values)
 {
   int nlocal = atom->nlocal;
   if (nlocal == nmax) grow(0);
@@ -1924,9 +1924,10 @@ void AtomVecMolecularKokkos::data_atom(double *coord, imageint imagetmp,
    initialize other atom quantities for this sub-style
 ------------------------------------------------------------------------- */
 
-int AtomVecMolecularKokkos::data_atom_hybrid(int nlocal, char **values)
+int AtomVecMolecularKokkos::data_atom_hybrid(int nlocal, const std::vector<std::string> &values,
+                                             int offset)
 {
-  h_molecule(nlocal) = utils::inumeric(FLERR,values[0],true,lmp);
+  h_molecule(nlocal) = utils::inumeric(FLERR,values[offset],true,lmp);
   h_num_bond(nlocal) = 0;
   h_num_angle(nlocal) = 0;
   h_num_dihedral(nlocal) = 0;

src/KOKKOS/atom_vec_molecular_kokkos.h

@@ -51,8 +51,8 @@ class AtomVecMolecularKokkos : public AtomVecKokkos {
   int pack_restart(int, double *);
   int unpack_restart(double *);
   void create_atom(int, double *);
-  void data_atom(double *, tagint, char **);
-  int data_atom_hybrid(int, char **);
+  void data_atom(double *, imageint, const std::vector<std::string> &);
+  int data_atom_hybrid(int, const std::vector<std::string> &, int);
   void pack_data(double **);
   int pack_data_hybrid(int, double *);
   void write_data(FILE *, int, double **);

src/KOKKOS/atom_vec_sphere_kokkos.cpp

@@ -2543,7 +2543,8 @@ void AtomVecSphereKokkos::create_atom(int itype, double *coord)
    initialize other atom quantities
 ------------------------------------------------------------------------- */
 
-void AtomVecSphereKokkos::data_atom(double *coord, imageint imagetmp, char **values)
+void AtomVecSphereKokkos::data_atom(double *coord, imageint imagetmp,
+                                    const std::vector<std::string> &values)
 {
   int nlocal = atom->nlocal;
   if (nlocal == nmax) grow(0);
@@ -2590,13 +2591,14 @@ void AtomVecSphereKokkos::data_atom(double *coord, imageint imagetmp, char **val
    initialize other atom quantities for this sub-style
 ------------------------------------------------------------------------- */
 
-int AtomVecSphereKokkos::data_atom_hybrid(int nlocal, char **values)
+int AtomVecSphereKokkos::data_atom_hybrid(int nlocal, const std::vector<std::string> &values,
+                                          int offset)
 {
-  radius[nlocal] = 0.5 * utils::numeric(FLERR,values[0],true,lmp);
+  radius[nlocal] = 0.5 * utils::numeric(FLERR,values[offset],true,lmp);
   if (radius[nlocal] < 0.0)
     error->one(FLERR,"Invalid radius in Atoms section of data file");
 
-  double density = utils::numeric(FLERR,values[1],true,lmp);
+  double density = utils::numeric(FLERR,values[offset+1],true,lmp);
   if (density <= 0.0)
     error->one(FLERR,"Invalid density in Atoms section of data file");
 
@@ -2615,15 +2617,16 @@ int AtomVecSphereKokkos::data_atom_hybrid(int nlocal, char **values)
    unpack one line from Velocities section of data file
 ------------------------------------------------------------------------- */
 
-void AtomVecSphereKokkos::data_vel(int m, char **values)
+void AtomVecSphereKokkos::data_vel(int m, const std::vector<std::string> &values)
 {
+  int ivalue = 1;
   atomKK->sync(Host,V_MASK|OMEGA_MASK);
-  h_v(m,0) = utils::numeric(FLERR,values[0],true,lmp);
-  h_v(m,1) = utils::numeric(FLERR,values[1],true,lmp);
-  h_v(m,2) = utils::numeric(FLERR,values[2],true,lmp);
-  h_omega(m,0) = utils::numeric(FLERR,values[3],true,lmp);
-  h_omega(m,1) = utils::numeric(FLERR,values[4],true,lmp);
-  h_omega(m,2) = utils::numeric(FLERR,values[5],true,lmp);
+  h_v(m,0) = utils::numeric(FLERR,values[ivalue++],true,lmp);
+  h_v(m,1) = utils::numeric(FLERR,values[ivalue++],true,lmp);
+  h_v(m,2) = utils::numeric(FLERR,values[ivalue++],true,lmp);
+  h_omega(m,0) = utils::numeric(FLERR,values[ivalue++],true,lmp);
+  h_omega(m,1) = utils::numeric(FLERR,values[ivalue++],true,lmp);
+  h_omega(m,2) = utils::numeric(FLERR,values[ivalue++],true,lmp);
   atomKK->modified(Host,V_MASK|OMEGA_MASK);
 }
 
@@ -2631,12 +2634,13 @@ void AtomVecSphereKokkos::data_vel(int m, char **values)
    unpack hybrid quantities from one line in Velocities section of data file
 ------------------------------------------------------------------------- */
 
-int AtomVecSphereKokkos::data_vel_hybrid(int m, char **values)
+int AtomVecSphereKokkos::data_vel_hybrid(int m, const std::vector<std::string> &values,
+                                         int offset)
 {
   atomKK->sync(Host,OMEGA_MASK);
-  omega[m][0] = utils::numeric(FLERR,values[0],true,lmp);
-  omega[m][1] = utils::numeric(FLERR,values[1],true,lmp);
-  omega[m][2] = utils::numeric(FLERR,values[2],true,lmp);
+  omega[m][0] = utils::numeric(FLERR,values[offset],true,lmp);
+  omega[m][1] = utils::numeric(FLERR,values[offset+1],true,lmp);
+  omega[m][2] = utils::numeric(FLERR,values[offset+2],true,lmp);
   atomKK->modified(Host,OMEGA_MASK);
   return 3;
 }

src/KOKKOS/atom_vec_sphere_kokkos.h

@@ -58,10 +58,10 @@ class AtomVecSphereKokkos : public AtomVecKokkos {
   int pack_restart(int, double *);
   int unpack_restart(double *);
   void create_atom(int, double *);
-  void data_atom(double *, imageint, char **);
-  int data_atom_hybrid(int, char **);
-  void data_vel(int, char **);
-  int data_vel_hybrid(int, char **);
+  void data_atom(double *, imageint, const std::vector<std::string> &);
+  int data_atom_hybrid(int, const std::vector<std::string> &, int);
+  void data_vel(int, const std::vector<std::string> &);
+  int data_vel_hybrid(int, const std::vector<std::string> &, int);
   void pack_data(double **);
   int pack_data_hybrid(int, double *);
   void write_data(FILE *, int, double **);

src/KOKKOS/atom_vec_spin_kokkos.cpp

@@ -1056,7 +1056,7 @@ void AtomVecSpinKokkos::create_atom(int itype, double *coord)
 ------------------------------------------------------------------------- */
 
 void AtomVecSpinKokkos::data_atom(double *coord, imageint imagetmp,
-                                    char **values)
+                                  const std::vector<std::string> &values)
 {
   int nlocal = atom->nlocal;
   if (nlocal == nmax) grow(0);
@@ -1098,12 +1098,13 @@ void AtomVecSpinKokkos::data_atom(double *coord, imageint imagetmp,
    initialize other atom quantities for this sub-style
 ------------------------------------------------------------------------- */
 
-int AtomVecSpinKokkos::data_atom_hybrid(int nlocal, char **values)
+int AtomVecSpinKokkos::data_atom_hybrid(int nlocal, const std::vector<std::string> &values,
+                                        int offset)
 {
-  h_sp(nlocal,3) = utils::numeric(FLERR,values[0],true,lmp);
-  h_sp(nlocal,0) = utils::numeric(FLERR,values[1],true,lmp);
-  h_sp(nlocal,1) = utils::numeric(FLERR,values[2],true,lmp);
-  h_sp(nlocal,2) = utils::numeric(FLERR,values[3],true,lmp);
+  h_sp(nlocal,3) = utils::numeric(FLERR,values[offset],true,lmp);
+  h_sp(nlocal,0) = utils::numeric(FLERR,values[offset+1],true,lmp);
+  h_sp(nlocal,1) = utils::numeric(FLERR,values[offset+2],true,lmp);
+  h_sp(nlocal,2) = utils::numeric(FLERR,values[offset+3],true,lmp);
   double inorm = 1.0/sqrt(sp[nlocal][0]*sp[nlocal][0] +
                           sp[nlocal][1]*sp[nlocal][1] +
                           sp[nlocal][2]*sp[nlocal][2]);

src/KOKKOS/atom_vec_spin_kokkos.h

@@ -45,8 +45,8 @@ class AtomVecSpinKokkos : public AtomVecKokkos {
   int pack_restart(int, double *);
   int unpack_restart(double *);
   void create_atom(int, double *);
-  void data_atom(double *, imageint, char **);
-  int data_atom_hybrid(int, char **);
+  void data_atom(double *, imageint, const std::vector<std::string> &);
+  int data_atom_hybrid(int, const std::vector<std::string> &, int);
   void pack_data(double **);
   int pack_data_hybrid(int, double *);
   void write_data(FILE *, int, double **);

src/KOKKOS/compute_ave_sphere_atom_kokkos.cpp

@@ -0,0 +1,209 @@
+/* ----------------------------------------------------------------------
+   LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
+   https://www.lammps.org/, Sandia National Laboratories
+   Steve Plimpton, sjplimp@sandia.gov
+
+   Copyright (2003) Sandia Corporation.  Under the terms of Contract
+   DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains
+   certain rights in this software.  This software is distributed under
+   the GNU General Public License.
+
+   See the README file in the top-level LAMMPS directory.
+------------------------------------------------------------------------- */
+
+#include "compute_ave_sphere_atom_kokkos.h"
+
+#include "atom_kokkos.h"
+#include "atom_masks.h"
+#include "comm.h"
+#include "error.h"
+#include "force.h"
+#include "memory_kokkos.h"
+#include "modify.h"
+#include "neigh_list.h"
+#include "neigh_request.h"
+#include "neighbor_kokkos.h"
+#include "pair.h"
+#include "update.h"
+#include "math_const.h"
+
+#include <cmath>
+#include <cstring>
+
+using namespace LAMMPS_NS;
+using namespace MathConst;
+
+
+/* ---------------------------------------------------------------------- */
+
+template<class DeviceType>
+ComputeAveSphereAtomKokkos<DeviceType>::ComputeAveSphereAtomKokkos(LAMMPS *lmp, int narg, char **arg) :
+  ComputeAveSphereAtom(lmp, narg, arg)
+{
+  kokkosable = 1;
+  atomKK = (AtomKokkos *) atom;
+  execution_space = ExecutionSpaceFromDevice<DeviceType>::space;
+  datamask_read = EMPTY_MASK;
+  datamask_modify = EMPTY_MASK;
+}
+
+/* ---------------------------------------------------------------------- */
+
+template<class DeviceType>
+ComputeAveSphereAtomKokkos<DeviceType>::~ComputeAveSphereAtomKokkos()
+{
+  if (copymode) return;
+
+  memoryKK->destroy_kokkos(k_result,result);
+}
+
+/* ---------------------------------------------------------------------- */
+
+template<class DeviceType>
+void ComputeAveSphereAtomKokkos<DeviceType>::init()
+{
+  ComputeAveSphereAtom::init();
+
+  // need an occasional full neighbor list
+
+  // irequest = neigh request made by parent class
+
+  int irequest = neighbor->nrequest - 1;
+
+  neighbor->requests[irequest]->
+    kokkos_host = std::is_same<DeviceType,LMPHostType>::value &&
+    !std::is_same<DeviceType,LMPDeviceType>::value;
+  neighbor->requests[irequest]->
+    kokkos_device = std::is_same<DeviceType,LMPDeviceType>::value;
+}
+
+/* ---------------------------------------------------------------------- */
+
+template<class DeviceType>
+void ComputeAveSphereAtomKokkos<DeviceType>::compute_peratom()
+{
+  invoked_peratom = update->ntimestep;
+
+  // grow result array if necessary
+
+  if (atom->nmax > nmax) {
+    memoryKK->destroy_kokkos(k_result,result);
+    nmax = atom->nmax;
+    memoryKK->create_kokkos(k_result,result,nmax,2,"ave/sphere/atom:result");
+    d_result = k_result.view<DeviceType>();
+    array_atom = result;
+  }
+
+  // need velocities of ghost atoms
+
+  atomKK->sync(Host,V_MASK);
+  comm->forward_comm_compute(this);
+  atomKK->modified(Host,V_MASK);
+
+  // invoke full neighbor list (will copy or build if necessary)
+
+  neighbor->build_one(list);
+  int inum = list->inum;
+
+  NeighListKokkos<DeviceType>* k_list = static_cast<NeighListKokkos<DeviceType>*>(list);
+  d_numneigh = k_list->d_numneigh;
+  d_neighbors = k_list->d_neighbors;
+  d_ilist = k_list->d_ilist;
+
+  // compute properties for each atom in group
+  // use full neighbor list to count atoms less than cutoff
+
+  atomKK->sync(execution_space,X_MASK|V_MASK|TYPE_MASK|MASK_MASK);
+  x = atomKK->k_x.view<DeviceType>();
+  v = atomKK->k_v.view<DeviceType>();
+  mask = atomKK->k_mask.view<DeviceType>();
+
+  Kokkos::deep_copy(d_result,0.0);
+
+  copymode = 1;
+  typename Kokkos::RangePolicy<DeviceType, TagComputeAveSphereAtom> policy(0,inum);
+  Kokkos::parallel_for("ComputeAveSphereAtom",policy,*this);
+  copymode = 0;
+
+  k_result.modify<DeviceType>();
+  k_result.sync_host();
+}
+
+template<class DeviceType>
+KOKKOS_INLINE_FUNCTION
+void ComputeAveSphereAtomKokkos<DeviceType>::operator()(TagComputeAveSphereAtom, const int &ii) const
+{
+  const int i = d_ilist[ii];
+  if (mask[i] & groupbit) {
+    const X_FLOAT xtmp = x(i,0);
+    const X_FLOAT ytmp = x(i,1);
+    const X_FLOAT ztmp = x(i,2);
+    const int jnum = d_numneigh[i];
+
+    // i atom contribution
+
+    int count = 1;
+    double vsum[3];
+    vsum[0] = v(i,0);
+    vsum[1] = v(i,1);
+    vsum[2] = v(i,2);
+
+    for (int jj = 0; jj < jnum; jj++) {
+      int j = d_neighbors(i,jj);
+      j &= NEIGHMASK;
+
+      const F_FLOAT delx = x(j,0) - xtmp;
+      const F_FLOAT dely = x(j,1) - ytmp;
+      const F_FLOAT delz = x(j,2) - ztmp;
+      const F_FLOAT rsq = delx*delx + dely*dely + delz*delz;
+      if (rsq < cutsq) {
+        count++;
+        vsum[0] += v(j,0);
+        vsum[1] += v(j,1);
+        vsum[2] += v(j,2);
+      }
+    }
+
+    double vavg[3];
+    vavg[0] = vsum[0]/count;
+    vavg[1] = vsum[1]/count;
+    vavg[2] = vsum[2]/count;
+
+    // i atom contribution
+
+    count = 1;
+    double vnet[3];
+    vnet[0] = v(i,0) - vavg[0];
+    vnet[1] = v(i,1) - vavg[1];
+    vnet[2] = v(i,2) - vavg[2];
+    double ke_sum = vnet[0]*vnet[0] + vnet[1]*vnet[1] + vnet[2]*vnet[2];
+
+    for (int jj = 0; jj < jnum; jj++) {
+      int j = d_neighbors(i,jj);
+      j &= NEIGHMASK;
+
+      const F_FLOAT delx = x(j,0) - xtmp;
+      const F_FLOAT dely = x(j,1) - ytmp;
+      const F_FLOAT delz = x(j,2) - ztmp;
+      const F_FLOAT rsq = delx*delx + dely*dely + delz*delz;
+      if (rsq < cutsq) {
+        count++;
+        vnet[0] = v(j,0) - vavg[0];
+        vnet[1] = v(j,1) - vavg[1];
+        vnet[2] = v(j,2) - vavg[2];
+        ke_sum += vnet[0]*vnet[0] + vnet[1]*vnet[1] + vnet[2]*vnet[2];
+      }
+    }
+    double density = count/sphere_vol;
+    double temp = ke_sum/3.0/count;
+    d_result(i,0) = density;
+    d_result(i,1) = temp;
+  }
+}
+
+namespace LAMMPS_NS {
+template class ComputeAveSphereAtomKokkos<LMPDeviceType>;
+#ifdef LMP_KOKKOS_GPU
+template class ComputeAveSphereAtomKokkos<LMPHostType>;
+#endif
+}

src/KOKKOS/compute_ave_sphere_atom_kokkos.h

@@ -0,0 +1,66 @@
+/* -*- c++ -*- ----------------------------------------------------------
+   LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
+   https://www.lammps.org/, Sandia National Laboratories
+   Steve Plimpton, sjplimp@sandia.gov
+
+   Copyright (2003) Sandia Corporation.  Under the terms of Contract
+   DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains
+   certain rights in this software.  This software is distributed under
+   the GNU General Public License.
+
+   See the README file in the top-level LAMMPS directory.
+------------------------------------------------------------------------- */
+
+#ifdef COMPUTE_CLASS
+
+ComputeStyle(ave/sphere/atom/kk,ComputeAveSphereAtomKokkos<LMPDeviceType>)
+ComputeStyle(ave/sphere/atom/kk/device,ComputeAveSphereAtomKokkos<LMPDeviceType>)
+ComputeStyle(ave/sphere/atom/kk/host,ComputeAveSphereAtomKokkos<LMPHostType>)
+
+#else
+
+#ifndef LMP_COMPUTE_AVE_SPHERE_ATOM_KOKKOS_H
+#define LMP_COMPUTE_AVE_SPHERE_ATOM_KOKKOS_H
+
+#include "compute_ave_sphere_atom.h"
+#include "kokkos_type.h"
+
+namespace LAMMPS_NS {
+
+struct TagComputeAveSphereAtom{};
+
+template<class DeviceType>
+class ComputeAveSphereAtomKokkos : public ComputeAveSphereAtom {
+ public:
+  typedef DeviceType device_type;
+  typedef ArrayTypes<DeviceType> AT;
+
+  ComputeAveSphereAtomKokkos(class LAMMPS *, int, char **);
+  virtual ~ComputeAveSphereAtomKokkos();
+  void init();
+  void compute_peratom();
+
+  KOKKOS_INLINE_FUNCTION
+  void operator()(TagComputeAveSphereAtom, const int&) const;
+
+ private:
+  typename AT::t_x_array_randomread x;
+  typename AT::t_v_array_randomread v;
+  typename ArrayTypes<DeviceType>::t_int_1d mask;
+
+  typename AT::t_neighbors_2d d_neighbors;
+  typename AT::t_int_1d_randomread d_ilist;
+  typename AT::t_int_1d_randomread d_numneigh;
+
+  DAT::tdual_float_2d k_result;
+  typename AT::t_float_2d d_result;
+};
+
+}
+
+#endif
+#endif
+
+/* ERROR/WARNING messages:
+
+*/

src/MANYBODY/pair_eam_cd.cpp

@@ -20,9 +20,6 @@
 
 #include "pair_eam_cd.h"
 
-#include <cmath>
-
-#include <cstring>
 #include "atom.h"
 #include "force.h"
 #include "comm.h"
@@ -31,11 +28,11 @@
 #include "error.h"
 #include "tokenizer.h"
 
-
+#include <cmath>
+#include <cstring>
 
 using namespace LAMMPS_NS;
 
-#define ASSERT(cond)
 #define MAXLINE 1024        // This sets the maximum line length in EAM input files.
 
 PairEAMCD::PairEAMCD(LAMMPS *lmp, int _cdeamVersion)
@@ -298,7 +295,7 @@ void PairEAMCD::compute(int eflag, int vflag)
     // It will be replaced by the concentration at site i if atom i is either A or B.
 
     double x_i = -1.0;
-    double D_i, h_prime_i;
+    double D_i = 0.0, h_prime_i;
 
     // This if-clause is only required for ternary alloys.
 
@@ -307,7 +304,6 @@ void PairEAMCD::compute(int eflag, int vflag)
       // Compute local concentration at site i.
 
       x_i = rhoB[i]/rho[i];
-      ASSERT(x_i >= 0 && x_i<=1.0);
 
       if (cdeamVersion == 1) {
 
@@ -317,8 +313,6 @@ void PairEAMCD::compute(int eflag, int vflag)
         D_i = D_values[i] * h_prime_i / (2.0 * rho[i] * rho[i]);
       } else if (cdeamVersion == 2) {
         D_i = D_values[i];
-      } else {
-        ASSERT(false);
       }
     }
 
@@ -354,14 +348,11 @@ void PairEAMCD::compute(int eflag, int vflag)
 
         // This code line is required for ternary alloy.
 
-        if (jtype == speciesA || jtype == speciesB) {
-          ASSERT(rho[i] != 0.0);
-          ASSERT(rho[j] != 0.0);
+        if ((jtype == speciesA || jtype == speciesB) && rho[j] != 0.0) {
 
           // Compute local concentration at site j.
 
           x_j = rhoB[j]/rho[j];
-          ASSERT(x_j >= 0 && x_j<=1.0);
 
           double D_j=0.0;
           if (cdeamVersion == 1) {
@@ -372,8 +363,6 @@ void PairEAMCD::compute(int eflag, int vflag)
             D_j = D_values[j] * h_prime_j / (2.0 * rho[j] * rho[j]);
           } else if (cdeamVersion == 2) {
             D_j = D_values[j];
-          } else {
-            ASSERT(false);
           }
           double t2 = -rhoB[j];
           if (itype == speciesB) t2 += rho[j];
@@ -422,8 +411,6 @@ void PairEAMCD::compute(int eflag, int vflag)
             // Calculate h(x_ij) polynomial function.
 
             h = evalH(x_ij);
-          } else {
-            ASSERT(false);
           }
           fpair += h * phip;
           phi *= h;
@@ -460,7 +447,8 @@ void PairEAMCD::coeff(int narg, char **arg)
   // Make sure the EAM file is a CD-EAM binary alloy.
 
   if (setfl->nelements < 2)
-    error->all(FLERR,"The EAM file must contain at least 2 elements to be used with the eam/cd pair style.");
+    error->all(FLERR,"The EAM file must contain at least 2 elements to be "
+                    "used with the eam/cd pair style.");
 
   // Read in the coefficients of the h polynomial from the end of the EAM file.
 
@@ -502,22 +490,28 @@ void PairEAMCD::read_h_coeff(char *filename)
     // Open potential file
 
     FILE *fptr;
-    char line[MAXLINE];
-    char nextline[MAXLINE];
     int convert_flag = unit_convert_flag;
     fptr = utils::open_potential(filename, lmp, &convert_flag);
     if (fptr == nullptr)
-      error->one(FLERR,"Cannot open EAMCD potential file {}",
-                                   filename);
+      error->one(FLERR,"Cannot open EAMCD potential file {}", filename);
 
     // h coefficients are stored at the end of the file.
-    // Skip to last line of file.
+    // Seek to end of file, read last part into a buffer and
+    // then skip over lines in buffer until reaching the end.
 
-    while (fgets(nextline, MAXLINE, fptr) != nullptr) {
-      strcpy(line, nextline);
-    }
+    platform::fseek(fptr, platform::END_OF_FILE);
+    platform::fseek(fptr, platform::ftell(fptr) - MAXLINE);
+    char *buf = new char[MAXLINE+1];
+    fread(buf, 1, MAXLINE, fptr);
+    buf[MAXLINE] = '\0';        // must 0-terminate buffer for string processing
+    Tokenizer lines(buf, "\n");
+    delete[] buf;
 
-    ValueTokenizer values(line);
+    std::string lastline;
+    while (lines.has_next())
+      lastline = lines.next();
+
+    ValueTokenizer values(lastline);
     int degree = values.next_int();
     nhcoeff = degree+1;
 
@@ -527,10 +521,8 @@ void PairEAMCD::read_h_coeff(char *filename)
     delete[] hcoeff;
     hcoeff = new double[nhcoeff];
 
-    int i = 0;
-    while (values.has_next()) {
-      hcoeff[i++] = values.next_double();
-    }
+    for (int i = 0; i < nhcoeff; ++i)
+      hcoeff[i] = values.next_double();
 
     // Close the potential file.
 
@@ -545,7 +537,6 @@ void PairEAMCD::read_h_coeff(char *filename)
   MPI_Bcast(hcoeff, nhcoeff, MPI_DOUBLE, 0, world);
 }
 
-
 /* ---------------------------------------------------------------------- */
 
 int PairEAMCD::pack_forward_comm(int n, int *list, double *buf,
@@ -572,7 +563,7 @@ int PairEAMCD::pack_forward_comm(int n, int *list, double *buf,
         buf[m++] = rhoB[j];
       }
       return m;
-    } else { ASSERT(false); return 0; }
+    } else return 0;
   } else if (communicationStage == 4) {
     for (i = 0; i < n; i++) {
       j = list[i];
@@ -604,8 +595,6 @@ void PairEAMCD::unpack_forward_comm(int n, int first, double *buf)
         rho[i] = buf[m++];
         rhoB[i] = buf[m++];
       }
-    } else {
-      ASSERT(false);
     }
   } else if (communicationStage == 4) {
     for (i = first; i < last; i++) {
@@ -636,7 +625,7 @@ int PairEAMCD::pack_reverse_comm(int n, int first, double *buf)
         buf[m++] = rhoB[i];
       }
       return m;
-    } else { ASSERT(false); return 0; }
+    } else return 0;
   } else if (communicationStage == 3) {
     for (i = first; i < last; i++) {
       buf[m++] = D_values[i];
@@ -666,8 +655,6 @@ void PairEAMCD::unpack_reverse_comm(int n, int *list, double *buf)
         rho[j] += buf[m++];
         rhoB[j] += buf[m++];
       }
-    } else {
-      ASSERT(false);
     }
   } else if (communicationStage == 3) {
     for (i = 0; i < n; i++) {

src/MANYBODY/pair_eam_cd.h

@@ -120,6 +120,7 @@ class PairEAMCD : public PairEAMAlloy {
     index.p = r * rdr + 1.0;
     index.m = static_cast<int>(index.p);
     index.m = index.m <= (nr - 1) ? index.m : (nr - 1);
+    index.m = index.m > 1 ? index.m : 1;
     index.p -= index.m;
     index.p = index.p <= 1.0 ? index.p : 1.0;
     return index;
@@ -132,6 +133,7 @@ class PairEAMCD : public PairEAMAlloy {
     index.p = rho * rdrho + 1.0;
     index.m = static_cast<int>(index.p);
     index.m = index.m <= (nrho - 1) ? index.m : (nrho - 1);
+    index.m = index.m > 1 ? index.m : 1;
     index.p -= index.m;
     index.p = index.p <= 1.0 ? index.p : 1.0;
     return index;

src/MC/fix_gcmc.cpp

@@ -504,8 +504,7 @@ void FixGCMC::init()
     int flagall;
     MPI_Allreduce(&flag,&flagall,1,MPI_INT,MPI_SUM,world);
     if (flagall && comm->me == 0)
-      error->all(FLERR,
-       "Fix gcmc cannot exchange individual atoms belonging to a molecule");
+      error->all(FLERR, "Fix gcmc cannot exchange individual atoms belonging to a molecule");
   }
 
   // if molecules are exchanged or moved, check for unset mol IDs
@@ -520,16 +519,13 @@ void FixGCMC::init()
     int flagall;
     MPI_Allreduce(&flag,&flagall,1,MPI_INT,MPI_SUM,world);
     if (flagall && comm->me == 0)
-      error->all(FLERR,
-       "All mol IDs should be set for fix gcmc group atoms");
+      error->all(FLERR, "All mol IDs should be set for fix gcmc group atoms");
   }
 
   if (exchmode == EXCHMOL || movemode == MOVEMOL)
     if (atom->molecule_flag == 0 || !atom->tag_enable
         || (atom->map_style == Atom::MAP_NONE))
-      error->all(FLERR,
-       "Fix gcmc molecule command requires that "
-       "atoms have molecule attributes");
+      error->all(FLERR, "Fix gcmc molecule command requires that atoms have molecule attributes");
 
   // if rigidflag defined, check for rigid/small fix
   // its molecule template must be same as this one
@@ -541,9 +537,7 @@ void FixGCMC::init()
     fixrigid = modify->fix[ifix];
     int tmp;
     if (&onemols[imol] != (Molecule **) fixrigid->extract("onemol",tmp))
-      error->all(FLERR,
-                 "Fix gcmc and fix rigid/small not using "
-                 "same molecule template ID");
+      error->all(FLERR, "Fix gcmc and fix rigid/small not using same molecule template ID");
   }
 
   // if shakeflag defined, check for SHAKE fix

src/MC/fix_widom.cpp

@@ -310,16 +310,13 @@ void FixWidom::init()
     int flagall;
     MPI_Allreduce(&flag,&flagall,1,MPI_INT,MPI_SUM,world);
     if (flagall && comm->me == 0)
-      error->all(FLERR,
-       "All mol IDs should be set for fix widom group atoms");
+      error->all(FLERR, "All mol IDs should be set for fix widom group atoms");
   }
 
   if (exchmode == EXCHMOL)
     if (atom->molecule_flag == 0 || !atom->tag_enable
         || (atom->map_style == Atom::MAP_NONE))
-      error->all(FLERR,
-       "Fix widom molecule command requires that "
-       "atoms have molecule attributes");
+      error->all(FLERR, "Fix widom molecule command requires that atoms have molecule attributes");
 
   if (domain->dimension == 2)
     error->all(FLERR,"Cannot use fix widom in a 2d simulation");

src/ML-IAP/mliap_descriptor_so3.cpp

@@ -125,28 +125,27 @@ void MLIAPDescriptorSO3::read_paramfile(char *paramfilename)
 
     // check for keywords with one value per element
 
-    if (strcmp(skeywd.c_str(), "elems") == 0 || strcmp(skeywd.c_str(), "radelems") == 0 ||
-        strcmp(skeywd.c_str(), "welems") == 0) {
+    if ((skeywd == "elems") || (skeywd == "radelems") || (skeywd == "welems"))  {
 
       if (nelementsflag == 0 || nwords != nelements + 1)
         error->all(FLERR, "Incorrect SO3 parameter file");
 
-      if (strcmp(skeywd.c_str(), "elems") == 0) {
+      if (skeywd == "elems") {
         for (int ielem = 0; ielem < nelements; ielem++) {
           elements[ielem] = utils::strdup(skeyval);
           if (ielem < nelements - 1) skeyval = p.next();
         }
 
         elementsflag = 1;
-      } else if (strcmp(skeywd.c_str(), "radelems") == 0) {
+      } else if (skeywd == "radelems")  {
         for (int ielem = 0; ielem < nelements; ielem++) {
-          radelem[ielem] = utils::numeric(FLERR, skeyval.c_str(), false, lmp);
+          radelem[ielem] = utils::numeric(FLERR, skeyval, false, lmp);
           if (ielem < nelements - 1) skeyval = p.next();
         }
         radelemflag = 1;
-      } else if (strcmp(skeywd.c_str(), "welems") == 0) {
+      } else if (skeywd == "welems") {
         for (int ielem = 0; ielem < nelements; ielem++) {
-          wjelem[ielem] = utils::numeric(FLERR, skeyval.c_str(), false, lmp);
+          wjelem[ielem] = utils::numeric(FLERR, skeyval, false, lmp);
           if (ielem < nelements - 1) skeyval = p.next();
         }
         wjelemflag = 1;
@@ -158,23 +157,23 @@ void MLIAPDescriptorSO3::read_paramfile(char *paramfilename)
 
       if (nwords != 2) error->all(FLERR, "Incorrect SO3 parameter file");
 
-      if (strcmp(skeywd.c_str(), "nelems") == 0) {
-        nelements = utils::inumeric(FLERR, skeyval.c_str(), false, lmp);
+      if (skeywd == "nelems") {
+        nelements = utils::inumeric(FLERR, skeyval, false, lmp);
         elements = new char *[nelements];
         memory->create(radelem, nelements, "mliap_so3_descriptor:radelem");
         memory->create(wjelem, nelements, "mliap_so3_descriptor:wjelem");
         nelementsflag = 1;
-      } else if (strcmp(skeywd.c_str(), "rcutfac") == 0) {
-        rcutfac = utils::numeric(FLERR, skeyval.c_str(), false, lmp);
+      } else if (skeywd == "rcutfac") {
+        rcutfac = utils::numeric(FLERR, skeyval, false, lmp);
         rcutfacflag = 1;
-      } else if (strcmp(skeywd.c_str(), "nmax") == 0) {
-        nmax = utils::inumeric(FLERR, skeyval.c_str(), false, lmp);
+      } else if (skeywd == "nmax") {
+        nmax = utils::inumeric(FLERR, skeyval, false, lmp);
         nmaxflag = 1;
-      } else if (strcmp(skeywd.c_str(), "lmax") == 0) {
-        lmax = utils::inumeric(FLERR, skeyval.c_str(), false, lmp);
+      } else if (skeywd == "lmax") {
+        lmax = utils::inumeric(FLERR, skeyval, false, lmp);
         lmaxflag = 1;
-      } else if (strcmp(skeywd.c_str(), "alpha") == 0) {
-        alpha = utils::numeric(FLERR, skeyval.c_str(), false, lmp);
+      } else if (skeywd == "alpha") {
+        alpha = utils::numeric(FLERR, skeyval, false, lmp);
         alphaflag = 1;
       } else
         error->all(FLERR, "Incorrect SO3 parameter file");

src/ML-RANN/pair_rann.cpp

@@ -440,7 +440,7 @@ void PairRANN::read_mass(const std::vector<std::string> &line1, const std::vecto
   if (nelements == -1)error->one(filename,linenum-1,"atom types must be defined before mass in potential file.");
   for (int i=0;i<nelements;i++) {
     if (line1[1].compare(elements[i])==0) {
-      mass[i]=utils::numeric(filename,linenum,line2[0].c_str(),true,lmp);
+      mass[i]=utils::numeric(filename,linenum,line2[0],true,lmp);
       return;
     }
   }
@@ -452,7 +452,7 @@ void PairRANN::read_fpe(std::vector<std::string> line,std::vector<std::string> l
   if (nelements == -1)error->one(filename,linenum-1,"atom types must be defined before fingerprints per element in potential file.");
   for (i=0;i<nelementsp;i++) {
     if (line[1].compare(elementsp[i])==0) {
-      fingerprintperelement[i] = utils::inumeric(filename,linenum,line1[0].c_str(),true,lmp);
+      fingerprintperelement[i] = utils::inumeric(filename,linenum,line1[0],true,lmp);
       fingerprints[i] = new RANN::Fingerprint *[fingerprintperelement[i]];
       for (int j=0;j<fingerprintperelement[i];j++) {
         fingerprints[i][j]=new RANN::Fingerprint(this);
@@ -491,7 +491,7 @@ void PairRANN::read_fingerprints(std::vector<std::string> line,std::vector<std::
     fingerprints[i][i1] = create_fingerprint(line1[k].c_str());
     if (fingerprints[i][i1]->n_body_type!=nwords-1) {error->one(filename,linenum,"invalid fingerprint for element combination");}
     k++;
-    fingerprints[i][i1]->init(atomtypes,utils::inumeric(filename,linenum,line1[k++].c_str(),true,lmp));
+    fingerprints[i][i1]->init(atomtypes,utils::inumeric(filename,linenum,line1[k++],true,lmp));
     fingerprintcount[i]++;
   }
   delete[] atomtypes;
@@ -523,7 +523,7 @@ void PairRANN::read_fingerprint_constants(std::vector<std::string> line,std::vec
     for (j=0;j<n_body_type;j++) {
       if (fingerprints[i][k]->atomtypes[j]!=atomtypes[j]) {break;}
       if (j==n_body_type-1) {
-        if (line[nwords-3].compare(fingerprints[i][k]->style)==0 && utils::inumeric(filename,linenum,line[nwords-2].c_str(),true,lmp)==fingerprints[i][k]->id) {
+        if (line[nwords-3].compare(fingerprints[i][k]->style)==0 && utils::inumeric(filename,linenum,line[nwords-2],true,lmp)==fingerprints[i][k]->id) {
           found=true;
           i1 = k;
           break;
@@ -542,7 +542,7 @@ void PairRANN::read_network_layers(std::vector<std::string> line,std::vector<std
   if (nelements == -1)error->one(filename,linenum-1,"atom types must be defined before network layers in potential file.");
   for (i=0;i<nelements;i++) {
     if (line[1].compare(elements[i])==0) {
-      net[i].layers = utils::inumeric(filename,linenum,line1[0].c_str(),true,lmp);
+      net[i].layers = utils::inumeric(filename,linenum,line1[0],true,lmp);
       if (net[i].layers < 1)error->one(filename,linenum,"invalid number of network layers");
       delete[] net[i].dimensions;
       weightdefined[i] = new bool [net[i].layers];
@@ -570,9 +570,9 @@ void PairRANN::read_layer_size(std::vector<std::string> line,std::vector<std::st
   for (i=0;i<nelements;i++) {
     if (line[1].compare(elements[i])==0) {
       if (net[i].layers==0)error->one(filename,linenum-1,"networklayers for each atom type must be defined before the corresponding layer sizes.");
-      int j = utils::inumeric(filename,linenum,line[2].c_str(),true,lmp);
+      int j = utils::inumeric(filename,linenum,line[2],true,lmp);
       if (j>=net[i].layers || j<0) {error->one(filename,linenum,"invalid layer in layer size definition");};
-      net[i].dimensions[j]= utils::inumeric(filename,linenum,line1[0].c_str(),true,lmp);
+      net[i].dimensions[j]= utils::inumeric(filename,linenum,line1[0],true,lmp);
       return;
     }
   }
@@ -587,7 +587,7 @@ void PairRANN::read_weight(std::vector<std::string> line,std::vector<std::string
   for (l=0;l<nelements;l++) {
     if (line[1].compare(elements[l])==0) {
       if (net[l].layers==0)error->one(filename,*linenum-1,"networklayers must be defined before weights.");
-      i=utils::inumeric(filename,*linenum,line[2].c_str(),true,lmp);
+      i=utils::inumeric(filename,*linenum,line[2],true,lmp);
       if (i>=net[l].layers || i<0)error->one(filename,*linenum-1,"invalid weight layer");
       if (net[l].dimensions[i]==0 || net[l].dimensions[i+1]==0) error->one(filename,*linenum-1,"network layer sizes must be defined before corresponding weight");
       net[l].Weights[i] = new double[net[l].dimensions[i]*net[l].dimensions[i+1]];
@@ -595,7 +595,7 @@ void PairRANN::read_weight(std::vector<std::string> line,std::vector<std::string
       nwords = line1.size();
       if (nwords != net[l].dimensions[i])error->one(filename,*linenum,"invalid weights per line");
       for (k=0;k<net[l].dimensions[i];k++) {
-        net[l].Weights[i][k] = utils::numeric(filename,*linenum,line1[k].c_str(),true,lmp);
+        net[l].Weights[i][k] = utils::numeric(filename,*linenum,line1[k],true,lmp);
       }
       for (j=1;j<net[l].dimensions[i+1];j++) {
         ptr = fgets(linetemp,longline,fp);
@@ -606,7 +606,7 @@ void PairRANN::read_weight(std::vector<std::string> line,std::vector<std::string
         nwords = line1.size();
         if (nwords != net[l].dimensions[i])error->one(filename,*linenum,"invalid weights per line");
         for (k=0;k<net[l].dimensions[i];k++) {
-          net[l].Weights[i][j*net[l].dimensions[i]+k] = utils::numeric(filename,*linenum,line1[k].c_str(),true,lmp);
+          net[l].Weights[i][j*net[l].dimensions[i]+k] = utils::numeric(filename,*linenum,line1[k],true,lmp);
         }
       }
       return;
@@ -621,19 +621,19 @@ void PairRANN::read_bias(std::vector<std::string> line,std::vector<std::string>
   for (l=0;l<nelements;l++) {
     if (line[1].compare(elements[l])==0) {
       if (net[l].layers==0)error->one(filename,*linenum-1,"networklayers must be defined before biases.");
-      i=utils::inumeric(filename,*linenum,line[2].c_str(),true,lmp);
+      i=utils::inumeric(filename,*linenum,line[2],true,lmp);
       if (i>=net[l].layers || i<0)error->one(filename,*linenum-1,"invalid bias layer");
       if (net[l].dimensions[i]==0) error->one(filename,*linenum-1,"network layer sizes must be defined before corresponding bias");
       biasdefined[l][i] = true;
       net[l].Biases[i] = new double[net[l].dimensions[i+1]];
-      net[l].Biases[i][0] = utils::numeric(filename,*linenum,line1[0].c_str(),true,lmp);
+      net[l].Biases[i][0] = utils::numeric(filename,*linenum,line1[0],true,lmp);
       for (j=1;j<net[l].dimensions[i+1];j++) {
         ptr=fgets(linetemp,MAXLINE,fp);
         if (ptr==nullptr)error->one(filename,*linenum,"unexpected end of potential file!");
         (*linenum)++;
         Tokenizer values1 = Tokenizer(linetemp,": ,\t_\n");
         line1 = values1.as_vector();
-        net[l].Biases[i][j] = utils::numeric(filename,*linenum,line1[0].c_str(),true,lmp);
+        net[l].Biases[i][j] = utils::numeric(filename,*linenum,line1[0],true,lmp);
       }
       return;
     }
@@ -680,10 +680,10 @@ void PairRANN::read_screening(std::vector<std::string> line,std::vector<std::str
   k = atomtypes[2];
   int index = i*nelements*nelements+j*nelements+k;
   if (line[4].compare("Cmin")==0)  {
-    screening_min[index] = utils::numeric(filename,linenum,line1[0].c_str(),true,lmp);
+    screening_min[index] = utils::numeric(filename,linenum,line1[0],true,lmp);
   }
   else if (line[4].compare("Cmax")==0) {
-    screening_max[index] = utils::numeric(filename,linenum,line1[0].c_str(),true,lmp);
+    screening_max[index] = utils::numeric(filename,linenum,line1[0],true,lmp);
   }
   else error->one(filename,linenum-1,"unrecognized screening keyword");
   delete[] atomtypes;

src/ML-SNAP/pair_snap.cpp

@@ -570,8 +570,8 @@ void PairSNAP::read_files(char *coefffilename, char *paramfilename)
     else
       elementflags[jelem] = 1;
 
-    radelem[jelem] = utils::numeric(FLERR,words[1].c_str(),false,lmp);
-    wjelem[jelem] = utils::numeric(FLERR,words[2].c_str(),false,lmp);
+    radelem[jelem] = utils::numeric(FLERR,words[1],false,lmp);
+    wjelem[jelem] = utils::numeric(FLERR,words[2],false,lmp);
 
     if (comm->me == 0)
       utils::logmesg(lmp,"SNAP Element = {}, Radius {}, Weight {}\n",
@@ -672,34 +672,33 @@ void PairSNAP::read_files(char *coefffilename, char *paramfilename)
       utils::logmesg(lmp,"SNAP keyword {} {}\n",keywd,keyval);
 
     if (keywd == "rcutfac") {
-      rcutfac = utils::numeric(FLERR,keyval.c_str(),false,lmp);
+      rcutfac = utils::numeric(FLERR,keyval,false,lmp);
       rcutfacflag = 1;
     } else if (keywd == "twojmax") {
-      twojmax = utils::inumeric(FLERR,keyval.c_str(),false,lmp);
+      twojmax = utils::inumeric(FLERR,keyval,false,lmp);
       twojmaxflag = 1;
     } else if (keywd == "rfac0")
-      rfac0 = utils::numeric(FLERR,keyval.c_str(),false,lmp);
+      rfac0 = utils::numeric(FLERR,keyval,false,lmp);
     else if (keywd == "rmin0")
-      rmin0 = utils::numeric(FLERR,keyval.c_str(),false,lmp);
+      rmin0 = utils::numeric(FLERR,keyval,false,lmp);
     else if (keywd == "switchflag")
-      switchflag = utils::inumeric(FLERR,keyval.c_str(),false,lmp);
+      switchflag = utils::inumeric(FLERR,keyval,false,lmp);
     else if (keywd == "bzeroflag")
-      bzeroflag = utils::inumeric(FLERR,keyval.c_str(),false,lmp);
+      bzeroflag = utils::inumeric(FLERR,keyval,false,lmp);
     else if (keywd == "quadraticflag")
-      quadraticflag = utils::inumeric(FLERR,keyval.c_str(),false,lmp);
+      quadraticflag = utils::inumeric(FLERR,keyval,false,lmp);
     else if (keywd == "chemflag")
-      chemflag = utils::inumeric(FLERR,keyval.c_str(),false,lmp);
+      chemflag = utils::inumeric(FLERR,keyval,false,lmp);
     else if (keywd == "bnormflag")
-      bnormflag = utils::inumeric(FLERR,keyval.c_str(),false,lmp);
+      bnormflag = utils::inumeric(FLERR,keyval,false,lmp);
     else if (keywd == "wselfallflag")
-      wselfallflag = utils::inumeric(FLERR,keyval.c_str(),false,lmp);
+      wselfallflag = utils::inumeric(FLERR,keyval,false,lmp);
     else if (keywd == "chunksize")
-      chunksize = utils::inumeric(FLERR,keyval.c_str(),false,lmp);
+      chunksize = utils::inumeric(FLERR,keyval,false,lmp);
     else if (keywd == "parallelthresh")
-      parallel_thresh = utils::inumeric(FLERR,keyval.c_str(),false,lmp);
+      parallel_thresh = utils::inumeric(FLERR,keyval,false,lmp);
     else
-      error->all(FLERR,"Unknown parameter '{}' in SNAP "
-                                   "parameter file", keywd);
+      error->all(FLERR,"Unknown parameter '{}' in SNAP parameter file", keywd);
   }
 
   if (rcutfacflag == 0 || twojmaxflag == 0)

src/MOFFF/angle_class2_p6.cpp

@@ -174,6 +174,8 @@ void AngleClass2P6::compute(int eflag, int vflag)
 
     // force & energy for bond-angle term
 
+    dr1 = r1 - ba_r1[type];
+    dr2 = r2 - ba_r2[type];
     aa1 = s * dr1 * ba_k1[type];
     aa2 = s * dr2 * ba_k2[type];
 
@@ -479,6 +481,9 @@ double AngleClass2P6::single(int type, int i1, int i2, int i3)
   double dr2 = r2 - bb_r2[type];
   energy += bb_k[type]*dr1*dr2;
 
+  dr1 = r1 - ba_r1[type];
+  dr2 = r2 - ba_r2[type];
   energy += ba_k1[type]*dr1*dtheta + ba_k2[type]*dr2*dtheta;
+
   return energy;
 }

src/MOLECULE/angle_table.cpp

@@ -470,9 +470,9 @@ void AngleTable::compute_table(Table *tb)
 
   memory->create(tb->ang,tablength,"angle:ang");
   memory->create(tb->e,tablength,"angle:e");
-  memory->create(tb->de,tlm1,"angle:de");
+  memory->create(tb->de,tablength,"angle:de");
   memory->create(tb->f,tablength,"angle:f");
-  memory->create(tb->df,tlm1,"angle:df");
+  memory->create(tb->df,tablength,"angle:df");
   memory->create(tb->e2,tablength,"angle:e2");
   memory->create(tb->f2,tablength,"angle:f2");
 
@@ -488,6 +488,9 @@ void AngleTable::compute_table(Table *tb)
     tb->de[i] = tb->e[i+1] - tb->e[i];
     tb->df[i] = tb->f[i+1] - tb->f[i];
   }
+  // get final elements from linear extrapolation
+  tb->de[tlm1] = 2.0*tb->de[tlm1-1] - tb->de[tlm1-2];
+  tb->df[tlm1] = 2.0*tb->df[tlm1-1] - tb->df[tlm1-2];
 
   double ep0 = - tb->f[0];
   double epn = - tb->f[tlm1];
@@ -575,7 +578,7 @@ void AngleTable::spline(double *x, double *y, int n,
   double p,qn,sig,un;
   double *u = new double[n];
 
-  if (yp1 > 0.99e30) y2[0] = u[0] = 0.0;
+  if (yp1 > 0.99e300) y2[0] = u[0] = 0.0;
   else {
     y2[0] = -0.5;
     u[0] = (3.0/(x[1]-x[0])) * ((y[1]-y[0]) / (x[1]-x[0]) - yp1);
@@ -587,7 +590,7 @@ void AngleTable::spline(double *x, double *y, int n,
     u[i] = (y[i+1]-y[i]) / (x[i+1]-x[i]) - (y[i]-y[i-1]) / (x[i]-x[i-1]);
     u[i] = (6.0*u[i] / (x[i+1]-x[i-1]) - sig*u[i-1]) / p;
   }
-  if (ypn > 0.99e30) qn = un = 0.0;
+  if (ypn > 0.99e300) qn = un = 0.0;
   else {
     qn = 0.5;
     un = (3.0/(x[n-1]-x[n-2])) * (ypn - (y[n-1]-y[n-2]) / (x[n-1]-x[n-2]));
@@ -615,8 +618,7 @@ double AngleTable::splint(double *xa, double *ya, double *y2a, int n, double x)
   h = xa[khi]-xa[klo];
   a = (xa[khi]-x) / h;
   b = (x-xa[klo]) / h;
-  y = a*ya[klo] + b*ya[khi] +
-    ((a*a*a-a)*y2a[klo] + (b*b*b-b)*y2a[khi]) * (h*h)/6.0;
+  y = a*ya[klo] + b*ya[khi] + ((a*a*a-a)*y2a[klo] + (b*b*b-b)*y2a[khi]) * (h*h)/6.0;
   return y;
 }
 
@@ -632,8 +634,9 @@ void AngleTable::uf_lookup(int type, double x, double &u, double &f)
 
   double fraction,a,b;
   const Table *tb = &tables[tabindex[type]];
-  int itable = static_cast<int> (x * tb->invdelta);
 
+  // invdelta is based on tablength-1
+  int itable = static_cast<int> (x * tb->invdelta);
   if (itable < 0) itable = 0;
   if (itable >= tablength) itable = tablength-1;
 
@@ -647,11 +650,9 @@ void AngleTable::uf_lookup(int type, double x, double &u, double &f)
     b = (x - tb->ang[itable]) * tb->invdelta;
     a = 1.0 - b;
     u = a * tb->e[itable] + b * tb->e[itable+1] +
-      ((a*a*a-a)*tb->e2[itable] + (b*b*b-b)*tb->e2[itable+1]) *
-      tb->deltasq6;
+      ((a*a*a-a)*tb->e2[itable] + (b*b*b-b)*tb->e2[itable+1]) * tb->deltasq6;
     f = a * tb->f[itable] + b * tb->f[itable+1] +
-      ((a*a*a-a)*tb->f2[itable] + (b*b*b-b)*tb->f2[itable+1]) *
-      tb->deltasq6;
+      ((a*a*a-a)*tb->f2[itable] + (b*b*b-b)*tb->f2[itable+1]) * tb->deltasq6;
   }
 }
 
@@ -681,7 +682,6 @@ void AngleTable::u_lookup(int type, double x, double &u)
     b = (x - tb->ang[itable]) * tb->invdelta;
     a = 1.0 - b;
     u = a * tb->e[itable] + b * tb->e[itable+1] +
-      ((a*a*a-a)*tb->e2[itable] + (b*b*b-b)*tb->e2[itable+1]) *
-      tb->deltasq6;
+      ((a*a*a-a)*tb->e2[itable] + (b*b*b-b)*tb->e2[itable+1]) * tb->deltasq6;
   }
 }

src/MOLECULE/bond_table.cpp

@@ -435,9 +435,9 @@ void BondTable::compute_table(Table *tb)
 
   memory->create(tb->r,tablength,"bond:r");
   memory->create(tb->e,tablength,"bond:e");
-  memory->create(tb->de,tlm1,"bond:de");
+  memory->create(tb->de,tablength,"bond:de");
   memory->create(tb->f,tablength,"bond:f");
-  memory->create(tb->df,tlm1,"bond:df");
+  memory->create(tb->df,tablength,"bond:df");
   memory->create(tb->e2,tablength,"bond:e2");
   memory->create(tb->f2,tablength,"bond:f2");
 
@@ -453,6 +453,9 @@ void BondTable::compute_table(Table *tb)
     tb->de[i] = tb->e[i+1] - tb->e[i];
     tb->df[i] = tb->f[i+1] - tb->f[i];
   }
+  // get final elements from linear extrapolation
+  tb->de[tlm1] = 2.0*tb->de[tlm1-1] - tb->de[tlm1-2];
+  tb->df[tlm1] = 2.0*tb->df[tlm1-1] - tb->df[tlm1-2];
 
   double ep0 = - tb->f[0];
   double epn = - tb->f[tlm1];
@@ -538,7 +541,7 @@ void BondTable::spline(double *x, double *y, int n,
   double p,qn,sig,un;
   double *u = new double[n];
 
-  if (yp1 > 0.99e30) y2[0] = u[0] = 0.0;
+  if (yp1 > 0.99e300) y2[0] = u[0] = 0.0;
   else {
     y2[0] = -0.5;
     u[0] = (3.0/(x[1]-x[0])) * ((y[1]-y[0]) / (x[1]-x[0]) - yp1);
@@ -550,7 +553,7 @@ void BondTable::spline(double *x, double *y, int n,
     u[i] = (y[i+1]-y[i]) / (x[i+1]-x[i]) - (y[i]-y[i-1]) / (x[i]-x[i-1]);
     u[i] = (6.0*u[i] / (x[i+1]-x[i-1]) - sig*u[i-1]) / p;
   }
-  if (ypn > 0.99e30) qn = un = 0.0;
+  if (ypn > 0.99e300) qn = un = 0.0;
   else {
     qn = 0.5;
     un = (3.0/(x[n-1]-x[n-2])) * (ypn - (y[n-1]-y[n-2]) / (x[n-1]-x[n-2]));
@@ -578,8 +581,7 @@ double BondTable::splint(double *xa, double *ya, double *y2a, int n, double x)
   h = xa[khi]-xa[klo];
   a = (xa[khi]-x) / h;
   b = (x-xa[klo]) / h;
-  y = a*ya[klo] + b*ya[khi] +
-    ((a*a*a-a)*y2a[klo] + (b*b*b-b)*y2a[khi]) * (h*h)/6.0;
+  y = a*ya[klo] + b*ya[khi] + ((a*a*a-a)*y2a[klo] + (b*b*b-b)*y2a[khi]) * (h*h)/6.0;
   return y;
 }
 
@@ -598,11 +600,9 @@ void BondTable::uf_lookup(int type, double x, double &u, double &f)
   const Table *tb = &tables[tabindex[type]];
   const int itable = static_cast<int> ((x - tb->lo) * tb->invdelta);
   if (itable < 0)
-    error->one(FLERR,"Bond length < table inner cutoff: "
-               "type {} length {:.8}",type,x);
+    error->one(FLERR,"Bond length < table inner cutoff: type {} length {:.8}",type,x);
   else if (itable >= tablength)
-    error->one(FLERR,"Bond length > table outer cutoff: "
-               "type {} length {:.8}",type,x);
+    error->one(FLERR,"Bond length > table outer cutoff: type {} length {:.8}",type,x);
 
   if (tabstyle == LINEAR) {
     fraction = (x - tb->r[itable]) * tb->invdelta;
@@ -614,10 +614,8 @@ void BondTable::uf_lookup(int type, double x, double &u, double &f)
     b = (x - tb->r[itable]) * tb->invdelta;
     a = 1.0 - b;
     u = a * tb->e[itable] + b * tb->e[itable+1] +
-      ((a*a*a-a)*tb->e2[itable] + (b*b*b-b)*tb->e2[itable+1]) *
-      tb->deltasq6;
+      ((a*a*a-a)*tb->e2[itable] + (b*b*b-b)*tb->e2[itable+1]) * tb->deltasq6;
     f = a * tb->f[itable] + b * tb->f[itable+1] +
-      ((a*a*a-a)*tb->f2[itable] + (b*b*b-b)*tb->f2[itable+1]) *
-      tb->deltasq6;
+      ((a*a*a-a)*tb->f2[itable] + (b*b*b-b)*tb->f2[itable+1]) * tb->deltasq6;
   }
 }

src/MOLECULE/dihedral_table.cpp

@@ -189,11 +189,8 @@ static int solve_cyc_tridiag( const double diag[], size_t d_stride,
    spline and splint routines modified from Numerical Recipes
 ------------------------------------------------------------------------- */
 
-static int cyc_spline(double const *xa,
-                      double const *ya,
-                      int n,
-                      double period,
-                      double *y2a, bool warn)
+static int cyc_spline(double const *xa, double const *ya, int n,
+                      double period, double *y2a, bool warn)
 {
 
   double *diag    = new double[n];
@@ -264,12 +261,8 @@ static int cyc_spline(double const *xa,
 //           range should not exceed period (ie xa[n-1] < xa[0] + period).
 //           x must lie in the range:  [(xa[n-1]-period), (xa[0]+period)]
 //           "period" is typically 2*PI.
-static double cyc_splint(double const *xa,
-                         double const *ya,
-                         double const *y2a,
-                         int n,
-                         double period,
-                         double x)
+static double cyc_splint(double const *xa, double const *ya, double const *y2a,
+                         int n, double period, double x)
 {
   int klo = -1;
   int khi = n;
@@ -302,11 +295,8 @@ static double cyc_splint(double const *xa,
 } // cyc_splint()
 
 
-static double cyc_lin(double const *xa,
-                      double const *ya,
-                      int n,
-                      double period,
-                      double x)
+static double cyc_lin(double const *xa, double const *ya,
+                      int n, double period, double x)
 {
   int klo = -1;
   int khi = n;
@@ -337,21 +327,14 @@ static double cyc_lin(double const *xa,
 
 } // cyc_lin()
 
-
-
-
 // cyc_splintD(): Evaluate the deriviative of a cyclic spline at position x,
 //           with n control points at xa[], ya[], with parameters y2a[].
 //           The xa[] must be monotonically increasing and their
 //           range should not exceed period (ie xa[n-1] < xa[0] + period).
 //           x must lie in the range:  [(xa[n-1]-period), (xa[0]+period)]
 //           "period" is typically 2*PI.
-static double cyc_splintD(double const *xa,
-                          double const *ya,
-                          double const *y2a,
-                          int n,
-                          double period,
-                          double x)
+static double cyc_splintD(double const *xa, double const *ya, double const *y2a,
+                          int n, double period, double x)
 {
   int klo = -1;
   int khi = n; // (not n-1)
@@ -829,9 +812,9 @@ void DihedralTable::coeff(int narg, char **arg)
   // We also want the angles to be sorted in increasing order.
   // This messy code fixes these problems with the user's data:
   {
-    double *phifile_tmp = new double [tb->ninput];  //temporary arrays
-    double *ffile_tmp = new double [tb->ninput];  //used for sorting
-    double *efile_tmp = new double [tb->ninput];
+    double *phifile_tmp = new double[tb->ninput];  //temporary arrays
+    double *ffile_tmp = new double[tb->ninput];  //used for sorting
+    double *efile_tmp = new double[tb->ninput];
 
     // After re-imaging, does the range of angles cross the 0 or 2*PI boundary?
     // If so, find the discontinuity:
@@ -1184,8 +1167,7 @@ void DihedralTable::compute_table(Table *tb)
       if (! tb->f_unspecified)
         tb->f[i] = cyc_splint(tb->phifile,tb->ffile,tb->f2file,tb->ninput,MY_2PI,phi);
     }
-  } // if (tabstyle == SPLINE)
-  else if (tabstyle == LINEAR) {
+  } else if (tabstyle == LINEAR) {
     if (! tb->f_unspecified) {
       for (int i = 0; i < tablength; i++) {
         double phi = i*tb->delta;
@@ -1193,8 +1175,7 @@ void DihedralTable::compute_table(Table *tb)
         tb->e[i]= cyc_lin(tb->phifile,tb->efile,tb->ninput,MY_2PI,phi);
         tb->f[i]= cyc_lin(tb->phifile,tb->ffile,tb->ninput,MY_2PI,phi);
       }
-    }
-    else {
+    } else {
       for (int i = 0; i < tablength; i++) {
         double phi = i*tb->delta;
         tb->phi[i] = phi;
@@ -1269,8 +1250,7 @@ void DihedralTable::param_extract(Table *tb, char *line)
       //else if (word == "EQ") {
       //  tb->theta0 = values.next_double();
       //}
-      else error->one(FLERR,"Invalid keyword in dihedral angle "
-                                        "table parameters ({})", word);
+      else error->one(FLERR,"Invalid keyword in dihedral angle table parameters ({})", word);
     }
   } catch (TokenizerException &e) {
     error->one(FLERR, e.what());

src/MOLECULE/fix_cmap.cpp

@@ -1072,10 +1072,10 @@ void FixCMAP::read_data_header(char *line)
    store CMAP interactions as if newton_bond = OFF, even if actually ON
 ------------------------------------------------------------------------- */
 
-void FixCMAP::read_data_section(char *keyword, int n, char *buf,
+void FixCMAP::read_data_section(char * /*keyword*/, int /*n*/, char *buf,
                                  tagint id_offset)
 {
-  int m,tmp,itype;
+  int m,itype;
   tagint atom1,atom2,atom3,atom4,atom5;
 
   auto lines = utils::split_lines(buf);

src/MOLFILE/molfile_interface.cpp

@@ -75,25 +75,25 @@ extern "C" {
 
   /* corresponding table of masses. */
   static const float pte_mass[] = {
-    /* X  */ 0.00000, 1.00794, 4.00260, 6.941, 9.012182, 10.811,
-    /* C  */ 12.0107, 14.0067, 15.9994, 18.9984032, 20.1797,
-    /* Na */ 22.989770, 24.3050, 26.981538, 28.0855, 30.973761,
-    /* S  */ 32.065, 35.453, 39.948, 39.0983, 40.078, 44.955910,
-    /* Ti */ 47.867, 50.9415, 51.9961, 54.938049, 55.845, 58.9332,
-    /* Ni */ 58.6934, 63.546, 65.409, 69.723, 72.64, 74.92160,
-    /* Se */ 78.96, 79.904, 83.798, 85.4678, 87.62, 88.90585,
-    /* Zr */ 91.224, 92.90638, 95.94, 98.0, 101.07, 102.90550,
-    /* Pd */ 106.42, 107.8682, 112.411, 114.818, 118.710, 121.760,
-    /* Te */ 127.60, 126.90447, 131.293, 132.90545, 137.327,
-    /* La */ 138.9055, 140.116, 140.90765, 144.24, 145.0, 150.36,
-    /* Eu */ 151.964, 157.25, 158.92534, 162.500, 164.93032,
-    /* Er */ 167.259, 168.93421, 173.04, 174.967, 178.49, 180.9479,
-    /* W  */ 183.84, 186.207, 190.23, 192.217, 195.078, 196.96655,
-    /* Hg */ 200.59, 204.3833, 207.2, 208.98038, 209.0, 210.0, 222.0,
-    /* Fr */ 223.0, 226.0, 227.0, 232.0381, 231.03588, 238.02891,
-    /* Np */ 237.0, 244.0, 243.0, 247.0, 247.0, 251.0, 252.0, 257.0,
-    /* Md */ 258.0, 259.0, 262.0, 261.0, 262.0, 266.0, 264.0, 269.0,
-    /* Mt */ 268.0, 271.0, 272.0
+    /* X  */ 0.00000f, 1.00794f, 4.00260f, 6.941f, 9.012182f, 10.811f,
+    /* C  */ 12.0107f, 14.0067f, 15.9994f, 18.9984032f, 20.1797f,
+    /* Na */ 22.989770f, 24.3050f, 26.981538f, 28.0855f, 30.973761f,
+    /* S  */ 32.065f, 35.453f, 39.948f, 39.0983f, 40.078f, 44.955910f,
+    /* Ti */ 47.867f, 50.9415f, 51.9961f, 54.938049f, 55.845f, 58.9332f,
+    /* Ni */ 58.6934f, 63.546f, 65.409f, 69.723f, 72.64f, 74.92160f,
+    /* Se */ 78.96f, 79.904f, 83.798f, 85.4678f, 87.62f, 88.90585f,
+    /* Zr */ 91.224f, 92.90638f, 95.94f, 98.0f, 101.07f, 102.90550f,
+    /* Pd */ 106.42f, 107.8682f, 112.411f, 114.818f, 118.710f, 121.760f,
+    /* Te */ 127.60f, 126.90447f, 131.293f, 132.90545f, 137.327f,
+    /* La */ 138.9055f, 140.116f, 140.90765f, 144.24f, 145.0f, 150.36f,
+    /* Eu */ 151.964f, 157.25f, 158.92534f, 162.500f, 164.93032f,
+    /* Er */ 167.259f, 168.93421f, 173.04f, 174.967f, 178.49f, 180.9479f,
+    /* W  */ 183.84f, 186.207f, 190.23f, 192.217f, 195.078f, 196.96655f,
+    /* Hg */ 200.59f, 204.3833f, 207.2f, 208.98038f, 209.0f, 210.0f, 222.0f,
+    /* Fr */ 223.0f, 226.0f, 227.0f, 232.0381f, 231.03588f, 238.02891f,
+    /* Np */ 237.0f, 244.0f, 243.0f, 247.0f, 247.0f, 251.0f, 252.0f, 257.0f,
+    /* Md */ 258.0f, 259.0f, 262.0f, 261.0f, 262.0f, 266.0f, 264.0f, 269.0f,
+    /* Mt */ 268.0f, 271.0f, 272.0f
   };
 
   /*
@@ -107,25 +107,25 @@ extern "C" {
    * Rmin/2 parameters for (SOD, POT, CLA, CAL, MG, CES) by default.
    */
   static const float pte_vdw_radius[] = {
-    /* X  */ 1.5, 1.2, 1.4, 1.82, 2.0, 2.0,
-    /* C  */ 1.7, 1.55, 1.52, 1.47, 1.54,
-    /* Na */ 1.36, 1.18, 2.0, 2.1, 1.8,
-    /* S  */ 1.8, 2.27, 1.88, 1.76, 1.37, 2.0,
-    /* Ti */ 2.0, 2.0, 2.0, 2.0, 2.0, 2.0,
-    /* Ni */ 1.63, 1.4, 1.39, 1.07, 2.0, 1.85,
-    /* Se */ 1.9, 1.85, 2.02, 2.0, 2.0, 2.0,
-    /* Zr */ 2.0, 2.0, 2.0, 2.0, 2.0, 2.0,
-    /* Pd */ 1.63, 1.72, 1.58, 1.93, 2.17, 2.0,
-    /* Te */ 2.06, 1.98, 2.16, 2.1, 2.0,
-    /* La */ 2.0, 2.0, 2.0, 2.0, 2.0, 2.0,
-    /* Eu */ 2.0, 2.0, 2.0, 2.0, 2.0,
-    /* Er */ 2.0, 2.0, 2.0, 2.0, 2.0, 2.0,
-    /* W  */ 2.0, 2.0, 2.0, 2.0, 1.72, 1.66,
-    /* Hg */ 1.55, 1.96, 2.02, 2.0, 2.0, 2.0, 2.0,
-    /* Fr */ 2.0, 2.0, 2.0, 2.0, 2.0, 1.86,
-    /* Np */ 2.0, 2.0, 2.0, 2.0, 2.0, 2.0, 2.0, 2.0,
-    /* Md */ 2.0, 2.0, 2.0, 2.0, 2.0, 2.0, 2.0, 2.0,
-    /* Mt */ 2.0, 2.0, 2.0
+    /* X  */ 1.5f, 1.2f, 1.4f, 1.82f, 2.0f, 2.0f,
+    /* C  */ 1.7f, 1.55f, 1.52f, 1.47f, 1.54f,
+    /* Na */ 1.36f, 1.18f, 2.0f, 2.1f, 1.8f,
+    /* S  */ 1.8f, 2.27f, 1.88f, 1.76f, 1.37f, 2.0f,
+    /* Ti */ 2.0f, 2.0f, 2.0f, 2.0f, 2.0f, 2.0f,
+    /* Ni */ 1.63f, 1.4f, 1.39f, 1.07f, 2.0f, 1.85f,
+    /* Se */ 1.9f, 1.85f, 2.02f, 2.0f, 2.0f, 2.0f,
+    /* Zr */ 2.0f, 2.0f, 2.0f, 2.0f, 2.0f, 2.0f,
+    /* Pd */ 1.63f, 1.72f, 1.58f, 1.93f, 2.17f, 2.0f,
+    /* Te */ 2.06f, 1.98f, 2.16f, 2.1f, 2.0f,
+    /* La */ 2.0f, 2.0f, 2.0f, 2.0f, 2.0f, 2.0f,
+    /* Eu */ 2.0f, 2.0f, 2.0f, 2.0f, 2.0f,
+    /* Er */ 2.0f, 2.0f, 2.0f, 2.0f, 2.0f, 2.0f,
+    /* W  */ 2.0f, 2.0f, 2.0f, 2.0f, 1.72f, 1.66f,
+    /* Hg */ 1.55f, 1.96f, 2.02f, 2.0f, 2.0f, 2.0f, 2.0f,
+    /* Fr */ 2.0f, 2.0f, 2.0f, 2.0f, 2.0f, 1.86f,
+    /* Np */ 2.0f, 2.0f, 2.0f, 2.0f, 2.0f, 2.0f, 2.0f, 2.0f,
+    /* Md */ 2.0f, 2.0f, 2.0f, 2.0f, 2.0f, 2.0f, 2.0f, 2.0f,
+    /* Mt */ 2.0f, 2.0f, 2.0f
   };
 
   /* lookup functions */

src/MPIIO/dump_atom_mpiio.cpp

@@ -38,7 +38,12 @@ using namespace LAMMPS_NS;
 
 /* ---------------------------------------------------------------------- */
 
-DumpAtomMPIIO::DumpAtomMPIIO(LAMMPS *lmp, int narg, char **arg) : DumpAtom(lmp, narg, arg) {}
+DumpAtomMPIIO::DumpAtomMPIIO(LAMMPS *lmp, int narg, char **arg)
+  : DumpAtom(lmp, narg, arg)
+{
+  if (me == 0)
+    error->warning(FLERR,"MPI-IO output is unmaintained and unreliable. Use with caution.");
+}
 
 /* ---------------------------------------------------------------------- */
 

src/MPIIO/dump_cfg_mpiio.cpp

@@ -51,7 +51,11 @@ using namespace LAMMPS_NS;
 /* ---------------------------------------------------------------------- */
 
 DumpCFGMPIIO::DumpCFGMPIIO(LAMMPS *lmp, int narg, char **arg) :
-  DumpCFG(lmp, narg, arg) {}
+  DumpCFG(lmp, narg, arg)
+{
+  if (me == 0)
+    error->warning(FLERR,"MPI-IO output is unmaintained and unreliable. Use with caution.");
+}
 
 /* ---------------------------------------------------------------------- */
 

src/MPIIO/dump_custom_mpiio.cpp

@@ -53,7 +53,12 @@ enum{ LT, LE, GT, GE, EQ, NEQ };
 // clang-format on
 /* ---------------------------------------------------------------------- */
 
-DumpCustomMPIIO::DumpCustomMPIIO(LAMMPS *lmp, int narg, char **arg) : DumpCustom(lmp, narg, arg) {}
+DumpCustomMPIIO::DumpCustomMPIIO(LAMMPS *lmp, int narg, char **arg)
+  : DumpCustom(lmp, narg, arg)
+{
+  if (me == 0)
+    error->warning(FLERR,"MPI-IO output is unmaintained and unreliable. Use with caution.");
+}
 
 /* ---------------------------------------------------------------------- */
 

src/MPIIO/dump_xyz_mpiio.cpp

@@ -52,7 +52,10 @@ enum{LT,LE,GT,GE,EQ,NEQ};
 /* ---------------------------------------------------------------------- */
 
 DumpXYZMPIIO::DumpXYZMPIIO(LAMMPS *lmp, int narg, char **arg) :
-  DumpXYZ(lmp, narg, arg) {}
+  DumpXYZ(lmp, narg, arg) {
+  if (me == 0)
+    error->warning(FLERR,"MPI-IO output is unmaintained and unreliable. Use with caution.");
+}
 
 /* ---------------------------------------------------------------------- */
 

src/NETCDF/dump_netcdf.cpp

@@ -19,6 +19,7 @@
 #if defined(LMP_HAS_NETCDF)
 
 #include "dump_netcdf.h"
+#include "netcdf_units.h"
 
 #include "atom.h"
 #include "comm.h"
@@ -43,6 +44,9 @@
 
 using namespace LAMMPS_NS;
 using namespace MathConst;
+using NetCDFUnits::Quantity;
+using NetCDFUnits::get_unit_for;
+using NetCDFUnits::LMP_MAX_VAR_DIMS;
 
 static const char NC_FRAME_STR[]         = "frame";
 static const char NC_SPATIAL_STR[]       = "spatial";
@@ -63,7 +67,6 @@ static const char NC_SCALE_FACTOR_STR[]  = "scale_factor";
 static constexpr int THIS_IS_A_FIX      = -1;
 static constexpr int THIS_IS_A_COMPUTE  = -2;
 static constexpr int THIS_IS_A_VARIABLE = -3;
-static constexpr int THIS_IS_A_BIGINT   = -4;
 
 /* ---------------------------------------------------------------------- */
 
@@ -102,6 +105,7 @@ DumpNetCDF::DumpNetCDF(LAMMPS *lmp, int narg, char **arg) :
     int ndims = 1;
     std::string mangled = earg[i];
     bool constant = false;
+    int quantity = Quantity::UNKNOWN;
 
     // name mangling
     // in the AMBER specification
@@ -109,26 +113,32 @@ DumpNetCDF::DumpNetCDF(LAMMPS *lmp, int narg, char **arg) :
       idim = mangled[0] - 'x';
       ndims = 3;
       mangled = "coordinates";
+      quantity = Quantity::DISTANCE;
     } else if ((mangled == "vx") || (mangled == "vy") || (mangled == "vz")) {
       idim = mangled[1] - 'x';
       ndims = 3;
       mangled = "velocities";
+      quantity = Quantity::VELOCITY;
     } else if ((mangled == "xs") || (mangled == "ys") || (mangled == "zs")) {
       idim = mangled[0] - 'x';
       ndims = 3;
       mangled = "scaled_coordinates";
+      // no unit for scaled coordinates
     } else if ((mangled == "xu") || (mangled == "yu") || (mangled == "zu")) {
       idim = mangled[0] - 'x';
       ndims = 3;
       mangled = "unwrapped_coordinates";
+      quantity = Quantity::DISTANCE;
     } else if ((mangled == "fx") || (mangled == "fy") || (mangled == "fz")) {
       idim = mangled[1] - 'x';
       ndims = 3;
       mangled = "forces";
+      quantity = Quantity::FORCE;
     } else if ((mangled == "mux") || (mangled == "muy") || (mangled == "muz")) {
       idim = mangled[2] - 'x';
       ndims = 3;
       mangled = "mu";
+      quantity = Quantity::DIPOLE_MOMENT;
     } else if (utils::strmatch(mangled, "^c_")) {
       std::size_t found = mangled.find('[');
       if (found != std::string::npos) {
@@ -175,13 +185,14 @@ DumpNetCDF::DumpNetCDF(LAMMPS *lmp, int narg, char **arg) :
     perat[inc].constant = constant;
     perat[inc].ndumped = 0;
     perat[inc].field[idim] = i;
+    perat[inc].quantity = quantity;
   }
 
   n_buffer = 0;
   int_buffer = nullptr;
   double_buffer = nullptr;
 
-  double_precision = false;
+  type_nc_real = NC_FLOAT;
 
   thermo = false;
   thermovar = nullptr;
@@ -196,7 +207,7 @@ DumpNetCDF::~DumpNetCDF()
   closefile();
 
   delete[] perat;
-  if (thermovar) delete[] thermovar;
+  delete[] thermovar;
 
   if (int_buffer) memory->sfree(int_buffer);
   if (double_buffer) memory->sfree(double_buffer);
@@ -224,7 +235,7 @@ void DumpNetCDF::openfile()
   }
 
   if (thermo && !singlefile_opened) {
-    if (thermovar)  delete[] thermovar;
+    delete[] thermovar;
     thermovar = new int[output->thermo->nfield];
   }
 
@@ -290,18 +301,18 @@ void DumpNetCDF::openfile()
       NCERRX( nc_inq_dimid(ncid, NC_LABEL_STR, &label_dim), NC_LABEL_STR );
 
       for (int i = 0; i < n_perat; i++) {
-        int dims = perat[i].dims;
-        if (vector_dim[dims] < 0) {
+        int dim = perat[i].dims;
+        if (vector_dim[dim] < 0) {
           char dimstr[1024];
-          if (dims == 3) {
+          if (dim == 3) {
             strcpy(dimstr, NC_SPATIAL_STR);
-          } else if (dims == 6) {
+          } else if (dim == 6) {
             strcpy(dimstr, NC_VOIGT_STR);
           } else {
-            sprintf(dimstr, "vec%i", dims);
+            sprintf(dimstr, "vec%i", dim);
           }
-          if (dims != 1) {
-            NCERRX( nc_inq_dimid(ncid, dimstr, &vector_dim[dims]), dimstr );
+          if (dim != 1) {
+            NCERRX( nc_inq_dimid(ncid, dimstr, &vector_dim[dim]), dimstr );
           }
         }
       }
@@ -339,9 +350,8 @@ void DumpNetCDF::openfile()
       if (framei != 0 && !multifile)
         error->all(FLERR,"at keyword requires use of 'append yes'");
 
-      int dims[NC_MAX_VAR_DIMS];
-      size_t index[NC_MAX_VAR_DIMS], count[NC_MAX_VAR_DIMS];
-      double d[1];
+      int dims[LMP_MAX_VAR_DIMS];
+      size_t index[LMP_MAX_VAR_DIMS], count[LMP_MAX_VAR_DIMS];
 
       if (singlefile_opened) return;
       singlefile_opened = 1;
@@ -373,22 +383,22 @@ void DumpNetCDF::openfile()
       }
 
       // default variables
-      dims[0] = 0;
+      dims[0] = vector_dim[3];
       NCERRX( nc_def_var(ncid, NC_SPATIAL_STR, NC_CHAR, 1, dims, &spatial_var), NC_SPATIAL_STR );
       NCERRX( nc_def_var(ncid, NC_CELL_SPATIAL_STR, NC_CHAR, 1, dims, &cell_spatial_var), NC_CELL_SPATIAL_STR );
-      dims[0] = 0;
+      dims[0] = vector_dim[3];
       dims[1] = label_dim;
       NCERRX( nc_def_var(ncid, NC_CELL_ANGULAR_STR, NC_CHAR, 2, dims, &cell_angular_var), NC_CELL_ANGULAR_STR );
 
       dims[0] = frame_dim;
-      NCERRX( nc_def_var(ncid, NC_TIME_STR, NC_DOUBLE, 1, dims, &time_var), NC_TIME_STR);
+      NCERRX( nc_def_var(ncid, NC_TIME_STR, type_nc_real, 1, dims, &time_var), NC_TIME_STR);
       dims[0] = frame_dim;
       dims[1] = cell_spatial_dim;
-      NCERRX( nc_def_var(ncid, NC_CELL_ORIGIN_STR, NC_DOUBLE, 2, dims, &cell_origin_var), NC_CELL_ORIGIN_STR );
-      NCERRX( nc_def_var(ncid, NC_CELL_LENGTHS_STR, NC_DOUBLE, 2, dims, &cell_lengths_var), NC_CELL_LENGTHS_STR );
+      NCERRX( nc_def_var(ncid, NC_CELL_ORIGIN_STR, type_nc_real, 2, dims, &cell_origin_var), NC_CELL_ORIGIN_STR );
+      NCERRX( nc_def_var(ncid, NC_CELL_LENGTHS_STR, type_nc_real, 2, dims, &cell_lengths_var), NC_CELL_LENGTHS_STR );
       dims[0] = frame_dim;
       dims[1] = cell_angular_dim;
-      NCERRX( nc_def_var(ncid, NC_CELL_ANGLES_STR, NC_DOUBLE, 2, dims, &cell_angles_var), NC_CELL_ANGLES_STR );
+      NCERRX( nc_def_var(ncid, NC_CELL_ANGLES_STR, type_nc_real, 2, dims, &cell_angles_var), NC_CELL_ANGLES_STR );
 
       // variables specified in the input file
       dims[0] = frame_dim;
@@ -397,7 +407,6 @@ void DumpNetCDF::openfile()
 
       for (int i = 0; i < n_perat; i++) {
         nc_type xtype;
-
         // Type mangling
         if (vtype[perat[i].field[0]] == Dump::INT) {
           xtype = NC_INT;
@@ -406,10 +415,7 @@ void DumpNetCDF::openfile()
         } else if (vtype[perat[i].field[0]] == Dump::STRING) {
           error->all(FLERR,"Dump netcdf currently does not support dumping string properties");
         } else {
-          if (double_precision)
-            xtype = NC_DOUBLE;
-          else
-            xtype = NC_FLOAT;
+          xtype = type_nc_real;
         }
 
         if (perat[i].constant) {
@@ -430,6 +436,11 @@ void DumpNetCDF::openfile()
             NCERRX( nc_def_var(ncid, perat[i].name, xtype, 3, dims, &perat[i].var), perat[i].name );
           }
         }
+
+        std::string unit = get_unit_for(update->unit_style, perat[i].quantity, error);
+        if (!unit.empty()) {
+          NCERR( nc_put_att_text(ncid, perat[i].var, NC_UNITS_STR, unit.size(), unit.c_str()) );
+        }
       }
 
       // perframe variables
@@ -437,7 +448,7 @@ void DumpNetCDF::openfile()
         Thermo *th = output->thermo;
         for (int i = 0; i < th->nfield; i++) {
           if (th->vtype[i] == Thermo::FLOAT) {
-            NCERRX( nc_def_var(ncid, th->keyword[i], NC_DOUBLE, 1, dims,
+            NCERRX( nc_def_var(ncid, th->keyword[i], type_nc_real, 1, dims,
                                &thermovar[i]), th->keyword[i] );
           } else if (th->vtype[i] == Thermo::INT) {
             NCERRX( nc_def_var(ncid, th->keyword[i], NC_INT, 1, dims,
@@ -461,43 +472,18 @@ void DumpNetCDF::openfile()
       NCERR( nc_put_att_text(ncid, NC_GLOBAL, "program", 6, "LAMMPS") );
       NCERR( nc_put_att_text(ncid, NC_GLOBAL, "programVersion",strlen(lmp->version), lmp->version) );
 
-      // units
-      if (!strcmp(update->unit_style, "lj")) {
-        NCERR( nc_put_att_text(ncid, time_var, NC_UNITS_STR, 2, "lj") );
-        NCERR( nc_put_att_text(ncid, cell_origin_var, NC_UNITS_STR, 2, "lj") );
-        NCERR( nc_put_att_text(ncid, cell_lengths_var, NC_UNITS_STR, 2, "lj") );
-      } else if (!strcmp(update->unit_style, "real")) {
-        NCERR( nc_put_att_text(ncid, time_var, NC_UNITS_STR, 11, "femtosecond") );
-        NCERR( nc_put_att_text(ncid, cell_origin_var, NC_UNITS_STR, 8, "Angstrom") );
-        NCERR( nc_put_att_text(ncid, cell_lengths_var, NC_UNITS_STR, 8, "Angstrom") );
-      } else if (!strcmp(update->unit_style, "metal")) {
-        NCERR( nc_put_att_text(ncid, time_var, NC_UNITS_STR, 10, "picosecond") );
-        NCERR( nc_put_att_text(ncid, cell_origin_var, NC_UNITS_STR, 8, "Angstrom") );
-        NCERR( nc_put_att_text(ncid, cell_lengths_var, NC_UNITS_STR, 8, "Angstrom") );
-      } else if (!strcmp(update->unit_style, "si")) {
-        NCERR( nc_put_att_text(ncid, time_var, NC_UNITS_STR, 6, "second") );
-        NCERR( nc_put_att_text(ncid, cell_origin_var, NC_UNITS_STR, 5, "meter") );
-        NCERR( nc_put_att_text(ncid, cell_lengths_var, NC_UNITS_STR, 5, "meter") );
-      } else if (!strcmp(update->unit_style, "cgs")) {
-        NCERR( nc_put_att_text(ncid, time_var, NC_UNITS_STR, 6, "second") );
-        NCERR( nc_put_att_text(ncid, cell_origin_var, NC_UNITS_STR, 10, "centimeter") );
-        NCERR( nc_put_att_text(ncid, cell_lengths_var, NC_UNITS_STR, 10, "centimeter") );
-      } else if (!strcmp(update->unit_style, "electron")) {
-        NCERR( nc_put_att_text(ncid, time_var, NC_UNITS_STR, 11, "femtosecond") );
-        NCERR( nc_put_att_text(ncid, cell_origin_var, NC_UNITS_STR, 4, "Bohr") );
-        NCERR( nc_put_att_text(ncid, cell_lengths_var, NC_UNITS_STR, 4, "Bohr") );
-      } else {
-        error->all(FLERR,"Unsupported unit style: {}", update->unit_style);
-      }
+      // units & scale
+      std::string unit = get_unit_for(update->unit_style, Quantity::TIME, error);
+      NCERR( nc_put_att_text(ncid, time_var, NC_UNITS_STR, unit.size(), unit.c_str()) );
 
-      NCERR( nc_put_att_text(ncid, cell_angles_var, NC_UNITS_STR,6, "degree") );
+      unit = get_unit_for(update->unit_style, Quantity::DISTANCE, error);
+      NCERR( nc_put_att_text(ncid, cell_origin_var, NC_UNITS_STR, unit.size(), unit.c_str()) );
+      NCERR( nc_put_att_text(ncid, cell_lengths_var, NC_UNITS_STR, unit.size(), unit.c_str()) );
 
-      d[0] = update->dt;
-      NCERR( nc_put_att_double(ncid, time_var, NC_SCALE_FACTOR_STR,NC_DOUBLE, 1, d) );
-      d[0] = 1.0;
-      NCERR( nc_put_att_double(ncid, cell_origin_var, NC_SCALE_FACTOR_STR,NC_DOUBLE, 1, d) );
-      d[0] = 1.0;
-      NCERR( nc_put_att_double(ncid, cell_lengths_var, NC_SCALE_FACTOR_STR,NC_DOUBLE, 1, d) );
+      NCERR( nc_put_att_text(ncid, cell_angles_var, NC_UNITS_STR, 6, "degree") );
+
+      float scale[1] = {static_cast<float>(update->dt)};
+      NCERR( nc_put_att_float(ncid, time_var, NC_SCALE_FACTOR_STR, NC_FLOAT, 1, scale) );
 
       /*
        * Finished with definition
@@ -735,8 +721,8 @@ void DumpNetCDF::write_header(bigint n)
 
 void DumpNetCDF::write_data(int n, double *mybuf)
 {
-  size_t start[NC_MAX_VAR_DIMS], count[NC_MAX_VAR_DIMS];
-  ptrdiff_t stride[NC_MAX_VAR_DIMS];
+  size_t start[LMP_MAX_VAR_DIMS], count[LMP_MAX_VAR_DIMS];
+  ptrdiff_t stride[LMP_MAX_VAR_DIMS];
 
   if (!int_buffer) {
     n_buffer = n;
@@ -872,7 +858,12 @@ int DumpNetCDF::modify_param(int narg, char **arg)
   if (strcmp(arg[iarg],"double") == 0) {
     iarg++;
     if (iarg >= narg) error->all(FLERR,"expected 'yes' or 'no' after 'double' keyword.");
-    double_precision = utils::logical(FLERR,arg[iarg],false,lmp) == 1;
+
+    if (utils::logical(FLERR,arg[iarg],false,lmp) == 1)
+      type_nc_real = NC_DOUBLE;
+    else
+      type_nc_real = NC_FLOAT;
+
     iarg++;
     return 2;
   } else if (strcmp(arg[iarg],"at") == 0) {
@@ -897,10 +888,10 @@ int DumpNetCDF::modify_param(int narg, char **arg)
 void DumpNetCDF::ncerr(int err, const char *descr, int line)
 {
   if (err != NC_NOERR) {
-    if (descr) error->one(FLERR,"NetCDF failed with error '{}' (while accessing '{}') "
-                          " in line {} of {}.", nc_strerror(err), descr, line, __FILE__);
-    else error->one(FLERR,"NetCDF failed with error '{}' in line {} of {}.",
-                    nc_strerror(err), line, __FILE__);
+    if (descr) error->one(__FILE__, line, "NetCDF failed with error '{}' (while accessing '{}') ",
+                          nc_strerror(err), descr);
+    else error->one(__FILE__, line,"NetCDF failed with error '{}' in line {} of {}.",
+                    nc_strerror(err));
   }
 }
 

src/NETCDF/dump_netcdf.h

@@ -24,15 +24,12 @@ DumpStyle(netcdf,DumpNetCDF);
 #else
 
 #ifndef LMP_DUMP_NETCDF_H
-#define LMP_DUMP_NETCDFC_H
+#define LMP_DUMP_NETCDF_H
 
 #include "dump_custom.h"
 
 namespace LAMMPS_NS {
 
-const int NC_FIELD_NAME_MAX = 100;
-const int DUMP_NC_MAX_DIMS = 100;
-
 class DumpNetCDF : public DumpCustom {
  public:
   DumpNetCDF(class LAMMPS *, int, char **);
@@ -40,12 +37,16 @@ class DumpNetCDF : public DumpCustom {
   virtual void write();
 
  private:
+  static constexpr int NC_FIELD_NAME_MAX = 100;
+  static constexpr int DUMP_NC_MAX_DIMS = 100;
+
   // per-atoms quantities (positions, velocities, etc.)
   struct nc_perat_t {
     int dims;                        // number of dimensions
     int field[DUMP_NC_MAX_DIMS];     // field indices corresponding to the dim.
     char name[NC_FIELD_NAME_MAX];    // field name
     int var;                         // NetCDF variable
+    int quantity;                    // type of the quantity
 
     bool constant;    // is this property per file (not per frame)
     int ndumped;      // number of enties written for this prop.
@@ -62,8 +63,8 @@ class DumpNetCDF : public DumpCustom {
 
   int *thermovar;    // NetCDF variables for thermo output
 
-  bool double_precision;    // write everything as double precision
-  bool thermo;              // write thermo output to netcdf file
+  int type_nc_real;    // netcdf type to use for real variables: float or double
+  bool thermo;         // write thermo output to netcdf file
 
   bigint n_buffer;          // size of buffer
   bigint *int_buffer;       // buffer for passing data to netcdf

src/NETCDF/dump_netcdf_mpiio.cpp

@@ -19,6 +19,7 @@
 #if defined(LMP_HAS_PNETCDF)
 
 #include "dump_netcdf_mpiio.h"
+#include "netcdf_units.h"
 
 #include "atom.h"
 #include "comm.h"
@@ -43,6 +44,9 @@
 
 using namespace LAMMPS_NS;
 using namespace MathConst;
+using NetCDFUnits::Quantity;
+using NetCDFUnits::get_unit_for;
+using NetCDFUnits::LMP_MAX_VAR_DIMS;
 
 static const char NC_FRAME_STR[]         = "frame";
 static const char NC_SPATIAL_STR[]       = "spatial";
@@ -63,7 +67,6 @@ static const char NC_SCALE_FACTOR_STR[]  = "scale_factor";
 static constexpr int THIS_IS_A_FIX      = -1;
 static constexpr int THIS_IS_A_COMPUTE  = -2;
 static constexpr int THIS_IS_A_VARIABLE = -3;
-static constexpr int THIS_IS_A_BIGINT   = -4;
 
 /* ---------------------------------------------------------------------- */
 
@@ -101,7 +104,7 @@ DumpNetCDFMPIIO::DumpNetCDFMPIIO(LAMMPS *lmp, int narg, char **arg) :
     int idim = 0;
     int ndims = 1;
     std::string mangled = earg[i];
-    bool constant = false;
+    int quantity = Quantity::UNKNOWN;
 
     // name mangling
     // in the AMBER specification
@@ -109,26 +112,32 @@ DumpNetCDFMPIIO::DumpNetCDFMPIIO(LAMMPS *lmp, int narg, char **arg) :
       idim = mangled[0] - 'x';
       ndims = 3;
       mangled = "coordinates";
+      quantity = Quantity::DISTANCE;
     } else if ((mangled == "vx") || (mangled == "vy") || (mangled == "vz")) {
       idim = mangled[1] - 'x';
       ndims = 3;
       mangled = "velocities";
+      quantity = Quantity::VELOCITY;
     } else if ((mangled == "xs") || (mangled == "ys") || (mangled == "zs")) {
       idim = mangled[0] - 'x';
       ndims = 3;
       mangled = "scaled_coordinates";
+      // no unit for scaled coordinates
     } else if ((mangled == "xu") || (mangled == "yu") || (mangled == "zu")) {
       idim = mangled[0] - 'x';
       ndims = 3;
       mangled = "unwrapped_coordinates";
+      quantity = Quantity::DISTANCE;
     } else if ((mangled == "fx") || (mangled == "fy") || (mangled == "fz")) {
       idim = mangled[1] - 'x';
       ndims = 3;
       mangled = "forces";
+      quantity = Quantity::FORCE;
     } else if ((mangled == "mux") || (mangled == "muy") || (mangled == "muz")) {
       idim = mangled[2] - 'x';
       ndims = 3;
       mangled = "mu";
+      quantity = Quantity::DIPOLE_MOMENT;
     } else if (utils::strmatch(mangled, "^c_")) {
       std::size_t found = mangled.find('[');
       if (found != std::string::npos) {
@@ -173,13 +182,14 @@ DumpNetCDFMPIIO::DumpNetCDFMPIIO(LAMMPS *lmp, int narg, char **arg) :
     }
 
     perat[inc].field[idim] = i;
+    perat[inc].quantity = quantity;
   }
 
   n_buffer = 0;
   int_buffer = nullptr;
   double_buffer = nullptr;
 
-  double_precision = false;
+  type_nc_real = NC_FLOAT;
 
   thermo = false;
   thermovar = nullptr;
@@ -194,7 +204,7 @@ DumpNetCDFMPIIO::~DumpNetCDFMPIIO()
   closefile();
 
   delete[] perat;
-  if (thermovar) delete[] thermovar;
+  delete[] thermovar;
 
   if (int_buffer) memory->sfree(int_buffer);
   if (double_buffer) memory->sfree(double_buffer);
@@ -211,8 +221,7 @@ void DumpNetCDFMPIIO::openfile()
     char *ptr = strchr(filestar,'*');
     *ptr = '\0';
     if (padflag == 0)
-      sprintf(filecurrent,"%s" BIGINT_FORMAT "%s",
-              filestar,update->ntimestep,ptr+1);
+      sprintf(filecurrent,"%s" BIGINT_FORMAT "%s", filestar,update->ntimestep,ptr+1);
     else {
       char bif[8],pad[16];
       strcpy(bif,BIGINT_FORMAT);
@@ -223,7 +232,7 @@ void DumpNetCDFMPIIO::openfile()
   }
 
   if (thermo && !singlefile_opened) {
-    if (thermovar)  delete[] thermovar;
+    delete[] thermovar;
     thermovar = new int[output->thermo->nfield];
   }
 
@@ -275,9 +284,6 @@ void DumpNetCDFMPIIO::openfile()
     if (!platform::file_is_readable(filecurrent))
       error->all(FLERR, "cannot append to non-existent file {}", filecurrent);
 
-    MPI_Offset index[NC_MAX_VAR_DIMS], count[NC_MAX_VAR_DIMS];
-    double d[1];
-
     if (singlefile_opened) return;
     singlefile_opened = 1;
 
@@ -291,18 +297,18 @@ void DumpNetCDFMPIIO::openfile()
     NCERRX( ncmpi_inq_dimid(ncid, NC_LABEL_STR, &label_dim), NC_LABEL_STR );
 
     for (int i = 0; i < n_perat; i++) {
-      int dims = perat[i].dims;
-      if (vector_dim[dims] < 0) {
+      int dim = perat[i].dims;
+      if (vector_dim[dim] < 0) {
         char dimstr[1024];
-        if (dims == 3) {
+        if (dim == 3) {
           strcpy(dimstr, NC_SPATIAL_STR);
-        } else if (dims == 6) {
+        } else if (dim == 6) {
           strcpy(dimstr, NC_VOIGT_STR);
         } else {
-          sprintf(dimstr, "vec%i", dims);
+          sprintf(dimstr, "vec%i", dim);
         }
-        if (dims != 1) {
-          NCERRX( ncmpi_inq_dimid(ncid, dimstr, &vector_dim[dims]), dimstr );
+        if (dim != 1) {
+          NCERRX( ncmpi_inq_dimid(ncid, dimstr, &vector_dim[dim]), dimstr );
         }
       }
     }
@@ -340,9 +346,8 @@ void DumpNetCDFMPIIO::openfile()
     if (framei != 0 && !multifile)
       error->all(FLERR,"at keyword requires use of 'append yes'");
 
-    int dims[NC_MAX_VAR_DIMS];
-    MPI_Offset index[NC_MAX_VAR_DIMS], count[NC_MAX_VAR_DIMS];
-    double d[1];
+    int dims[LMP_MAX_VAR_DIMS];
+    MPI_Offset index[LMP_MAX_VAR_DIMS], count[LMP_MAX_VAR_DIMS];
 
     if (singlefile_opened) return;
     singlefile_opened = 1;
@@ -356,19 +361,24 @@ void DumpNetCDFMPIIO::openfile()
     NCERRX( ncmpi_def_dim(ncid, NC_CELL_ANGULAR_STR, 3, &cell_angular_dim), NC_CELL_ANGULAR_STR );
     NCERRX( ncmpi_def_dim(ncid, NC_LABEL_STR, 10, &label_dim), NC_LABEL_STR );
 
+    if (vector_dim[3] < 0)
+      NCERRX( ncmpi_def_dim(ncid, NC_SPATIAL_STR, 3, &vector_dim[3]), NC_SPATIAL_STR );
+    if (vector_dim[6] < 0)
+      NCERRX( ncmpi_def_dim(ncid, NC_VOIGT_STR, 6, &vector_dim[6]), NC_VOIGT_STR );
+
     for (int i = 0; i < n_perat; i++) {
-      int dims = perat[i].dims;
-      if (vector_dim[dims] < 0) {
+      int dim = perat[i].dims;
+      if (vector_dim[dim] < 0) {
         char dimstr[1024];
-        if (dims == 3) {
+        if (dim == 3) {
           strcpy(dimstr, NC_SPATIAL_STR);
-        } else if (dims == 6) {
+        } else if (dim == 6) {
           strcpy(dimstr, NC_VOIGT_STR);
         } else {
-          sprintf(dimstr, "vec%i", dims);
+          sprintf(dimstr, "vec%i", dim);
         }
-        if (dims != 1) {
-          NCERRX( ncmpi_def_dim(ncid, dimstr, dims, &vector_dim[dims]), dimstr );
+        if (dim != 1) {
+          NCERRX( ncmpi_def_dim(ncid, dimstr, dim, &vector_dim[dim]), dimstr );
         }
       }
     }
@@ -380,16 +390,15 @@ void DumpNetCDFMPIIO::openfile()
     dims[0] = vector_dim[3];
     dims[1] = label_dim;
     NCERRX( ncmpi_def_var(ncid, NC_CELL_ANGULAR_STR, NC_CHAR, 2, dims, &cell_angular_var), NC_CELL_ANGULAR_STR );
-
     dims[0] = frame_dim;
-    NCERRX( ncmpi_def_var(ncid, NC_TIME_STR, NC_DOUBLE, 1, dims, &time_var), NC_TIME_STR);
+    NCERRX( ncmpi_def_var(ncid, NC_TIME_STR, type_nc_real, 1, dims, &time_var), NC_TIME_STR);
     dims[0] = frame_dim;
     dims[1] = cell_spatial_dim;
-    NCERRX( ncmpi_def_var(ncid, NC_CELL_ORIGIN_STR, NC_DOUBLE, 2, dims, &cell_origin_var), NC_CELL_ORIGIN_STR );
-    NCERRX( ncmpi_def_var(ncid, NC_CELL_LENGTHS_STR, NC_DOUBLE, 2, dims, &cell_lengths_var), NC_CELL_LENGTHS_STR );
+    NCERRX( ncmpi_def_var(ncid, NC_CELL_ORIGIN_STR, type_nc_real, 2, dims, &cell_origin_var), NC_CELL_ORIGIN_STR );
+    NCERRX( ncmpi_def_var(ncid, NC_CELL_LENGTHS_STR, type_nc_real, 2, dims, &cell_lengths_var), NC_CELL_LENGTHS_STR );
     dims[0] = frame_dim;
     dims[1] = cell_angular_dim;
-    NCERRX( ncmpi_def_var(ncid, NC_CELL_ANGLES_STR, NC_DOUBLE, 2, dims, &cell_angles_var), NC_CELL_ANGLES_STR );
+    NCERRX( ncmpi_def_var(ncid, NC_CELL_ANGLES_STR, type_nc_real, 2, dims, &cell_angles_var), NC_CELL_ANGLES_STR );
 
     // variables specified in the input file
     dims[0] = frame_dim;
@@ -405,10 +414,7 @@ void DumpNetCDFMPIIO::openfile()
       } else if (vtype[perat[i].field[0]] == Dump::BIGINT) {
         xtype = NC_INT64;
       } else {
-        if (double_precision)
-          xtype = NC_DOUBLE;
-        else
-          xtype = NC_FLOAT;
+        xtype = type_nc_real;
       }
 
       if (perat[i].dims == 1) {
@@ -418,6 +424,11 @@ void DumpNetCDFMPIIO::openfile()
         dims[2] = vector_dim[perat[i].dims];
         NCERRX( ncmpi_def_var(ncid, perat[i].name, xtype, 3, dims, &perat[i].var), perat[i].name );
       }
+
+      std::string unit = get_unit_for(update->unit_style, perat[i].quantity, error);
+      if (!unit.empty()) {
+        NCERR( ncmpi_put_att_text(ncid, perat[i].var, NC_UNITS_STR, unit.size(), unit.c_str()) );
+      }
     }
 
     // perframe variables
@@ -425,7 +436,7 @@ void DumpNetCDFMPIIO::openfile()
       Thermo *th = output->thermo;
       for (int i = 0; i < th->nfield; i++) {
         if (th->vtype[i] == Thermo::FLOAT) {
-          NCERRX( ncmpi_def_var(ncid, th->keyword[i], NC_DOUBLE, 1, dims, &thermovar[i]), th->keyword[i] );
+          NCERRX( ncmpi_def_var(ncid, th->keyword[i], type_nc_real, 1, dims, &thermovar[i]), th->keyword[i] );
         } else if (th->vtype[i] == Thermo::INT) {
           NCERRX( ncmpi_def_var(ncid, th->keyword[i], NC_INT, 1, dims, &thermovar[i]), th->keyword[i] );
         } else if (th->vtype[i] == Thermo::BIGINT) {
@@ -445,43 +456,18 @@ void DumpNetCDFMPIIO::openfile()
     NCERR( ncmpi_put_att_text(ncid, NC_GLOBAL, "program", 6, "LAMMPS") );
     NCERR( ncmpi_put_att_text(ncid, NC_GLOBAL, "programVersion", strlen(lmp->version), lmp->version) );
 
-    // units
-    if (!strcmp(update->unit_style, "lj")) {
-      NCERR( ncmpi_put_att_text(ncid, time_var, NC_UNITS_STR, 2, "lj") );
-      NCERR( ncmpi_put_att_text(ncid, cell_origin_var, NC_UNITS_STR, 2, "lj") );
-      NCERR( ncmpi_put_att_text(ncid, cell_lengths_var, NC_UNITS_STR, 2, "lj") );
-    } else if (!strcmp(update->unit_style, "real")) {
-      NCERR( ncmpi_put_att_text(ncid, time_var, NC_UNITS_STR, 11, "femtosecond") );
-      NCERR( ncmpi_put_att_text(ncid, cell_origin_var, NC_UNITS_STR, 8, "Angstrom") );
-      NCERR( ncmpi_put_att_text(ncid, cell_lengths_var, NC_UNITS_STR, 8, "Angstrom") );
-    } else if (!strcmp(update->unit_style, "metal")) {
-      NCERR( ncmpi_put_att_text(ncid, time_var, NC_UNITS_STR, 10, "picosecond") );
-      NCERR( ncmpi_put_att_text(ncid, cell_origin_var, NC_UNITS_STR, 8, "Angstrom") );
-      NCERR( ncmpi_put_att_text(ncid, cell_lengths_var, NC_UNITS_STR, 8, "Angstrom") );
-    } else if (!strcmp(update->unit_style, "si")) {
-      NCERR( ncmpi_put_att_text(ncid, time_var, NC_UNITS_STR, 6, "second") );
-      NCERR( ncmpi_put_att_text(ncid, cell_origin_var, NC_UNITS_STR, 5, "meter") );
-      NCERR( ncmpi_put_att_text(ncid, cell_lengths_var, NC_UNITS_STR, 5, "meter") );
-    } else if (!strcmp(update->unit_style, "cgs")) {
-      NCERR( ncmpi_put_att_text(ncid, time_var, NC_UNITS_STR, 6, "second") );
-      NCERR( ncmpi_put_att_text(ncid, cell_origin_var, NC_UNITS_STR, 10, "centimeter") );
-      NCERR( ncmpi_put_att_text(ncid, cell_lengths_var, NC_UNITS_STR, 10, "centimeter") );
-    } else if (!strcmp(update->unit_style, "electron")) {
-      NCERR( ncmpi_put_att_text(ncid, time_var, NC_UNITS_STR, 11, "femtosecond") );
-      NCERR( ncmpi_put_att_text(ncid, cell_origin_var, NC_UNITS_STR, 4, "Bohr") );
-      NCERR( ncmpi_put_att_text(ncid, cell_lengths_var, NC_UNITS_STR, 4, "Bohr") );
-    } else {
-      error->all(FLERR,"Unsupported unit style: {}", update->unit_style);
-    }
+    // units & scale
+    std::string unit = get_unit_for(update->unit_style, Quantity::TIME, error);
+    NCERR( ncmpi_put_att_text(ncid, time_var, NC_UNITS_STR, unit.size(), unit.c_str()) );
+
+    unit = get_unit_for(update->unit_style, Quantity::DISTANCE, error);
+    NCERR( ncmpi_put_att_text(ncid, cell_origin_var, NC_UNITS_STR, unit.size(), unit.c_str()) );
+    NCERR( ncmpi_put_att_text(ncid, cell_lengths_var, NC_UNITS_STR, unit.size(), unit.c_str()) );
 
     NCERR( ncmpi_put_att_text(ncid, cell_angles_var, NC_UNITS_STR, 6, "degree") );
 
-    d[0] = update->dt;
-    NCERR( ncmpi_put_att_double(ncid, time_var, NC_SCALE_FACTOR_STR, NC_DOUBLE, 1, d) );
-    d[0] = 1.0;
-    NCERR( ncmpi_put_att_double(ncid, cell_origin_var, NC_SCALE_FACTOR_STR, NC_DOUBLE, 1, d) );
-    d[0] = 1.0;
-    NCERR( ncmpi_put_att_double(ncid, cell_lengths_var, NC_SCALE_FACTOR_STR, NC_DOUBLE, 1, d) );
+    float scale[1] = {static_cast<float>(update->dt)};
+    NCERR( ncmpi_put_att_float(ncid, time_var, NC_SCALE_FACTOR_STR, NC_FLOAT, 1, scale) );
 
     /*
      * Finished with definition
@@ -502,16 +488,13 @@ void DumpNetCDFMPIIO::openfile()
       index[1] = 0;
       count[0] = 1;
       count[1] = 5;
-      NCERR( ncmpi_put_vara_text(ncid, cell_angular_var, index, count,
-                                 "alpha") );
+      NCERR( ncmpi_put_vara_text(ncid, cell_angular_var, index, count, "alpha") );
       index[0] = 1;
       count[1] = 4;
-      NCERR( ncmpi_put_vara_text(ncid, cell_angular_var, index, count,
-                                 "beta") );
+      NCERR( ncmpi_put_vara_text(ncid, cell_angular_var, index, count, "beta") );
       index[0] = 2;
       count[1] = 5;
-      NCERR( ncmpi_put_vara_text(ncid, cell_angular_var, index, count,
-                                 "gamma") );
+      NCERR( ncmpi_put_vara_text(ncid, cell_angular_var, index, count, "gamma") );
     }
 
     NCERR( ncmpi_end_indep_data(ncid) );
@@ -753,8 +736,7 @@ void DumpNetCDFMPIIO::write_time_and_cell()
 
 void DumpNetCDFMPIIO::write_data(int n, double *mybuf)
 {
-  MPI_Offset start[NC_MAX_VAR_DIMS], count[NC_MAX_VAR_DIMS];
-  MPI_Offset stride[NC_MAX_VAR_DIMS];
+  MPI_Offset start[LMP_MAX_VAR_DIMS], count[LMP_MAX_VAR_DIMS], stride[LMP_MAX_VAR_DIMS];
 
   if (!int_buffer) {
     n_buffer = std::max(1, n);
@@ -867,7 +849,12 @@ int DumpNetCDFMPIIO::modify_param(int narg, char **arg)
   if (strcmp(arg[iarg],"double") == 0) {
     iarg++;
     if (iarg >= narg) error->all(FLERR,"expected 'yes' or 'no' after 'double' keyword.");
-    double_precision = utils::logical(FLERR,arg[iarg],false,lmp) == 1;
+
+    if (utils::logical(FLERR,arg[iarg],false,lmp) == 1)
+      type_nc_real = NC_DOUBLE;
+    else
+      type_nc_real = NC_FLOAT;
+
     iarg++;
     return 2;
   } else if (strcmp(arg[iarg],"at") == 0) {
@@ -892,10 +879,9 @@ int DumpNetCDFMPIIO::modify_param(int narg, char **arg)
 void DumpNetCDFMPIIO::ncerr(int err, const char *descr, int line)
 {
   if (err != NC_NOERR) {
-    if (descr) error->one(FLERR,"NetCDF failed with error '{}' (while accessing '{}') "
-                          " in line {} of {}.", ncmpi_strerror(err), descr, line, __FILE__);
-    else error->one(FLERR,"NetCDF failed with error '{}' in line {} of {}.",
-                    ncmpi_strerror(err), line, __FILE__);
+    if (descr) error->one(__FILE__, line, "NetCDF failed with error '{}' (while accessing '{}') ",
+                          ncmpi_strerror(err), descr);
+    else error->one(__FILE__, line,"NetCDF failed with error '{}'.", ncmpi_strerror(err));
   }
 }
 

src/NETCDF/dump_netcdf_mpiio.h

@@ -30,9 +30,6 @@ DumpStyle(netcdf/mpiio,DumpNetCDFMPIIO);
 
 namespace LAMMPS_NS {
 
-const int NC_MPIIO_FIELD_NAME_MAX = 100;
-const int DUMP_NC_MPIIO_MAX_DIMS = 100;
-
 class DumpNetCDFMPIIO : public DumpCustom {
  public:
   DumpNetCDFMPIIO(class LAMMPS *, int, char **);
@@ -40,16 +37,18 @@ class DumpNetCDFMPIIO : public DumpCustom {
   virtual void write();
 
  private:
+  static constexpr int NC_MPIIO_FIELD_NAME_MAX = 100;
+  static constexpr int DUMP_NC_MPIIO_MAX_DIMS = 100;
+
   // per-atoms quantities (positions, velocities, etc.)
   struct nc_perat_t {
     int dims;                              // number of dimensions
     int field[DUMP_NC_MPIIO_MAX_DIMS];     // field indices corresponding to the dim.
     char name[NC_MPIIO_FIELD_NAME_MAX];    // field name
     int var;                               // NetCDF variable
+    int quantity;                          // type of the quantity
   };
 
-  typedef void (DumpNetCDFMPIIO::*funcptr_t)(void *);
-
   int framei;    // current frame index
   int blocki;    // current block index
   int ndata;     // number of data blocks to expect
@@ -61,8 +60,8 @@ class DumpNetCDFMPIIO : public DumpCustom {
 
   int *thermovar;    // NetCDF variables for thermo output
 
-  bool double_precision;    // write everything as double precision
-  bool thermo;              // write thermo output to netcdf file
+  int type_nc_real;    // netcdf type to use for real variables: float or double
+  bool thermo;         // write thermo output to netcdf file
 
   bigint n_buffer;          // size of buffer
   bigint *int_buffer;       // buffer for passing data to netcdf

src/NETCDF/netcdf_units.cpp

@@ -0,0 +1,145 @@
+/* ----------------------------------------------------------------------
+   LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
+   https://www.lammps.org/, Sandia National Laboratories
+   Steve Plimpton, sjplimp@sandia.gov
+
+   Copyright (2003) Sandia Corporation.  Under the terms of Contract
+   DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains
+   certain rights in this software.  This software is distributed under
+   the GNU General Public License.
+
+   See the README file in the top-level LAMMPS directory.
+------------------------------------------------------------------------- */
+
+/* ----------------------------------------------------------------------
+   Contributing author: Lars Pastewka (University of Freiburg)
+------------------------------------------------------------------------- */
+
+#if defined(LMP_HAS_NETCDF) || defined(LMP_HAS_PNETCDF)
+
+#include "netcdf_units.h"
+
+#include "error.h"
+
+using namespace LAMMPS_NS;
+
+std::string NetCDFUnits::get_unit_for(const char *unit_style, int quantity, Error *error)
+{
+  if (!strcmp(unit_style, "lj")) {
+    if (quantity == Quantity::UNKNOWN) {
+      return "";
+    } else {
+      return "lj";
+    }
+  } else if (!strcmp(unit_style, "real")) {
+    switch (quantity) {
+      case Quantity::UNKNOWN:
+        return "";
+      case Quantity::TIME:
+        return "femtosecond";
+      case Quantity::DISTANCE:
+        return "angstrom";
+      case Quantity::VELOCITY:
+        return "angstrom/femtosecond";
+      case Quantity::FORCE:
+        return "(Kcal/mol)/angstrom)";
+      case Quantity::DIPOLE_MOMENT:
+        return "e * angstrom";
+    }
+  } else if (!strcmp(unit_style, "metal")) {
+    switch (quantity) {
+      case Quantity::UNKNOWN:
+        return "";
+      case Quantity::TIME:
+        return "picosecond";
+      case Quantity::DISTANCE:
+        return "angstrom";
+      case Quantity::VELOCITY:
+        return "angstrom/picosecond";
+      case Quantity::FORCE:
+        return "eV/angstrom";
+      case Quantity::DIPOLE_MOMENT:
+        return "e * angstrom";
+    }
+  } else if (!strcmp(unit_style, "si")) {
+    switch (quantity) {
+      case Quantity::UNKNOWN:
+        return "";
+      case Quantity::TIME:
+        return "second";
+      case Quantity::DISTANCE:
+        return "meter";
+      case Quantity::VELOCITY:
+        return "meter/second";
+      case Quantity::FORCE:
+        return "Newton";
+      case Quantity::DIPOLE_MOMENT:
+        return "Coulomb * meter";
+    }
+  } else if (!strcmp(unit_style, "cgs")) {
+    switch (quantity) {
+      case Quantity::UNKNOWN:
+        return "";
+      case Quantity::TIME:
+        return "second";
+      case Quantity::DISTANCE:
+        return "centimeter";
+      case Quantity::VELOCITY:
+        return "centimeter/second";
+      case Quantity::FORCE:
+        return "dynes";
+      case Quantity::DIPOLE_MOMENT:
+        return "statcoul * cm";
+    }
+  } else if (!strcmp(unit_style, "electron")) {
+    switch (quantity) {
+      case Quantity::UNKNOWN:
+        return "";
+      case Quantity::TIME:
+        return "femtoseconds";
+      case Quantity::DISTANCE:
+        return "Bohr";
+      case Quantity::VELOCITY:
+        return "Bohr/atomic time units";
+      case Quantity::FORCE:
+        return "Hartree/Bohr";
+      case Quantity::DIPOLE_MOMENT:
+        return "Debye";
+    }
+  } else if (!strcmp(unit_style, "micro")) {
+    switch (quantity) {
+      case Quantity::UNKNOWN:
+        return "";
+      case Quantity::TIME:
+        return "microseconds";
+      case Quantity::DISTANCE:
+        return "micrometers";
+      case Quantity::VELOCITY:
+        return "micrometers/microsecond";
+      case Quantity::FORCE:
+        return "picogram * micrometer/microsecond^2";
+      case Quantity::DIPOLE_MOMENT:
+        return "picocoulomb * micrometer";
+    }
+  } else if (!strcmp(unit_style, "nano")) {
+    switch (quantity) {
+      case Quantity::UNKNOWN:
+        return "";
+      case Quantity::TIME:
+        return "nanoseconds";
+      case Quantity::DISTANCE:
+        return "nanometers";
+      case Quantity::VELOCITY:
+        return "nanometers/nanosecond";
+      case Quantity::FORCE:
+        return "attogram * nanometer/nanosecond^2";
+      case Quantity::DIPOLE_MOMENT:
+        return "e * nanometer";
+    }
+  }
+
+  error->all(FLERR, "Unsupported unit style: {}", unit_style);
+  return "";
+}
+
+#endif

src/NETCDF/netcdf_units.h

@@ -0,0 +1,49 @@
+/* -*- c++ -*- ----------------------------------------------------------
+   LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
+   https://www.lammps.org/, Sandia National Laboratories
+   Steve Plimpton, sjplimp@sandia.gov
+
+   Copyright (2003) Sandia Corporation.  Under the terms of Contract
+   DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains
+   certain rights in this software.  This software is distributed under
+   the GNU General Public License.
+
+   See the README file in the top-level LAMMPS directory.
+------------------------------------------------------------------------- */
+
+/* ----------------------------------------------------------------------
+   Contributing authors: Lars Pastewka (University of Freiburg), Guillaume Fraux (EPFL)
+------------------------------------------------------------------------- */
+
+#ifndef LMP_NETCDF_UNITS_H
+#define LMP_NETCDF_UNITS_H
+
+#if defined(LMP_HAS_NETCDF) || defined(LMP_HAS_PNETCDF)
+
+#include <string>
+
+namespace LAMMPS_NS {
+class Error;
+
+namespace NetCDFUnits {
+  // type of quantity for per-atom values (used to get the unit)
+  enum Quantity {
+    UNKNOWN = 0,
+    TIME,
+    DISTANCE,
+    VELOCITY,
+    FORCE,
+    DIPOLE_MOMENT,
+  };
+
+  // for compatibility with older NetCDF versions
+  static constexpr int LMP_MAX_VAR_DIMS = 1024;
+
+  // get the name of the unit for the given `quantity` in the given LAMMPS
+  // `unit_style` any error will be reported through `error`
+  std::string get_unit_for(const char *unit_style, int quantity, Error *error);
+}    // namespace NetCDFUnits
+}    // namespace LAMMPS_NS
+
+#endif
+#endif

src/OPENMP/angle_class2_omp.cpp

@@ -175,6 +175,8 @@ void AngleClass2OMP::eval(int nfrom, int nto, ThrData * const thr)
 
     // force & energy for bond-angle term
 
+    dr1 = r1 - ba_r1[type];
+    dr2 = r2 - ba_r2[type];
     aa1 = s * dr1 * ba_k1[type];
     aa2 = s * dr2 * ba_k2[type];
 

src/OPENMP/angle_table_omp.cpp

@@ -16,15 +16,16 @@
    Contributing author: Axel Kohlmeyer (Temple U)
 ------------------------------------------------------------------------- */
 
-#include "omp_compat.h"
 #include "angle_table_omp.h"
-#include <cmath>
+
 #include "atom.h"
 #include "comm.h"
 #include "force.h"
 #include "neighbor.h"
 
+#include <cmath>
 
+#include "omp_compat.h"
 #include "suffix.h"
 using namespace LAMMPS_NS;
 

src/OPENMP/bond_table_omp.cpp

@@ -16,16 +16,16 @@
    Contributing author: Axel Kohlmeyer (Temple U)
 ------------------------------------------------------------------------- */
 
-#include "omp_compat.h"
 #include "bond_table_omp.h"
+
 #include "atom.h"
 #include "comm.h"
 #include "force.h"
 #include "neighbor.h"
 
-
 #include <cmath>
 
+#include "omp_compat.h"
 #include "suffix.h"
 using namespace LAMMPS_NS;
 

src/OPENMP/fix_rigid_small_omp.h

@@ -21,12 +21,16 @@ FixStyle(rigid/small/omp,FixRigidSmallOMP);
 #define LMP_FIX_RIGID_SMALL_OMP_H
 
 #include "fix_rigid_small.h"
+#include "force.h"
 
 namespace LAMMPS_NS {
 
 class FixRigidSmallOMP : public FixRigidSmall {
  public:
-  FixRigidSmallOMP(class LAMMPS *lmp, int narg, char **args) : FixRigidSmall(lmp, narg, args){};
+  FixRigidSmallOMP(class LAMMPS *lmp, int narg, char **args) : FixRigidSmall(lmp, narg, args)
+  {
+    centroidstressflag = CENTROID_NOTAVAIL;
+  }
   virtual ~FixRigidSmallOMP(){};
 
   virtual void initial_integrate(int);

src/PHONON/Install.sh

@@ -26,16 +26,6 @@ action () {
   fi
 }
 
-# PHONON uses the parallel FFT wrapper used in PPPM,
-# so we must require the KSPACE package to be installed.
-
-if (test $1 = 1) then
-  if (test ! -e ../fft3d_wrap.h) then
-    echo "Must install KSPACE package with PHONON"
-    exit 1
-  fi
-fi
-
 # list of files with optional dependcies
 
 action fix_phonon.cpp fft3d_wrap.h

src/REAXFF/fix_reaxff_species.cpp

@@ -1,4 +1,3 @@
-// clang-format off
 /* ----------------------------------------------------------------------
    LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
    https://www.lammps.org/, Sandia National Laboratories
@@ -247,8 +246,8 @@ FixReaxFFSpecies::~FixReaxFFSpecies()
     if (posflag && multipos_opened) fclose(pos);
   }
 
-  modify->delete_compute("SPECATOM");
-  modify->delete_fix("SPECBOND");
+  modify->delete_compute(fmt::format("SPECATOM_{}",id));
+  modify->delete_fix(fmt::format("SPECBOND_{}",id));
 }
 
 /* ---------------------------------------------------------------------- */
@@ -288,22 +287,16 @@ void FixReaxFFSpecies::init()
   if (nvalid != update->ntimestep)
     nvalid = update->ntimestep+nfreq;
 
-  // check if this fix has been called twice
-  int count = 0;
-  for (int i = 0; i < modify->nfix; i++)
-    if (strcmp(modify->fix[i]->style,"reaxff/species") == 0) count++;
-  if (count > 1 && comm->me == 0)
-    error->warning(FLERR,"More than one fix reaxff/species");
-
   if (!setupflag) {
     // create a compute to store properties
-    modify->add_compute("SPECATOM all SPEC/ATOM q x y z vx vy vz abo01 abo02 abo03 abo04 "
-                        "abo05 abo06 abo07 abo08 abo09 abo10 abo11 abo12 abo13 abo14 "
-                        "abo15 abo16 abo17 abo18 abo19 abo20 abo21 abo22 abo23 abo24");
+    modify->add_compute(fmt::format("SPECATOM_{} all SPEC/ATOM q x y z vx vy vz abo01 abo02 "
+                                    "abo03 abo04 abo05 abo06 abo07 abo08 abo09 abo10 abo11 "
+                                    "abo12 abo13 abo14 abo15 abo16 abo17 abo18 abo19 abo20 "
+                                    "abo21 abo22 abo23 abo24",id));
 
     // create a fix to point to fix_ave_atom for averaging stored properties
-    auto fixcmd = fmt::format("SPECBOND all ave/atom {} {} {}",nevery,nrepeat,nfreq);
-    for (int i = 1; i < 32; ++i) fixcmd += " c_SPECATOM[" + std::to_string(i) + "]";
+    auto fixcmd = fmt::format("SPECBOND_{} all ave/atom {} {} {}",id,nevery,nrepeat,nfreq);
+    for (int i = 1; i < 32; ++i) fixcmd += fmt::format(" c_SPECATOM_{}[{}]",id,i);
     f_SPECBOND = (FixAveAtom *) modify->add_fix(fixcmd);
     setupflag = 1;
   }

src/REAXFF/reaxff_ffield.cpp

@@ -583,6 +583,7 @@ namespace ReaxFF {
       } catch (std::exception &e) {
         error->one(FLERR,e.what());
       }
+      fclose(fp);
     }
 
     // broadcast global parameters and allocate list on ranks != 0

src/RIGID/fix_rigid.cpp

@@ -68,6 +68,7 @@ FixRigid::FixRigid(LAMMPS *lmp, int narg, char **arg) :
   create_attribute = 1;
   dof_flag = 1;
   enforce2d_flag = 1;
+  centroidstressflag = CENTROID_NOTAVAIL;
 
   MPI_Comm_rank(world,&me);
   MPI_Comm_size(world,&nprocs);

src/RIGID/fix_rigid_small.cpp

@@ -73,6 +73,7 @@ FixRigidSmall::FixRigidSmall(LAMMPS *lmp, int narg, char **arg) :
   dof_flag = 1;
   enforce2d_flag = 1;
   stores_ids = 1;
+  centroidstressflag = CENTROID_AVAIL;
 
   MPI_Comm_rank(world,&me);
   MPI_Comm_size(world,&nprocs);
@@ -785,7 +786,7 @@ void FixRigidSmall::initial_integrate(int vflag)
   // forward communicate updated info of all bodies
 
   commflag = INITIAL;
-  comm->forward_comm_fix(this,26);
+  comm->forward_comm_fix(this,29);
 
   // set coords/orient and velocity/rotation of atoms in rigid bodies
 
@@ -879,6 +880,7 @@ void FixRigidSmall::enforce2d()
     b->xcm[2] = 0.0;
     b->vcm[2] = 0.0;
     b->fcm[2] = 0.0;
+    b->xgc[2] = 0.0;
     b->torque[0] = 0.0;
     b->torque[1] = 0.0;
     b->angmom[0] = 0.0;
@@ -1349,10 +1351,22 @@ void FixRigidSmall::set_xv()
       vr[4] = 0.5*x0*fc2;
       vr[5] = 0.5*x1*fc2;
 
-      v_tally(1,&i,1.0,vr);
+      double rlist[][3] = {x0, x1, x2};
+      double flist[][3] = {0.5*fc0, 0.5*fc1, 0.5*fc2};
+      v_tally(1,&i,1.0,vr,rlist,flist,b->xgc);
     }
   }
 
+  // update the position of geometric center
+  for (int ibody = 0; ibody < nlocal_body + nghost_body; ibody++) {
+    Body *b = &body[ibody];
+    MathExtra::matvec(b->ex_space,b->ey_space,b->ez_space,
+                      b->xgc_body,b->xgc);
+    b->xgc[0] += b->xcm[0];
+    b->xgc[1] += b->xcm[1];
+    b->xgc[2] += b->xcm[2];
+  }
+
   // set orientation, omega, angmom of each extended particle
 
   if (extended) {
@@ -1499,7 +1513,9 @@ void FixRigidSmall::set_v()
       vr[4] = 0.5*x0*fc2;
       vr[5] = 0.5*x1*fc2;
 
-      v_tally(1,&i,1.0,vr);
+      double rlist[][3] = {x0, x1, x2};
+      double flist[][3] = {0.5*fc0, 0.5*fc1, 0.5*fc2};
+      v_tally(1,&i,1.0,vr,rlist,flist,b->xgc);
     }
   }
 
@@ -1905,11 +1921,15 @@ void FixRigidSmall::setup_bodies_static()
   double **x = atom->x;
 
   double *xcm;
+  double *xgc;
 
   for (ibody = 0; ibody < nlocal_body+nghost_body; ibody++) {
     xcm = body[ibody].xcm;
+    xgc = body[ibody].xgc;
     xcm[0] = xcm[1] = xcm[2] = 0.0;
+    xgc[0] = xgc[1] = xgc[2] = 0.0;
     body[ibody].mass = 0.0;
+    body[ibody].natoms = 0;
   }
 
   double unwrap[3];
@@ -1924,22 +1944,31 @@ void FixRigidSmall::setup_bodies_static()
 
     domain->unmap(x[i],xcmimage[i],unwrap);
     xcm = b->xcm;
+    xgc = b->xgc;
     xcm[0] += unwrap[0] * massone;
     xcm[1] += unwrap[1] * massone;
     xcm[2] += unwrap[2] * massone;
+    xgc[0] += unwrap[0];
+    xgc[1] += unwrap[1];
+    xgc[2] += unwrap[2];
     b->mass += massone;
+    b->natoms++;
   }
 
   // reverse communicate xcm, mass of all bodies
 
   commflag = XCM_MASS;
-  comm->reverse_comm_fix(this,4);
+  comm->reverse_comm_fix(this,8);
 
   for (ibody = 0; ibody < nlocal_body; ibody++) {
     xcm = body[ibody].xcm;
+    xgc = body[ibody].xgc;
     xcm[0] /= body[ibody].mass;
     xcm[1] /= body[ibody].mass;
     xcm[2] /= body[ibody].mass;
+    xgc[0] /= body[ibody].natoms;
+    xgc[1] /= body[ibody].natoms;
+    xgc[2] /= body[ibody].natoms;
   }
 
   // set vcm, angmom = 0.0 in case inpfile is used
@@ -2124,12 +2153,22 @@ void FixRigidSmall::setup_bodies_static()
     // create initial quaternion
 
     MathExtra::exyz_to_q(ex,ey,ez,body[ibody].quat);
+
+    // convert geometric center position to principal axis coordinates
+    // xcm is wrapped, but xgc is not initially
+    xcm = body[ibody].xcm;
+    xgc = body[ibody].xgc;
+    double delta[3];
+    MathExtra::sub3(xgc,xcm,delta);
+    domain->minimum_image(delta);
+    MathExtra::transpose_matvec(ex,ey,ez,delta,body[ibody].xgc_body);
+    MathExtra::add3(xcm,delta,xgc);
   }
 
   // forward communicate updated info of all bodies
 
   commflag = INITIAL;
-  comm->forward_comm_fix(this,26);
+  comm->forward_comm_fix(this,29);
 
   // displace = initial atom coords in basis of principal axes
   // set displace = 0.0 for atoms not in any rigid body
@@ -2807,6 +2846,10 @@ void FixRigidSmall::set_molecule(int nlocalprev, tagint tagprev, int imol,
       if (nlocal_body == nmax_body) grow_body();
       Body *b = &body[nlocal_body];
       b->mass = onemols[imol]->masstotal;
+      b->natoms = onemols[imol]->natoms;
+      b->xgc[0] = xgeom[0];
+      b->xgc[1] = xgeom[1];
+      b->xgc[2] = xgeom[2];
 
       // new COM = Q (onemols[imol]->xcm - onemols[imol]->center) + xgeom
       // Q = rotation matrix associated with quat
@@ -2829,6 +2872,12 @@ void FixRigidSmall::set_molecule(int nlocalprev, tagint tagprev, int imol,
       MathExtra::quatquat(quat,onemols[imol]->quat,b->quat);
       MathExtra::q_to_exyz(b->quat,b->ex_space,b->ey_space,b->ez_space);
 
+      MathExtra::transpose_matvec(b->ex_space,b->ey_space,b->ez_space,
+                                  ctr2com_rotate,b->xgc_body);
+      b->xgc_body[0] *= -1;
+      b->xgc_body[1] *= -1;
+      b->xgc_body[2] *= -1;
+
       b->angmom[0] = b->angmom[1] = b->angmom[2] = 0.0;
       b->omega[0] = b->omega[1] = b->omega[2] = 0.0;
       b->conjqm[0] = b->conjqm[1] = b->conjqm[2] = b->conjqm[3] = 0.0;
@@ -2961,7 +3010,7 @@ int FixRigidSmall::pack_forward_comm(int n, int *list, double *buf,
                                      int /*pbc_flag*/, int * /*pbc*/)
 {
   int i,j;
-  double *xcm,*vcm,*quat,*omega,*ex_space,*ey_space,*ez_space,*conjqm;
+  double *xcm,*xgc,*vcm,*quat,*omega,*ex_space,*ey_space,*ez_space,*conjqm;
 
   int m = 0;
 
@@ -2973,6 +3022,10 @@ int FixRigidSmall::pack_forward_comm(int n, int *list, double *buf,
       buf[m++] = xcm[0];
       buf[m++] = xcm[1];
       buf[m++] = xcm[2];
+      xgc = body[bodyown[j]].xgc;
+      buf[m++] = xgc[0];
+      buf[m++] = xgc[1];
+      buf[m++] = xgc[2];
       vcm = body[bodyown[j]].vcm;
       buf[m++] = vcm[0];
       buf[m++] = vcm[1];
@@ -3048,7 +3101,7 @@ int FixRigidSmall::pack_forward_comm(int n, int *list, double *buf,
 void FixRigidSmall::unpack_forward_comm(int n, int first, double *buf)
 {
   int i,j,last;
-  double *xcm,*vcm,*quat,*omega,*ex_space,*ey_space,*ez_space,*conjqm;
+  double *xcm,*xgc,*vcm,*quat,*omega,*ex_space,*ey_space,*ez_space,*conjqm;
 
   int m = 0;
   last = first + n;
@@ -3060,6 +3113,10 @@ void FixRigidSmall::unpack_forward_comm(int n, int first, double *buf)
       xcm[0] = buf[m++];
       xcm[1] = buf[m++];
       xcm[2] = buf[m++];
+      xgc = body[bodyown[i]].xgc;
+      xgc[0] = buf[m++];
+      xgc[1] = buf[m++];
+      xgc[2] = buf[m++];
       vcm = body[bodyown[i]].vcm;
       vcm[0] = buf[m++];
       vcm[1] = buf[m++];
@@ -3135,7 +3192,7 @@ void FixRigidSmall::unpack_forward_comm(int n, int first, double *buf)
 int FixRigidSmall::pack_reverse_comm(int n, int first, double *buf)
 {
   int i,j,m,last;
-  double *fcm,*torque,*vcm,*angmom,*xcm;
+  double *fcm,*torque,*vcm,*angmom,*xcm, *xgc;
 
   m = 0;
   last = first + n;
@@ -3170,10 +3227,15 @@ int FixRigidSmall::pack_reverse_comm(int n, int first, double *buf)
     for (i = first; i < last; i++) {
       if (bodyown[i] < 0) continue;
       xcm = body[bodyown[i]].xcm;
+      xgc = body[bodyown[i]].xgc;
       buf[m++] = xcm[0];
       buf[m++] = xcm[1];
       buf[m++] = xcm[2];
+      buf[m++] = xgc[0];
+      buf[m++] = xgc[1];
+      buf[m++] = xgc[2];
       buf[m++] = body[bodyown[i]].mass;
+      buf[m++] = static_cast<double>(body[bodyown[i]].natoms);
     }
 
   } else if (commflag == ITENSOR) {
@@ -3208,7 +3270,7 @@ int FixRigidSmall::pack_reverse_comm(int n, int first, double *buf)
 void FixRigidSmall::unpack_reverse_comm(int n, int *list, double *buf)
 {
   int i,j,k;
-  double *fcm,*torque,*vcm,*angmom,*xcm;
+  double *fcm,*torque,*vcm,*angmom,*xcm, *xgc;
 
   int m = 0;
 
@@ -3245,10 +3307,15 @@ void FixRigidSmall::unpack_reverse_comm(int n, int *list, double *buf)
       j = list[i];
       if (bodyown[j] < 0) continue;
       xcm = body[bodyown[j]].xcm;
+      xgc = body[bodyown[j]].xgc;
       xcm[0] += buf[m++];
       xcm[1] += buf[m++];
       xcm[2] += buf[m++];
+      xgc[0] += buf[m++];
+      xgc[1] += buf[m++];
+      xgc[2] += buf[m++];
       body[bodyown[j]].mass += buf[m++];
+      body[bodyown[j]].natoms += static_cast<int>(buf[m++]);
     }
 
   } else if (commflag == ITENSOR) {