Merge branch 'master' of ssh://github.com/lammps/lammps into atomvec-custom

2020-05-06 15:53:21 -06:00
parent 7d9091cacc 14fb49c1cb
commit 8301c7a9c6
81 changed files with 7857 additions and 184 deletions
--- a/cmake/CMakeLists.txt
+++ b/cmake/CMakeLists.txt
@ -425,7 +425,12 @@ if(PKG_USER-ATC)
  if(LAMMPS_SIZES STREQUAL BIGBIG)
    message(FATAL_ERROR "The USER-ATC Package is not compatible with -DLAMMPS_BIGBIG")
  endif()
-  target_link_libraries(atc PRIVATE ${LAPACK_LIBRARIES} MPI::MPI_CXX)
+  target_link_libraries(atc PRIVATE ${LAPACK_LIBRARIES})
  if(BUILD_MPI)
    target_link_libraries(atc PRIVATE MPI::MPI_CXX)
  else()
    target_link_libraries(atc PRIVATE mpi_stubs)
  endif()
  target_include_directories(atc PRIVATE ${LAMMPS_SOURCE_DIR})
  target_compile_definitions(atc PRIVATE -DLAMMPS_${LAMMPS_SIZES})
 endif()
@ -549,10 +554,16 @@ if(BUILD_TOOLS)
  include(CheckGeneratorSupport)
  if(CMAKE_GENERATOR_SUPPORT_FORTRAN)
    include(CheckLanguage)
    check_language(Fortran)
    if(CMAKE_Fortran_COMPILER)
      enable_language(Fortran)
      add_executable(chain.x ${LAMMPS_TOOLS_DIR}/chain.f)
      target_link_libraries(chain.x PRIVATE ${CMAKE_Fortran_IMPLICIT_LINK_LIBRARIES})
      install(TARGETS chain.x DESTINATION ${CMAKE_INSTALL_BINDIR})
    else()
      message(WARNING "No suitable Fortran compiler found, skipping building 'chain.x'")
    endif()
  else()
    message(WARNING "CMake build doesn't support fortran, skipping building 'chain.x'")
  endif()
--- a/cmake/Modules/FindTBB_MALLOC.cmake
+++ b/cmake/Modules/FindTBB_MALLOC.cmake
@ -11,7 +11,8 @@
 # TBB Malloc
 find_path(TBB_MALLOC_INCLUDE_DIR NAMES tbb/tbb.h PATHS $ENV{TBBROOT}/include)
-find_library(TBB_MALLOC_LIBRARY NAMES tbbmalloc PATHS $ENV{TBBROOT}/lib/intel64/gcc4.7
+find_library(TBB_MALLOC_LIBRARY NAMES tbbmalloc PATHS $ENV{TBBROOT}/lib/intel64/gcc4.8
                                                      $ENV{TBBROOT}/lib/intel64/gcc4.7
                                                      $ENV{TBBROOT}/lib/intel64/gcc4.4
                                                      $ENV{TBBROOT}/lib/intel64/gcc4.1)
--- a/cmake/Modules/Packages/USER-H5MD.cmake
+++ b/cmake/Modules/Packages/USER-H5MD.cmake
@ -1,5 +1,8 @@
 enable_language(C)
 if(CMAKE_VERSION VERSION_GREATER_EQUAL 3.12)
  cmake_policy(SET CMP0074 NEW)
 endif()
 find_package(HDF5 REQUIRED)
 target_link_libraries(h5md PRIVATE ${HDF5_LIBRARIES})
 target_include_directories(h5md PUBLIC ${HDF5_INCLUDE_DIRS})
--- a/cmake/Modules/Packages/USER-NETCDF.cmake
+++ b/cmake/Modules/Packages/USER-NETCDF.cmake
@ -1,11 +1,18 @@
 # USER-NETCDF can use NetCDF, Parallel NetCDF (PNetCDF), or both. At least one necessary.
 # NetCDF library enables dump style "netcdf", while PNetCDF enables dump style "netcdf/mpiio"
 if(CMAKE_VERSION VERSION_GREATER_EQUAL 3.12)
  cmake_policy(SET CMP0074 NEW)
 endif()
 if(NOT BUILD_MPI)
  find_package(NetCDF REQUIRED)
 else()
  find_package(NetCDF)
  if(NETCDF_FOUND)
    find_package(PNetCDF)
-else(NETCDF_FOUND)
+  else()
    find_package(PNetCDF REQUIRED)
-endif(NETCDF_FOUND)
+  endif()
 endif()
 if(NETCDF_FOUND)
  target_link_libraries(lammps PRIVATE NetCDF::NetCDF)
--- a/cmake/Modules/Packages/VORONOI.cmake
+++ b/cmake/Modules/Packages/VORONOI.cmake
@ -24,7 +24,11 @@ if(DOWNLOAD_VORO)
  ExternalProject_Add(voro_build
    URL https://download.lammps.org/thirdparty/voro++-0.4.6.tar.gz
    URL_MD5 2338b824c3b7b25590e18e8df5d68af9
-    CONFIGURE_COMMAND "" BUILD_COMMAND make ${VORO_BUILD_OPTIONS} BUILD_IN_SOURCE 1 INSTALL_COMMAND ""
+    PATCH_COMMAND patch -b -p0 < ${LAMMPS_LIB_SOURCE_DIR}/voronoi/voro-make.patch
    CONFIGURE_COMMAND ""
    BUILD_COMMAND make ${VORO_BUILD_OPTIONS}
    BUILD_IN_SOURCE 1
    INSTALL_COMMAND ""
    BUILD_BYPRODUCTS <SOURCE_DIR>/src/libvoro++.a
    )
  ExternalProject_get_property(voro_build SOURCE_DIR)
--- a/cmake/presets/most.cmake
+++ b/cmake/presets/most.cmake
@ -5,10 +5,10 @@
 set(ALL_PACKAGES ASPHERE BODY CLASS2 COLLOID COMPRESS CORESHELL
        DIPOLE GRANULAR KSPACE MANYBODY MC MISC MOLECULE OPT PERI
        POEMS PYTHON QEQ REPLICA RIGID SHOCK SNAP SPIN SRD VORONOI
-        USER-CGDNA USER-CGSDK USER-COLVARS USER-DIFFRACTION
+        USER-BOCS USER-CGDNA USER-CGSDK USER-COLVARS USER-DIFFRACTION
-        USER-DPD USER-DRUDE USER-FEP USER-MEAMC USER-MESODPD
+        USER-DPD USER-DRUDE USER-EFF USER-FEP USER-MEAMC USER-MESODPD
        USER-MISC USER-MOFFF USER-OMP USER-PHONON USER-REACTION
-        USER-REAXC USER-SPH USER-SMD USER-UEF USER-YAFF)
+        USER-REAXC USER-SDPD USER-SPH USER-SMD USER-UEF USER-YAFF)
 foreach(PKG ${ALL_PACKAGES})
  set(PKG_${PKG} ON CACHE BOOL "" FORCE)
--- a/doc/lammps.1
+++ b/doc/lammps.1
@ -1,4 +1,4 @@
-.TH LAMMPS "15 April 2020" "2020-04-15"
+.TH LAMMPS "5 May 2020" "2020-05-5"
 .SH NAME
 .B LAMMPS
 \- Molecular Dynamics Simulator.
--- a/doc/src/Commands_compute.rst
+++ b/doc/src/Commands_compute.rst
@ -42,7 +42,7 @@ KOKKOS, o = USER-OMP, t = OPT.
   * :doc:`com <compute_com>`
   * :doc:`com/chunk <compute_com_chunk>`
   * :doc:`contact/atom <compute_contact_atom>`
-   * :doc:`coord/atom <compute_coord_atom>`
+   * :doc:`coord/atom (k) <compute_coord_atom>`
   * :doc:`damage/atom <compute_damage_atom>`
   * :doc:`dihedral <compute_dihedral>`
   * :doc:`dihedral/local <compute_dihedral_local>`
@ -161,5 +161,6 @@ KOKKOS, o = USER-OMP, t = OPT.
   * :doc:`torque/chunk <compute_torque_chunk>`
   * :doc:`vacf <compute_vacf>`
   * :doc:`vcm/chunk <compute_vcm_chunk>`
   * :doc:`viscosity/cos <compute_viscosity_cos>`
   * :doc:`voronoi/atom <compute_voronoi_atom>`
   * :doc:`xrd <compute_xrd>`
--- a/doc/src/Commands_fix.rst
+++ b/doc/src/Commands_fix.rst
@ -22,6 +22,7 @@ OPT.
 .. table_from_list::
   :columns: 5
   * :doc:`accelerate/cos <fix_accelerate_cos>`
   * :doc:`adapt <fix_adapt>`
   * :doc:`adapt/fep <fix_adapt_fep>`
   * :doc:`addforce <fix_addforce>`
--- a/doc/src/Howto_viscosity.rst
+++ b/doc/src/Howto_viscosity.rst
@ -1,10 +1,11 @@
 Calculate viscosity
 ===================
-The shear viscosity eta of a fluid can be measured in at least 5 ways
+The shear viscosity eta of a fluid can be measured in at least 6 ways
 using various options in LAMMPS.  See the examples/VISCOSITY directory
 for scripts that implement the 5 methods discussed here for a simple
-Lennard-Jones fluid model.  Also, see the :doc:`Howto kappa <Howto_kappa>` doc page for an analogous discussion for
+Lennard-Jones fluid model and 1 method for SPC/E water model.
 Also, see the :doc:`Howto kappa <Howto_kappa>` doc page for an analogous discussion for
 thermal conductivity.
 Eta is a measure of the propensity of a fluid to transmit momentum in
@ -130,9 +131,25 @@ time-integrated momentum fluxes play the role of Cartesian
 coordinates, whose mean-square displacement increases linearly
 with time at sufficiently long times.
 The sixth is periodic perturbation method. It is also a non-equilibrium MD method.
 However, instead of measure the momentum flux in response of applied velocity gradient,
 it measures the velocity profile in response of applied stress.
 A cosine-shaped periodic acceleration is added to the system via the
 :doc:`fix accelerate/cos <fix_accelerate_cos>` command,
 and the :doc:`compute viscosity/cos<compute_viscosity_cos>` command is used to monitor the
 generated velocity profile and remove the velocity bias before thermostatting.
 .. note::
    An article by :ref:`(Hess) <Hess3>` discussed the accuracy and efficiency of these methods.
 ----------
 .. _Daivis-viscosity:
 **(Daivis and Todd)** Daivis and Todd, Nonequilibrium Molecular Dynamics (book),
 Cambridge University Press, https://doi.org/10.1017/9781139017848, (2017).
 .. _Hess3:
 **(Hess)** Hess, B. The Journal of Chemical Physics 2002, 116 (1), 209-217.
--- a/doc/src/compute.rst
+++ b/doc/src/compute.rst
@ -307,6 +307,7 @@ The individual style names on the :doc:`Commands compute <Commands_compute>` doc
 * :doc:`torque/chunk <compute_torque_chunk>` - torque applied on each chunk
 * :doc:`vacf <compute_vacf>` - velocity auto-correlation function of group of atoms
 * :doc:`vcm/chunk <compute_vcm_chunk>` - velocity of center-of-mass for each chunk
 * :doc:`viscosity/cos <compute_viscosity_cos>` - velocity profile under cosine-shaped acceleration
 * :doc:`voronoi/atom <compute_voronoi_atom>` - Voronoi volume and neighbors for each atom
 * :doc:`xrd <compute_xrd>` - x-ray diffraction intensity on a mesh of reciprocal lattice nodes
--- a/doc/src/compute_coord_atom.rst
+++ b/doc/src/compute_coord_atom.rst
@ -3,6 +3,9 @@
 compute coord/atom command
 ==========================
 compute coord/atom/kk command
 ===================================
 Syntax
 """"""
@ -109,6 +112,30 @@ too frequently.
   :doc:`special_bonds <special_bonds>` command that includes all pairs in
   the neighbor list.
 ----------
 Styles with a *gpu*\ , *intel*\ , *kk*\ , *omp*\ , or *opt* suffix are
 functionally the same as the corresponding style without the suffix.
 They have been optimized to run faster, depending on your available
 hardware, as discussed on the :doc:`Speed packages <Speed_packages>` doc
 page.  The accelerated styles take the same arguments and should
 produce the same results, except for round-off and precision issues.
 These accelerated styles are part of the GPU, USER-INTEL, KOKKOS,
 USER-OMP and OPT packages, respectively.  They are only enabled if
 LAMMPS was built with those packages.  See the :doc:`Build package <Build_package>` doc page for more info.
 You can specify the accelerated styles explicitly in your input script
 by including their suffix, or you can use the :doc:`-suffix command-line switch <Run_options>` when you invoke LAMMPS, or you can use the
 :doc:`suffix <suffix>` command in your input script.
 See the :doc:`Speed packages <Speed_packages>` doc page for more
 instructions on how to use the accelerated styles effectively.
 ----------
 **Output info:**
 For *cstyle* cutoff, this compute can calculate a per-atom vector or
--- a/doc/src/compute_orientorder_atom.rst
+++ b/doc/src/compute_orientorder_atom.rst
@ -19,13 +19,14 @@ Syntax
  .. parsed-literal::
-     keyword = *cutoff* or *nnn* or *degrees* or *components*
+     keyword = *cutoff* or *nnn* or *degrees* or *components* or *chunksize*
       *cutoff* value = distance cutoff
       *nnn* value = number of nearest neighbors
       *degrees* values = nlvalues, l1, l2,...
       *wl* value = yes or no
       *wl/hat* value = yes or no
       *components* value = ldegree
       *chunksize* value = number of atoms in each pass 
 Examples
 """"""""
@ -107,6 +108,14 @@ in conjunction with :doc:`compute coord_atom <compute_coord_atom>` to
 calculate the ten Wolde's criterion to identify crystal-like
 particles, as discussed in :ref:`ten Wolde <tenWolde2>`.
 The optional keyword *chunksize* is only applicable when using the
 the KOKKOS package and is ignored otherwise. This keyword controls
 the number of atoms in each pass used to compute the bond-orientational
 order parameters and is used to avoid running out of memory. For example
 if there are 4000 atoms in the simulation and the *chunksize*
 is set to 2000, the parameter calculation will be broken up
 into two passes.
 The value of :math:`Q_l` is set to zero for atoms not in the
 specified compute group, as well as for atoms that have less than
 *nnn* neighbors within the distance cutoff, unless *nnn* is NULL.
@ -192,7 +201,7 @@ Default
 The option defaults are *cutoff* = pair style cutoff, *nnn* = 12,
 *degrees* = 5 4 6 8 10 12 i.e. :math:`Q_4`, :math:`Q_6`, :math:`Q_8`, :math:`Q_{10}`, and :math:`Q_{12}`,
-*wl* = no, *wl/hat* = no, and *components* off
+*wl* = no, *wl/hat* = no, *components* off, and *chunksize* = 2000
 ----------
--- a/doc/src/compute_viscosity_cos.rst
+++ b/doc/src/compute_viscosity_cos.rst
@ -0,0 +1,156 @@
 .. index:: compute viscosity/cos
 compute viscosity/cos command
 =============================
 Syntax
 """"""
 .. parsed-literal::
   compute ID group-ID viscosity/cos
 * ID, group-ID are documented in :doc:`compute <compute>` command
 * viscosity/cos = style name of this compute command
 Examples
 """"""""
 .. code-block:: LAMMPS
   compute  cos all viscosity/cos
   variable V equal c_cos[7]
   variable A equal 0.02E-5
   variable density equal density
   variable lz equal lz
   variable reciprocalViscosity equal v_V/${A}/v_density*39.4784/v_lz/v_lz*100
 Description
 """""""""""
 Define a computation that calculates the velocity amplitude of a group of atoms
 with an cosine-shaped velocity profile and the temperature of them
 after subtracting out the velocity profile before computing the kinetic energy.
 A compute of this style can be used by any command that computes a temperature,
 e.g. :doc:`thermo_modify <thermo_modify>`, :doc:`fix npt <fix_nh>`, etc.
 This command together with :doc:`fix_accelerate/cos<fix_accelerate_cos>`
 enables viscosity calculation with periodic perturbation method,
 as described by :ref:`Hess<Hess1>`.
 An acceleration along the x-direction is applied to the simulation system
 by using :doc:`fix_accelerate/cos<fix_accelerate_cos>` command.
 The acceleration is a periodic function along the z-direction:
 .. math::
   a_{x}(z) = A \cos \left(\frac{2 \pi z}{l_{z}}\right)
 where :math:`A` is the acceleration amplitude, :math:`l_z` is the z-length
 of the simulation box. At steady state, the acceleration generates
 a velocity profile:
 .. math::
   v_{x}(z) = V \cos \left(\frac{2 \pi z}{l_{z}}\right)
 The generated velocity amplitude :math:`V` is related to the
 shear viscosity :math:`\eta` by:
 .. math::
   V = \frac{A \rho}{\eta}\left(\frac{l_{z}}{2 \pi}\right)^{2}
 and it can be obtained from ensemble average of the velocity profile:
 .. math::
   V = \frac{\sum_i 2 m_{i} v_{i, x} \cos \left(\frac{2 \pi z_i}{l_{z}}\right)}{\sum_i m_{i}}
 where :math:`m_i`, :math:`v_{i,x}` and :math:`z_i` are the mass,
 x-component velocity and z coordinate of a particle.
 After the cosine-shaped collective velocity in :math:`x` direction
 has been subtracted for each atom, the temperature is calculated by the formula
 KE = dim/2 N k T, where KE = total kinetic energy of the group of
 atoms (sum of 1/2 m v\^2), dim = 2 or 3 = dimensionality of the
 simulation, N = number of atoms in the group, k = Boltzmann constant,
 and T = temperature.
 A kinetic energy tensor, stored as a 6-element vector, is also
 calculated by this compute for use in the computation of a pressure
 tensor. The formula for the components of the tensor is the same as
 the above formula, except that v\^2 is replaced by vx\*vy for the xy
 component, etc. The 6 components of the vector are ordered xx, yy,
 zz, xy, xz, yz.
 The number of atoms contributing to the temperature is assumed to be
 constant for the duration of the run; use the *dynamic* option of the
 :doc:`compute_modify <compute_modify>` command if this is not the case.
 However, in order to get meaningful result, the group ID of this compute should be all.
 The removal of the cosine-shaped velocity component by this command is
 essentially computing the temperature after a "bias" has been removed
 from the velocity of the atoms.  If this compute is used with a fix
 command that performs thermostatting then this bias will be subtracted
 from each atom, thermostatting of the remaining thermal velocity will
 be performed, and the bias will be added back in.  Thermostatting
 fixes that work in this way include :doc:`fix nvt <fix_nh>`, :doc:`fix temp/rescale <fix_temp_rescale>`, :doc:`fix temp/berendsen <fix_temp_berendsen>`, and :doc:`fix langevin <fix_langevin>`.
 This compute subtracts out degrees-of-freedom due to fixes that
 constrain molecular motion, such as :doc:`fix shake <fix_shake>` and
 :doc:`fix rigid <fix_rigid>`.  This means the temperature of groups of
 atoms that include these constraints will be computed correctly.  If
 needed, the subtracted degrees-of-freedom can be altered using the
 *extra* option of the :doc:`compute_modify <compute_modify>` command.
 See the :doc:`Howto thermostat <Howto_thermostat>` doc page for a
 discussion of different ways to compute temperature and perform
 thermostatting.
 ----------
 **Output info:**
 This compute calculates a global scalar (the temperature) and a global
 vector of length 7, which can be accessed by indices 1-7.
 The first 6 elements of the vector are the KE tensor,
 and the 7-th is the cosine-shaped velocity amplitude :math:`V`,
 which can be used to calculate the reciprocal viscosity, as shown in the example.
 These values can be used by any command that uses global scalar or
 vector values from a compute as input.
 See the :doc:`Howto output <Howto_output>` doc page for an overview of LAMMPS output options.
 The scalar value calculated by this compute is "intensive".  The
 first 6 elements of vector values are "extensive",
 and the 7-th element of vector values is "intensive".
 The scalar value will be in temperature :doc:`units <units>`.  The
 first 6 elements of vector values will be in energy :doc:`units <units>`.
 The 7-th element of vector value will be in velocity :doc:`units <units>`.
 Restrictions
 """"""""""""
 This command is only available when LAMMPS was built with the USER-MISC package.
 Since this compute depends on :doc:`fix accelerate/cos <fix_accelerate_cos>` which can
 only work for 3d systems, it cannot be used for 2d systems.
 Related commands
 """"""""""""""""
 :doc:`fix accelerate/cos <fix_accelerate_cos>`
 Default
 """""""
 none
 ----------
 .. _Hess1:
 **(Hess)** Hess, B. The Journal of Chemical Physics 2002, 116 (1), 209-217.
--- a/doc/src/fix.rst
+++ b/doc/src/fix.rst
@ -165,6 +165,7 @@ The individual style names on the :doc:`Commands fix <Commands_fix>` doc
 page are followed by one or more of (g,i,k,o,t) to indicate which
 accelerated styles exist.
 * :doc:`accelerate/cos <fix_accelerate_cos>` - apply cosine-shaped acceleration to atoms
 * :doc:`adapt <fix_adapt>` - change a simulation parameter over time
 * :doc:`adapt/fep <fix_adapt_fep>` - enhanced version of fix adapt
 * :doc:`addforce <fix_addforce>` - add a force to each atom
--- a/doc/src/fix_accelerate_cos.rst
+++ b/doc/src/fix_accelerate_cos.rst
@ -0,0 +1,104 @@
 .. index:: fix accelerate/cos
 fix accelerate/cos command
 ==========================
 Syntax
 """"""
 .. parsed-literal::
   fix ID group-ID accelerate value
 * ID, group-ID are documented in :doc:`fix <fix>` command
 * accelerate/cos = style name of this fix command
 * value = amplitude of acceleration (in unit of force/mass)
 Examples
 """"""""
 .. code-block:: LAMMPS
   fix 1 all accelerate/cos 0.02e-5
 Description
 """""""""""
 Give each atom a acceleration in x-direction based on its z coordinate.
 The acceleration is a periodic function along the z-direction:
 .. math::
   a_{x}(z) = A \cos \left(\frac{2 \pi z}{l_{z}}\right)
 where :math:`A` is the acceleration amplitude, :math:`l_z` is the z-length
 of the simulation box.
 At steady state, the acceleration generates a velocity profile:
 .. math::
   v_{x}(z) = V \cos \left(\frac{2 \pi z}{l_{z}}\right)
 The generated velocity amplitude :math:`V` is related to the
 shear viscosity :math:`\eta` by:
 .. math::
   V = \frac{A \rho}{\eta}\left(\frac{l_{z}}{2 \pi}\right)^{2}
 and it can be obtained from ensemble average of the velocity profile:
 .. math::
   V = \frac{\sum_i 2 m_{i} v_{i, x} \cos \left(\frac{2 \pi z_i}{l_{z}}\right)}{\sum_i m_{i}}
 where :math:`m_i`, :math:`v_{i,x}` and :math:`z_i` are the mass,
 x-component velocity and z coordinate of a particle.
 The velocity amplitude :math:`V` can be calculated with :doc:`compute viscosity/cos <compute_viscosity_cos>`,
 which enables viscosity calculation with periodic perturbation method,
 as described by :ref:`Hess<Hess2>`.
 Because the applied acceleration drives the system away from equilibration,
 the calculated shear viscosity is lower than the intrinsic viscosity
 due to the shear-thinning effect.
 Extrapolation to zero acceleration should generally be performed to
 predict the zero-shear viscosity.
 As the shear stress decreases, the signal-noise ratio decreases rapidly,
 the simulation time must be extended accordingly to get converged result.
 In order to get meaningful result, the group ID of this fix should be all.
 ----------
 **Restart, fix_modify, output, run start/stop, minimize info:**
 No information about this fix is written to binary restart files.
 None of the fix_modify options are relevant to this fix.
 No global or per-atom quantities are stored by this fix for access by various output commands.
 No parameter of this fix can be used with the start/stop keywords of the run command.
 This fix is not invoked during energy minimization.
 Restrictions
 """"""""""""
 This command is only available when LAMMPS was built with the USER-MISC package.
 Since this fix depends on the z-coordinate of atoms, it cannot be used in 2d simulations.
 Related commands
 """"""""""""""""
 :doc:`compute viscosity/cos <compute_viscosity_cos>`
 Default
 """""""
 none
 ----------
 .. _Hess2:
 **(Hess)** Hess, B. The Journal of Chemical Physics 2002, 116 (1), 209-217.
--- a/doc/src/fix_ave_correlate_long.rst
+++ b/doc/src/fix_ave_correlate_long.rst
@ -26,7 +26,7 @@ Syntax
       v_name = global value calculated by an equal-style variable with name
 * zero or more keyword/arg pairs may be appended
-* keyword = *type* or *start* or *file* or *overwrite* or *title1* or *title2* or *ncorr* or *p* or *m*
+* keyword = *type* or *start* or *file* or *overwrite* or *title1* or *title2* or *ncorr* or *nlen* or *ncount*
  .. parsed-literal::
--- a/doc/src/fix_halt.rst
+++ b/doc/src/fix_halt.rst
@ -160,7 +160,7 @@ the :doc:`run <run>` command.
 Restrictions
 """"""""""""
-The *diskfree* attribute is currently only supported on Linux and MacOS.
+The *diskfree* attribute is currently only supported on Linux, MacOSX, and BSD.
 Related commands
 """"""""""""""""
--- a/doc/src/pair_reaxc.rst
+++ b/doc/src/pair_reaxc.rst
@ -21,12 +21,13 @@ Syntax
  .. parsed-literal::
-     keyword = *checkqeq* or *lgvdw* or *safezone* or *mincap*
+     keyword = *checkqeq* or *lgvdw* or *safezone* or *mincap* or *minhbonds*
       *checkqeq* value = *yes* or *no* = whether or not to require qeq/reax fix
       *enobonds* value = *yes* or *no* = whether or not to tally energy of atoms with no bonds
       *lgvdw* value = *yes* or *no* = whether or not to use a low gradient vdW correction
       *safezone* = factor used for array allocation
       *mincap* = minimum size for array allocation
       *minhbonds* = minimum size use for storing hydrogen bonds
 Examples
 """"""""
@ -146,11 +147,11 @@ zero.  The latter behavior is usual not desired, as it causes
 discontinuities in the potential energy when the bonding of an atom
 drops to zero.
-Optional keywords *safezone* and *mincap* are used for allocating
+Optional keywords *safezone*\ , *mincap*\ , and *minhbonds* are used
-reax/c arrays.  Increasing these values can avoid memory problems,
+for allocating reax/c arrays.  Increasing these values can avoid memory
-such as segmentation faults and bondchk failed errors, that could
+problems, such as segmentation faults and bondchk failed errors, that
-occur under certain conditions. These keywords are not used by the
+could occur under certain conditions. These keywords are not used by
-Kokkos version, which instead uses a more robust memory allocation
+the Kokkos version, which instead uses a more robust memory allocation
 scheme that checks if the sizes of the arrays have been exceeded and
 automatically allocates more memory.
@ -352,7 +353,7 @@ Default
 """""""
 The keyword defaults are checkqeq = yes, enobonds = yes, lgvdw = no,
-safezone = 1.2, mincap = 50.
+safezone = 1.2, mincap = 50, minhbonds = 25.
 ----------
--- a/doc/utils/sphinx-config/false_positives.txt
+++ b/doc/utils/sphinx-config/false_positives.txt
@ -1789,6 +1789,7 @@ Militzer
 Minary
 mincap
 Mindlin
 minhbonds
 mingw
 minima
 minimizations
--- a/examples/VISCOSITY/README
+++ b/examples/VISCOSITY/README
@ -1,5 +1,6 @@
-This directory has 5 scripts that compute the viscosity (eta) of a
+This directory has 6 scripts that compute the viscosity (eta) of fluid
-Lennard-Jones fluid using 5 different methods.  See the discussion in
+using 6 different methods. 5 of them are for a Lennard-Jones fluid
 and the last one is for SPC/E water model. See the discussion in
 Section 6.21 of the manual for an overview of the methods and pointers
 to doc pages for the commands which implement them.  Citations for the
 various methods can also be found in the manual.
@ -10,7 +11,7 @@ enough to generate good statistics and highly accurate results.
 -------------
-These are the 5 methods for computing viscosity.  The first 3 are
+These are the 5 methods for computing viscosity of a LJ fluid.  The first 3 are
 non-equilibrium methods; the last 2 are equilibrium methods.
 in.wall = move a wall to shear the fluid between two walls
@ -89,3 +90,18 @@ heat/flux doc page - the resulting value prints at the end of the run
 and is in the log file
 eta = 1.07
 -------------
 in.cos.1000SPCE is an example script of using cosine periodic perturbation method
 to calculate the viscosity of SPC/E water model.
 The reciprocal of eta is computed within the script, and printed out as v_invVis
 in thermo_style command. The result will converge after hundreds of picoseconds.
 Then eta is obtained from the reciprocal of time average of v_invVis.
 eta = 0.75 mPa*s
 Note that the calculated viscosity by this method decreases with increased acceleration.
 It is therefore generally necessary to perform calculation at different accelerations
 and extrapolate the viscosity to zero shear.
--- a/examples/VISCOSITY/data.cos.1000SPCE
+++ b/examples/VISCOSITY/data.cos.1000SPCE
--- a/examples/VISCOSITY/in.cos.1000SPCE
+++ b/examples/VISCOSITY/in.cos.1000SPCE
@ -0,0 +1,53 @@
 # DFF generated Lammps input file
 units          real
 atom_style     full
 boundary       p p p
 pair_style     lj/cut/coul/long     10.0
 pair_modify    mix arithmetic
 pair_modify    tail yes
 kspace_style   pppm 1.0e-4
 dielectric     1.0
 special_bonds  amber
 bond_style     harmonic
 angle_style    harmonic
 dihedral_style none
 improper_style none
 read_data      data.cos.1000SPCE
 variable T equal 300
 variable P equal 1.0
 velocity       all create ${T} 12345 mom yes rot yes dist gaussian
 timestep       1.0
 # Constraint ##################################
 fix com        all momentum 100 linear 1 1 1
 fix rigid      all shake 1e-4 20 0 b 1 a 1
 # Viscosity  ##################################
 variable       A equal 0.05e-5 # angstrom/fs^2
 fix            cos all accelerate/cos ${A}
 compute        cos all viscosity/cos
 variable       density equal density
 variable       lz equal lz
 variable       vMax equal c_cos[7] # velocity of atoms at z=0
 variable       invVis equal v_vMax/${A}/v_density*39.4784/v_lz/v_lz*100 # reciprocal of viscosity 1/Pa/s
 fix            npt all npt temp ${T} ${T} 100 iso ${P} ${P} 1000
 fix_modify     npt temp cos
 thermo_style   custom step cpu temp press pe density v_vMax v_invVis
 thermo_modify  temp cos
 thermo         100
 ################################################
 dump           1 all custom 10000 dump.lammpstrj id mol type element q xu yu zu
 dump_modify    1 sort id element O H
 run            2000
--- a/examples/VISCOSITY/log.30Apr20.cos.1000SPCE.g++.1
+++ b/examples/VISCOSITY/log.30Apr20.cos.1000SPCE.g++.1
@ -0,0 +1,160 @@
 LAMMPS (3 Mar 2020)
 OMP_NUM_THREADS environment is not set. Defaulting to 1 thread. (src/comm.cpp:94)
  using 1 OpenMP thread(s) per MPI task
 # DFF generated Lammps input file
 units          real
 atom_style     full
 boundary       p p p
 pair_style     lj/cut/coul/long     10.0
 pair_modify    mix arithmetic
 pair_modify    tail yes
 kspace_style   pppm 1.0e-4
 dielectric     1.0
 special_bonds  amber
 bond_style     harmonic
 angle_style    harmonic
 dihedral_style none
 improper_style none
 read_data      data.1000SPCE.lmp
  orthogonal box = (0 0 0) to (31.043 31.043 31.043)
  2 by 2 by 2 MPI processor grid
  reading atoms ...
  3000 atoms
  scanning bonds ...
  2 = max bonds/atom
  scanning angles ...
  1 = max angles/atom
  reading bonds ...
  2000 bonds
  reading angles ...
  1000 angles
  2 = max # of 1-2 neighbors
  1 = max # of 1-3 neighbors
  1 = max # of 1-4 neighbors
  2 = max # of special neighbors
  special bonds CPU = 0.00114917 secs
  read_data CPU = 0.00953543 secs
 variable T equal 300
 variable P equal 1.0
 velocity       all create ${T} 12345 mom yes rot yes dist gaussian
 velocity       all create 300 12345 mom yes rot yes dist gaussian
 timestep       1.0
 # Constraint ##################################
 fix com        all momentum 100 linear 1 1 1
 fix rigid      all shake 1e-4 20 0 b 1 a 1
  0 = # of size 2 clusters
  0 = # of size 3 clusters
  0 = # of size 4 clusters
  1000 = # of frozen angles
  find clusters CPU = 0.000442737 secs
 # Viscosity  ##################################
 variable       A equal 0.02e-5 # angstrom/fs^2
 fix            cos all accelerate/cos ${A}
 fix            cos all accelerate/cos 2e-07
 compute        cos all viscosity/cos
 variable       density equal density
 variable       lz equal lz
 variable       vMax equal c_cos[7] # velocity of atoms at z=0
 variable       invVis equal v_vMax/${A}/v_density*39.4784/v_lz/v_lz*100 # reciprocal of viscosity 1/Pa/s
 variable       invVis equal v_vMax/2e-07/v_density*39.4784/v_lz/v_lz*100 
 fix            npt all npt temp ${T} ${T} 100 iso ${P} ${P} 1000
 fix            npt all npt temp 300 ${T} 100 iso ${P} ${P} 1000
 fix            npt all npt temp 300 300 100 iso ${P} ${P} 1000
 fix            npt all npt temp 300 300 100 iso 1 ${P} 1000
 fix            npt all npt temp 300 300 100 iso 1 1 1000
 fix_modify     npt temp cos
 thermo_style   custom step cpu temp press pe density v_vMax v_invVis
 thermo_modify  temp cos
 thermo         100
 ################################################
 dump           1 all custom 10000 dump.lammpstrj id mol type element q xu yu zu
 dump_modify    1 sort id element O H
 run            2000
 PPPM initialization ...
  using 12-bit tables for long-range coulomb (src/kspace.cpp:332)
  G vector (1/distance) = 0.263539
  grid = 16 16 16
  stencil order = 5
  estimated absolute RMS force accuracy = 0.0325342
  estimated relative force accuracy = 9.79757e-05
  using double precision MKL FFT
  3d grid and FFT values/proc = 3375 512
 Neighbor list info ...
  update every 1 steps, delay 10 steps, check yes
  max neighbors/atom: 2000, page size: 100000
  master list distance cutoff = 12
  ghost atom cutoff = 12
  binsize = 6, bins = 6 6 6
  1 neighbor lists, perpetual/occasional/extra = 1 0 0
  (1) pair lj/cut/coul/long, perpetual
      attributes: half, newton on
      pair build: half/bin/newton
      stencil: half/bin/3d/newton
      bin: standard
 Per MPI rank memory allocation (min/avg/max) = 10.14 | 10.14 | 10.14 Mbytes
 Step CPU Temp Press PotEng Density v_vMax v_invVis 
       0            0    450.04468    9838.6886    -7651.736   0.99999331 0.0001293705    2649.9663 
     100   0.33736925    497.65155    2024.4827   -8342.2499   0.98595028 0.00019602427    4034.2452 
     200   0.67116638    446.48518    27.075082   -8908.9684    0.9652009 0.00027615455    5723.7907 
     300   0.99760895    401.79875   -776.99871   -9381.8986   0.95205822 0.00019626685    4086.6103 
     400    1.3239019    369.65373    -510.5336   -9803.6463   0.94817309 0.00022998514    4795.2142 
     500    1.6488092    343.35807   -936.31982   -10146.023   0.94840581 0.0001434423    2990.5423 
     600    1.9826063    319.86131   -1381.3302   -10405.812   0.95459202 6.6411532e-05    1381.5767 
     700    2.3172637    307.74606   -98.775733     -10643.5   0.96669652 0.00010446317    2164.0664 
     800    2.6562841    305.14214   -540.57804   -10865.742   0.97808367 5.4381233e-05    1122.1765 
     900    2.9938415    288.01316    639.00486    -10925.39   0.98686357 0.00010878474    2238.1355 
    1000     3.327893    295.07773   -226.06503   -11033.826   0.99128496 0.00011935058    2451.8608 
    1100    3.6618862    299.21578    306.34231   -11049.152   0.99552203 8.9538943e-05    1836.8166 
    1200    3.9984287    301.82462    85.804646   -11013.564   0.99713434 0.00015912276      3262.51 
    1300    4.3320735     308.6009    268.08897   -11009.836   0.99695358 0.00026212596      5374.72 
    1400     4.668875    298.36903   -258.75495   -10962.299   0.99503447 0.00033087355    6788.7027 
    1500    5.0003694    299.96073    99.512082   -10980.551   0.99315631 0.00033996557    6979.6425 
    1600    5.3367337    304.18018   -500.65441   -11002.054    0.9914558 0.00039075642    8026.9849 
    1700    5.6780828    301.63978   -499.07458    -10992.88   0.99234354 0.00038101175    7824.4738 
    1800    6.0140638    303.25858    640.03432   -11053.335   0.99553958 0.00041336203    8479.7267 
    1900    6.3532521    301.40882    208.28331   -11119.481   0.99534534 0.00032474734    6662.3144 
    2000    6.6938104     298.0462   -236.47954   -11162.212   0.99421846 0.00023869721    4898.8129 
 Loop time of 6.69387 on 8 procs for 2000 steps with 3000 atoms
 Performance: 25.815 ns/day, 0.930 hours/ns, 298.781 timesteps/s
 99.7% CPU use with 8 MPI tasks x 1 OpenMP threads
 MPI task timing breakdown:
 Section |  min time  |  avg time  |  max time  |%varavg| %total
 ---------------------------------------------------------------
 Pair    | 4.36       | 4.4981     | 4.6781     |   5.1 | 67.20
 Bond    | 0.00074545 | 0.00088463 | 0.0012464  |   0.0 |  0.01
 Kspace  | 0.86696    | 1.0476     | 1.1863     |  10.5 | 15.65
 Neigh   | 0.37733    | 0.37785    | 0.3784     |   0.1 |  5.64
 Comm    | 0.19874    | 0.20168    | 0.20729    |   0.6 |  3.01
 Output  | 0.0015529  | 0.0015803  | 0.0017546  |   0.2 |  0.02
 Modify  | 0.54083    | 0.55143    | 0.55445    |   0.6 |  8.24
 Other   |            | 0.01483    |            |       |  0.22
 Nlocal:    375 ave 385 max 361 min
 Histogram: 1 1 0 0 1 0 2 0 1 2
 Nghost:    5772.25 ave 5789 max 5757 min
 Histogram: 1 1 2 0 0 0 2 0 0 2
 Neighs:    135285 ave 144189 max 127550 min
 Histogram: 1 2 1 1 0 0 0 0 1 2
 Total # of neighbors = 1082280
 Ave neighs/atom = 360.76
 Ave special neighs/atom = 2
 Neighbor list builds = 101
 Dangerous builds = 1
 Total wall time: 0:00:06
--- a/lib/atc/Function.cpp
+++ b/lib/atc/Function.cpp
@ -1,5 +1,11 @@
-#ifndef _WIN32
+#if !defined(_WIN32)
 #if defined(__FreeBSD__) || defined(__DragonFly__) || defined(__OpenBSD__) || defined(__NetBSD__)
 #include <cstdlib>
 #else
 #include <alloca.h>
 #endif
 #endif
 #include "Function.h"
 #include "ATC_Error.h"
--- a/lib/voronoi/voro-make.patch
+++ b/lib/voronoi/voro-make.patch
@ -0,0 +1,30 @@
 --- src/Makefile.orig	2020-05-03 03:50:23.501557199 -0400
 +++ src/Makefile	2020-05-03 03:53:32.147681674 -0400
@@ -10,10 +10,10 @@
 # List of the common source files
 objs=cell.o common.o container.o unitcell.o v_compute.o c_loops.o \
      v_base.o wall.o pre_container.o container_prd.o
 -src=$(patsubst %.o,%.cc,$(objs))
 +src=$(objs:.o=.cc)
 # Makefile rules
 -all: libvoro++.a voro++
 +all: depend libvoro++.a voro++
 depend:
 	$(CXX) -MM $(src) >Makefile.dep
@@ -22,12 +22,12 @@
 libvoro++.a: $(objs)
 	rm -f libvoro++.a
 -	ar rs libvoro++.a $^
 +	$(AR) rs libvoro++.a $(objs)
 voro++: libvoro++.a cmd_line.cc
 	$(CXX) $(CFLAGS) -L. -o voro++ cmd_line.cc -lvoro++
 -%.o: %.cc
 +.cc.o:
 	$(CXX) $(CFLAGS) -c $<
 help: Doxyfile $(SOURCE)
--- a/src/KIM/pair_kim.cpp
+++ b/src/KIM/pair_kim.cpp
@ -926,7 +926,7 @@ void PairKIM::set_argument_pointers()
        kimerror = kimerror ||
        KIM_ComputeArguments_SetArgumentPointerDouble(
          pargs,KIM_COMPUTE_ARGUMENT_NAME_partialEnergy,
-          reinterpret_cast<double * const>(NULL));
+          static_cast<double * const>(NULL));
      }
  }
@ -947,7 +947,7 @@ void PairKIM::set_argument_pointers()
    kimerror = kimerror || KIM_ComputeArguments_SetArgumentPointerDouble(
      pargs,
      KIM_COMPUTE_ARGUMENT_NAME_partialParticleEnergy,
-      reinterpret_cast<double * const>(NULL));
+      static_cast<double * const>(NULL));
  } else if (KIM_SupportStatus_NotEqual(kim_model_support_for_particleEnergy,
                                      KIM_SUPPORT_STATUS_notSupported)) {
    kimerror = kimerror || KIM_ComputeArguments_SetArgumentPointerDouble(
@ -960,7 +960,7 @@ void PairKIM::set_argument_pointers()
    kimerror = kimerror || KIM_ComputeArguments_SetArgumentPointerDouble(
      pargs,
      KIM_COMPUTE_ARGUMENT_NAME_partialForces,
-      reinterpret_cast<double * const>(NULL));
+      static_cast<double * const>(NULL));
  } else {
    kimerror = kimerror || KIM_ComputeArguments_SetArgumentPointerDouble(
        pargs, KIM_COMPUTE_ARGUMENT_NAME_partialForces, &(atom->f[0][0]));
@ -983,7 +983,7 @@ void PairKIM::set_argument_pointers()
    kimerror = kimerror || KIM_ComputeArguments_SetArgumentPointerDouble(
      pargs,
      KIM_COMPUTE_ARGUMENT_NAME_partialParticleVirial,
-      reinterpret_cast<double * const>(NULL));
+      static_cast<double * const>(NULL));
  } else if (KIM_SupportStatus_NotEqual(kim_model_support_for_particleVirial,
                                        KIM_SUPPORT_STATUS_notSupported)) {
    kimerror = kimerror || KIM_ComputeArguments_SetArgumentPointerDouble(
--- a/src/KOKKOS/Install.sh
+++ b/src/KOKKOS/Install.sh
@ -85,6 +85,8 @@ action comm_kokkos.cpp
 action comm_kokkos.h
 action comm_tiled_kokkos.cpp
 action comm_tiled_kokkos.h
 action compute_coord_atom_kokkos.cpp
 action compute_coord_atom_kokkos.h
 action compute_orientorder_atom_kokkos.cpp
 action compute_orientorder_atom_kokkos.h
 action compute_temp_kokkos.cpp
--- a/src/KOKKOS/compute_coord_atom_kokkos.cpp
+++ b/src/KOKKOS/compute_coord_atom_kokkos.cpp
@ -0,0 +1,249 @@
 /* ----------------------------------------------------------------------
   LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
   http://lammps.sandia.gov, 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_coord_atom_kokkos.h"
 #include <cmath>
 #include <cstring>
 #include "compute_orientorder_atom_kokkos.h"
 #include "atom_kokkos.h"
 #include "update.h"
 #include "modify.h"
 #include "neighbor_kokkos.h"
 #include "neigh_list.h"
 #include "neigh_request.h"
 #include "force.h"
 #include "pair.h"
 #include "comm.h"
 #include "group.h"
 #include "memory_kokkos.h"
 #include "error.h"
 #include "atom_masks.h"
 using namespace LAMMPS_NS;
 #define INVOKED_PERATOM 8
 /* ---------------------------------------------------------------------- */
 template<class DeviceType>
 ComputeCoordAtomKokkos<DeviceType>::ComputeCoordAtomKokkos(LAMMPS *lmp, int narg, char **arg) :
  ComputeCoordAtom(lmp, narg, arg)
 {
  kokkosable = 1;
  atomKK = (AtomKokkos *) atom;
  execution_space = ExecutionSpaceFromDevice<DeviceType>::space;
  datamask_read = EMPTY_MASK;
  datamask_modify = EMPTY_MASK;
  d_typelo = typename AT::t_int_1d("coord/atom:typelo",ncol);
  d_typehi = typename AT::t_int_1d("coord/atom:typehi",ncol);
  auto h_typelo = Kokkos::create_mirror_view(d_typelo);
  auto h_typehi = Kokkos::create_mirror_view(d_typehi);
  for (int i = 0; i < ncol; i++) {
    h_typelo(i) = typelo[i];
    h_typehi(i) = typehi[i];
  }
  Kokkos::deep_copy(d_typelo,h_typelo);
  Kokkos::deep_copy(d_typehi,h_typehi);
 }
 /* ---------------------------------------------------------------------- */
 template<class DeviceType>
 ComputeCoordAtomKokkos<DeviceType>::~ComputeCoordAtomKokkos<DeviceType>()
 {
  if (copymode) return;
  memoryKK->destroy_kokkos(k_cvec,cvec);
  memoryKK->destroy_kokkos(k_carray,carray);
 }
 /* ---------------------------------------------------------------------- */
 template<class DeviceType>
 void ComputeCoordAtomKokkos<DeviceType>::init()
 {
  ComputeCoordAtom::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 ComputeCoordAtomKokkos<DeviceType>::compute_peratom()
 {
  invoked_peratom = update->ntimestep;
  // grow coordination array if necessary
  if (atom->nmax > nmax) {
    if (ncol == 1) {
      memoryKK->destroy_kokkos(k_cvec,cvec);
      nmax = atom->nmax;
      memoryKK->create_kokkos(k_cvec,cvec,nmax,"coord/atom:cvec");
      vector_atom = cvec;
      d_cvec = k_cvec.template view<DeviceType>();
    } else {
      memoryKK->destroy_kokkos(k_carray,carray);
      nmax = atom->nmax;
      memoryKK->create_kokkos(k_carray,carray,nmax,ncol,"coord/atom:carray");
      array_atom = carray;
      d_carray = k_carray.template view<DeviceType>();
    }
  }
  if (cstyle == ORIENT) {
    if (!(c_orientorder->invoked_flag & INVOKED_PERATOM)) {
      c_orientorder->compute_peratom();
      c_orientorder->invoked_flag |= INVOKED_PERATOM;
    }
    nqlist = c_orientorder->nqlist;
    normv = c_orientorder->array_atom;
    comm->forward_comm_compute(this);
    if (!c_orientorder->kokkosable)
      error->all(FLERR,"Must use compute orientorder/atom/kk with compute coord/atom/kk");
    if (c_orientorder->execution_space == Host) {
      ComputeOrientOrderAtomKokkos<LMPHostType>* c_orientorder_kk;
      c_orientorder_kk = (ComputeOrientOrderAtomKokkos<LMPHostType>*) c_orientorder;
      c_orientorder_kk->k_qnarray.modify<LMPHostType>();
      c_orientorder_kk->k_qnarray.sync<DeviceType>();
      d_normv = c_orientorder_kk->k_qnarray.view<DeviceType>();
    } else {
      ComputeOrientOrderAtomKokkos<LMPDeviceType>* c_orientorder_kk;
      c_orientorder_kk = (ComputeOrientOrderAtomKokkos<LMPDeviceType>*) c_orientorder;
      c_orientorder_kk->k_qnarray.modify<LMPHostType>();
      c_orientorder_kk->k_qnarray.sync<DeviceType>();
      d_normv = c_orientorder_kk->k_qnarray.view<DeviceType>();
    }
  }
  // invoke full neighbor list (will copy or build if necessary)
  neighbor->build_one(list);
  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 coordination number(s) for each atom in group
  // use full neighbor list to count atoms less than cutoff
  atomKK->sync(execution_space,X_MASK|TYPE_MASK|MASK_MASK);
  x = atomKK->k_x.view<DeviceType>();
  type = atomKK->k_type.view<DeviceType>();
  mask = atomKK->k_mask.view<DeviceType>();
  copymode = 1;
  if (cstyle == CUTOFF) {
    if (ncol == 1) {
      typename Kokkos::RangePolicy<DeviceType, TagComputeCoordAtom<CUTOFF,1> > policy(0,inum);
      Kokkos::parallel_for("ComputeCoordAtom",policy,*this);
    } else {
      typename Kokkos::RangePolicy<DeviceType, TagComputeCoordAtom<CUTOFF,0> > policy(0,inum);
      Kokkos::parallel_for("ComputeCoordAtom",policy,*this);
    }
  } else if (cstyle == ORIENT) {
    typename Kokkos::RangePolicy<DeviceType, TagComputeCoordAtom<ORIENT,1> > policy(0,inum);
    Kokkos::parallel_for("ComputeCoordAtom",policy,*this);
  }
  copymode = 0;
  if (ncol == 1 || cstyle == ORIENT) {
    k_cvec.modify<DeviceType>();
    k_cvec.sync<LMPHostType>();
  } else {
    k_carray.modify<DeviceType>();
    k_carray.sync<LMPHostType>();
  }
 }
 template<class DeviceType>
 template<int CSTYLE, int NCOL>
 KOKKOS_INLINE_FUNCTION
 void ComputeCoordAtomKokkos<DeviceType>::operator()(TagComputeCoordAtom<CSTYLE,NCOL>, const int &ii) const
 {
  const int i = d_ilist[ii];
  if (NCOL == 1 || CSTYLE == ORIENT)
    d_cvec(i) = 0.0;
  else
    for (int m = 0; m < ncol; m++) d_carray(i,m) = 0.0;
  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];
    int n = 0;
    for (int jj = 0; jj < jnum; jj++) {
      int j = d_neighbors(i,jj);
      j &= NEIGHMASK;
      if (NCOL == 1)
        if (!(mask[j] & jgroupbit)) continue;
      const int jtype = type[j];
      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) {
        if (CSTYLE == CUTOFF) {
          if (NCOL == 1) {
            if (jtype >= d_typelo[0] && jtype <= d_typehi[0])
              n++;
          } else {
              for (int m = 0; m < ncol; m++)
                if (jtype >= d_typelo[m] && jtype <= d_typehi[m])
                  d_carray(i,m) += 1.0;
          }
        } else if (CSTYLE == ORIENT) {
            double dot_product = 0.0;
            for (int m=0; m < 2*(2*l+1); m++) {
              dot_product += d_normv(i,nqlist+m)*d_normv(j,nqlist+m);
            }
            if (dot_product > threshold) n++;
        }
      }
    }
    if (NCOL == 1 || CSTYLE == ORIENT)
      d_cvec[i] = n;
  }
 }
 namespace LAMMPS_NS {
 template class ComputeCoordAtomKokkos<LMPDeviceType>;
 #ifdef KOKKOS_ENABLE_CUDA
 template class ComputeCoordAtomKokkos<LMPHostType>;
 #endif
 }
--- a/src/KOKKOS/compute_coord_atom_kokkos.h
+++ b/src/KOKKOS/compute_coord_atom_kokkos.h
@ -0,0 +1,78 @@
 /* -*- c++ -*- ----------------------------------------------------------
   LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
   http://lammps.sandia.gov, 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(coord/atom/kk,ComputeCoordAtomKokkos<LMPDeviceType>)
 ComputeStyle(coord/atom/kk/device,ComputeCoordAtomKokkos<LMPDeviceType>)
 ComputeStyle(coord/atom/kk/host,ComputeCoordAtomKokkos<LMPHostType>)
 #else
 #ifndef LMP_COMPUTE_COORD_ATOM_KOKKOS_H
 #define LMP_COMPUTE_COORD_ATOM_KOKKOS_H
 #include "compute_coord_atom.h"
 #include "kokkos_type.h"
 namespace LAMMPS_NS {
 template<int CSTYLE, int NCOL>
 struct TagComputeCoordAtom{};
 template<class DeviceType>
 class ComputeCoordAtomKokkos : public ComputeCoordAtom {
 public:
  typedef DeviceType device_type;
  typedef ArrayTypes<DeviceType> AT;
  ComputeCoordAtomKokkos(class LAMMPS *, int, char **);
  virtual ~ComputeCoordAtomKokkos();
  void init();
  void compute_peratom();
  enum {NONE,CUTOFF,ORIENT};
  template<int CSTYLE, int NCOL>
  KOKKOS_INLINE_FUNCTION
  void operator()(TagComputeCoordAtom<CSTYLE,NCOL>, const int&) const;
 private:
  int inum;
  typename AT::t_x_array_randomread x;
  typename ArrayTypes<DeviceType>::t_int_1d_randomread type;
  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;
  typename AT::t_int_1d d_typelo;
  typename AT::t_int_1d d_typehi;
  DAT::tdual_float_1d k_cvec;
  typename AT::t_float_1d d_cvec;
  DAT::tdual_float_2d k_carray;
  typename AT::t_float_2d d_carray;
  typename AT::t_float_2d d_normv;
 };
 }
 #endif
 #endif
 /* ERROR/WARNING messages:
 */
--- a/src/KOKKOS/fft3d_kokkos.cpp
+++ b/src/KOKKOS/fft3d_kokkos.cpp
@ -38,7 +38,8 @@ FFT3dKokkos<DeviceType>::FFT3dKokkos(LAMMPS *lmp, MPI_Comm comm, int nfast, int
             int in_klo, int in_khi,
             int out_ilo, int out_ihi, int out_jlo, int out_jhi,
             int out_klo, int out_khi,
-             int scaled, int permute, int *nbuf, int usecollective) :
+             int scaled, int permute, int *nbuf, int usecollective,
             int usecuda_aware) :
  Pointers(lmp)
 {
  int nthreads = lmp->kokkos->nthreads;
@ -70,7 +71,7 @@ FFT3dKokkos<DeviceType>::FFT3dKokkos(LAMMPS *lmp, MPI_Comm comm, int nfast, int
  plan = fft_3d_create_plan_kokkos(comm,nfast,nmid,nslow,
                            in_ilo,in_ihi,in_jlo,in_jhi,in_klo,in_khi,
                            out_ilo,out_ihi,out_jlo,out_jhi,out_klo,out_khi,
-                            scaled,permute,nbuf,usecollective,nthreads);
+                            scaled,permute,nbuf,usecollective,nthreads,usecuda_aware);
  if (plan == NULL) error->one(FLERR,"Could not create 3d FFT plan");
 }
@ -368,6 +369,7 @@ void FFT3dKokkos<DeviceType>::fft_3d_kokkos(typename FFT_AT::t_FFT_DATA_1d d_in,
                          2 = permute twice = slow->fast, fast->mid, mid->slow
   nbuf                 returns size of internal storage buffers used by FFT
   usecollective        use collective MPI operations for remapping data
   usecuda_aware        use CUDA-Aware MPI or not
 ------------------------------------------------------------------------- */
 template<class DeviceType>
@ -378,7 +380,7 @@ struct fft_plan_3d_kokkos<DeviceType>* FFT3dKokkos<DeviceType>::fft_3d_create_pl
       int out_ilo, int out_ihi, int out_jlo, int out_jhi,
       int out_klo, int out_khi,
       int scaled, int permute, int *nbuf, int usecollective,
-       int nthreads)
+       int nthreads, int usecuda_aware)
 {
  struct fft_plan_3d_kokkos<DeviceType> *plan;
  int me,nprocs;
@ -435,7 +437,8 @@ struct fft_plan_3d_kokkos<DeviceType>* FFT3dKokkos<DeviceType>::fft_3d_create_pl
    plan->pre_plan =
      remapKK->remap_3d_create_plan_kokkos(comm,in_ilo,in_ihi,in_jlo,in_jhi,in_klo,in_khi,
                           first_ilo,first_ihi,first_jlo,first_jhi,
-                           first_klo,first_khi,2,0,0,FFT_PRECISION,0);
+                           first_klo,first_khi,2,0,0,FFT_PRECISION,
                           usecollective,usecuda_aware);
    if (plan->pre_plan == NULL) return NULL;
  }
@ -460,7 +463,7 @@ struct fft_plan_3d_kokkos<DeviceType>* FFT3dKokkos<DeviceType>::fft_3d_create_pl
                           first_klo,first_khi,
                           second_ilo,second_ihi,second_jlo,second_jhi,
                           second_klo,second_khi,2,1,0,FFT_PRECISION,
-                           usecollective);
+                           usecollective,usecuda_aware);
  if (plan->mid1_plan == NULL) return NULL;
  // 1d FFTs along mid axis
@ -500,7 +503,8 @@ struct fft_plan_3d_kokkos<DeviceType>* FFT3dKokkos<DeviceType>::fft_3d_create_pl
                         second_jlo,second_jhi,second_klo,second_khi,
                         second_ilo,second_ihi,
                         third_jlo,third_jhi,third_klo,third_khi,
-                         third_ilo,third_ihi,2,1,0,FFT_PRECISION,usecollective);
+                         third_ilo,third_ihi,2,1,0,FFT_PRECISION,
                         usecollective,usecuda_aware);
  if (plan->mid2_plan == NULL) return NULL;
  // 1d FFTs along slow axis
@ -527,7 +531,8 @@ struct fft_plan_3d_kokkos<DeviceType>* FFT3dKokkos<DeviceType>::fft_3d_create_pl
                           third_klo,third_khi,third_ilo,third_ihi,
                           third_jlo,third_jhi,
                           out_klo,out_khi,out_ilo,out_ihi,
-                           out_jlo,out_jhi,2,(permute+1)%3,0,FFT_PRECISION,0);
+                           out_jlo,out_jhi,2,(permute+1)%3,0,FFT_PRECISION,
                           usecollective,usecuda_aware);
    if (plan->post_plan == NULL) return NULL;
  }
--- a/src/KOKKOS/fft3d_kokkos.h
+++ b/src/KOKKOS/fft3d_kokkos.h
@ -77,7 +77,7 @@ class FFT3dKokkos : protected Pointers {
  FFT3dKokkos(class LAMMPS *, MPI_Comm,
        int,int,int,int,int,int,int,int,int,int,int,int,int,int,int,
-        int,int,int *,int);
+        int,int,int *,int,int);
  ~FFT3dKokkos();
  void compute(typename FFT_AT::t_FFT_SCALAR_1d, typename FFT_AT::t_FFT_SCALAR_1d, int);
  void timing1d(typename FFT_AT::t_FFT_SCALAR_1d, int, int);
@ -95,7 +95,7 @@ class FFT3dKokkos : protected Pointers {
  struct fft_plan_3d_kokkos<DeviceType> *fft_3d_create_plan_kokkos(MPI_Comm, int, int, int,
                                         int, int, int, int, int,
                                         int, int, int, int, int, int, int,
-                                         int, int, int *, int, int);
+                                         int, int, int *, int, int, int);
  void fft_3d_destroy_plan_kokkos(struct fft_plan_3d_kokkos<DeviceType> *);
--- a/src/KOKKOS/npair_kokkos.cpp
+++ b/src/KOKKOS/npair_kokkos.cpp
@ -289,7 +289,10 @@ void NPairKokkos<DeviceType,HALF_NEIGH,GHOST,TRI,SIZE>::build(NeighList *list_)
            if (team_size <= team_size_max) {
              Kokkos::TeamPolicy<DeviceType> config((mbins+factor-1)/factor,team_size);
              Kokkos::parallel_for(config, f);
-            } else Kokkos::parallel_for(nall, f); // fall back to flat method
+            } else { // fall back to flat method
              f.sharedsize = 0;
              Kokkos::parallel_for(nall, f);
            }
          } else
            Kokkos::parallel_for(nall, f);
 #else
--- a/src/KOKKOS/pppm_kokkos.cpp
+++ b/src/KOKKOS/pppm_kokkos.cpp
@ -840,21 +840,23 @@ void PPPMKokkos<DeviceType>::allocate()
  // 2nd FFT returns data in 3d brick decomposition
  // remap takes data from 3d brick to FFT decomposition
  int collective_flag = 0; // not yet supported in Kokkos version
  int cuda_aware_flag = lmp->kokkos->cuda_aware_flag;
  int tmp;
  fft1 = new FFT3dKokkos<DeviceType>(lmp,world,nx_pppm,ny_pppm,nz_pppm,
                         nxlo_fft,nxhi_fft,nylo_fft,nyhi_fft,nzlo_fft,nzhi_fft,
                         nxlo_fft,nxhi_fft,nylo_fft,nyhi_fft,nzlo_fft,nzhi_fft,
-                         0,0,&tmp,collective_flag);
+                         0,0,&tmp,collective_flag,cuda_aware_flag);
  fft2 = new FFT3dKokkos<DeviceType>(lmp,world,nx_pppm,ny_pppm,nz_pppm,
                         nxlo_fft,nxhi_fft,nylo_fft,nyhi_fft,nzlo_fft,nzhi_fft,
                         nxlo_in,nxhi_in,nylo_in,nyhi_in,nzlo_in,nzhi_in,
-                         0,0,&tmp,collective_flag);
+                         0,0,&tmp,collective_flag,cuda_aware_flag);
  remap = new RemapKokkos<DeviceType>(lmp,world,
                          nxlo_in,nxhi_in,nylo_in,nyhi_in,nzlo_in,nzhi_in,
                          nxlo_fft,nxhi_fft,nylo_fft,nyhi_fft,nzlo_fft,nzhi_fft,
-                          1,0,0,FFT_PRECISION,collective_flag);
+                          1,0,0,FFT_PRECISION,collective_flag,cuda_aware_flag);
  // create ghost grid object for rho and electric field communication
--- a/src/KOKKOS/remap_kokkos.cpp
+++ b/src/KOKKOS/remap_kokkos.cpp
@ -38,12 +38,14 @@ RemapKokkos<DeviceType>::RemapKokkos(LAMMPS *lmp, MPI_Comm comm,
             int out_ilo, int out_ihi, int out_jlo, int out_jhi,
             int out_klo, int out_khi,
             int nqty, int permute, int memory,
-             int precision, int usecollective) : Pointers(lmp)
+             int precision, int usecollective,
             int usecuda_aware) : Pointers(lmp)
 {
  plan = remap_3d_create_plan_kokkos(comm,
                              in_ilo,in_ihi,in_jlo,in_jhi,in_klo,in_khi,
                              out_ilo,out_ihi,out_jlo,out_jhi,out_klo,out_khi,
-                              nqty,permute,memory,precision,usecollective);
+                              nqty,permute,memory,precision,usecollective,
                              usecuda_aware);
  if (plan == NULL) error->one(FLERR,"Could not create 3d remap plan");
 }
@ -119,20 +121,36 @@ void RemapKokkos<DeviceType>::remap_3d_kokkos(typename FFT_AT::t_FFT_SCALAR_1d d
  // post all recvs into scratch space
  FFT_SCALAR* v_scratch = d_scratch.data();
  if (!plan->usecuda_aware) {
    plan->h_scratch = Kokkos::create_mirror_view(d_scratch);
    v_scratch = plan->h_scratch.data();
  }
  for (irecv = 0; irecv < plan->nrecv; irecv++) {
-    FFT_SCALAR* scratch = d_scratch.data() + plan->recv_bufloc[irecv];
+    FFT_SCALAR* scratch = v_scratch + plan->recv_bufloc[irecv];
    MPI_Irecv(scratch,plan->recv_size[irecv],
              MPI_FFT_SCALAR,plan->recv_proc[irecv],0,
              plan->comm,&plan->request[irecv]);
  }
  FFT_SCALAR* v_sendbuf = plan->d_sendbuf.data();
  if (!plan->usecuda_aware) {
    plan->h_sendbuf = Kokkos::create_mirror_view(plan->d_sendbuf);
    v_sendbuf = plan->h_sendbuf.data();
  }
  // send all messages to other procs
  for (isend = 0; isend < plan->nsend; isend++) {
    int in_offset = plan->send_offset[isend];
    plan->pack(d_in,in_offset,
               plan->d_sendbuf,0,&plan->packplan[isend]);
-    MPI_Send(plan->d_sendbuf.data(),plan->send_size[isend],MPI_FFT_SCALAR,
+
    if (!plan->usecuda_aware)
      Kokkos::deep_copy(plan->h_sendbuf,plan->d_sendbuf);
    MPI_Send(v_sendbuf,plan->send_size[isend],MPI_FFT_SCALAR,
             plan->send_proc[isend],0,plan->comm);
  }
@ -161,6 +179,9 @@ void RemapKokkos<DeviceType>::remap_3d_kokkos(typename FFT_AT::t_FFT_SCALAR_1d d
    int scratch_offset = plan->recv_bufloc[irecv];
    int out_offset = plan->recv_offset[irecv];
    if (!plan->usecuda_aware)
      Kokkos::deep_copy(d_scratch,plan->h_scratch);
    plan->unpack(d_scratch,scratch_offset,
                 d_out,out_offset,&plan->unpackplan[irecv]);
  }
@ -189,6 +210,7 @@ void RemapKokkos<DeviceType>::remap_3d_kokkos(typename FFT_AT::t_FFT_SCALAR_1d d
                          1 = single precision (4 bytes per datum)
                          2 = double precision (8 bytes per datum)
   usecollective        whether to use collective MPI or point-to-point
   usecuda_aware        whether to use CUDA-Aware MPI or not
 ------------------------------------------------------------------------- */
 template<class DeviceType>
@ -198,7 +220,8 @@ struct remap_plan_3d_kokkos<DeviceType>* RemapKokkos<DeviceType>::remap_3d_creat
  int in_klo, int in_khi,
  int out_ilo, int out_ihi, int out_jlo, int out_jhi,
  int out_klo, int out_khi,
-  int nqty, int permute, int memory, int precision, int usecollective)
+  int nqty, int permute, int memory, int precision,
  int usecollective, int usecuda_aware)
 {
  struct remap_plan_3d_kokkos<DeviceType> *plan;
@ -216,6 +239,7 @@ struct remap_plan_3d_kokkos<DeviceType>* RemapKokkos<DeviceType>::remap_3d_creat
  plan = new struct remap_plan_3d_kokkos<DeviceType>;
  if (plan == NULL) return NULL;
  plan->usecollective = usecollective;
  plan->usecuda_aware = usecuda_aware;
  // store parameters in local data structs
--- a/src/KOKKOS/remap_kokkos.h
+++ b/src/KOKKOS/remap_kokkos.h
@ -28,7 +28,9 @@ struct remap_plan_3d_kokkos {
  typedef DeviceType device_type;
  typedef FFTArrayTypes<DeviceType> FFT_AT;
  typename FFT_AT::t_FFT_SCALAR_1d d_sendbuf;                  // buffer for MPI sends
  FFT_HAT::t_FFT_SCALAR_1d h_sendbuf;                          // host buffer for MPI sends
  typename FFT_AT::t_FFT_SCALAR_1d d_scratch;                  // scratch buffer for MPI recvs
  FFT_HAT::t_FFT_SCALAR_1d h_scratch;                          // host scratch buffer for MPI recvs
  void (*pack)(typename FFT_AT::t_FFT_SCALAR_1d_um, int, typename FFT_AT::t_FFT_SCALAR_1d_um, int, struct pack_plan_3d *);
                                    // which pack function to use
  void (*unpack)(typename FFT_AT::t_FFT_SCALAR_1d_um, int, typename FFT_AT::t_FFT_SCALAR_1d_um, int, struct pack_plan_3d *);
@ -51,6 +53,7 @@ struct remap_plan_3d_kokkos {
  int usecollective;                // use collective or point-to-point MPI
  int commringlen;                  // length of commringlist
  int *commringlist;                // ranks on communication ring of this plan
  int usecuda_aware;                // use CUDA-Aware MPI or not
 };
 template<class DeviceType>
@ -60,7 +63,7 @@ class RemapKokkos : protected Pointers {
  typedef FFTArrayTypes<DeviceType> FFT_AT;
  RemapKokkos(class LAMMPS *);
  RemapKokkos(class LAMMPS *, MPI_Comm,int,int,int,int,int,int,
-        int,int,int,int,int,int,int,int,int,int,int);
+        int,int,int,int,int,int,int,int,int,int,int,int);
  ~RemapKokkos();
  void perform(typename FFT_AT::t_FFT_SCALAR_1d, typename FFT_AT::t_FFT_SCALAR_1d, typename FFT_AT::t_FFT_SCALAR_1d);
@ -70,7 +73,7 @@ class RemapKokkos : protected Pointers {
  struct remap_plan_3d_kokkos<DeviceType> *remap_3d_create_plan_kokkos(MPI_Comm,
                                             int, int, int, int, int, int,
                                             int, int, int, int, int, int,
-                                             int, int, int, int, int);
+                                             int, int, int, int, int, int);
  void remap_3d_destroy_plan_kokkos(struct remap_plan_3d_kokkos<DeviceType> *);
 };
--- a/src/MISC/xdr_compat.h
+++ b/src/MISC/xdr_compat.h
@ -59,7 +59,7 @@ extern "C" {
 typedef int bool_t;
-#if defined(__MINGW32__) || defined(__APPLE__)
+#if defined(__MINGW32__) || defined(__APPLE__) || defined(__FreeBSD__) || defined(__DragonFly__) || defined(__OpenBSD__) || defined(__NetBSD__)
 typedef char * caddr_t;
 typedef unsigned int u_int;
 #endif
--- a/src/USER-INTEL/verlet_lrt_intel.cpp
+++ b/src/USER-INTEL/verlet_lrt_intel.cpp
@ -102,7 +102,7 @@ void VerletLRTIntel::setup(int flag)
  }
  #if defined(_LMP_INTEL_LRT_PTHREAD)
-  #if defined(__linux)
+  #if defined(__linux__)
  if (comm->me == 0) {
    cpu_set_t cpuset;
    sched_getaffinity(0, sizeof(cpuset), &cpuset);
--- a/src/USER-MEAMC/meam_setup_global.cpp
+++ b/src/USER-MEAMC/meam_setup_global.cpp
@ -67,7 +67,14 @@ MEAM::meam_setup_global(int nelt, lattice_t* lat, int* ielement, double* /*atwt*
      case DIA3:
        this->re_meam[i][i] = tmplat[i] * sqrt(3.0) / 4.0;
        break;
-      //default:
+      case B1:
      case B2:
      case C11:
      case L12:
        // do nothing
        break;
      default:
        ;
      //  error
    }
  }
--- a/src/USER-MISC/README
+++ b/src/USER-MISC/README
@ -38,6 +38,7 @@ compute momentum, Rupert Nash (University of Edinburgh), r.nash at epcc.ed.ac.uk
 compute stress/mop, Romain Vermorel (U Pau) & Laurent Joly (U Lyon), romain.vermorel at univ-pau.fr & ljoly.ulyon at gmail.com, 5 Sep 18
 compute stress/mop/profile, Romain Vermorel (U Pau) & Laurent Joly (U Lyon), romain.vermorel at univ-pau.fr & ljoly.ulyon at gmail.com, 5 Sep 18
 compute temp/rotate, Laurent Joly (U Lyon), ljoly.ulyon at gmail.com, 8 Aug 11
 compute viscosity/cos, Zheng Gong (ENS de Lyon), z.gong@outlook.com, 24 Apr 20
 compute PRESSURE/GREM, David Stelter, dstelter@bu.edu, 22 Nov 16
 dihedral_style cosine/shift/exp, Carsten Svaneborg, science at zqex.dk, 8 Aug 11
 dihedral_style fourier, Loukas Peristeras, loukas.peristeras at scienomics.com, 27 Oct 12
@ -46,6 +47,7 @@ dihedral_style quadratic, Loukas Peristeras, loukas.peristeras at scienomics.com
 dihedral_style spherical, Andrew Jewett, jewett.aij@gmail.com, 15 Jul 16
 dihedral_style table, Andrew Jewett, jewett.aij@gmail.com, 10 Jan 12
 dihedral_style table/cut, Mike Salerno, ksalerno@pha.jhu.edu, 11 May 18
 fix accelerate/cos, Zheng Gong (ENS de Lyon), z.gong@outlook.com, 24 Apr 20
 fix addtorque, Laurent Joly (U Lyon), ljoly.ulyon at gmail.com, 8 Aug 11
 fix ave/correlate/long, Jorge Ramirez (UPM Madrid), jorge.ramirez at upm.es, 21 Oct 2015
 fix electron/stopping, Konstantin Avchaciov, k.avchachov at gmail.com, 26 Feb 2019
--- a/src/USER-MISC/compute_viscosity_cos.cpp
+++ b/src/USER-MISC/compute_viscosity_cos.cpp
@ -0,0 +1,293 @@
 /* ----------------------------------------------------------------------
   LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
   http://lammps.sandia.gov, 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: Zheng GONG (ENS de Lyon, z.gong@outlook.com)
 ------------------------------------------------------------------------- */
 #include <mpi.h>
 #include "compute_viscosity_cos.h"
 #include "atom.h"
 #include "update.h"
 #include "force.h"
 #include "domain.h"
 #include "comm.h"
 #include "group.h"
 #include "error.h"
 #include "math_const.h"
 using namespace LAMMPS_NS;
 using namespace MathConst;
 /* ---------------------------------------------------------------------- */
 ComputeViscosityCos::ComputeViscosityCos(LAMMPS *lmp, int narg, char **arg) :
    Compute(lmp, narg, arg) {
  if (narg != 3) error->all(FLERR, "Illegal compute viscosity/cos command");
  scalar_flag = vector_flag = 1;
  size_vector = 7;
  extscalar = 0;
  extvector = -1;
  extlist = new int[7]{1,1,1,1,1,1,0};
  tempflag = 1;
  tempbias = 1;
  maxbias = 0;
  vbiasall = NULL;
  vector = new double[7];
 }
 /* ---------------------------------------------------------------------- */
 ComputeViscosityCos::~ComputeViscosityCos() {
  if (!copymode)
    delete[] vector;
 }
 /* ---------------------------------------------------------------------- */
 void ComputeViscosityCos::setup() {
  dynamic = 0;
  if (dynamic_user || group->dynamic[igroup]) dynamic = 1;
  dof_compute();
 }
 /* ---------------------------------------------------------------------- */
 void ComputeViscosityCos::dof_compute() {
  adjust_dof_fix();
  natoms_temp = group->count(igroup);
  dof = domain->dimension * natoms_temp;
  dof -= extra_dof + fix_dof;
  if (dof > 0.0) tfactor = force->mvv2e / (dof * force->boltz);
  else tfactor = 0.0;
 }
 /* ---------------------------------------------------------------------- */
 void ComputeViscosityCos::calc_V() {
  double zlo = domain->boxlo[2];
  double zhi = domain->boxhi[2];
  double **x = atom->x;
  double **v = atom->v;
  double *mass = atom->mass;
  double *rmass = atom->rmass;
  int *type = atom->type;
  int *mask = atom->mask;
  int nlocal = atom->nlocal;
  double massone;
  double V_m[2];
  double V_m_local[2] = {0, 0};
  for (int i = 0; i < nlocal; i++)
    if (mask[i] & groupbit) {
      if (rmass) massone = rmass[i];
      else massone = mass[type[i]];
      V_m_local[0] +=
          2 * massone * v[i][0] * cos(MY_2PI * (x[i][2] - zlo) / (zhi - zlo));
      V_m_local[1] += massone;
    }
  MPI_Allreduce(V_m_local, V_m, 2, MPI_DOUBLE, MPI_SUM, world);
  V = V_m[0] / V_m[1];
 }
 /* ---------------------------------------------------------------------- */
 double ComputeViscosityCos::compute_scalar() {
  invoked_scalar = update->ntimestep;
  double zlo = domain->boxlo[2];
  double zhi = domain->boxhi[2];
  double **x = atom->x;
  double **v = atom->v;
  double *mass = atom->mass;
  double *rmass = atom->rmass;
  int *type = atom->type;
  int *mask = atom->mask;
  int nlocal = atom->nlocal;
  double t = 0.0;
  double vx_acc;
  double massone;
  calc_V();
  for (int i = 0; i < nlocal; i++)
    if (mask[i] & groupbit) {
      if (rmass) massone = rmass[i];
      else massone = mass[type[i]];
      vx_acc = V * cos(MY_2PI * (x[i][2] - zlo) / (zhi - zlo));
      t += ((v[i][0] - vx_acc) * (v[i][0] - vx_acc) + v[i][1] * v[i][1] +
            v[i][2] * v[i][2]) * massone;
    }
  MPI_Allreduce(&t, &scalar, 1, MPI_DOUBLE, MPI_SUM, world);
  if (dynamic) dof_compute();
  if (dof < 0.0 && natoms_temp > 0.0)
    error->all(FLERR, "Temperature compute degrees of freedom < 0");
  scalar *= tfactor;
  return scalar;
 }
 /* ---------------------------------------------------------------------- */
 void ComputeViscosityCos::compute_vector() {
  int i;
  invoked_vector = update->ntimestep;
  double zlo = domain->boxlo[2];
  double zhi = domain->boxhi[2];
  double **x = atom->x;
  double **v = atom->v;
  double *mass = atom->mass;
  double *rmass = atom->rmass;
  int *type = atom->type;
  int *mask = atom->mask;
  int nlocal = atom->nlocal;
  double vx_acc;
  double massone, t[6];
  for (i = 0; i < 6; i++) t[i] = 0.0;
  for (i = 0; i < nlocal; i++)
    if (mask[i] & groupbit) {
      if (rmass) massone = rmass[i];
      else massone = mass[type[i]];
      vx_acc = V * cos(MY_2PI * (x[i][2] - zlo) / (zhi - zlo));
      t[0] += massone * (v[i][0] - vx_acc) * (v[i][0] - vx_acc);
      t[1] += massone * v[i][1] * v[i][1];
      t[2] += massone * v[i][2] * v[i][2];
      t[3] += massone * (v[i][0] - vx_acc) * v[i][1];
      t[4] += massone * (v[i][0] - vx_acc) * v[i][2];
      t[5] += massone * v[i][1] * v[i][2];
    }
  MPI_Allreduce(t, vector, 6, MPI_DOUBLE, MPI_SUM, world);
  for (i = 0; i < 6; i++) vector[i] *= force->mvv2e;
  vector[6] = V;
 }
 /* ----------------------------------------------------------------------
   remove velocity bias from atom I to leave thermal velocity
 ------------------------------------------------------------------------- */
 void ComputeViscosityCos::remove_bias(int i, double *v) {
  double zlo = domain->boxlo[2];
  double zhi = domain->boxhi[2];
  double **x = atom->x;
  vbias[0] = V * cos(MY_2PI * (x[i][2] - zlo) / (zhi - zlo));
  vbias[1] = 0;
  vbias[2] = 0;
  v[0] -= vbias[0];
 //  v[1] -= vbias[1];
 //  v[2] -= vbias[2];
 }
 /* ----------------------------------------------------------------------
   remove velocity bias from atom I to leave thermal velocity
 ------------------------------------------------------------------------- */
 void ComputeViscosityCos::remove_bias_thr(int i, double *v, double *b) {
  double zlo = domain->boxlo[2];
  double zhi = domain->boxhi[2];
  double **x = atom->x;
  b[0] = V * cos(MY_2PI * (x[i][2] - zlo) / (zhi - zlo));
  b[1] = 0;
  b[2] = 0;
  v[0] -= b[0];
 //  v[1] -= b[1];
 //  v[2] -= b[2];
 }
 /* ----------------------------------------------------------------------
   remove velocity bias from all atoms to leave thermal velocity
 ------------------------------------------------------------------------- */
 void ComputeViscosityCos::remove_bias_all() {
  double zlo = domain->boxlo[2];
  double zhi = domain->boxhi[2];
  double **x = atom->x;
  double **v = atom->v;
  int *mask = atom->mask;
  int nlocal = atom->nlocal;
  for (int i = 0; i < nlocal; i++)
    if (mask[i] & groupbit) {
      vbiasall[i][0] = V * cos(MY_2PI * (x[i][2] - zlo) / (zhi - zlo));
      vbiasall[i][1] = 0;
      vbiasall[i][2] = 0;
      v[i][0] -= vbiasall[i][0];
 //      v[i][1] -= vbiasall[i][1];
 //      v[i][2] -= vbiasall[i][2];
    }
 }
 /* ----------------------------------------------------------------------
   add back in velocity bias to atom I removed by remove_bias()
   assume remove_bias() was previously called
 ------------------------------------------------------------------------- */
 void ComputeViscosityCos::restore_bias(int i, double *v) {
  v[0] += vbias[0];
  v[1] += vbias[1];
  v[2] += vbias[2];
 }
 /* ----------------------------------------------------------------------
   add back in velocity bias to atom I removed by remove_bias_thr()
   assume remove_bias_thr() was previously called with the same buffer b
 ------------------------------------------------------------------------- */
 void ComputeViscosityCos::restore_bias_thr(int i, double *v, double *b) {
  v[0] += b[0];
  v[1] += b[1];
  v[2] += b[2];
 }
 /* ----------------------------------------------------------------------
   add back in velocity bias to all atoms removed by remove_bias_all()
   assume remove_bias_all() was previously called
 ------------------------------------------------------------------------- */
 void ComputeViscosityCos::restore_bias_all() {
  double **v = atom->v;
  int *mask = atom->mask;
  int nlocal = atom->nlocal;
  for (int i = 0; i < nlocal; i++)
    if (mask[i] & groupbit) {
      v[i][0] += vbiasall[i][0];
      v[i][1] += vbiasall[i][1];
      v[i][2] += vbiasall[i][2];
    }
 }
--- a/src/USER-MISC/compute_viscosity_cos.h
+++ b/src/USER-MISC/compute_viscosity_cos.h
@ -0,0 +1,73 @@
 /* -*- c++ -*- ----------------------------------------------------------
   LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
   http://lammps.sandia.gov, 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: Zheng GONG (ENS de Lyon, z.gong@outlook.com)
 ------------------------------------------------------------------------- */
 #ifdef COMPUTE_CLASS
 ComputeStyle(viscosity/cos,ComputeViscosityCos)
 #else
 #ifndef LMP_COMPUTE_VISCOSITY_COS_H
 #define LMP_COMPUTE_VISCOSITY_COS_H
 #include "compute.h"
 namespace LAMMPS_NS {
 class ComputeViscosityCos : public Compute {
 public:
  ComputeViscosityCos(class LAMMPS *, int, char **);
  virtual ~ComputeViscosityCos();
  void init() {}
  void setup();
  virtual double compute_scalar();
  virtual void compute_vector();
  void remove_bias(int, double *);
  void remove_bias_thr(int, double *, double *);
  void remove_bias_all();
  void restore_bias(int, double *);
  void restore_bias_thr(int, double *, double *);
  void restore_bias_all();
 protected:
  double tfactor;
  double V;
  void dof_compute();
  void calc_V();
 };
 }
 #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: Temperature compute degrees of freedom < 0
 This should not happen if you are calculating the temperature
 on a valid set of atoms.
 */
--- a/src/USER-MISC/fix_accelerate_cos.cpp
+++ b/src/USER-MISC/fix_accelerate_cos.cpp
@ -0,0 +1,85 @@
 /* ----------------------------------------------------------------------
   LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
   http://lammps.sandia.gov, 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: Zheng GONG (ENS de Lyon, z.gong@outlook.com)
 ------------------------------------------------------------------------- */
 #include "fix_accelerate_cos.h"
 #include "atom.h"
 #include "update.h"
 #include "modify.h"
 #include "domain.h"
 #include "region.h"
 #include "respa.h"
 #include "input.h"
 #include "variable.h"
 #include "memory.h"
 #include "error.h"
 #include "force.h"
 #include "math_const.h"
 using namespace LAMMPS_NS;
 using namespace FixConst;
 /* ---------------------------------------------------------------------- */
 FixAccelerateCos::FixAccelerateCos(LAMMPS *lmp, int narg, char **arg) :
    Fix(lmp, narg, arg) {
  if (narg < 4) error->all(FLERR, "Illegal fix accelerate/cos command");
  acceleration = force->numeric(FLERR, arg[3]);
  if (domain->dimension == 2)
    error->all(FLERR,"Fix accelerate/cos cannot be used with 2d systems");
 }
 /* ---------------------------------------------------------------------- */
 int FixAccelerateCos::setmask() {
  int mask = 0;
  mask |= POST_FORCE;
  return mask;
 }
 /* ---------------------------------------------------------------------- */
 void FixAccelerateCos::setup(int vflag) {
  post_force(vflag);
 }
 /* ---------------------------------------------------------------------- */
 void FixAccelerateCos::post_force(int vflag) {
  double **x = atom->x;
  double **f = atom->f;
  int *type = atom->type;
  double *mass = atom->mass;
  double *rmass = atom->rmass;
  int *mask = atom->mask;
  int nlocal = atom->nlocal;
  double massone, force_x, acc_x;
  double zlo = domain->boxlo[2];
  double zhi = domain->boxhi[2];
  for (int i = 0; i < nlocal; i++)
    if (mask[i] & groupbit) {
      if (rmass) massone = rmass[i];
      else massone = mass[type[i]];
      acc_x = acceleration *
              cos(MathConst::MY_2PI * (x[i][2] - zlo) / (zhi - zlo));
      force_x = acc_x * massone * force->mvv2e;
      f[i][0] += force_x;
    }
 }
--- a/src/USER-MISC/fix_accelerate_cos.h
+++ b/src/USER-MISC/fix_accelerate_cos.h
@ -0,0 +1,61 @@
 /* -*- c++ -*- ----------------------------------------------------------
   LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
   http://lammps.sandia.gov, 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: Zheng GONG (ENS de Lyon, z.gong@outlook.com)
 ------------------------------------------------------------------------- */
 #ifdef FIX_CLASS
 FixStyle(accelerate/cos,FixAccelerateCos)
 #else
 #ifndef LMP_FIX_ACCELERATE_COS_H
 #define LMP_FIX_ACCELERATE_COS_H
 #include "fix.h"
 namespace LAMMPS_NS {
 class FixAccelerateCos: public Fix {
 public:
  FixAccelerateCos(class LAMMPS *, int, char **);
  virtual ~FixAccelerateCos() {};
  int setmask();
  virtual void init() {};
  void setup(int);
  virtual void post_force(int);
 protected:
  double acceleration;
 };
 }
 #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: Fix accelerate/cos cannot be used with 2d systems
 Self-explanatory.
 */
--- a/src/USER-MISC/fix_ave_correlate_long.cpp
+++ b/src/USER-MISC/fix_ave_correlate_long.cpp
@ -495,7 +495,7 @@ void FixAveCorrelateLong::end_of_step()
    if(overwrite) fseek(fp,filepos,SEEK_SET);
    fprintf(fp,"# Timestep: " BIGINT_FORMAT "\n", ntimestep);
    for (unsigned int i=0;i<npcorr;++i) {
-      fprintf(fp, "%lg ", t[i]*update->dt);
+      fprintf(fp, "%lg ", t[i]*update->dt*nevery);
      for (int j=0;j<npair;++j) {
        fprintf(fp, "%lg ", f[j][i]);
      }
--- a/src/USER-MISC/fix_imd.cpp
+++ b/src/USER-MISC/fix_imd.cpp
@ -1237,7 +1237,7 @@ void *imdsock_accept(void * v) {
 #elif defined(SOCKLEN_T)
  SOCKLEN_T len;
 #define _SOCKLEN_TYPE SOCKLEN_T
-#elif defined(_POSIX_SOURCE) || (defined(__APPLE__) && defined(__MACH__)) || defined(__linux)
+#elif defined(_POSIX_SOURCE) || (defined(__APPLE__) && defined(__MACH__)) || defined(__linux__) || defined(__FreeBSD__) || defined(__DragonFly__) || defined(__OpenBSD__) || defined(__NetBSD__)
  socklen_t len;
 #define _SOCKLEN_TYPE socklen_t
 #else
--- a/src/USER-MISC/pair_meam_spline.cpp
+++ b/src/USER-MISC/pair_meam_spline.cpp
@ -717,6 +717,7 @@ void PairMEAMSpline::SplineFunction::prepareSpline(Error* error)
    Y2[i] /= h*6.0;
 #endif
  }
  inv_h = (1/h);
  xmax_shifted = xmax - xmin;
 }
@ -732,6 +733,7 @@ void PairMEAMSpline::SplineFunction::communicate(MPI_Comm& world, int me)
  MPI_Bcast(&isGridSpline, 1, MPI_INT, 0, world);
  MPI_Bcast(&h, 1, MPI_DOUBLE, 0, world);
  MPI_Bcast(&hsq, 1, MPI_DOUBLE, 0, world);
  MPI_Bcast(&inv_h, 1, MPI_DOUBLE, 0, world);
  if(me != 0) {
    X = new double[N];
    Xs = new double[N];
--- a/src/USER-MISC/pair_meam_spline.h
+++ b/src/USER-MISC/pair_meam_spline.h
@ -142,7 +142,7 @@ protected:
          ((a*a*a - a) * Y2[klo] + (b*b*b - b) * Y2[khi])*(h*h)/6.0;
 #else
        // For a spline with regular grid, we directly calculate the interval X is in.
-        int klo = (int)(x / h);
+        int klo = (int)(x*inv_h);
        int khi = klo + 1;
        double a = Xs[khi] - x;
        double b = h - a;
@ -186,7 +186,7 @@ protected:
           (b*b*b - b) * Y2[khi]) * (h*h) / 6.0;
 #else
        // For a spline with regular grid, we directly calculate the interval X is in.
-        int klo = (int)(x / h);
+        int klo = (int)(x*inv_h);
        int khi = klo + 1;
        double a = Xs[khi] - x;
        double b = h - a;
@ -224,6 +224,7 @@ protected:
    int isGridSpline;// Indicates that all spline knots are on a regular grid.
    double h;        // The distance between knots if this is a grid spline with equidistant knots.
    double hsq;      // The squared distance between knots if this is a grid spline with equidistant knots.
    double inv_h;    // (1/h), used to avoid numerical errors in binnning for grid spline with equidistant knots.
    double xmax_shifted; // The end of the spline interval after it has been shifted to begin at X=0.
  };
--- a/src/USER-MISC/pair_meam_sw_spline.cpp
+++ b/src/USER-MISC/pair_meam_sw_spline.cpp
@ -674,6 +674,7 @@ void PairMEAMSWSpline::SplineFunction::prepareSpline(Error* error)
                Y2[i] /= h*6.0;
 #endif
        }
        inv_h = 1/h;
        xmax_shifted = xmax - xmin;
 }
@ -689,6 +690,7 @@ void PairMEAMSWSpline::SplineFunction::communicate(MPI_Comm& world, int me)
        MPI_Bcast(&isGridSpline, 1, MPI_INT, 0, world);
        MPI_Bcast(&h, 1, MPI_DOUBLE, 0, world);
        MPI_Bcast(&hsq, 1, MPI_DOUBLE, 0, world);
        MPI_Bcast(&inv_h, 1, MPI_DOUBLE, 0, world);
        if(me != 0) {
                X = new double[N];
                Xs = new double[N];
--- a/src/USER-MISC/pair_meam_sw_spline.h
+++ b/src/USER-MISC/pair_meam_sw_spline.h
@ -130,7 +130,7 @@ protected:
 #else
                                // For a spline with grid points, we can directly calculate the interval X is in.
                                //
-                                int klo = (int)(x / h);
+                                int klo = (int)(x*inv_h);
                                if ( klo > N - 2 ) klo = N - 2;
                                int khi = klo + 1;
                                double a = Xs[khi] - x;
@ -170,7 +170,7 @@ protected:
                                return a * Y[klo] + b * Y[khi] + ((a*a*a - a) * Y2[klo] + (b*b*b - b) * Y2[khi]) * (h*h) / 6.0;
 #else
                                // For a spline with grid points, we can directly calculate the interval X is in.
-                                int klo = (int)(x / h);
+                                int klo = (int)(x*inv_h);
                                if ( klo > N - 2 ) klo = N - 2;
                                int khi = klo + 1;
                                double a = Xs[khi] - x;
@ -207,6 +207,7 @@ protected:
                int isGridSpline;                // Indicates that all spline knots are on a regular grid.
                double h;                                // The distance between knots if this is a grid spline with equidistant knots.
                double hsq;                                // The squared distance between knots if this is a grid spline with equidistant knots.
                double inv_h;    // (1/h), used to avoid numerical errors in binnning for grid spline with equidistant knots.
                double xmax_shifted;        // The end of the spline interval after it has been shifted to begin at X=0.
        };
--- a/src/USER-OMP/pair_reaxc_omp.cpp
+++ b/src/USER-OMP/pair_reaxc_omp.cpp
@ -508,7 +508,7 @@ int PairReaxCOMP::estimate_reax_lists()
    num_nbrs += numneigh[i];
  }
-  int new_estimate = MAX (num_nbrs, mincap*MIN_NBRS);
+  int new_estimate = MAX(num_nbrs, mincap*REAX_MIN_NBRS);
  return new_estimate;
 }
--- a/src/USER-OMP/reaxc_forces_omp.cpp
+++ b/src/USER-OMP/reaxc_forces_omp.cpp
@ -310,7 +310,7 @@ void Validate_ListsOMP(reax_system *system, storage * /*workspace*/, reax_list *
      Hindex = system->my_atoms[i].Hindex;
      if (Hindex > -1) {
        system->my_atoms[i].num_hbonds =
-          (int)(MAX( Num_Entries(Hindex, hbonds)*saferzone, MIN_HBONDS ));
+          (int)(MAX(Num_Entries(Hindex,hbonds)*saferzone,system->minhbonds));
        if (Hindex < numH-1)
          comp = Start_Index(Hindex+1, hbonds);
--- a/src/USER-OMP/reaxc_init_md_omp.cpp
+++ b/src/USER-OMP/reaxc_init_md_omp.cpp
@ -68,7 +68,7 @@ int Init_ListsOMP(reax_system *system, control_params *control,
      system->my_atoms[i].num_hbonds = hb_top[i];
      total_hbonds += hb_top[i];
    }
-    total_hbonds = (int)(MAX( total_hbonds*saferzone, mincap*MIN_HBONDS ));
+    total_hbonds = (int)(MAX(total_hbonds*saferzone,mincap*system->minhbonds));
    if( !Make_List( system->Hcap, total_hbonds, TYP_HBOND,
                    *lists+HBONDS ) ) {
--- a/src/USER-REAXC/fix_qeq_reax.cpp
+++ b/src/USER-REAXC/fix_qeq_reax.cpp
@ -43,11 +43,8 @@ using namespace LAMMPS_NS;
 using namespace FixConst;
 #define EV_TO_KCAL_PER_MOL 14.4
 //#define DANGER_ZONE     0.95
 //#define LOOSE_ZONE      0.7
 #define SQR(x) ((x)*(x))
 #define CUBE(x) ((x)*(x)*(x))
 #define MIN_NBRS 100
 static const char cite_fix_qeq_reax[] =
  "fix qeq/reax command:\n\n"
@ -300,8 +297,8 @@ void FixQEqReax::allocate_matrix()
    mincap = reaxc->system->mincap;
    safezone = reaxc->system->safezone;
  } else {
-    mincap = MIN_CAP;
+    mincap = REAX_MIN_CAP;
-    safezone = SAFE_ZONE;
+    safezone = REAX_SAFE_ZONE;
  }
  n = atom->nlocal;
@ -324,7 +321,7 @@ void FixQEqReax::allocate_matrix()
    i = ilist[ii];
    m += numneigh[i];
  }
-  m_cap = MAX( (int)(m * safezone), mincap * MIN_NBRS);
+  m_cap = MAX( (int)(m * safezone), mincap * REAX_MIN_NBRS);
  H.n = n_cap;
  H.m = m_cap;
--- a/src/USER-REAXC/pair_reaxc.cpp
+++ b/src/USER-REAXC/pair_reaxc.cpp
@ -244,9 +244,10 @@ void PairReaxC::settings(int narg, char **arg)
  qeqflag = 1;
  control->lgflag = 0;
  control->enobondsflag = 1;
-  system->mincap = MIN_CAP;
+  system->mincap = REAX_MIN_CAP;
-  system->safezone = SAFE_ZONE;
+  system->minhbonds = REAX_MIN_HBONDS;
-  system->saferzone = SAFER_ZONE;
+  system->safezone = REAX_SAFE_ZONE;
  system->saferzone = REAX_SAFER_ZONE;
  // process optional keywords
@ -284,6 +285,12 @@ void PairReaxC::settings(int narg, char **arg)
      if (system->mincap < 0)
        error->all(FLERR,"Illegal pair_style reax/c mincap command");
      iarg += 2;
    } else if (strcmp(arg[iarg],"minhbonds") == 0) {
      if (iarg+2 > narg) error->all(FLERR,"Illegal pair_style reax/c command");
      system->minhbonds = force->inumeric(FLERR,arg[iarg+1]);
      if (system->minhbonds < 0)
        error->all(FLERR,"Illegal pair_style reax/c minhbonds command");
      iarg += 2;
    } else error->all(FLERR,"Illegal pair_style reax/c command");
  }
@ -712,7 +719,7 @@ int PairReaxC::estimate_reax_lists()
  free( marked );
-  return static_cast<int> (MAX( num_nbrs*safezone, mincap*MIN_NBRS ));
+  return static_cast<int> (MAX(num_nbrs*safezone, mincap*REAX_MIN_NBRS));
 }
 /* ---------------------------------------------------------------------- */
--- a/src/USER-REAXC/reaxc_allocate.cpp
+++ b/src/USER-REAXC/reaxc_allocate.cpp
@ -340,7 +340,7 @@ static int Reallocate_HBonds_List( reax_system *system, reax_list *hbonds )
    if ((system->my_atoms[i].Hindex) >= 0) {
      total_hbonds += system->my_atoms[i].num_hbonds;
    }
-  total_hbonds = (int)(MAX( total_hbonds*saferzone, mincap*MIN_HBONDS ));
+  total_hbonds = (int)(MAX(total_hbonds*saferzone, mincap*system->minhbonds));
  Delete_List( hbonds);
  if (!Make_List( system->Hcap, total_hbonds, TYP_HBOND, hbonds )) {
@ -456,7 +456,7 @@ void ReAllocate( reax_system *system, control_params *control,
      }
      newsize = static_cast<int>
-        (MAX( realloc->num_far*safezone, mincap*MIN_NBRS ));
+        (MAX( realloc->num_far*safezone, mincap*REAX_MIN_NBRS));
      Reallocate_Neighbor_List( far_nbrs, system->total_cap, newsize);
      realloc->num_far = 0;
--- a/src/USER-REAXC/reaxc_defs.h
+++ b/src/USER-REAXC/reaxc_defs.h
@ -75,28 +75,23 @@
 #define MAX_TOKENS          1024
 #define MAX_TOKEN_LEN       1024
 #define MAX_ATOM_ID         100000
 #define MAX_RESTRICT        15
 #define MAX_MOLECULE_SIZE   20
 #define MAX_ATOM_TYPES      25
 #define NUM_INTRS      10
 #define ALMOST_ZERO    1e-10
 #define NEG_INF       -1e10
 #define NO_BOND        1e-3  // 0.001
 #define HB_THRESHOLD   1e-2  // 0.01
-#define MIN_CAP        50
+#define REAX_MIN_CAP    50
-#define MIN_NBRS       100
+#define REAX_MIN_NBRS   100
 #define MIN_HENTRIES    100
 #define MAX_BONDS       30
 #define MIN_BONDS       25
-#define MIN_HBONDS     25
+#define REAX_MIN_HBONDS 25
 #define MIN_3BODIES     1000
 #define MIN_GCELL_POPL  50
 #define MIN_SEND        100
-#define SAFE_ZONE      1.2
+#define REAX_SAFE_ZONE  1.2
-#define SAFER_ZONE     1.4
+#define REAX_SAFER_ZONE 1.4
 #define DANGER_ZONE     0.90
 #define LOOSE_ZONE      0.75
 #define MAX_3BODY_PARAM 5
--- a/src/USER-REAXC/reaxc_forces.cpp
+++ b/src/USER-REAXC/reaxc_forces.cpp
@ -152,7 +152,7 @@ void Validate_Lists( reax_system *system, storage * /*workspace*/, reax_list **l
      Hindex = system->my_atoms[i].Hindex;
      if (Hindex > -1) {
        system->my_atoms[i].num_hbonds =
-          (int)(MAX( Num_Entries(Hindex, hbonds)*saferzone, MIN_HBONDS ));
+          (int)(MAX(Num_Entries(Hindex, hbonds)*saferzone, system->minhbonds));
        //if( Num_Entries(i, hbonds) >=
        //(Start_Index(i+1,hbonds)-Start_Index(i,hbonds))*0.90/*DANGER_ZONE*/){
@ -423,7 +423,7 @@ void Estimate_Storages( reax_system *system, control_params *control,
  *Htop = (int)(MAX( *Htop * safezone, mincap * MIN_HENTRIES ));
  for( i = 0; i < system->n; ++i )
-    hb_top[i] = (int)(MAX( hb_top[i] * saferzone, MIN_HBONDS ));
+    hb_top[i] = (int)(MAX(hb_top[i] * saferzone, system->minhbonds));
  for( i = 0; i < system->N; ++i ) {
    *num_3body += SQR(bond_top[i]);
--- a/src/USER-REAXC/reaxc_init_md.cpp
+++ b/src/USER-REAXC/reaxc_init_md.cpp
@ -179,7 +179,7 @@ int  Init_Lists( reax_system *system, control_params *control,
      system->my_atoms[i].num_hbonds = hb_top[i];
      total_hbonds += hb_top[i];
    }
-    total_hbonds = (int)(MAX( total_hbonds*saferzone, mincap*MIN_HBONDS ));
+    total_hbonds = (int)(MAX(total_hbonds*saferzone,mincap*system->minhbonds));
    if( !Make_List( system->Hcap, total_hbonds, TYP_HBOND,
                    *lists+HBONDS ) ) {
--- a/src/USER-REAXC/reaxc_lookup.cpp
+++ b/src/USER-REAXC/reaxc_lookup.cpp
@ -27,7 +27,6 @@
 #include "reaxc_lookup.h"
 #include <mpi.h>
 #include <cstdlib>
 #include "reaxc_defs.h"
 #include "reaxc_nonbonded.h"
 #include "reaxc_tool_box.h"
@ -153,7 +152,7 @@ int Init_Lookup_Tables( reax_system *system, control_params *control,
 {
  int i, j, r;
  int num_atom_types;
-  int existing_types[MAX_ATOM_TYPES], aggregated[MAX_ATOM_TYPES];
+  int existing_types[REAX_MAX_ATOM_TYPES], aggregated[REAX_MAX_ATOM_TYPES];
  double dr;
  double *h, *fh, *fvdw, *fele, *fCEvd, *fCEclmb;
  double v0_vdw, v0_ele, vlast_vdw, vlast_ele;
@ -186,12 +185,12 @@ int Init_Lookup_Tables( reax_system *system, control_params *control,
    LR[i] = (LR_lookup_table*)
      scalloc(system->error_ptr,  num_atom_types, sizeof(LR_lookup_table), "lookup:LR[i]");
-  for( i = 0; i < MAX_ATOM_TYPES; ++i )
+  for( i = 0; i < REAX_MAX_ATOM_TYPES; ++i )
    existing_types[i] = 0;
  for( i = 0; i < system->n; ++i )
    existing_types[ system->my_atoms[i].type ] = 1;
-  MPI_Allreduce( existing_types, aggregated, MAX_ATOM_TYPES,
+  MPI_Allreduce( existing_types, aggregated, REAX_MAX_ATOM_TYPES,
                 MPI_INT, MPI_SUM, mpi_data->world );
  for( i = 0; i < num_atom_types; ++i ) {
--- a/src/USER-REAXC/reaxc_traj.cpp
+++ b/src/USER-REAXC/reaxc_traj.cpp
@ -39,11 +39,11 @@ int Reallocate_Output_Buffer( LAMMPS_NS::Error *error_ptr, output_controls *out_
  if (out_control->buffer_len > 0)
    free( out_control->buffer );
-  out_control->buffer_len = (int)(req_space*SAFE_ZONE);
+  out_control->buffer_len = (int)(req_space*REAX_SAFE_ZONE);
  out_control->buffer = (char*) malloc(out_control->buffer_len*sizeof(char));
  if (out_control->buffer == NULL) {
    char errmsg[256];
-    snprintf(errmsg, 256, "Insufficient memory for required buffer size %d", (int) (req_space*SAFE_ZONE));
+    snprintf(errmsg, 256, "Insufficient memory for required buffer size %d", (int) (req_space*REAX_SAFE_ZONE));
    error_ptr->one(FLERR,errmsg);
  }
--- a/src/USER-REAXC/reaxc_types.h
+++ b/src/USER-REAXC/reaxc_types.h
@ -409,7 +409,7 @@ struct _reax_system
  class LAMMPS_NS::Error *error_ptr;
  class LAMMPS_NS::Pair *pair_ptr;
  int my_bonds;
-  int mincap;
+  int mincap,minhbonds;
  double safezone, saferzone;
  _LR_lookup_table **LR;
--- a/src/compute_coord_atom.cpp
+++ b/src/compute_coord_atom.cpp
@ -27,6 +27,7 @@
 #include "group.h"
 #include "memory.h"
 #include "error.h"
 #include "utils.h"
 using namespace LAMMPS_NS;
@ -93,7 +94,7 @@ ComputeCoordAtom::ComputeCoordAtom(LAMMPS *lmp, int narg, char **arg) :
    int iorientorder = modify->find_compute(id_orientorder);
    if (iorientorder < 0)
      error->all(FLERR,"Could not find compute coord/atom compute ID");
-    if (strcmp(modify->compute[iorientorder]->style,"orientorder/atom") != 0)
+    if (!utils::strmatch(modify->compute[iorientorder]->style,"^orientorder/atom"))
      error->all(FLERR,"Compute coord/atom compute ID is not orientorder/atom");
    threshold = force->numeric(FLERR,arg[5]);
@ -119,6 +120,8 @@ ComputeCoordAtom::ComputeCoordAtom(LAMMPS *lmp, int narg, char **arg) :
 ComputeCoordAtom::~ComputeCoordAtom()
 {
  if (copymode) return;
  delete [] group2;
  delete [] typelo;
  delete [] typehi;
@ -138,7 +141,7 @@ void ComputeCoordAtom::init()
    l = c_orientorder->qlcomp;
    //  communicate real and imaginary 2*l+1 components of the normalized vector
    comm_forward = 2*(2*l+1);
-    if (c_orientorder->iqlcomp < 0)
+    if (!(c_orientorder->qlcompflag))
      error->all(FLERR,"Compute coord/atom requires components "
                 "option in compute orientorder/atom");
  }
--- a/src/compute_coord_atom.h
+++ b/src/compute_coord_atom.h
@ -27,16 +27,16 @@ namespace LAMMPS_NS {
 class ComputeCoordAtom : public Compute {
 public:
  ComputeCoordAtom(class LAMMPS *, int, char **);
-  ~ComputeCoordAtom();
+  virtual ~ComputeCoordAtom();
-  void init();
+  virtual void init();
  void init_list(int, class NeighList *);
-  void compute_peratom();
+  virtual void compute_peratom();
  int pack_forward_comm(int, int *, double *, int, int *);
  void unpack_forward_comm(int, int, double *);
  double memory_usage();
  enum {NONE,CUTOFF,ORIENT};
- private:
+ protected:
  int nmax,ncol;
  double cutsq;
  class NeighList *list;
--- a/src/fix_halt.cpp
+++ b/src/fix_halt.cpp
@ -306,12 +306,12 @@ double FixHalt::tlimit()
 /* ----------------------------------------------------------------------
   determine available disk space, if supported. Return -1 if not.
 ------------------------------------------------------------------------- */
-#if defined(__linux) || defined(__APPLE__)
+#if defined(__linux__) || defined(__APPLE__) || defined(__FreeBSD__) || defined(__DragonFly__) || defined(__OpenBSD__) || defined(__NetBSD__)
 #include <sys/statvfs.h>
 #endif
 double FixHalt::diskfree()
 {
-#if defined(__linux) || defined(__APPLE__)
+#if defined(__linux__) || defined(__APPLE__) || defined(__FreeBSD__) || defined(__DragonFly__) || defined(__OpenBSD__) || defined(__NetBSD__)
  struct statvfs fs;
  double disk_free = -1.0;
@ -319,9 +319,9 @@ double FixHalt::diskfree()
    disk_free = 1.0e100;
    int rv = statvfs(dlimit_path,&fs);
    if (rv == 0) {
-#if defined(__linux)
+#if defined(__linux__)
      disk_free = fs.f_bavail*fs.f_bsize/1048576.0;
-#elif defined(__APPLE__)
+#elif defined(__APPLE__) || defined(__FreeBSD__) || defined(__DragonFly__) || defined(__OpenBSD__) || defined(__NetBSD__)
      disk_free = fs.f_bavail*fs.f_frsize/1048576.0;
 #endif
    } else
--- a/src/fix_restrain.cpp
+++ b/src/fix_restrain.cpp
@ -44,8 +44,9 @@ enum{BOND,LBOUND,ANGLE,DIHEDRAL};
 FixRestrain::FixRestrain(LAMMPS *lmp, int narg, char **arg) :
  Fix(lmp, narg, arg),
-  rstyle(NULL), mult(NULL), ids(NULL), kstart(NULL), kstop(NULL), target(NULL),
+  rstyle(NULL), mult(NULL), ids(NULL), kstart(NULL), kstop(NULL),
-  deqstart(NULL), deqstop(NULL), cos_target(NULL), sin_target(NULL)
+  deqstart(NULL), deqstop(NULL), target(NULL), cos_target(NULL),
  sin_target(NULL)
 {
  if (narg < 4) error->all(FLERR,"Illegal fix restrain command");
--- a/src/info.cpp
+++ b/src/info.cpp
@ -60,7 +60,7 @@
 #include <sys/utsname.h>
 #endif
-#if defined(__linux)
+#if defined(__linux__)
 #include <malloc.h>
 #endif
@ -345,7 +345,7 @@ void Info::command(int narg, char **arg)
    fprintf(out,"Maximum working set size: %.4g Mbyte\n",
            (double)pmc.PeakWorkingSetSize/1048576.0);
 #else
-#if defined(__linux)
+#if defined(__linux__)
    struct mallinfo mi;
    mi = mallinfo();
    fprintf(out,"Current reserved memory pool size: %.4g Mbyte\n",
--- a/src/random_mars.cpp
+++ b/src/random_mars.cpp
@ -201,7 +201,7 @@ void RanMars::select_subset(bigint ntarget, int nmine, int *mark, int *next)
  for (int i = 0; i < nmine; i++) mark[i] = 0;
  for (int i = 0; i < nmine; i++) next[i] = i+1;
-  next[nmine-1] = -1;
+  if (nmine > 0) next[nmine-1] = -1;
  nmark = 0;
  niter = 0;
@ -241,7 +241,7 @@ void RanMars::select_subset(bigint ntarget, int nmine, int *mark, int *next)
    // flip each value based on RN < thresh
    nflip = 0;
-    while (index >= 0) {
+    while ((nmine > 0) && (index >= 0)) {
      if (uniform() < thresh) {
        mark[index] = newvalue;
        nflip++;
--- a/src/utils.cpp
+++ b/src/utils.cpp
@ -17,7 +17,7 @@
 #include "lammps.h"
 #include "error.h"
-#if defined(__linux)
+#if defined(__linux__)
 #include <unistd.h>  // for readlink
 #endif
@ -95,7 +95,7 @@ static const char *guesspath(char *buf, int len, FILE *fp)
 {
  memset(buf,0,len);
-#if defined(__linux)
+#if defined(__linux__)
  char procpath[32];
  int fd = fileno(fp);
  snprintf(procpath,32,"/proc/self/fd/%d",fd);
--- a/src/version.h
+++ b/src/version.h
@ -1 +1 @@
-#define LAMMPS_VERSION "15 Apr 2020"
+#define LAMMPS_VERSION "5 May 2020"
--- a/tools/singularity/centos7.def
+++ b/tools/singularity/centos7.def
@ -7,16 +7,35 @@ From: centos:7
        yum -y install vim-enhanced \
            ccache gcc-c++ gcc-gfortran clang gdb valgrind-openmpi \
            make cmake cmake3 ninja-build patch which file git Lmod \
-            libpng-devel libjpeg-devel openmpi-devel mpich-devel python-devel \
+            libpng-devel libjpeg-devel openmpi-devel mpich-devel python-devel python36-devel \
-            python-virtualenv fftw-devel voro++-devel eigen3-devel gsl-devel openblas-devel enchant
+            hdf5-devel python36-virtualenv python36-pip python-pip \
            netcdf-devel netcdf-cxx-devel netcdf-mpich-devel netcdf-openmpi-devel \
            python-virtualenv fftw-devel voro++-devel eigen3-devel gsl-devel openblas-devel enchant \
            blas-devel lapack-devel
        yum clean all
 %environment
        LC_ALL=C
        export LC_ALL
        # we need to reset any module variables
        # inherited from the host.
        unset __LMOD_REF_COUNT__LMFILES_
        unset __LMOD_REF_COUNT_PATH
        unset __LMOD_REF_COUNT_LD_LIBRARY_PATH
        unset __LMOD_REF_COUNT_MANPATH
        unset __LMOD_REF_COUNT_MODULEPATH
        unset __LMOD_REF_COUNT_LOADEDMODULES
        unset _LMFILES_
        unset MODULEPATH
        unset MODULESHOME
        unset MODULEPATH_ROOT
        unset LOADEDMODULES
-        . /etc/profile.d/z00_lmod.sh
+        unset LMOD_SYSTEM_DEFAULT_MODULES
-        module purge
+
        # load MPI by default
        . /etc/profile
        module load mpi
 %labels
-        Author akohlmey
+        Author akohlmey, rbberger
--- a/tools/singularity/centos8.def
+++ b/tools/singularity/centos8.def
@ -8,21 +8,39 @@ From: centos:8
        dnf -y install vim-enhanced git file make cmake patch which file ninja-build \
               ccache gcc-c++ gcc-gfortran clang gdb valgrind libubsan libasan libtsan \
               eigen3-devel openblas-devel libpng-devel libjpeg-devel platform-python-devel \
-               openmpi-devel mpich-devel fftw-devel voro++-devel gsl-devel \
+               openmpi-devel mpich-devel fftw-devel voro++-devel gsl-devel hdf5-devel \
               netcdf-devel netcdf-cxx-devel netcdf-mpich-devel netcdf-openmpi-devel \
               enchant python3-virtualenv doxygen \
               texlive-latex-fonts texlive-pslatex texlive-collection-latexrecommended \
               texlive-latex texlive-latexconfig doxygen-latex texlive-collection-latex \
               texlive-latex-bin texlive-lualatex-math texlive-fncychap texlive-tabulary \
               texlive-framed texlive-wrapfig texlive-upquote texlive-capt-of \
-               texlive-needspace texlive-titlesec texlive-anysize texlive-dvipng
+               texlive-needspace texlive-titlesec texlive-anysize texlive-dvipng \
               blas-devel lapack-devel
        dnf clean all
 %environment
        LC_ALL=C
        export LC_ALL
        # we need to reset any module variables
        # inherited from the host.
        unset __LMOD_REF_COUNT__LMFILES_
        unset __LMOD_REF_COUNT_PATH
        unset __LMOD_REF_COUNT_LD_LIBRARY_PATH
        unset __LMOD_REF_COUNT_MANPATH
        unset __LMOD_REF_COUNT_MODULEPATH
        unset __LMOD_REF_COUNT_LOADEDMODULES
        unset _LMFILES_
        unset MODULEPATH
        unset MODULESHOME
        unset MODULEPATH_ROOT
        unset LOADEDMODULES
-        . /etc/profile.d/modules.sh
+        unset LMOD_SYSTEM_DEFAULT_MODULES
-        module purge
+
        # load MPI by default
        . /etc/profile
        module load mpi
 %labels
-        Author akohlmey
+        Author akohlmey, rbberger
--- a/tools/singularity/fedora30_mingw.def
+++ b/tools/singularity/fedora30_mingw.def
@ -7,7 +7,8 @@ From: fedora:30
               ninja-build clang libomp-devel libubsan libasan libtsan \
               dos2unix findutils rsync python-devel libjpeg-devel libpng-devel \
               ccache gcc-c++ gcc-gfortran gdb valgrind eigen3-devel openblas-devel \
-               openmpi-devel mpich-devel fftw-devel voro++-devel gsl-devel \
+               openmpi-devel mpich-devel fftw-devel voro++-devel gsl-devel hdf5-devel \
               netcdf-devel netcdf-cxx-devel netcdf-mpich-devel netcdf-openmpi-devel \
               mingw-filesystem-base mingw32-nsis mingw-binutils-generic \
               mingw32-filesystem mingw32-pkg-config \
               mingw64-filesystem mingw64-pkg-config \
@ -33,17 +34,32 @@ From: fedora:30
               texlive-latex texlive-latexconfig doxygen-latex texlive-collection-latex \
               texlive-latex-bin texlive-lualatex-math texlive-fncychap texlive-tabulary \
               texlive-framed texlive-wrapfig texlive-upquote texlive-capt-of \
-               texlive-needspace texlive-titlesec texlive-anysize texlive-dvipng
+               texlive-needspace texlive-titlesec texlive-anysize texlive-dvipng \
               blas-devel lapack-devel
        dnf clean all
 %environment
        LC_ALL=C
        export LC_ALL
        # we need to reset any module variables
        # inherited from the host.
        unset __LMOD_REF_COUNT__LMFILES_
        unset __LMOD_REF_COUNT_PATH
        unset __LMOD_REF_COUNT_LD_LIBRARY_PATH
        unset __LMOD_REF_COUNT_MANPATH
        unset __LMOD_REF_COUNT_MODULEPATH
        unset __LMOD_REF_COUNT_LOADEDMODULES
        unset _LMFILES_
        unset MODULEPATH
        unset MODULESHOME
        unset MODULEPATH_ROOT
        unset LOADEDMODULES
-        . /etc/profile.d/modules.sh
+        unset LMOD_SYSTEM_DEFAULT_MODULES
-        module purge
+
        # load MPI by default
        . /etc/profile
        module load mpi
 %labels
-        Author akohlmey
+        Author akohlmey, rbberger
--- a/tools/singularity/fedora32_mingw.def
+++ b/tools/singularity/fedora32_mingw.def
@ -7,7 +7,8 @@ From: fedora:32
               ninja-build clang libomp-devel libubsan libasan libtsan \
               dos2unix findutils rsync python-devel libjpeg-devel libpng-devel \
               ccache gcc-c++ gcc-gfortran gdb valgrind eigen3-devel openblas-devel \
-               openmpi-devel mpich-devel fftw-devel voro++-devel gsl-devel \
+               openmpi-devel mpich-devel fftw-devel voro++-devel gsl-devel hdf5-devel \
               netcdf-devel netcdf-cxx-devel netcdf-mpich-devel netcdf-openmpi-devel \
               mingw-filesystem-base mingw32-nsis mingw-binutils-generic \
               mingw32-filesystem mingw32-pkg-config \
               mingw64-filesystem mingw64-pkg-config \
@ -33,17 +34,29 @@ From: fedora:32
               texlive-latex texlive-latexconfig doxygen-latex texlive-collection-latex \
               texlive-latex-bin texlive-lualatex-math texlive-fncychap texlive-tabulary \
               texlive-framed texlive-wrapfig texlive-upquote texlive-capt-of \
-               texlive-needspace texlive-titlesec texlive-anysize texlive-dvipng
+               texlive-needspace texlive-titlesec texlive-anysize texlive-dvipng \
               blas-devel lapack-devel diffutils
        dnf clean all
 %environment
        LC_ALL=C
        export LC_ALL
        # we need to reset any module variables
        # inherited from the host.
        unset __LMOD_REF_COUNT__LMFILES_
        unset __LMOD_REF_COUNT_PATH
        unset __LMOD_REF_COUNT_LD_LIBRARY_PATH
        unset __LMOD_REF_COUNT_MANPATH
        unset __LMOD_REF_COUNT_MODULEPATH
        unset __LMOD_REF_COUNT_LOADEDMODULES
        unset _LMFILES_
        unset MODULEPATH
        unset MODULESHOME
        unset MODULEPATH_ROOT
        unset LOADEDMODULES
-        . /etc/profile.d/modules.sh
+        unset LMOD_SYSTEM_DEFAULT_MODULES
-        module purge
+
        # load MPI by default
        . /etc/profile
        module load mpi
 %labels
-        Author akohlmey
+        Author akohlmey, rbberger
--- a/tools/singularity/ubuntu18.04.def
+++ b/tools/singularity/ubuntu18.04.def
@ -30,6 +30,7 @@ From: ubuntu:18.04
        libhwloc-dev \
        libjpeg-dev \
        liblapack-dev \
        libnetcdf-dev \
        libomp-dev \
        libopenblas-dev \
        libnuma-dev \
@ -61,11 +62,17 @@ From: ubuntu:18.04
        virtualenv \
        voro++-dev \
        wget \
-        xxd
+        xxd \
        valgrind \
        gdb
    # clean cache
    rm -rf /var/lib/apt/lists/*
 %environment
    LC_ALL=C
    export LC_ALL
    export PATH=/usr/lib/ccache:$PATH
 %labels
    Author akohlmey, rbberger
--- a/tools/singularity/ubuntu18.04_amd_rocm.def
+++ b/tools/singularity/ubuntu18.04_amd_rocm.def
@ -61,7 +61,9 @@ From: rocm/dev-ubuntu-18.04
        virtualenv \
        voro++-dev \
        wget \
-        xxd
+        xxd \
        valgrind \
        gdb
    export PATH=$PATH:/opt/rocm/bin:/opt/rocm/profiler/bin:/opt/rocm/opencl/bin/x86_64
    git clone -b master-rocm-3.3 https://github.com/ROCmSoftwarePlatform/hipCUB.git
@ -72,6 +74,9 @@ From: rocm/dev-ubuntu-18.04
    make package
    make install
    # clean cache
    rm -rf /var/lib/apt/lists/*
 %environment
    LC_ALL=C
    export LC_ALL
--- a/tools/singularity/ubuntu18.04_amd_rocm_cuda.def
+++ b/tools/singularity/ubuntu18.04_amd_rocm_cuda.def
@ -2,11 +2,11 @@ BootStrap: library
 From: lammps/default/lammps_development:ubuntu18.04_amd_rocm
 %environment
-    export PATH=/usr/local/nvidia/bin:/usr/local/cuda/bin:${PATH}
+    export PATH=/usr/local/cuda-10.2/bin:${PATH}:/opt/rocm/bin:/opt/rocm/profiler/bin:/opt/rocm/opencl/bin/x86_64
-    export CUDADIR=/usr/local/cuda
+    export CUDADIR=/usr/local/cuda-10.2
-    export CUDA_PATH=/usr/local/cuda
+    export CUDA_PATH=/usr/local/cuda-10.2
-    export LD_LIBRARY_PATH=$LD_LIBRARY_PATH:/usr/local/nvidia/lib:/usr/local/nvidia/lib64
+    export LD_LIBRARY_PATH=$LD_LIBRARY_PATH:/usr/local/cuda-10.2/lib64
-    export LIBRARY_PATH=/usr/local/cuda/lib64/stubs
+    export LIBRARY_PATH=/usr/local/cuda-10.2/lib64/stubs
 %post
    export DEBIAN_FRONTEND=noninteractive
@ -15,11 +15,24 @@ From: lammps/default/lammps_development:ubuntu18.04_amd_rocm
    apt-key adv --fetch-keys https://developer.download.nvidia.com/compute/cuda/repos/ubuntu1804/x86_64/7fa2af80.pub
    add-apt-repository "deb http://developer.download.nvidia.com/compute/cuda/repos/ubuntu1804/x86_64/ /"
    apt-get update
-    apt-get -y install cuda cuda-compiler-10-2
+
    export CUDA_PKG_VERSION=10.2
    apt-get install -y --no-install-recommends \
        cuda-libraries-$CUDA_PKG_VERSION \
        cuda-command-line-tools-$CUDA_PKG_VERSION \
        cuda-libraries-dev-$CUDA_PKG_VERSION \
        cuda-minimal-build-$CUDA_PKG_VERSION \
        cuda-compat-$CUDA_PKG_VERSION \
        libcublas10 \
        libcublas-dev
    echo "/usr/local/nvidia/lib" >> /etc/ld.so.conf.d/nvidia.conf
    echo "/usr/local/nvidia/lib64" >> /etc/ld.so.conf.d/nvidia.conf
    # clean cache
    rm -rf /var/lib/apt/lists/*
 %labels
    Author rbberger
--- a/tools/singularity/ubuntu18.04_gpu.def
+++ b/tools/singularity/ubuntu18.04_gpu.def
@ -2,11 +2,11 @@ BootStrap: docker
 From: rocm/dev-ubuntu-18.04
 %environment
-    export PATH=/usr/local/nvidia/bin:/usr/local/cuda/bin:${PATH}:/opt/rocm/bin:/opt/rocm/profiler/bin:/opt/rocm/opencl/bin/x86_64
+    export PATH=/usr/lib/ccache:/usr/local/cuda-10.2/bin:${PATH}:/opt/rocm/bin:/opt/rocm/profiler/bin:/opt/rocm/opencl/bin/x86_64
-    export CUDADIR=/usr/local/cuda
+    export CUDADIR=/usr/local/cuda-10.2
-    export CUDA_PATH=/usr/local/cuda
+    export CUDA_PATH=/usr/local/cuda-10.2
-    export LD_LIBRARY_PATH=$LD_LIBRARY_PATH:/usr/local/nvidia/lib:/usr/local/nvidia/lib64
+    export LD_LIBRARY_PATH=$LD_LIBRARY_PATH:/usr/local/cuda-10.2/lib64
-    export LIBRARY_PATH=/usr/local/cuda/lib64/stubs
+    export LIBRARY_PATH=/usr/local/cuda-10.2/lib64/stubs
 %post
    export DEBIAN_FRONTEND=noninteractive
    apt-get update
@ -64,15 +64,26 @@ From: rocm/dev-ubuntu-18.04
        virtualenv \
        voro++-dev \
        wget \
-        xxd
+        xxd \
        valgrind \
        gdb
    wget https://developer.download.nvidia.com/compute/cuda/repos/ubuntu1804/x86_64/cuda-ubuntu1804.pin
    mv cuda-ubuntu1804.pin /etc/apt/preferences.d/cuda-repository-pin-600
    apt-key adv --fetch-keys https://developer.download.nvidia.com/compute/cuda/repos/ubuntu1804/x86_64/7fa2af80.pub
    add-apt-repository "deb http://developer.download.nvidia.com/compute/cuda/repos/ubuntu1804/x86_64/ /"
    apt-get update
-    apt-get install --no-install-recommends -y \
+
-        cuda \
+    export CUDA_PKG_VERSION=10.2
-        cuda-compiler-10-2 \
+
    apt-get install -y --no-install-recommends \
        cuda-libraries-$CUDA_PKG_VERSION \
        cuda-command-line-tools-$CUDA_PKG_VERSION \
        cuda-libraries-dev-$CUDA_PKG_VERSION \
        cuda-minimal-build-$CUDA_PKG_VERSION \
        cuda-compat-$CUDA_PKG_VERSION \
        libcublas10 \
        libcublas-dev
    export PATH=$PATH:/opt/rocm/bin:/opt/rocm/profiler/bin:/opt/rocm/opencl/bin/x86_64
    git clone -b master-rocm-3.3 https://github.com/ROCmSoftwarePlatform/hipCUB.git
@ -86,6 +97,13 @@ From: rocm/dev-ubuntu-18.04
    echo "/usr/local/nvidia/lib" >> /etc/ld.so.conf.d/nvidia.conf
    echo "/usr/local/nvidia/lib64" >> /etc/ld.so.conf.d/nvidia.conf
    # clean cache
    rm -rf /var/lib/apt/lists/*
    # add missing symlink
    ln -s /usr/local/cuda-10.2 /usr/local/cuda
    ln -s /usr/local/cuda-10.2/lib64/stubs/libcuda.so /usr/local/cuda-10.2/lib64/stubs/libcuda.so.1
 %environment
    LC_ALL=C
    export LC_ALL
--- a/tools/singularity/ubuntu18.04_intel_opencl.def
+++ b/tools/singularity/ubuntu18.04_intel_opencl.def
@ -59,12 +59,17 @@ From: ubuntu:18.04
        virtualenv \
        voro++-dev \
        wget \
-        xxd
+        xxd \
        valgrind \
        gdb
    add-apt-repository ppa:intel-opencl/intel-opencl
    apt-get update
    apt-get install -y intel-opencl-icd
    # clean cache
    rm -rf /var/lib/apt/lists/*
 %environment
        LC_ALL=C
        export LC_ALL
--- a/tools/singularity/ubuntu18.04_nvidia.def
+++ b/tools/singularity/ubuntu18.04_nvidia.def
@ -59,11 +59,16 @@ From: nvidia/cuda:10.2-devel-ubuntu18.04
        virtualenv \
        voro++-dev \
        wget \
-        xxd
+        xxd \
        valgrind \
        gdb
    mkdir -p /etc/OpenCL/vendors
    echo "libnvidia-opencl.so.1" > /etc/OpenCL/vendors/nvidia.icd
    # clean cache
    rm -rf /var/lib/apt/lists/*
 %environment
        LC_ALL=C
        export LC_ALL
--- a/tools/singularity/ubuntu20.04.def
+++ b/tools/singularity/ubuntu20.04.def
@ -30,6 +30,7 @@ From: ubuntu:20.04
        libhwloc-dev \
        libjpeg-dev \
        liblapack-dev \
        libnetcdf-dev \
        libomp-dev \
        libopenblas-dev \
        libnuma-dev \
@ -57,7 +58,12 @@ From: ubuntu:20.04
        virtualenv \
        voro++-dev \
        wget \
-        xxd
+        xxd \
        valgrind \
        gdb
    # clean cache
    rm -rf /var/lib/apt/lists/*
 %environment
    LC_ALL=C
`@ -1 +1 @@`
	`#define LAMMPS_VERSION "15 Apr 2020"`	`#define LAMMPS_VERSION "5 May 2020"`