Merge pull request #4044 from lammps/maintenance

Third Set of Collected Bug Fixes and Maintenance Updates for 2 August 2023 Stable release

.github/CODEOWNERS

@@ -151,12 +151,12 @@ tools/vim/*           @hammondkd
 unittest/*            @akohlmey
 
 # cmake
-cmake/*               @rbberger
+cmake/*               @akohlmey
 cmake/Modules/LAMMPSInterfacePlugin.cmake @akohlmey
 cmake/Modules/MPI4WIN.cmake @akohlmey
 cmake/Modules/OpenCLLoader.cmake @akohlmey
-cmake/Modules/Packages/COLVARS.cmake @rbberger @giacomofiorin
-cmake/Modules/Packages/KIM.cmake @rbberger @ellio167
+cmake/Modules/Packages/COLVARS.cmake @giacomofiorin
+cmake/Modules/Packages/KIM.cmake @ellio167
 cmake/presets/*.cmake @akohlmey
 
 # python

cmake/CMakeLists.txt

@@ -2,7 +2,6 @@
 ########################################
 # CMake build system
 # This file is part of LAMMPS
-# Created by Christoph Junghans and Richard Berger
 cmake_minimum_required(VERSION 3.10)
 ########################################
 # set policy to silence warnings about ignoring <PackageName>_ROOT but use it
@@ -112,7 +111,7 @@ if(CMAKE_CXX_COMPILER_ID STREQUAL "Intel")
     if(CMAKE_CXX_COMPILER_VERSION VERSION_EQUAL 17.3 OR CMAKE_CXX_COMPILER_VERSION VERSION_EQUAL 17.4)
       set(CMAKE_TUNE_DEFAULT "-xCOMMON-AVX512")
     else()
-      set(CMAKE_TUNE_DEFAULT "-xHost -fp-model fast=2 -no-prec-div -qoverride-limits -diag-disable=10441 -diag-disable=2196")
+      set(CMAKE_TUNE_DEFAULT "-xHost -fp-model fast=2 -no-prec-div -qoverride-limits -diag-disable=10441 -diag-disable=11074 -diag-disable=11076 -diag-disable=2196")
     endif()
   endif()
 endif()
@@ -127,6 +126,19 @@ if((CMAKE_CXX_COMPILER_ID STREQUAL "NVHPC") OR (CMAKE_CXX_COMPILER_ID STREQUAL "
   set(CMAKE_TUNE_DEFAULT "-Minform=severe")
 endif()
 
+# this hack is required to compile fmt lib with CrayClang version 15.0.2
+# CrayClang is only directly recognized by version 3.28 and later
+if(CMAKE_VERSION VERSION_LESS 3.28)
+  get_filename_component(_exe "${CMAKE_CXX_COMPILER}" NAME)
+  if((CMAKE_CXX_COMPILER_ID STREQUAL "Clang") AND (_exe STREQUAL "crayCC"))
+    set(CMAKE_TUNE_DEFAULT "-DFMT_STATIC_THOUSANDS_SEPARATOR")
+  endif()
+else()
+  if(CMAKE_CXX_COMPILER_ID STREQUAL "CrayClang")
+    set(CMAKE_TUNE_DEFAULT "-DFMT_STATIC_THOUSANDS_SEPARATOR")
+  endif()
+endif()
+
 # silence nvcc warnings
 if((PKG_KOKKOS) AND (Kokkos_ENABLE_CUDA) AND NOT (CMAKE_CXX_COMPILER_ID STREQUAL "Clang"))
   set(CMAKE_TUNE_DEFAULT "${CMAKE_TUNE_DEFAULT} -Xcudafe --diag_suppress=unrecognized_pragma")
@@ -151,6 +163,7 @@ if(MSVC)
     add_compile_options(/Zc:__cplusplus)
     add_compile_options(/wd4244)
     add_compile_options(/wd4267)
+    add_compile_options(/wd4250)
     add_compile_options(/EHsc)
   endif()
   add_compile_definitions(_CRT_SECURE_NO_WARNINGS)
@@ -212,6 +225,10 @@ endif()
 add_executable(lmp ${MAIN_SOURCES})
 target_link_libraries(lmp PRIVATE lammps)
 set_target_properties(lmp PROPERTIES OUTPUT_NAME ${LAMMPS_BINARY})
+# re-export all symbols for plugins
+if(PKG_PLUGIN AND (NOT ((CMAKE_SYSTEM_NAME STREQUAL "Windows"))))
+  set_target_properties(lmp PROPERTIES ENABLE_EXPORTS TRUE)
+endif()
 install(TARGETS lmp EXPORT LAMMPS_Targets DESTINATION ${CMAKE_INSTALL_BINDIR})
 
 option(CMAKE_VERBOSE_MAKEFILE "Generate verbose Makefiles" OFF)
@@ -426,6 +443,7 @@ if(BUILD_OMP)
       (CMAKE_CXX_COMPILER_ID STREQUAL "IntelLLVM") OR (CMAKE_CXX_COMPILER_ID STREQUAL "XLClang") OR
       ((CMAKE_CXX_COMPILER_ID STREQUAL "AppleClang") AND (CMAKE_CXX_COMPILER_VERSION VERSION_GREATER_EQUAL 10.0)) OR
       ((CMAKE_CXX_COMPILER_ID STREQUAL "Clang") AND (CMAKE_CXX_COMPILER_VERSION VERSION_GREATER_EQUAL 10.0)) OR
+      ((CMAKE_CXX_COMPILER_ID STREQUAL "CrayClang") AND (CMAKE_CXX_COMPILER_VERSION VERSION_GREATER_EQUAL 10.0)) OR
       ((CMAKE_CXX_COMPILER_ID STREQUAL "Intel") AND (CMAKE_CXX_COMPILER_VERSION VERSION_GREATER_EQUAL 19.0)))
     # GCC 9.x and later plus Clang 10.x and later implement strict OpenMP 4.0 semantics for consts.
     # Intel 18.0 was tested to support both, so we switch to OpenMP 4+ from 19.x onward to be safe.
@@ -438,6 +456,18 @@ if(BUILD_OMP)
   target_link_libraries(lmp PRIVATE OpenMP::OpenMP_CXX)
 endif()
 
+# lower C++ standard for fmtlib sources when using Intel classic compiler
+if((CMAKE_CXX_COMPILER_ID STREQUAL "Intel") AND (CMAKE_CXX_STANDARD GREATER_EQUAL 17)
+    AND (CMAKE_CXX_COMPILER_VERSION VERSION_LESS 2021.10))
+  message(STATUS "Lowering C++ standard for compiling fmtlib sources with Intel Classic compiler")
+  get_filename_component(LMP_UTILS_SRC "${LAMMPS_SOURCE_DIR}/utils.cpp" ABSOLUTE)
+  get_filename_component(LMP_VARIABLE_SRC "${LAMMPS_SOURCE_DIR}/variable.cpp" ABSOLUTE)
+  get_filename_component(FMT_FORMAT_SRC "${LAMMPS_SOURCE_DIR}/fmtlib_format.cpp" ABSOLUTE)
+  get_filename_component(FMT_OS_SRC "${LAMMPS_SOURCE_DIR}/fmtlib_os.cpp" ABSOLUTE)
+  set_source_files_properties("${FMT_FORMAT_SRC}" "${FMT_OS_SRC}" "${LMP_VARIABLE_SRC}" "${LMP_UTILS_SRC}"
+          PROPERTIES COMPILE_OPTIONS "-std=c++14")
+endif()
+
 if(PKG_MSCG OR PKG_ATC OR PKG_AWPMD OR PKG_ML-QUIP OR PKG_ML-POD OR PKG_ELECTRODE OR BUILD_TOOLS)
   enable_language(C)
   if (NOT USE_INTERNAL_LINALG)

cmake/Modules/ExternalCMakeProject.cmake

@@ -43,5 +43,5 @@ function(ExternalCMakeProject target url hash basedir cmakedir cmakefile)
       "${CMAKE_BINARY_DIR}/_deps/${target}-src/${cmakedir}/CMakeLists.txt")
   endif()
   add_subdirectory("${CMAKE_BINARY_DIR}/_deps/${target}-src/${cmakedir}"
-    "${CMAKE_BINARY_DIR}/_deps/${target}-build")
+    "${CMAKE_BINARY_DIR}/_deps/${target}-build" EXCLUDE_FROM_ALL)
 endfunction(ExternalCMakeProject)

cmake/Modules/LAMMPSUtils.cmake

@@ -83,17 +83,17 @@ function(check_for_autogen_files source_dir)
     file(GLOB SRC_AUTOGEN_FILES ${CONFIGURE_DEPENDS} ${source_dir}/style_*.h)
     file(GLOB SRC_AUTOGEN_PACKAGES ${CONFIGURE_DEPENDS} ${source_dir}/packages_*.h)
     list(APPEND SRC_AUTOGEN_FILES ${SRC_AUTOGEN_PACKAGES} ${source_dir}/lmpinstalledpkgs.h ${source_dir}/lmpgitversion.h)
-    list(APPEND SRC_AUTOGEN_FILES ${SRC_AUTOGEN_PACKAGES} ${source_dir}/mliap_model_python_couple.h ${source_dir}/mliap_model_python_couple.cpp)
+    list(APPEND SRC_AUTOGEN_FILES ${source_dir}/mliap_model_python_couple.h ${source_dir}/mliap_model_python_couple.cpp)
     foreach(_SRC ${SRC_AUTOGEN_FILES})
       get_filename_component(FILENAME "${_SRC}" NAME)
       if(EXISTS ${source_dir}/${FILENAME})
         message(FATAL_ERROR "\n########################################################################\n"
-                              "Found header file(s) generated by the make-based build system\n"
-                              "\n"
-                              "Please run\n"
-                              "make -C ${source_dir} purge\n"
-                              "to remove\n"
-                              "########################################################################")
+                            "Found header file ${source_dir}/${FILENAME} generated by the make-based build system\n"
+                            "\n"
+                            "Please run\n"
+                            "make -C ${source_dir} purge\n"
+                            "to remove\n"
+                            "########################################################################")
       endif()
     endforeach()
 endfunction()

cmake/Modules/Packages/KOKKOS.cmake

@@ -132,8 +132,12 @@ if(PKG_KSPACE)
                                  ${KOKKOS_PKG_SOURCES_DIR}/remap_kokkos.cpp)
   if(Kokkos_ENABLE_CUDA)
     if(NOT (FFT STREQUAL "KISS"))
+      find_library(CUFFT_LIBRARY cufft)
+      if (CUFFT_LIBRARY STREQUAL "CUFFT_LIBRARY-NOTFOUND")
+        message(FATAL_ERROR "Required cuFFT library not found. Check your environment or set CUFFT_LIBRARY to its location")
+      endif()
       target_compile_definitions(lammps PRIVATE -DFFT_CUFFT)
-      target_link_libraries(lammps PRIVATE cufft)
+      target_link_libraries(lammps PRIVATE ${CUFFT_LIBRARY})
     endif()
   elseif(Kokkos_ENABLE_HIP)
     if(NOT (FFT STREQUAL "KISS"))

cmake/Modules/Packages/MDI.cmake

@@ -8,8 +8,8 @@ option(DOWNLOAD_MDI "Download and compile the MDI library instead of using an al
 
 if(DOWNLOAD_MDI)
   message(STATUS "MDI download requested - we will build our own")
-  set(MDI_URL "https://github.com/MolSSI-MDI/MDI_Library/archive/v1.4.16.tar.gz" CACHE STRING "URL for MDI tarball")
-  set(MDI_MD5 "407db44e2d79447ab5c1233af1965f65" CACHE STRING "MD5 checksum for MDI tarball")
+  set(MDI_URL "https://github.com/MolSSI-MDI/MDI_Library/archive/v1.4.26.tar.gz" CACHE STRING "URL for MDI tarball")
+  set(MDI_MD5 "3124bb85259471e2a53a891f04bf697a" CACHE STRING "MD5 checksum for MDI tarball")
   mark_as_advanced(MDI_URL)
   mark_as_advanced(MDI_MD5)
   GetFallbackURL(MDI_URL MDI_FALLBACK)

cmake/Modules/Packages/ML-PACE.cmake

@@ -1,6 +1,6 @@
-set(PACELIB_URL "https://github.com/ICAMS/lammps-user-pace/archive/refs/tags/v.2023.01.3.fix.tar.gz" CACHE STRING "URL for PACE evaluator library sources")
+set(PACELIB_URL "https://github.com/ICAMS/lammps-user-pace/archive/refs/tags/v.2023.10.04.tar.gz" CACHE STRING "URL for PACE evaluator library sources")
 
-set(PACELIB_MD5 "4f0b3b5b14456fe9a73b447de3765caa" CACHE STRING "MD5 checksum of PACE evaluator library tarball")
+set(PACELIB_MD5 "70ff79f4e59af175e55d24f3243ad1ff" CACHE STRING "MD5 checksum of PACE evaluator library tarball")
 mark_as_advanced(PACELIB_URL)
 mark_as_advanced(PACELIB_MD5)
 GetFallbackURL(PACELIB_URL PACELIB_FALLBACK)

cmake/Modules/Packages/ML-QUIP.cmake

@@ -18,7 +18,7 @@ if(DOWNLOAD_QUIP)
     set(temp "${temp}F77FLAGS += -fpp -fixed -fPIC\n")
     set(temp "${temp}F95_PRE_FILENAME_FLAG = -Tf\n")
   elseif(CMAKE_Fortran_COMPILER_ID STREQUAL GNU)
-    set(temp "${temp}FPP=${CMAKE_Fortran_COMPILER} -E -x f95-cpp-input\nOPTIM=${CMAKE_Fortran_FLAGS_${BTYPE}}\n")
+    set(temp "${temp}FPP=${CMAKE_Fortran_COMPILER} -E -x f95-cpp-input\nOPTIM=${CMAKE_Fortran_FLAGS_${BTYPE}} -fmax-stack-var-size=6553600\n")
     set(temp "${temp}DEFINES += -DGETARG_F2003 -DGETENV_F2003 -DGFORTRAN -DFORTRAN_UNDERSCORE\n")
     set(temp "${temp}F95FLAGS += -x f95-cpp-input -ffree-line-length-none -ffree-form -fno-second-underscore -fPIC\n")
     set(temp "${temp}F77FLAGS += -x f77-cpp-input -fno-second-underscore -fPIC\n")
@@ -41,6 +41,11 @@ if(DOWNLOAD_QUIP)
   set(temp "${temp}HAVE_TURBOGAP=0\nHAVE_QR=1\nHAVE_THIRDPARTY=0\nHAVE_FX=0\nHAVE_SCME=0\nHAVE_MTP=0\n")
   set(temp "${temp}HAVE_MBD=0\nHAVE_TTM_NF=0\nHAVE_CH4=0\nHAVE_NETCDF4=0\nHAVE_MDCORE=0\nHAVE_ASAP=0\n")
   set(temp "${temp}HAVE_CGAL=0\nHAVE_METIS=0\nHAVE_LMTO_TBE=0\nHAVE_SCALAPACK=0\n")
+  # for gfortran, the -std= flag, if present, *must* be -std=gnu or else the compilation will fail.
+  if(CMAKE_Fortran_COMPILER_ID STREQUAL GNU)
+    string(REGEX REPLACE -std=f[0-9]+ -std=gnu newtemp "${temp}")
+    set(temp "${newtemp}")
+  endif()
   file(WRITE ${CMAKE_BINARY_DIR}/quip.config "${temp}")
 
   message(STATUS "QUIP download via git requested - we will build our own")
@@ -56,7 +61,7 @@ if(DOWNLOAD_QUIP)
     GIT_SUBMODULES "src/fox;src/GAP"
     PATCH_COMMAND ${CMAKE_COMMAND} -E copy_if_different ${CMAKE_BINARY_DIR}/quip.config <SOURCE_DIR>/arch/Makefile.lammps
     CONFIGURE_COMMAND env QUIP_ARCH=lammps make config
-    BUILD_COMMAND env QUIP_ARCH=lammps make libquip
+    BUILD_COMMAND env QUIP_ARCH=lammps make -j1 libquip
     INSTALL_COMMAND ""
     BUILD_IN_SOURCE YES
     BUILD_BYPRODUCTS <SOURCE_DIR>/build/lammps/${CMAKE_STATIC_LIBRARY_PREFIX}quip${CMAKE_STATIC_LIBRARY_SUFFIX}

cmake/presets/kokkos-cuda.cmake

@@ -6,6 +6,8 @@ set(Kokkos_ENABLE_SERIAL ON CACHE BOOL "" FORCE)
 set(Kokkos_ENABLE_CUDA   ON CACHE BOOL "" FORCE)
 set(Kokkos_ARCH_PASCAL60 ON CACHE BOOL "" FORCE)
 set(BUILD_OMP ON CACHE BOOL "" FORCE)
+get_filename_component(NVCC_WRAPPER_CMD ${CMAKE_CURRENT_SOURCE_DIR}/../lib/kokkos/bin/nvcc_wrapper ABSOLUTE)
+set(CMAKE_CXX_COMPILER ${NVCC_WRAPPER_CMD} CACHE FILEPATH "" FORCE)
 
 # hide deprecation warnings temporarily for stable release
 set(Kokkos_ENABLE_DEPRECATION_WARNINGS OFF CACHE BOOL "" FORCE)

cmake/presets/oneapi.cmake

@@ -18,11 +18,11 @@ set(MPI_CXX_COMPILER "mpicxx" CACHE STRING "" FORCE)
 
 unset(HAVE_OMP_H_INCLUDE CACHE)
 set(OpenMP_C "icx" CACHE STRING "" FORCE)
-set(OpenMP_C_FLAGS "-qopenmp -qopenmp-simd" CACHE STRING "" FORCE)
+set(OpenMP_C_FLAGS "-qopenmp;-qopenmp-simd" CACHE STRING "" FORCE)
 set(OpenMP_C_LIB_NAMES "omp" CACHE STRING "" FORCE)
 set(OpenMP_CXX "icpx" CACHE STRING "" FORCE)
-set(OpenMP_CXX_FLAGS "-qopenmp -qopenmp-simd" CACHE STRING "" FORCE)
+set(OpenMP_CXX_FLAGS "-qopenmp;-qopenmp-simd" CACHE STRING "" FORCE)
 set(OpenMP_CXX_LIB_NAMES "omp" CACHE STRING "" FORCE)
-set(OpenMP_Fortran_FLAGS "-qopenmp -qopenmp-simd" CACHE STRING "" FORCE)
+set(OpenMP_Fortran_FLAGS "-qopenmp;-qopenmp-simd" CACHE STRING "" FORCE)
 set(OpenMP_omp_LIBRARY "libiomp5.so" CACHE PATH "" FORCE)
 

doc/Makefile

@@ -63,6 +63,7 @@ help:
 	@echo "  anchor_check  scan for duplicate anchor labels"
 	@echo "  style_check   check for complete and consistent style lists"
 	@echo "  package_check check for complete and consistent package lists"
+	@echo "  role_check    check for misformatted role keywords"
 	@echo "  spelling      spell-check the manual"
 
 # ------------------------------------------
@@ -98,6 +99,7 @@ html: xmlgen $(VENV) $(SPHINXCONFIG)/conf.py $(ANCHORCHECK) $(MATHJAX)
 		env LC_ALL=C grep -n '[^ -~]' $(RSTDIR)/*.rst ;\
 		env LC_ALL=C grep -n ' :[a-z]\+`' $(RSTDIR)/*.rst ;\
 		env LC_ALL=C grep -n ' `[^`]\+<[a-z][^`]\+`[^_]' $(RSTDIR)/*.rst ;\
+		env LC_ALL=C grep -n ':\(ref\|doc\):[^`]' $(RSTDIR)/*.rst ;\
 		$(PYTHON) $(BUILDDIR)/utils/check-styles.py -s ../src -d src ;\
 		echo "############################################" ;\
 		deactivate ;\
@@ -179,6 +181,7 @@ pdf: xmlgen $(VENV) $(SPHINXCONFIG)/conf.py $(ANCHORCHECK)
 		env LC_ALL=C grep -n '[^ -~]' $(RSTDIR)/*.rst ;\
 		env LC_ALL=C grep -n ' :[a-z]\+`' $(RSTDIR)/*.rst ;\
 		env LC_ALL=C grep -n ' `[^`]\+<[a-z][^`]\+`[^_]' $(RSTDIR)/*.rst ;\
+		env LC_ALL=C grep -n ':\(ref\|doc\):[^`]' $(RSTDIR)/*.rst ;\
 		$(PYTHON) utils/check-styles.py -s ../src -d src ;\
 		echo "############################################" ;\
 		deactivate ;\
@@ -227,6 +230,7 @@ char_check :
 role_check :
 	@( env LC_ALL=C grep -n ' :[a-z]\+`' $(RSTDIR)/*.rst && exit 1 || : )
 	@( env LC_ALL=C grep -n ' `[^`]\+<[a-z][^`]\+`[^_]' $(RSTDIR)/*.rst && exit 1 || : )
+	@( env LC_ALL=C grep -n ':\(ref\|doc\):[^`]' $(RSTDIR)/*.rst && exit 1 || : )
 
 link_check : $(VENV) html
 	@(\

doc/src/Howto_lammps_gui.rst

@@ -90,8 +90,8 @@ The run can be stopped cleanly by using either the ``Stop LAMMPS`` entry
 in the ``Run`` menu, the hotkey `Ctrl-/` (`Command-/` on macOS), or
 clicking on the red button in the status bar.  This will cause that the
 running LAMMPS process will complete the current iteration and then
-stop.  This is equivalent to the command `timer timeout 0 <timer>` and
-implemented by calling the :cpp:func:`lammps_force_timeout()` function
+stop.  This is equivalent to the command :doc:`timer timeout 0 <timer>`
+and implemented by calling the :cpp:func:`lammps_force_timeout()` function
 of the LAMMPS C-library interface.
 
 

doc/src/compute_saed.rst

@@ -68,7 +68,7 @@ reciprocal lattice nodes. The mesh spacing is defined either (a)  by
 the entire simulation domain or (b) manually using selected values as
 shown in the 2D diagram below.
 
-.. image:: img/saed_mesh.jpg
+.. image:: img/saed_mesh.png
    :scale: 75%
    :align: center
 

doc/src/compute_xrd.rst

@@ -72,7 +72,7 @@ reciprocal lattice nodes. The mesh spacing is defined either (a) by the entire
 simulation domain or (b) manually using selected values as
 shown in the 2D diagram below.
 
-.. image:: img/xrd_mesh.jpg
+.. image:: img/xrd_mesh.png
    :scale: 75%
    :align: center
 

doc/src/fix_adapt_fep.rst

@@ -307,7 +307,9 @@ the :doc:`run <run>` command.  This fix is not invoked during
 
 Restrictions
 """"""""""""
- none
+
+The keyword "scale yes" is not supported for scaling per-atom parameters
+diameter and change. You can use :doc:`fix adapt <fix_adapt>` for those.
 
 Related commands
 """"""""""""""""

doc/src/fix_atom_swap.rst

@@ -181,6 +181,12 @@ This fix is part of the MC package.  It is only enabled if LAMMPS was
 built with that package.  See the :doc:`Build package <Build_package>`
 doc page for more info.
 
+This fix cannot be used with systems that do not have per-type masses
+(e.g. atom style sphere) since the implemented algorithm pre-computes
+velocity rescaling factors from per-type masses and ignores any per-atom
+masses, if present.  In case both, per-type and per-atom masses are
+present, a warning is printed.
+
 Related commands
 """"""""""""""""
 

doc/src/fix_ave_chunk.rst

@@ -541,10 +541,10 @@ Restrictions
 Related commands
 """"""""""""""""
 
-:doc:`compute <compute>`, :doc:`fix ave/atom <fix_ave_atom>`, `fix
-:doc:ave/histo <fix_ave_histo>`, :doc:`fix ave/time <fix_ave_time>`,
-:doc:`variable <variable>`, :doc:`fix ave/correlate
-:doc:<fix_ave_correlate>`, `fix ave/atogrid <fix_ave_grid>`
+:doc:`compute <compute>`, :doc:`fix ave/atom <fix_ave_atom>`,
+:doc:`fix ave/histo <fix_ave_histo>`, :doc:`fix ave/time <fix_ave_time>`,
+:doc:`variable <variable>`, :doc:`fix ave/correlate <fix_ave_correlate>`,
+:doc:`fix ave/grid <fix_ave_grid>`
 
 
 Default

doc/src/fix_charge_regulation.rst

@@ -253,11 +253,11 @@ built with that package.  See the :doc:`Build package <Build_package>`
 page for more info.
 
 The :doc:`atom_style <atom_style>`, used must contain the charge
-property, for example, the style could be *charge* or *full*. Only
-usable for 3D simulations. Atoms specified as free ions cannot be part
-of rigid bodies or molecules and cannot have bonding interactions. The
-scheme is limited to integer charges, any atoms with non-integer charges
-will not be considered by the fix.
+property and have per atom type masses, for example, the style could be
+*charge* or *full*. Only usable for 3D simulations.  Atoms specified as
+free ions cannot be part of rigid bodies or molecules and cannot have
+bonding interactions.  The scheme is limited to integer charges, any
+atoms with non-integer charges will not be considered by the fix.
 
 All interaction potentials used must be continuous, otherwise the MD
 integration and the particle exchange MC moves do not correspond to the

doc/src/fix_gcmc.rst

@@ -440,8 +440,11 @@ This fix is part of the MC package.  It is only enabled if LAMMPS was
 built with that package.  See the :doc:`Build package <Build_package>`
 doc page for more info.
 
+This fix style requires an :doc:`atom style <atom_style>` with per atom
+type masses.
+
 Do not set "neigh_modify once yes" or else this fix will never be
-called.  Reneighboring is required.
+called.  Reneighboring is **required**.
 
 Only usable for 3D simulations.
 

doc/src/fix_nh_uef.rst

@@ -62,7 +62,7 @@ performed using the :doc:`fix deform <fix_deform>`, :doc:`fix nvt/sllod
 <fix_nvt_sllod>`, and :doc:`compute temp/deform <compute_temp_deform>`
 commands.
 
-The applied flow field is set by the *eps* keyword. The values
+The applied flow field is set by the *erate* keyword. The values
 *edot_x* and *edot_y* correspond to the strain rates in the xx and yy
 directions.  It is implicitly assumed that the flow field is
 traceless, and therefore the strain rate in the zz direction is eqal

doc/src/fix_qeq.rst

@@ -232,8 +232,6 @@ These fixes are part of the QEQ package.  They are only enabled if
 LAMMPS was built with that package.  See the :doc:`Build package
 <Build_package>` page for more info.
 
-These qeq fixes are not compatible with the GPU and USER-INTEL packages.
-
 These qeq fixes will ignore electric field contributions from
 :doc:`fix efield <fix_efield>`.
 

doc/src/fix_sgcmc.rst

@@ -155,6 +155,9 @@ This fix is part of the MC package. It is only enabled if LAMMPS was
 built with that package.  See the :doc:`Build package <Build_package>`
 page for more info.
 
+This fix style requires an :doc:`atom style <atom_style>` with per atom
+type masses.
+
 At present the fix provides optimized subroutines for EAM type
 potentials (see above) that calculate potential energy changes due to
 *local* atom type swaps very efficiently.  Other potentials are

doc/src/fix_widom.rst

@@ -195,8 +195,11 @@ doc page for more info.
 Do not set "neigh_modify once yes" or else this fix will never be
 called.  Reneighboring is **required**.
 
-Can be run in parallel, but aspects of the GCMC part will not scale well
-in parallel. Only usable for 3D simulations.
+This fix style requires an :doc:`atom style <atom_style>` with per atom
+type masses.
+
+Can be run in parallel, but some aspects of the insertion procedure
+will not scale well in parallel. Only usable for 3D simulations.
 
 
 Related commands

doc/src/improper_amoeba.rst

@@ -68,8 +68,8 @@ for more info.
 Related commands
 """"""""""""""""
 
-:doc:`improper_coeff <improper_coeff>`, `improper_harmonic
-:doc:<improper_harmonic>`
+:doc:`improper_coeff <improper_coeff>`,
+:doc:`improper_harmonic <improper_harmonic>`
 
 Default
 """""""

doc/src/pair_beck.rst

@@ -30,11 +30,11 @@ Description
 
 Style *beck* computes interactions based on the potential by
 :ref:`(Beck) <Beck>`, originally designed for simulation of Helium.  It
-includes truncation at a cutoff distance Rc.
+includes truncation at a cutoff distance :math:`r_c`.
 
 .. math::
 
-   E(r) &= A \exp\left[-\alpha r - \beta r^6\right] - \frac{B}{\left(r^2+a^2\right)^3} \left(1+\frac{2.709+3a^2}{r^2+a^2}\right) \qquad r < R_c \\
+   E(r) &= A \exp\left[-\alpha r - \beta r^6\right] - \frac{B}{\left(r^2+a^2\right)^3} \left(1+\frac{2.709+3a^2}{r^2+a^2}\right) \qquad r < r_c \\
 
 The following coefficients must be defined for each pair of atoms
 types via the :doc:`pair_coeff <pair_coeff>` command as in the examples
@@ -50,7 +50,7 @@ commands.
 * cutoff (distance units)
 
 The last coefficient is optional.  If not specified, the global cutoff
-:math:`R_c` is used.
+:math:`r_c` is used.
 
 ----------
 

doc/src/pair_dsmc.rst

@@ -138,8 +138,12 @@ This pair style can only be used via the *pair* keyword of the
 Restrictions
 """"""""""""
 
-This style is part of the MC package.  It is only enabled if LAMMPS
-was built with that package.  See the :doc:`Build package <Build_package>` page for more info.
+This pair style is part of the MC package.  It is only enabled if LAMMPS
+was built with that package.  See the :doc:`Build package
+<Build_package>` page for more info.
+
+This pair style requires an :doc:`atom style <atom_style>` with per
+atom type masses.
 
 Related commands
 """"""""""""""""

doc/src/pair_lj_smooth_linear.rst

@@ -31,13 +31,13 @@ Style *lj/smooth/linear* computes a truncated and force-shifted LJ
 interaction (aka Shifted Force Lennard-Jones) that combines the
 standard 12/6 Lennard-Jones function and subtracts a linear term based
 on the cutoff distance, so that both, the potential and the force, go
-continuously to zero at the cutoff Rc :ref:`(Toxvaerd) <Toxvaerd>`:
+continuously to zero at the cutoff :math:`r_c` :ref:`(Toxvaerd) <Toxvaerd>`:
 
 .. math::
 
    \phi\left(r\right) & =  4 \epsilon \left[ \left(\frac{\sigma}{r}\right)^{12} -
                        \left(\frac{\sigma}{r}\right)^6 \right] \\
-   E\left(r\right) & =  \phi\left(r\right)  - \phi\left(R_c\right) - \left(r - R_c\right) \left.\frac{d\phi}{d r} \right|_{r=R_c}       \qquad r < R_c
+   E\left(r\right) & =  \phi\left(r\right)  - \phi\left(r_c\right) - \left(r - r_c\right) \left.\frac{d\phi}{d r} \right|_{r=r_c}       \qquad r < r_c
 
 The following coefficients must be defined for each pair of atoms
 types via the :doc:`pair_coeff <pair_coeff>` command as in the examples
@@ -77,8 +77,9 @@ tail option for adding long-range tail corrections to energy and
 pressure, since the energy of the pair interaction is smoothed to 0.0
 at the cutoff.
 
-This pair style writes its information to :doc:`binary restart files <restart>`, so pair_style and pair_coeff commands do not need
-to be specified in an input script that reads a restart file.
+This pair style writes its information to :doc:`binary restart files <restart>`,
+so pair_style and pair_coeff commands do not need to be specified
+in an input script that reads a restart file.
 
 This pair style can only be used via the *pair* keyword of the
 :doc:`run_style respa <run_style>` command.  It does not support the

doc/src/pair_mie.rst

@@ -35,7 +35,7 @@ The *mie/cut* style computes the Mie potential, given by
    E =  C \epsilon \left[ \left(\frac{\sigma}{r}\right)^{\gamma_{rep}} - \left(\frac{\sigma}{r}\right)^{\gamma_{att}} \right]
                          \qquad r < r_c
 
-Rc is the cutoff and C is a function that depends on the repulsive and
+:math:`r_c` is the cutoff and C is a function that depends on the repulsive and
 attractive exponents, given by:
 
 .. math::

doc/src/pair_morse.rst

@@ -53,7 +53,7 @@ Style *morse* computes pairwise interactions with the formula
    E = D_0 \left[ e^{- 2 \alpha (r - r_0)} - 2 e^{- \alpha (r - r_0)} \right]
        \qquad r < r_c
 
-Rc is the cutoff.
+:math:`r_c` is the cutoff.
 
 The following coefficients must be defined for each pair of atoms
 types via the :doc:`pair_coeff <pair_coeff>` command as in the examples
@@ -78,7 +78,7 @@ so that both, potential energy and force, go to zero at the cut-off:
 .. math::
 
    \phi\left(r\right) & =  D_0 \left[ e^{- 2 \alpha (r - r_0)} - 2 e^{- \alpha (r - r_0)} \right] \qquad r < r_c \\
-   E\left(r\right) & =  \phi\left(r\right)  - \phi\left(R_c\right) - \left(r - R_c\right) \left.\frac{d\phi}{d r} \right|_{r=R_c}       \qquad r < R_c
+   E\left(r\right) & =  \phi\left(r\right)  - \phi\left(r_c\right) - \left(r - r_c\right) \left.\frac{d\phi}{d r} \right|_{r=r_c}       \qquad r < r_c
 
 The syntax of the pair_style and pair_coeff commands are the same for
 the *morse* and *morse/smooth/linear* styles.

doc/src/pair_soft.rst

@@ -44,8 +44,9 @@ It is useful for pushing apart overlapping atoms, since it does not
 blow up as r goes to 0.  A is a prefactor that can be made to vary in
 time from the start to the end of the run (see discussion below),
 e.g. to start with a very soft potential and slowly harden the
-interactions over time.  Rc is the cutoff.  See the :doc:`fix nve/limit <fix_nve_limit>` command for another way to push apart
-overlapping atoms.
+interactions over time.  :math:`r_c` is the cutoff.
+See the :doc:`fix nve/limit <fix_nve_limit>` command for another way
+to push apart overlapping atoms.
 
 The following coefficients must be defined for each pair of atom types
 via the :doc:`pair_coeff <pair_coeff>` command as in the examples above,

doc/src/pair_spica.rst

@@ -81,7 +81,7 @@ given by
 
 as required for the SPICA (formerly called SDK) and the pSPICA Coarse-grained MD parameterization discussed in
 :ref:`(Shinoda) <Shinoda3>`, :ref:`(DeVane) <DeVane>`, :ref:`(Seo) <Seo>`, and :ref:`(Miyazaki) <Miyazaki>`.
-Rc is the cutoff.
+:math:`r_c` is the cutoff.
 Summary information on these force fields can be found at https://www.spica-ff.org
 
 Style *lj/spica/coul/long* computes the adds Coulombic interactions

doc/utils/sphinx-config/LAMMPSLexer.py

@@ -76,12 +76,12 @@ class LAMMPSLexer(RegexLexer):
             include('conditionals'),
             include('keywords'),
             (r'#.*?\n', Comment),
-            ('"', String, 'string'),
-            ('\'', String, 'single_quote_string'),
+            (r'"', String, 'string'),
+            (r'\'', String, 'single_quote_string'),
             (r'[0-9]+:[0-9]+(:[0-9]+)?', Number),
             (r'[0-9]+(\.[0-9]+)?([eE]\-?[0-9]+)?', Number),
-            ('\$?\(', Name.Variable, 'expression'),
-            ('\$\{', Name.Variable, 'variable'),
+            (r'\$?\(', Name.Variable, 'expression'),
+            (r'\$\{', Name.Variable, 'variable'),
             (r'[\w_\.\[\]]+', Name),
             (r'\$[\w_]+', Name.Variable),
             (r'\s+', Whitespace),
@@ -97,21 +97,21 @@ class LAMMPSLexer(RegexLexer):
         ]
         ,
         'variable' : [
-            ('[^\}]+', Name.Variable),
-            ('\}', Name.Variable, '#pop'),
+            (r'[^\}]+', Name.Variable),
+            (r'\}', Name.Variable, '#pop'),
         ],
         'string' : [
-            ('[^"]+', String),
-            ('"', String, '#pop'),
+            (r'[^"]+', String),
+            (r'"', String, '#pop'),
         ],
         'single_quote_string' : [
-            ('[^\']+', String),
-            ('\'', String, '#pop'),
+            (r'[^\']+', String),
+            (r'\'', String, '#pop'),
         ],
         'expression' : [
-            ('[^\(\)]+', Name.Variable),
-            ('\(', Name.Variable, 'expression'),
-            ('\)', Name.Variable, '#pop'),
+            (r'[^\(\)]+', Name.Variable),
+            (r'\(', Name.Variable, 'expression'),
+            (r'\)', Name.Variable, '#pop'),
         ],
         'modify_cmd' : [
             (r'[\w_\-\.\[\]]+', Name.Variable.Identifier),

examples/PACKAGES/cgdna/util/generate.py

@@ -22,22 +22,26 @@
 """
 Import basic modules
 """
+
+# for python2/3 compatibility
+from __future__ import print_function
+
 import sys, os, timeit
 
 from timeit import default_timer as timer
 start_time = timer()
 """
-Try to import numpy; if failed, import a local version mynumpy 
+Try to import numpy; if failed, import a local version mynumpy
 which needs to be provided
 """
 try:
     import numpy as np
 except:
-    print >> sys.stderr, "numpy not found. Exiting."
+    print("numpy not found. Exiting.", file=sys.stderr)
     sys.exit(1)
 
 """
-Check that the required arguments (box offset and size in simulation units 
+Check that the required arguments (box offset and size in simulation units
 and the sequence file were provided
 """
 try:
@@ -45,8 +49,8 @@ try:
     box_length = float(sys.argv[2])
     infile = sys.argv[3]
 except:
-    print >> sys.stderr, "Usage: %s <%s> <%s> <%s>" % (sys.argv[0], \
-	"box offset", "box length", "file with sequences")
+    print( "Usage: %s <%s> <%s> <%s>" % (sys.argv[0], \
+	"box offset", "box length", "file with sequences"), file=sys.stderr)
     sys.exit(1)
 box = np.array ([box_length, box_length, box_length])
 
@@ -57,8 +61,7 @@ try:
     inp = open (infile, 'r')
     inp.close()
 except:
-    print >> sys.stderr, "Could not open file '%s' for reading. \
-					      Aborting." % infile
+    print( "Could not open file '%s' for reading. Aborting." % infile, file=sys.stderr)
     sys.exit(2)
 
 # return parts of a string
@@ -86,7 +89,7 @@ Define auxiliary variables for the construction of a helix
 # center of the double strand
 CM_CENTER_DS = POS_BASE + 0.2
 
-# ideal distance between base sites of two nucleotides 
+# ideal distance between base sites of two nucleotides
 # which are to be base paired in a duplex
 BASE_BASE = 0.3897628551303122
 
@@ -118,7 +121,7 @@ strandnum = []
 
 bonds = []
 
-""" 
+"""
 Convert local body frame to quaternion DOF
 """
 def exyz_to_quat (mya1, mya3):
@@ -135,25 +138,25 @@ def exyz_to_quat (mya1, mya3):
     # compute other components from it
 
     if q0sq >= 0.25:
-	myquat[0] = np.sqrt(q0sq)
-	myquat[1] = (mya2[2] - mya3[1]) / (4.0*myquat[0])
-	myquat[2] = (mya3[0] - mya1[2]) / (4.0*myquat[0])
-	myquat[3] = (mya1[1] - mya2[0]) / (4.0*myquat[0])
+        myquat[0] = np.sqrt(q0sq)
+        myquat[1] = (mya2[2] - mya3[1]) / (4.0*myquat[0])
+        myquat[2] = (mya3[0] - mya1[2]) / (4.0*myquat[0])
+        myquat[3] = (mya1[1] - mya2[0]) / (4.0*myquat[0])
     elif q1sq >= 0.25:
-	myquat[1] = np.sqrt(q1sq)
-	myquat[0] = (mya2[2] - mya3[1]) / (4.0*myquat[1])
-	myquat[2] = (mya2[0] + mya1[1]) / (4.0*myquat[1])
-	myquat[3] = (mya1[2] + mya3[0]) / (4.0*myquat[1])
+        myquat[1] = np.sqrt(q1sq)
+        myquat[0] = (mya2[2] - mya3[1]) / (4.0*myquat[1])
+        myquat[2] = (mya2[0] + mya1[1]) / (4.0*myquat[1])
+        myquat[3] = (mya1[2] + mya3[0]) / (4.0*myquat[1])
     elif q2sq >= 0.25:
-	myquat[2] = np.sqrt(q2sq)
-	myquat[0] = (mya3[0] - mya1[2]) / (4.0*myquat[2])
-	myquat[1] = (mya2[0] + mya1[1]) / (4.0*myquat[2])
-	myquat[3] = (mya3[1] + mya2[2]) / (4.0*myquat[2])
+        myquat[2] = np.sqrt(q2sq)
+        myquat[0] = (mya3[0] - mya1[2]) / (4.0*myquat[2])
+        myquat[1] = (mya2[0] + mya1[1]) / (4.0*myquat[2])
+        myquat[3] = (mya3[1] + mya2[2]) / (4.0*myquat[2])
     elif q3sq >= 0.25:
-	myquat[3] = np.sqrt(q3sq)
-	myquat[0] = (mya1[1] - mya2[0]) / (4.0*myquat[3])
-	myquat[1] = (mya3[0] + mya1[2]) / (4.0*myquat[3])
-	myquat[2] = (mya3[1] + mya2[2]) / (4.0*myquat[3])
+        myquat[3] = np.sqrt(q3sq)
+        myquat[0] = (mya1[1] - mya2[0]) / (4.0*myquat[3])
+        myquat[1] = (mya3[0] + mya1[2]) / (4.0*myquat[3])
+        myquat[2] = (mya3[1] + mya2[2]) / (4.0*myquat[3])
 
     norm = 1.0/np.sqrt(myquat[0]*myquat[0] + myquat[1]*myquat[1] + \
 			  myquat[2]*myquat[2] + myquat[3]*myquat[3])
@@ -169,62 +172,62 @@ Adds a strand to the system by appending it to the array of previous strands
 """
 def add_strands (mynewpositions, mynewa1s, mynewa3s):
     overlap = False
-	
-    # This is a simple check for each of the particles where for previously 
-    # placed particles i we check whether it overlaps with any of the 
+
+    # This is a simple check for each of the particles where for previously
+    # placed particles i we check whether it overlaps with any of the
     # newly created particles j
 
-    print >> sys.stdout, "## Checking for overlaps"
+    print( "## Checking for overlaps", file=sys.stdout)
 
-    for i in xrange(len(positions)):
+    for i in range(len(positions)):
 
-	p = positions[i]
-	pa1 = a1s[i]
+        p = positions[i]
+        pa1 = a1s[i]
 
-	for j in xrange (len(mynewpositions)):
+        for j in range (len(mynewpositions)):
 
-	    q = mynewpositions[j]
-	    qa1 = mynewa1s[j]
+            q = mynewpositions[j]
+            qa1 = mynewa1s[j]
 
-	    # skip particles that are anyway too far away
-	    dr = p - q
-	    dr -= box * np.rint (dr / box)
-	    if np.dot(dr, dr) > RC2:
-		continue
+            # skip particles that are anyway too far away
+            dr = p - q
+            dr -= box * np.rint(dr / box)
+            if np.dot(dr, dr) > RC2:
+                continue
 
-	    # base site and backbone site of the two particles
+            # base site and backbone site of the two particles
             p_pos_back = p + pa1 * POS_BACK
             p_pos_base = p + pa1 * POS_BASE
             q_pos_back = q + qa1 * POS_BACK
             q_pos_base = q + qa1 * POS_BASE
 
-	    # check for no overlap between the two backbone sites
+            # check for no overlap between the two backbone sites
             dr = p_pos_back - q_pos_back
-            dr -= box * np.rint (dr / box)
+            dr -= box * np.rint(dr / box)
             if np.dot(dr, dr) < RC2_BACK:
                 overlap = True
 
-	    # check for no overlap between the two base sites
+            # check for no overlap between the two base sites
             dr = p_pos_base -  q_pos_base
-            dr -= box * np.rint (dr / box)
+            dr -= box * np.rint(dr / box)
             if np.dot(dr, dr) < RC2_BASE:
                 overlap = True
 
-	    # check for no overlap between backbone site of particle p 
-	    # with base site of particle q
+            # check for no overlap between backbone site of particle p
+            # with base site of particle q
             dr = p_pos_back - q_pos_base
             dr -= box * np.rint (dr / box)
             if np.dot(dr, dr) < RC2_BACK_BASE:
                 overlap = True
 
-	    # check for no overlap between base site of particle p and 
-	    # backbone site of particle q
+            # check for no overlap between base site of particle p and
+            # backbone site of particle q
             dr = p_pos_base - q_pos_back
             dr -= box * np.rint (dr / box)
             if np.dot(dr, dr) < RC2_BACK_BASE:
                 overlap = True
 
-	    # exit if there is an overlap
+            # exit if there is an overlap
             if overlap:
                 return False
 
@@ -237,10 +240,10 @@ def add_strands (mynewpositions, mynewa1s, mynewa3s):
             a1s.append (p)
         for p in mynewa3s:
             a3s.append (p)
-	# calculate quaternion from local body frame and append
-	for ia in xrange(len(mynewpositions)):
-	    mynewquaternions = exyz_to_quat(mynewa1s[ia],mynewa3s[ia])
-	    quaternions.append(mynewquaternions)
+        # calculate quaternion from local body frame and append
+        for ia in range(len(mynewpositions)):
+            mynewquaternions = exyz_to_quat(mynewa1s[ia],mynewa3s[ia])
+            quaternions.append(mynewquaternions)
 
     return True
 
@@ -281,7 +284,7 @@ def get_rotation_matrix(axis, anglest):
                     [olc*x*z-st*y, olc*y*z+st*x, olc*z*z+ct]])
 
 """
-Generates the position and orientation vectors of a 
+Generates the position and orientation vectors of a
 (single or double) strand from a sequence string
 """
 def generate_strand(bp, sequence=None, start_pos=np.array([0, 0, 0]), \
@@ -295,76 +298,75 @@ def generate_strand(bp, sequence=None, start_pos=np.array([0, 0, 0]), \
     # overall direction of the helix
     dir = np.array(dir, dtype=float)
     if sequence == None:
-	sequence = np.random.randint(1, 5, bp)
+        sequence = np.random.randint(1, 5, bp)
 
-    # the elseif here is most likely redundant 
+    # the elseif here is most likely redundant
     elif len(sequence) != bp:
-	n = bp - len(sequence)
-	sequence += np.random.randint(1, 5, n)
-	print >> sys.stderr, "sequence is too short, adding %d random bases" % n
+        n = bp - len(sequence)
+        sequence += np.random.randint(1, 5, n)
+        print( "sequence is too short, adding %d random bases" % n, file=sys.stderr)
 
     # normalize direction
     dir_norm = np.sqrt(np.dot(dir,dir))
     if dir_norm < 1e-10:
-	print >> sys.stderr, "direction must be a valid vector, \
-			      defaulting to (0, 0, 1)"
-	dir = np.array([0, 0, 1])
+        print( "direction must be a valid vector, defaulting to (0, 0, 1)", file=sys.stderr)
+        dir = np.array([0, 0, 1])
     else: dir /= dir_norm
 
     # find a vector orthogonal to dir to act as helix direction,
     # if not provided switch off random orientation
     if perp is None or perp is False:
-	v1 = np.random.random_sample(3)
-	v1 -= dir * (np.dot(dir, v1))
-	v1 /= np.sqrt(sum(v1*v1))
+        v1 = np.random.random_sample(3)
+        v1 -= dir * (np.dot(dir, v1))
+        v1 /= np.sqrt(sum(v1*v1))
     else:
-	v1 = perp;
+        v1 = perp;
 
     # generate rotational matrix representing the overall rotation of the helix
     R0 = get_rotation_matrix(dir, rot)
-	    
+
     # rotation matrix corresponding to one step along the helix
     R = get_rotation_matrix(dir, [1, "bp"])
 
-    # set the vector a1 (backbone to base) to v1 
+    # set the vector a1 (backbone to base) to v1
     a1 = v1
-    
-    # apply the global rotation to a1 
+
+    # apply the global rotation to a1
     a1 = np.dot(R0, a1)
-    
+
     # set the position of the fist backbone site to start_pos
     rb = np.array(start_pos)
-	    
+
     # set a3 to the direction of the helix
     a3 = dir
     for i in range(bp):
     # work out the position of the centre of mass of the nucleotide
-	rcdm = rb - CM_CENTER_DS * a1
-	
-	# append to newpositions
-	mynewpositions.append(rcdm)
-	mynewa1s.append(a1)
-	mynewa3s.append(a3)
-	
-	# if we are not at the end of the helix, we work out a1 and rb for the 
-	# next nucleotide along the helix
-	if i != bp - 1:
-	    a1 = np.dot(R, a1)
-	    rb += a3 * BASE_BASE
+        rcdm = rb - CM_CENTER_DS * a1
 
-    # if we are working on a double strand, we do a cycle similar 
+        # append to newpositions
+        mynewpositions.append(rcdm)
+        mynewa1s.append(a1)
+        mynewa3s.append(a3)
+
+        # if we are not at the end of the helix, we work out a1 and rb for the
+        # next nucleotide along the helix
+        if i != bp - 1:
+            a1 = np.dot(R, a1)
+            rb += a3 * BASE_BASE
+
+    # if we are working on a double strand, we do a cycle similar
     # to the previous one but backwards
     if double == True:
-	a1 = -a1
-	a3 = -dir
-	R = R.transpose()
-	for i in range(bp):
-	    rcdm = rb - CM_CENTER_DS * a1
-	    mynewpositions.append (rcdm)
-	    mynewa1s.append (a1)
-	    mynewa3s.append (a3)
-	    a1 = np.dot(R, a1)
-	    rb += a3 * BASE_BASE
+        a1 = -a1
+        a3 = -dir
+        R = R.transpose()
+        for i in range(bp):
+            rcdm = rb - CM_CENTER_DS * a1
+            mynewpositions.append (rcdm)
+            mynewa1s.append (a1)
+            mynewa3s.append (a3)
+            a1 = np.dot(R, a1)
+            rb += a3 * BASE_BASE
 
     assert (len (mynewpositions) > 0)
 
@@ -391,10 +393,10 @@ def read_strands(filename):
     try:
         infile = open (filename)
     except:
-        print >> sys.stderr, "Could not open file '%s'. Aborting." % filename
+        print( "Could not open file '%s'. Aborting." % filename, file=sys.stderr )
         sys.exit(2)
 
-    # This block works out the number of nucleotides and strands by reading 
+    # This block works out the number of nucleotides and strands by reading
     # the number of non-empty lines in the input file and the number of letters,
     # taking the possible DOUBLE keyword into account.
     nstrands, nnucl, nbonds = 0, 0, 0
@@ -406,30 +408,29 @@ def read_strands(filename):
         if line[:6] == 'DOUBLE':
             line = line.split()[1]
             length = len(line)
-            print >> sys.stdout, "## Found duplex of %i base pairs" % length
+            print( "## Found duplex of %i base pairs" % length, file=sys.stdout)
             nnucl += 2*length
             nstrands += 2
-	    nbonds += (2*length-2)
+            nbonds += (2*length-2)
         else:
             line = line.split()[0]
             length = len(line)
-            print >> sys.stdout, \
-		    "## Found single strand of %i bases" % length
+            print( "## Found single strand of %i bases" % length, file=sys.stdout)
             nnucl += length
             nstrands += 1
-	    nbonds += length-1
+            nbonds += length-1
     # rewind the sequence input file
     infile.seek(0)
 
-    print >> sys.stdout, "## nstrands, nnucl = ", nstrands, nnucl
+    print( "## nstrands, nnucl = ", nstrands, nnucl, file=sys.stdout)
 
     # generate the data file in LAMMPS format
     try:
         out = open ("data.oxdna", "w")
     except:
-        print >> sys.stderr, "Could not open data file for writing. Aborting."
+        print( "Could not open data file for writing. Aborting.", file=sys.stderr)
         sys.exit(2)
-	
+
     lines = infile.readlines()
     nlines = len(lines)
     i = 1
@@ -440,115 +441,114 @@ def read_strands(filename):
         line = line.upper().strip()
 
         # skip empty lines
-        if len(line) == 0: 
-	    i += 1
-	    continue
+        if len(line) == 0:
+            i += 1
+            continue
 
-	# block for duplexes: last argument of the generate function 
-	# is set to 'True'
+        # block for duplexes: last argument of the generate function
+        # is set to 'True'
         if line[:6] == 'DOUBLE':
             line = line.split()[1]
             length = len(line)
             seq = [(base_to_number[x]) for x in line]
 
-	    myns += 1
-	    for b in xrange(length):
-		basetype.append(seq[b])
-		strandnum.append(myns)
+            myns += 1
+            for b in range(length):
+                basetype.append(seq[b])
+                strandnum.append(myns)
 
-	    for b in xrange(length-1):
-		bondpair = [noffset + b, noffset + b + 1]
-		bonds.append(bondpair)
-	    noffset += length
+            for b in range(length-1):
+                bondpair = [noffset + b, noffset + b + 1]
+                bonds.append(bondpair)
+            noffset += length
 
-	    # create the sequence of the second strand as made of 
-	    # complementary bases
-	    seq2 = [5-s for s in seq]
-	    seq2.reverse()
+            # create the sequence of the second strand as made of
+            # complementary bases
+            seq2 = [5-s for s in seq]
+            seq2.reverse()
 
-	    myns += 1
-	    for b in xrange(length):
-		basetype.append(seq2[b])
-		strandnum.append(myns)
+            myns += 1
+            for b in range(length):
+                basetype.append(seq2[b])
+                strandnum.append(myns)
 
-	    for b in xrange(length-1):
-		bondpair = [noffset + b, noffset + b + 1]
-		bonds.append(bondpair)
-	    noffset += length
- 
-            print >> sys.stdout, "## Created duplex of %i bases" % (2*length)
+            for b in range(length-1):
+                bondpair = [noffset + b, noffset + b + 1]
+                bonds.append(bondpair)
+            noffset += length
 
-	    # generate random position of the first nucleotide
+            print( "## Created duplex of %i bases" % (2*length), file=sys.stdout)
+
+            # generate random position of the first nucleotide
             cdm = box_offset + np.random.random_sample(3) * box
 
-            # generate the random direction of the helix 
+            # generate the random direction of the helix
             axis = np.random.random_sample(3)
             axis /= np.sqrt(np.dot(axis, axis))
 
-            # use the generate function defined above to create 
-	    # the position and orientation vector of the strand 
+            # use the generate function defined above to create
+            # the position and orientation vector of the strand
             newpositions, newa1s, newa3s = generate_strand(len(line), \
-		    sequence=seq, dir=axis, start_pos=cdm, double=True)
+                sequence=seq, dir=axis, start_pos=cdm, double=True)
 
             # generate a new position for the strand until it does not overlap
-	    # with anything already present
-	    start = timer()
+            # with anything already present
+            start = timer()
             while not add_strands(newpositions, newa1s, newa3s):
                 cdm = box_offset + np.random.random_sample(3) * box
                 axis = np.random.random_sample(3)
                 axis /= np.sqrt(np.dot(axis, axis))
                 newpositions, newa1s, newa3s = generate_strand(len(line), \
-		      sequence=seq, dir=axis, start_pos=cdm, double=True)
-                print >> sys.stdout, "## Trying %i" % i
-	    end = timer()
-            print >> sys.stdout, "## Added duplex of %i bases (line %i/%i) in %.2fs, now at %i/%i" % \
-				      (2*length, i, nlines, end-start, len(positions), nnucl)
+                    sequence=seq, dir=axis, start_pos=cdm, double=True)
+                print( "## Trying %i" % i, file=sys.stdout)
+            end = timer()
+            print( "## Added duplex of %i bases (line %i/%i) in %.2fs, now at %i/%i" % \
+				      (2*length, i, nlines, end-start, len(positions), nnucl), file=sys.stdout)
 
-	# block for single strands: last argument of the generate function 
-	# is set to 'False'
+        # block for single strands: last argument of the generate function
+        # is set to 'False'
         else:
             length = len(line)
             seq = [(base_to_number[x]) for x in line]
 
-	    myns += 1
-	    for b in xrange(length):
-		basetype.append(seq[b])
-		strandnum.append(myns)
+            myns += 1
+            for b in range(length):
+                basetype.append(seq[b])
+                strandnum.append(myns)
 
-	    for b in xrange(length-1):
-		bondpair = [noffset + b, noffset + b + 1]
-		bonds.append(bondpair)
-	    noffset += length
+            for b in range(length-1):
+                bondpair = [noffset + b, noffset + b + 1]
+                bonds.append(bondpair)
+            noffset += length
 
-	    # generate random position of the first nucleotide
+            # generate random position of the first nucleotide
             cdm = box_offset + np.random.random_sample(3) * box
 
-            # generate the random direction of the helix 
+            # generate the random direction of the helix
             axis = np.random.random_sample(3)
             axis /= np.sqrt(np.dot(axis, axis))
 
-            print >> sys.stdout, \
-		      "## Created single strand of %i bases" % length
+            print("## Created single strand of %i bases" % length, file=sys.stdout)
 
             newpositions, newa1s, newa3s = generate_strand(length, \
 		      sequence=seq, dir=axis, start_pos=cdm, double=False)
-	    start = timer()
+            start = timer()
             while not add_strands(newpositions, newa1s, newa3s):
                 cdm = box_offset + np.random.random_sample(3) * box
                 axis = np.random.random_sample(3)
-		axis /= np.sqrt(np.dot(axis, axis))
+                axis /= np.sqrt(np.dot(axis, axis))
                 newpositions, newa1s, newa3s = generate_strand(length, \
-			  sequence=seq, dir=axis, start_pos=cdm, double=False)
+                    sequence=seq, dir=axis, start_pos=cdm, double=False)
                 print >> sys.stdout, "## Trying  %i" % (i)
-	    end = timer()
-            print >> sys.stdout, "## Added single strand of %i bases (line %i/%i) in %.2fs, now at %i/%i" % \
-				      (length, i, nlines, end-start,len(positions), nnucl)
+            end = timer()
+            print( "## Added single strand of %i bases (line %i/%i) in %.2fs, now at %i/%i" % \
+				      (length, i, nlines, end-start,len(positions), nnucl), file=sys.stdout)
 
         i += 1
 
     # sanity check
     if not len(positions) == nnucl:
-        print len(positions), nnucl
+        print( len(positions), nnucl )
         raise AssertionError
 
     out.write('# LAMMPS data file\n')
@@ -580,44 +580,41 @@ def read_strands(filename):
     out.write('Atoms\n')
     out.write('\n')
 
-    for i in xrange(nnucl):
-	out.write('%d %d %22.15le %22.15le %22.15le %d 1 1\n' \
-		  % (i+1, basetype[i], \
-		     positions[i][0], positions[i][1], positions[i][2], \
-		     strandnum[i]))
+    for i in range(nnucl):
+        out.write('%d %d %22.15le %22.15le %22.15le %d 1 1\n' \
+            % (i+1, basetype[i], positions[i][0], positions[i][1], positions[i][2], strandnum[i]))
 
     out.write('\n')
     out.write('# Atom-ID, translational, rotational velocity\n')
     out.write('Velocities\n')
     out.write('\n')
 
-    for i in xrange(nnucl):
-	out.write("%d %22.15le %22.15le %22.15le %22.15le %22.15le %22.15le\n" \
-		  % (i+1,0.0,0.0,0.0,0.0,0.0,0.0))
+    for i in range(nnucl):
+        out.write("%d %22.15le %22.15le %22.15le %22.15le %22.15le %22.15le\n" \
+            % (i+1,0.0,0.0,0.0,0.0,0.0,0.0))
 
     out.write('\n')
     out.write('# Atom-ID, shape, quaternion\n')
     out.write('Ellipsoids\n')
     out.write('\n')
 
-    for i in xrange(nnucl):
-	out.write(\
-    "%d %22.15le %22.15le %22.15le %22.15le %22.15le %22.15le %22.15le\n"  \
-      % (i+1,1.1739845031423408,1.1739845031423408,1.1739845031423408, \
-	quaternions[i][0],quaternions[i][1], quaternions[i][2],quaternions[i][3]))
- 
+    for i in range(nnucl):
+        out.write("%d %22.15le %22.15le %22.15le %22.15le %22.15le %22.15le %22.15le\n"  \
+            % (i+1,1.1739845031423408,1.1739845031423408,1.1739845031423408, \
+            quaternions[i][0],quaternions[i][1], quaternions[i][2],quaternions[i][3]))
+
     out.write('\n')
     out.write('# Bond topology\n')
     out.write('Bonds\n')
     out.write('\n')
 
-    for i in xrange(nbonds):
-	out.write("%d  %d  %d  %d\n" % (i+1,1,bonds[i][0],bonds[i][1]))
+    for i in range(nbonds):
+        out.write("%d  %d  %d  %d\n" % (i+1,1,bonds[i][0],bonds[i][1]))
 
     out.close()
 
-    print >> sys.stdout, "## Wrote data to 'data.oxdna'"
-    print >> sys.stdout, "## DONE"
+    print("## Wrote data to 'data.oxdna'", file=sys.stdout)
+    print("## DONE", file=sys.stdout)
 
 # call the above main() function, which executes the program
 read_strands (infile)
@@ -627,4 +624,6 @@ runtime = end_time-start_time
 hours = runtime/3600
 minutes = (runtime-np.rint(hours)*3600)/60
 seconds = (runtime-np.rint(hours)*3600-np.rint(minutes)*60)%60
-print >> sys.stdout, "## Total runtime %ih:%im:%.2fs" % (hours,minutes,seconds)
+print( "## Total runtime %ih:%im:%.2fs" % (hours,minutes,seconds), file=sys.stdout)
+
+

examples/PACKAGES/cgdna/util/generate_simple.py

@@ -1,5 +1,8 @@
 # Setup tool for oxDNA input in LAMMPS format.
 
+# for python2/3 compatibility
+from __future__ import print_function
+
 import math,numpy as np,sys,os
 
 # system size
@@ -250,59 +253,59 @@ def duplex_array():
       qrot3=math.sin(0.5*twist)
 
       for letter in strand[2]:
-	temp1=[]
-	temp2=[]
+        temp1=[]
+        temp2=[]
 
-	temp1.append(nt2num[letter])
-	temp2.append(compnt2num[letter])
+        temp1.append(nt2num[letter])
+        temp2.append(compnt2num[letter])
 
-	temp1.append([posx1,posy1,posz1])
-	temp2.append([posx2,posy2,posz2])
+        temp1.append([posx1,posy1,posz1])
+        temp2.append([posx2,posy2,posz2])
 
-	vel=[0,0,0,0,0,0]
-	temp1.append(vel)
-	temp2.append(vel)
+        vel=[0,0,0,0,0,0]
+        temp1.append(vel)
+        temp2.append(vel)
 
-	temp1.append(shape)
-	temp2.append(shape)
+        temp1.append(shape)
+        temp2.append(shape)
 
-	temp1.append(quat1)
-	temp2.append(quat2)
+        temp1.append(quat1)
+        temp2.append(quat2)
 
-	quat1_0 = quat1[0]*qrot0 - quat1[1]*qrot1 - quat1[2]*qrot2 - quat1[3]*qrot3 
-	quat1_1 = quat1[0]*qrot1 + quat1[1]*qrot0 + quat1[2]*qrot3 - quat1[3]*qrot2 
-	quat1_2 = quat1[0]*qrot2 + quat1[2]*qrot0 + quat1[3]*qrot1 - quat1[1]*qrot3 
-	quat1_3 = quat1[0]*qrot3 + quat1[3]*qrot0 + quat1[1]*qrot2 + quat1[2]*qrot1 
+        quat1_0 = quat1[0]*qrot0 - quat1[1]*qrot1 - quat1[2]*qrot2 - quat1[3]*qrot3
+        quat1_1 = quat1[0]*qrot1 + quat1[1]*qrot0 + quat1[2]*qrot3 - quat1[3]*qrot2
+        quat1_2 = quat1[0]*qrot2 + quat1[2]*qrot0 + quat1[3]*qrot1 - quat1[1]*qrot3
+        quat1_3 = quat1[0]*qrot3 + quat1[3]*qrot0 + quat1[1]*qrot2 + quat1[2]*qrot1
 
-	quat1 = [quat1_0,quat1_1,quat1_2,quat1_3]
+        quat1 = [quat1_0,quat1_1,quat1_2,quat1_3]
 
-	posx1=axisx - dcomh*(quat1[0]**2+quat1[1]**2-quat1[2]**2-quat1[3]**2)
-	posy1=axisy - dcomh*(2*(quat1[1]*quat1[2]+quat1[0]*quat1[3]))
-	posz1=posz1+risez
+        posx1=axisx - dcomh*(quat1[0]**2+quat1[1]**2-quat1[2]**2-quat1[3]**2)
+        posy1=axisy - dcomh*(2*(quat1[1]*quat1[2]+quat1[0]*quat1[3]))
+        posz1=posz1+risez
 
-	quat2_0 = quat2[0]*qrot0 - quat2[1]*qrot1 - quat2[2]*qrot2 + quat2[3]*qrot3 
-	quat2_1 = quat2[0]*qrot1 + quat2[1]*qrot0 - quat2[2]*qrot3 - quat2[3]*qrot2 
-	quat2_2 = quat2[0]*qrot2 + quat2[2]*qrot0 + quat2[3]*qrot1 + quat2[1]*qrot3 
-	quat2_3 =-quat2[0]*qrot3 + quat2[3]*qrot0 + quat2[1]*qrot2 + quat2[2]*qrot1 
+        quat2_0 = quat2[0]*qrot0 - quat2[1]*qrot1 - quat2[2]*qrot2 + quat2[3]*qrot3
+        quat2_1 = quat2[0]*qrot1 + quat2[1]*qrot0 - quat2[2]*qrot3 - quat2[3]*qrot2
+        quat2_2 = quat2[0]*qrot2 + quat2[2]*qrot0 + quat2[3]*qrot1 + quat2[1]*qrot3
+        quat2_3 =-quat2[0]*qrot3 + quat2[3]*qrot0 + quat2[1]*qrot2 + quat2[2]*qrot1
 
-	quat2 = [quat2_0,quat2_1,quat2_2,quat2_3]
+        quat2 = [quat2_0,quat2_1,quat2_2,quat2_3]
 
-	posx2=axisx + dcomh*(quat1[0]**2+quat1[1]**2-quat1[2]**2-quat1[3]**2)
-	posy2=axisy + dcomh*(2*(quat1[1]*quat1[2]+quat1[0]*quat1[3]))
-	posz2=posz1
+        posx2=axisx + dcomh*(quat1[0]**2+quat1[1]**2-quat1[2]**2-quat1[3]**2)
+        posy2=axisy + dcomh*(2*(quat1[1]*quat1[2]+quat1[0]*quat1[3]))
+        posz2=posz1
 
-	if (len(nucleotide)+1 > strandstart):
-	  topology.append([1,len(nucleotide),len(nucleotide)+1])
-	  comptopo.append([1,len(nucleotide)+len(strand[2]),len(nucleotide)+len(strand[2])+1])
+        if (len(nucleotide)+1 > strandstart):
+          topology.append([1,len(nucleotide),len(nucleotide)+1])
+          comptopo.append([1,len(nucleotide)+len(strand[2]),len(nucleotide)+len(strand[2])+1])
 
-	nucleotide.append(temp1)
-	compstrand.append(temp2)
+        nucleotide.append(temp1)
+        compstrand.append(temp2)
 
       for ib in range(len(compstrand)):
-	nucleotide.append(compstrand[len(compstrand)-1-ib])
+        nucleotide.append(compstrand[len(compstrand)-1-ib])
 
       for ib in range(len(comptopo)):
-	topology.append(comptopo[ib])
+        topology.append(comptopo[ib])
 
   return
 

examples/meam/msmeam/HGa.meam

@@ -1,30 +0,0 @@
-bkgd_dyn        =       1
-emb_lin_neg = 1
-augt1=0 
-ialloy=1 
-rc	=	 5.9 
-#H
-attrac(1,1)=0.460 
-repuls(1,1)=0.460 
-Cmin(1,1,1)=1.3 # PuMS
-Cmax(1,1,1)= 2.80 
-nn2(1,1)=1
-#Ga
-rho0(2)         =       0.6
-attrac(2,2)=0.097 
-repuls(2,2)=0.097 
-nn2(2,2)=1
-#HGa
-attrac(1,2)=0.300 
-repuls(1,2)=0.300 
-lattce(1,2)=l12 
-re(1,2)=3.19 
-delta(1,2)=-0.48  
-alpha(1,2)=6.6 
-Cmin(1,1,2)=2.0 
-Cmin(2,1,2)= 2.0 
-Cmin(1,2,1)=2.0 
-Cmin(2,2,1)     =       1.4
-Cmin(1,2,2)     =       1.4
-Cmin(1,1,2)     =       1.4
-nn2(1,2)=1

examples/meam/msmeam/HGa.msmeam

@@ -0,0 +1 @@
+../../../potentials/HGa.msmeam

examples/meam/msmeam/data.msmeam.bu

@@ -1,25 +0,0 @@
-LAMMPS data file via write_data, version 16 Feb 2016, timestep = 1
-
-3 atoms
-2 atom types
-
--4.0000000000000000e+00 4.0000000000000000e+00 xlo xhi
--4.0000000000000000e+00 4.0000000000000000e+00 ylo yhi
--4.0000000000000000e+00 4.0000000000000000e+00 zlo zhi
-
-Masses
-
-1 1.0079
-2 69.723
-
-Atoms # atomic
-
-1 1 0.0000000000000000e+00 0.0000000000000000e+00 0.0000000000000000e+00 0 0 0
-2 2 2.2000000000000002e+00 0.0000000000000000e+00 0.0000000000000000e+00 0 0 0
-3 2 2.9999999999999999e-01 2.2999999999999998e+00 0.0000000000000000e+00 0 0 0
-
-Velocities
-
-1 0.0000000000000000e+00 0.0000000000000000e+00 0.0000000000000000e+00
-2 0.0000000000000000e+00 0.0000000000000000e+00 0.0000000000000000e+00
-3 0.0000000000000000e+00 0.0000000000000000e+00 0.0000000000000000e+00

examples/meam/msmeam/in.msmeam

@@ -1,5 +1,3 @@
-echo	both
-log	log.msmeam
 # Test of MEAM potential for HGa
 
 # ------------------------ INITIALIZATION ----------------------------
@@ -21,11 +19,11 @@ create_atoms    1 single 0 0 0  units box
 create_atoms    2 single 2.2 0 0  units box
 create_atoms    2 single 0.3 2.3 0  units box
 # ---------- Define Settings ---------------------
-variable	teng equal "c_eatoms"
+variable        teng equal "c_eatoms"
 compute pot_energy all pe/atom
 compute stress all stress/atom NULL
 dump 1 all custom 1 dump.msmeam id x y z fx fy fz c_pot_energy c_stress[1] c_stress[2] c_stress[3] c_stress[4] c_stress[5] c_stress[6]
-run	1
-write_data	data.msmeam
+run     1
+#write_data      data.msmeam
 
 print "All done!"

examples/meam/msmeam/library.msmeam

@@ -1,14 +0,0 @@
-# DATE: 2018-09-22 UNITS: metal CONTRIBUTOR: Steve Valone, smv@lanl.gov CITATION: Baskes, PRB 1992; smv, sr, mib, JNM 2010
-# ms-meam data format May 2010
-#  elt   lat    z     ielement      atwt
-#  alpha b0     b1    b2     b3     b1m     b2m   b3m      alat   esub   asub
-#    -   t0     t1    t2     t3     t1m     t2m   t3m      rozero ibar
-#  NOTE:  leading character cannot be a space
-
-'H'    'dim'  1.0   1      1.0079
-2.960  2.960  3.0   1.0    1.0    1.0    3.0  1.0       0.741  2.235  2.50
-1.0    0.44721 0.0  0.00   0.0    0.31623 0  6.70 0
-
-'Ga4'  'fcc'  12.0  31     69.723
-4.42   4.80   3.10  6.00   0.00   0.0    0.0  0.5       4.247  2.897  0.97
-1.0    1.649 1.435  0.00   0.0    0.0  2.0       0.70   0

examples/meam/msmeam/log.1Mar2024.msmeam.g++.1

@@ -0,0 +1,126 @@
+LAMMPS (7 Feb 2024 - Development - patch_7Feb2024_update1-182-g93942f2013-modified)
+OMP_NUM_THREADS environment is not set. Defaulting to 1 thread. (src/comm.cpp:98)
+  using 1 OpenMP thread(s) per MPI task
+# Test of MEAM potential for HGa
+
+# ------------------------ INITIALIZATION ----------------------------
+units           metal
+dimension       3
+boundary        p       p       p
+atom_style      atomic
+variable latparam equal 4.646
+variable ncell equal 3
+
+# ----------------------- ATOM DEFINITION ----------------------------
+region          box block -4 4 -4 4 -4 4
+create_box      2 box
+Created orthogonal box = (-4 -4 -4) to (4 4 4)
+  1 by 1 by 1 MPI processor grid
+
+#
+
+include potential.mod
+# NOTE: This script can be modified for different pair styles
+# See in.elastic for more info.
+
+variable Pu string H
+print "potential chosen ${Pu}"
+potential chosen H
+# Choose potential
+pair_style meam/ms
+print		"we just executed"
+we just executed
+
+pair_coeff      * * library.msmeam ${Pu} Ga4  HGa.msmeam ${Pu} Ga4
+pair_coeff      * * library.msmeam H Ga4  HGa.msmeam ${Pu} Ga4
+pair_coeff      * * library.msmeam H Ga4  HGa.msmeam H Ga4
+Reading MEAM library file library.msmeam with DATE: 2018-09-22
+# Setup neighbor style
+neighbor 1.0 bin
+neigh_modify once no every 1 delay 0 check yes
+
+# Setup minimization style
+variable dmax equal 1.0e-2
+min_style	     cg
+min_modify	     dmax ${dmax} line quadratic
+min_modify	     dmax 0.01 line quadratic
+compute eng all pe/atom
+compute eatoms all reduce sum c_eng
+
+# Setup output
+thermo		100
+thermo_style custom step temp etotal  press pxx pyy pzz pxy pxz pyz lx ly lz vol c_eatoms
+thermo_modify norm yes
+create_atoms    1 single 0 0 0  units box
+Created 1 atoms
+  using box units in orthogonal box = (-4 -4 -4) to (4 4 4)
+  create_atoms CPU = 0.000 seconds
+create_atoms    2 single 2.2 0 0  units box
+Created 1 atoms
+  using box units in orthogonal box = (-4 -4 -4) to (4 4 4)
+  create_atoms CPU = 0.000 seconds
+create_atoms    2 single 0.3 2.3 0  units box
+Created 1 atoms
+  using box units in orthogonal box = (-4 -4 -4) to (4 4 4)
+  create_atoms CPU = 0.000 seconds
+# ---------- Define Settings ---------------------
+variable        teng equal "c_eatoms"
+compute pot_energy all pe/atom
+compute stress all stress/atom NULL
+# dump 1 all custom 1 dump.msmeam id x y z fx fy fz c_pot_energy c_stress[1] c_stress[2] c_stress[3] c_stress[4] c_stress[5] c_stress[6]
+run     1
+WARNING: No fixes with time integration, atoms won't move (src/verlet.cpp:60)
+Neighbor list info ...
+  update: every = 1 steps, delay = 0 steps, check = yes
+  max neighbors/atom: 2000, page size: 100000
+  master list distance cutoff = 6.9
+  ghost atom cutoff = 6.9
+  binsize = 3.45, bins = 3 3 3
+  2 neighbor lists, perpetual/occasional/extra = 2 0 0
+  (1) pair meam/ms, perpetual
+      attributes: full, newton on
+      pair build: full/bin/atomonly
+      stencil: full/bin/3d
+      bin: standard
+  (2) pair meam/ms, perpetual, half/full from (1)
+      attributes: half, newton on
+      pair build: halffull/newton
+      stencil: none
+      bin: none
+Per MPI rank memory allocation (min/avg/max) = 8.587 | 8.587 | 8.587 Mbytes
+   Step          Temp          TotEng         Press           Pxx            Pyy            Pzz            Pxy            Pxz            Pyz             Lx             Ly             Lz           Volume        c_eatoms   
+         0   0              15.433079      491354.7       838670.96      635393.15      0              80195.797      0              0              8              8              8              512            15.433079    
+         1   0              15.433079      491354.7       838670.96      635393.15      0              80195.797      0              0              8              8              8              512            15.433079    
+Loop time of 4.4446e-05 on 1 procs for 1 steps with 3 atoms
+
+Performance: 1943.932 ns/day, 0.012 hours/ns, 22499.213 timesteps/s, 67.498 katom-step/s
+31.5% CPU use with 1 MPI tasks x 1 OpenMP threads
+
+MPI task timing breakdown:
+Section |  min time  |  avg time  |  max time  |%varavg| %total
+---------------------------------------------------------------
+Pair    | 2.9908e-05 | 2.9908e-05 | 2.9908e-05 |   0.0 | 67.29
+Neigh   | 0          | 0          | 0          |   0.0 |  0.00
+Comm    | 1.033e-06  | 1.033e-06  | 1.033e-06  |   0.0 |  2.32
+Output  | 9.347e-06  | 9.347e-06  | 9.347e-06  |   0.0 | 21.03
+Modify  | 2.02e-07   | 2.02e-07   | 2.02e-07   |   0.0 |  0.45
+Other   |            | 3.956e-06  |            |       |  8.90
+
+Nlocal:              3 ave           3 max           3 min
+Histogram: 1 0 0 0 0 0 0 0 0 0
+Nghost:             78 ave          78 max          78 min
+Histogram: 1 0 0 0 0 0 0 0 0 0
+Neighs:              7 ave           7 max           7 min
+Histogram: 1 0 0 0 0 0 0 0 0 0
+FullNghs:           14 ave          14 max          14 min
+Histogram: 1 0 0 0 0 0 0 0 0 0
+
+Total # of neighbors = 14
+Ave neighs/atom = 4.6666667
+Neighbor list builds = 0
+Dangerous builds = 0
+#write_data      data.msmeam
+
+print "All done!"
+All done!
+Total wall time: 0:00:00

examples/meam/msmeam/log.1Mar2024.msmeam.g++.4

@@ -0,0 +1,126 @@
+LAMMPS (7 Feb 2024 - Development - patch_7Feb2024_update1-182-g93942f2013-modified)
+OMP_NUM_THREADS environment is not set. Defaulting to 1 thread. (src/comm.cpp:98)
+  using 1 OpenMP thread(s) per MPI task
+# Test of MEAM potential for HGa
+
+# ------------------------ INITIALIZATION ----------------------------
+units           metal
+dimension       3
+boundary        p       p       p
+atom_style      atomic
+variable latparam equal 4.646
+variable ncell equal 3
+
+# ----------------------- ATOM DEFINITION ----------------------------
+region          box block -4 4 -4 4 -4 4
+create_box      2 box
+Created orthogonal box = (-4 -4 -4) to (4 4 4)
+  1 by 2 by 2 MPI processor grid
+
+#
+
+include potential.mod
+# NOTE: This script can be modified for different pair styles
+# See in.elastic for more info.
+
+variable Pu string H
+print "potential chosen ${Pu}"
+potential chosen H
+# Choose potential
+pair_style meam/ms
+print		"we just executed"
+we just executed
+
+pair_coeff      * * library.msmeam ${Pu} Ga4  HGa.msmeam ${Pu} Ga4
+pair_coeff      * * library.msmeam H Ga4  HGa.msmeam ${Pu} Ga4
+pair_coeff      * * library.msmeam H Ga4  HGa.msmeam H Ga4
+Reading MEAM library file library.msmeam with DATE: 2018-09-22
+# Setup neighbor style
+neighbor 1.0 bin
+neigh_modify once no every 1 delay 0 check yes
+
+# Setup minimization style
+variable dmax equal 1.0e-2
+min_style	     cg
+min_modify	     dmax ${dmax} line quadratic
+min_modify	     dmax 0.01 line quadratic
+compute eng all pe/atom
+compute eatoms all reduce sum c_eng
+
+# Setup output
+thermo		100
+thermo_style custom step temp etotal  press pxx pyy pzz pxy pxz pyz lx ly lz vol c_eatoms
+thermo_modify norm yes
+create_atoms    1 single 0 0 0  units box
+Created 1 atoms
+  using box units in orthogonal box = (-4 -4 -4) to (4 4 4)
+  create_atoms CPU = 0.000 seconds
+create_atoms    2 single 2.2 0 0  units box
+Created 1 atoms
+  using box units in orthogonal box = (-4 -4 -4) to (4 4 4)
+  create_atoms CPU = 0.000 seconds
+create_atoms    2 single 0.3 2.3 0  units box
+Created 1 atoms
+  using box units in orthogonal box = (-4 -4 -4) to (4 4 4)
+  create_atoms CPU = 0.000 seconds
+# ---------- Define Settings ---------------------
+variable        teng equal "c_eatoms"
+compute pot_energy all pe/atom
+compute stress all stress/atom NULL
+# dump 1 all custom 1 dump.msmeam id x y z fx fy fz c_pot_energy c_stress[1] c_stress[2] c_stress[3] c_stress[4] c_stress[5] c_stress[6]
+run     1
+WARNING: No fixes with time integration, atoms won't move (src/verlet.cpp:60)
+Neighbor list info ...
+  update: every = 1 steps, delay = 0 steps, check = yes
+  max neighbors/atom: 2000, page size: 100000
+  master list distance cutoff = 6.9
+  ghost atom cutoff = 6.9
+  binsize = 3.45, bins = 3 3 3
+  2 neighbor lists, perpetual/occasional/extra = 2 0 0
+  (1) pair meam/ms, perpetual
+      attributes: full, newton on
+      pair build: full/bin/atomonly
+      stencil: full/bin/3d
+      bin: standard
+  (2) pair meam/ms, perpetual, half/full from (1)
+      attributes: half, newton on
+      pair build: halffull/newton
+      stencil: none
+      bin: none
+Per MPI rank memory allocation (min/avg/max) = 7.965 | 8.123 | 8.594 Mbytes
+   Step          Temp          TotEng         Press           Pxx            Pyy            Pzz            Pxy            Pxz            Pyz             Lx             Ly             Lz           Volume        c_eatoms   
+         0   0              15.433079      491354.7       838670.96      635393.15      0              80195.797      0              0              8              8              8              512            15.433079    
+         1   0              15.433079      491354.7       838670.96      635393.15      0              80195.797      0              0              8              8              8              512            15.433079    
+Loop time of 8.70645e-05 on 4 procs for 1 steps with 3 atoms
+
+Performance: 992.368 ns/day, 0.024 hours/ns, 11485.738 timesteps/s, 34.457 katom-step/s
+29.0% CPU use with 4 MPI tasks x 1 OpenMP threads
+
+MPI task timing breakdown:
+Section |  min time  |  avg time  |  max time  |%varavg| %total
+---------------------------------------------------------------
+Pair    | 4.3957e-05 | 4.67e-05   | 5.1056e-05 |   0.0 | 53.64
+Neigh   | 0          | 0          | 0          |   0.0 |  0.00
+Comm    | 1.105e-05  | 1.3822e-05 | 1.7033e-05 |   0.0 | 15.88
+Output  | 1.5765e-05 | 1.9045e-05 | 2.5216e-05 |   0.0 | 21.87
+Modify  | 2.58e-07   | 3.465e-07  | 3.81e-07   |   0.0 |  0.40
+Other   |            | 7.151e-06  |            |       |  8.21
+
+Nlocal:           0.75 ave           3 max           0 min
+Histogram: 3 0 0 0 0 0 0 0 0 1
+Nghost:          38.25 ave          42 max          36 min
+Histogram: 2 0 0 0 0 1 0 0 0 1
+Neighs:           1.75 ave           7 max           0 min
+Histogram: 3 0 0 0 0 0 0 0 0 1
+FullNghs:          3.5 ave          14 max           0 min
+Histogram: 3 0 0 0 0 0 0 0 0 1
+
+Total # of neighbors = 14
+Ave neighs/atom = 4.6666667
+Neighbor list builds = 0
+Dangerous builds = 0
+#write_data      data.msmeam
+
+print "All done!"
+All done!
+Total wall time: 0:00:00

examples/meam/msmeam/log.msmeam.bu

@@ -1,107 +0,0 @@
-# Test of MEAM potential for HGa
-
-# ------------------------ INITIALIZATION ----------------------------
-units           metal
-dimension       3
-boundary        p       p       p
-atom_style      atomic
-variable latparam equal 4.646
-variable ncell equal 3
-
-# ----------------------- ATOM DEFINITION ----------------------------
-region          box block -4 4 -4 4 -4 4
-create_box      2 box
-Created orthogonal box = (-4 -4 -4) to (4 4 4)
-  1 by 1 by 1 MPI processor grid
-
-#
-
-include potential.mod
-# NOTE: This script can be modified for different pair styles
-# See in.elastic for more info.
-
-variable Pu string H
-print "potential chosen ${Pu}"
-potential chosen H
-# Choose potential
-pair_style      MSmeam
-print		"we just executed"
-we just executed
-
-pair_coeff      * * library.MSmeam ${Pu} Ga4  HGaMS.meam ${Pu} Ga4
-pair_coeff      * * library.MSmeam H Ga4  HGaMS.meam ${Pu} Ga4
-pair_coeff      * * library.MSmeam H Ga4  HGaMS.meam H Ga4
-Reading potential file library.MSmeam with DATE: 2018-09-22
-# Setup neighbor style
-neighbor 1.0 nsq
-neigh_modify once no every 1 delay 0 check yes
-
-# Setup minimization style
-variable dmax equal 1.0e-2
-min_style	     cg
-min_modify	     dmax ${dmax} line quadratic
-min_modify	     dmax 0.01 line quadratic
-compute eng all pe/atom
-compute eatoms all reduce sum c_eng
-
-# Setup output
-thermo		100
-thermo_style custom step temp etotal  press pxx pyy pzz pxy pxz pyz lx ly lz vol c_eatoms
-thermo_modify norm yes
-create_atoms    1 single 0 0 0  units box
-Created 1 atoms
-create_atoms    2 single 2.2 0 0  units box
-Created 1 atoms
-create_atoms    2 single 0.3 2.3 0  units box
-Created 1 atoms
-# ---------- Define Settings ---------------------
-variable	teng equal "c_eatoms"
-compute pot_energy all pe/atom
-compute stress all stress/atom NULL
-dump 1 all custom 1 dump.msmeam id x y z fx fy fz c_pot_energy c_stress[1] c_stress[2] c_stress[3] c_stress[4] c_stress[5] c_stress[6]
-run	1
-WARNING: No fixes defined, atoms won't move (../verlet.cpp:55)
-Neighbor list info ...
-  2 neighbor list requests
-  update every 1 steps, delay 0 steps, check yes
-  max neighbors/atom: 2000, page size: 100000
-  master list distance cutoff = 6.9
-  ghost atom cutoff = 6.9
-Memory usage per processor = 12.9295 Mbytes
-Step Temp TotEng Press Pxx Pyy Pzz Pxy Pxz Pyz Lx Ly Lz Volume eatoms 
-       0            0    15.433079    491354.68    838670.91    635393.13            0    80195.793            0            0            8            8            8          512    15.433079 
-       1            0    15.433079    491354.68    838670.91    635393.13            0    80195.793            0            0            8            8            8          512    15.433079 
-Loop time of 0.000172138 on 1 procs for 1 steps with 3 atoms
-
-Performance: 501.922 ns/day, 0.048 hours/ns, 5809.285 timesteps/s
-81.3% CPU use with 1 MPI tasks x no OpenMP threads
-
-MPI task timing breakdown:
-Section |  min time  |  avg time  |  max time  |%varavg| %total
----------------------------------------------------------------
-Pair    | 6.6996e-05 | 6.6996e-05 | 6.6996e-05 |   0.0 | 38.92
-Neigh   | 0          | 0          | 0          |   0.0 |  0.00
-Comm    | 1.9073e-06 | 1.9073e-06 | 1.9073e-06 |   0.0 |  1.11
-Output  | 9.7036e-05 | 9.7036e-05 | 9.7036e-05 |   0.0 | 56.37
-Modify  | 0          | 0          | 0          |   0.0 |  0.00
-Other   |            | 6.199e-06  |            |       |  3.60
-
-Nlocal:    3 ave 3 max 3 min
-Histogram: 1 0 0 0 0 0 0 0 0 0
-Nghost:    78 ave 78 max 78 min
-Histogram: 1 0 0 0 0 0 0 0 0 0
-Neighs:    7 ave 7 max 7 min
-Histogram: 1 0 0 0 0 0 0 0 0 0
-FullNghs:  14 ave 14 max 14 min
-Histogram: 1 0 0 0 0 0 0 0 0 0
-
-Total # of neighbors = 14
-Ave neighs/atom = 4.66667
-Neighbor list builds = 0
-Dangerous builds = 0
-write_data	data.msmeam
-
-print "All done!"
-All done!
-Total wall time: 0:00:00
-

examples/meam/msmeam/msmeam.dump.bu

@@ -1,24 +0,0 @@
-ITEM: TIMESTEP
-0
-ITEM: NUMBER OF ATOMS
-3
-ITEM: BOX BOUNDS pp pp pp
--4 4
--4 4
--4 4
-ITEM: ATOMS id x y z fx fy fz c_pot_energy c_stress[1] c_stress[2] c_stress[3] c_stress[4] c_stress[5] c_stress[6] 
-1 0 0 0 -131.925 -88.3005 0 22.9153 -2.147e+08 -1.62661e+08 -0 -2.05301e+07 -0 -0 
-2 2.2 0 0 120.809 -0.482171 0 14.7692 -2.12028e+08 -0 -0 403352 -0 -0 
-3 0.3 2.3 0 11.1159 88.7827 0 8.61478 -2.67145e+06 -1.62661e+08 -0 -2.09335e+07 -0 -0 
-ITEM: TIMESTEP
-1
-ITEM: NUMBER OF ATOMS
-3
-ITEM: BOX BOUNDS pp pp pp
--4 4
--4 4
--4 4
-ITEM: ATOMS id x y z fx fy fz c_pot_energy c_stress[1] c_stress[2] c_stress[3] c_stress[4] c_stress[5] c_stress[6] 
-1 0 0 0 -131.925 -88.3005 0 22.9153 -2.147e+08 -1.62661e+08 -0 -2.05301e+07 -0 -0 
-2 2.2 0 0 120.809 -0.482171 0 14.7692 -2.12028e+08 -0 -0 403352 -0 -0 
-3 0.3 2.3 0 11.1159 88.7827 0 8.61478 -2.67145e+06 -1.62661e+08 -0 -2.09335e+07 -0 -0 

examples/meam/msmeam/potential.mod

@@ -7,7 +7,7 @@ print "potential chosen ${Pu}"
 pair_style meam/ms
 print		"we just executed"
 
-pair_coeff      * * library.msmeam ${Pu} Ga4  HGa.meam ${Pu} Ga4
+pair_coeff      * * library.msmeam ${Pu} Ga4  HGa.msmeam ${Pu} Ga4
 # Setup neighbor style
 neighbor 1.0 bin
 neigh_modify once no every 1 delay 0 check yes

lib/gpu/geryon/nvd_device.h

@@ -138,7 +138,7 @@ class UCL_Device {
   /** \note You cannot delete the default stream **/
   inline void pop_command_queue() {
     if (_cq.size()<2) return;
-    CU_SAFE_CALL_NS(cuStreamDestroy(_cq.back()));
+    cuStreamDestroy(_cq.back());
     _cq.pop_back();
   }
 
@@ -426,8 +426,8 @@ void UCL_Device::clear() {
   if (_device>-1) {
     for (int i=1; i<num_queues(); i++) pop_command_queue();
 #if GERYON_NVD_PRIMARY_CONTEXT
-    CU_SAFE_CALL_NS(cuCtxSetCurrent(_old_context));
-    CU_SAFE_CALL_NS(cuDevicePrimaryCtxRelease(_cu_device));
+    cuCtxSetCurrent(_old_context);
+    cuDevicePrimaryCtxRelease(_cu_device);
 #else
     cuCtxDestroy(_context);
 #endif

lib/mdi/Install.py

@@ -32,7 +32,7 @@ make lib-mdi args="-m mpi" # build MDI lib with same settings as in the mpi Make
 
 # settings
 
-version = "1.4.16"
+version = "1.4.26"
 url = "https://github.com/MolSSI-MDI/MDI_Library/archive/v%s.tar.gz" % version
 
 # known checksums for different MDI versions. used to validate the download.
@@ -41,6 +41,7 @@ checksums = { \
               '1.4.12' : '7a222353ae8e03961d5365e6cd48baee', \
               '1.4.14' : '7a059bb12535360fdcb7de8402f9a0fc', \
               '1.4.16' : '407db44e2d79447ab5c1233af1965f65', \
+              '1.4.26' : '3124bb85259471e2a53a891f04bf697a', \
               }
 
 # print error message or help

lib/pace/Install.py

@@ -18,11 +18,11 @@ from install_helpers import fullpath, geturl, checkmd5sum, getfallback
 # settings
 
 thisdir = fullpath('.')
-version ='v.2023.01.3.fix'
+version ='v.2023.10.04'
 
 # known checksums for different PACE versions. used to validate the download.
 checksums = { \
-    'v.2023.01.3.fix': '4f0b3b5b14456fe9a73b447de3765caa'
+    'v.2023.10.04': '70ff79f4e59af175e55d24f3243ad1ff'
 }
 
 parser = ArgumentParser(prog='Install.py', description="LAMMPS library build wrapper script")

potentials/HGa.msmeam

@@ -1,29 +1,29 @@
 bkgd_dyn        =       1
 emb_lin_neg = 1
-augt1=0 
-ialloy=1 
-rc	=	 5.9 
+augt1=0
+ialloy=1
+rc      =        5.9
 #H
-attrac(1,1)=0.460 
-repuls(1,1)=0.460 
+attrac(1,1)=0.460
+repuls(1,1)=0.460
 Cmin(1,1,1)=1.3 # PuMS
-Cmax(1,1,1)= 2.80 
+Cmax(1,1,1)= 2.80
 nn2(1,1)=1
 #Ga
 rho0(2)         =       0.6
-attrac(2,2)=0.097 
-repuls(2,2)=0.097 
+attrac(2,2)=0.097
+repuls(2,2)=0.097
 nn2(2,2)=1
 #HGa
-attrac(1,2)=0.300 
-repuls(1,2)=0.300 
-lattce(1,2)=l12 
-re(1,2)=3.19 
-delta(1,2)=-0.48  
-alpha(1,2)=6.6 
-Cmin(1,1,2)=2.0 
-Cmin(2,1,2)= 2.0 
-Cmin(1,2,1)=2.0 
+attrac(1,2)=0.300
+repuls(1,2)=0.300
+lattce(1,2)=l12
+re(1,2)=3.19
+delta(1,2)=-0.48
+alpha(1,2)=6.6
+Cmin(1,1,2)=2.0
+Cmin(2,1,2)= 2.0
+Cmin(1,2,1)=2.0
 Cmin(2,2,1)     =       1.4
 Cmin(1,2,2)     =       1.4
 Cmin(1,1,2)     =       1.4

python/install.py

@@ -18,7 +18,7 @@ parser = ArgumentParser(prog='install.py',
                         description='LAMMPS python package installer script')
 
 parser.add_argument("-p", "--package", required=True,
-                    help="path to the LAMMPS Python package")
+                    help="path to the LAMMPS Python package folder")
 parser.add_argument("-l", "--lib", required=True,
                     help="path to the compiled LAMMPS shared library")
 parser.add_argument("-n", "--noinstall", action="store_true", default=False,
@@ -34,15 +34,21 @@ args = parser.parse_args()
 
 if args.package:
   if not os.path.exists(args.package):
-    print("ERROR: LAMMPS package folder %s does not exist" % args.package)
+    print("\nERROR: LAMMPS package folder %s does not exist\n" % args.package)
     parser.print_help()
     sys.exit(1)
   else:
     args.package = os.path.abspath(args.package)
+    if ((os.path.basename(args.package) != "lammps")
+        and ((os.path.basename(os.path.dirname(args.package)) != "python"))):
+             print("\nERROR: LAMMPS package folder path %s does not end in %s\n"
+                   % (args.package, os.path.join("python", "lammps")))
+             parser.print_help()
+             sys.exit(1)
 
 if args.lib:
   if not os.path.exists(args.lib):
-    print("ERROR: LAMMPS shared library %s does not exist" % args.lib)
+    print("\nERROR: LAMMPS shared library %s does not exist\n" % args.lib)
     parser.print_help()
     sys.exit(1)
   else:
@@ -50,7 +56,7 @@ if args.lib:
 
 if args.wheeldir:
   if not os.path.exists(args.wheeldir):
-    print("ERROR: directory %s to store the wheel does not exist" % args.wheeldir)
+    print("\nERROR: directory %s to store the wheel does not exist\n" % args.wheeldir)
     parser.print_help()
     sys.exit(1)
   else:
@@ -58,7 +64,7 @@ if args.wheeldir:
 
 if args.versionfile:
   if not os.path.exists(args.versionfile):
-    print("ERROR: LAMMPS version file at %s does not exist" % args.versionfile)
+    print("\nERROR: LAMMPS version file at %s does not exist\n" % args.versionfile)
     parser.print_help()
     sys.exit(1)
   else:

python/lammps/core.py

@@ -379,8 +379,9 @@ class lammps(object):
           for i in range(narg):
             if type(cmdargs[i]) is str:
               cmdargs[i] = cmdargs[i].encode()
-          cargs = (c_char_p*narg)(*cmdargs)
-          self.lib.lammps_open.argtypes = [c_int, c_char_p*narg, MPI_Comm, c_void_p]
+          cargs = (c_char_p*(narg+1))(*cmdargs)
+          cargs[narg] = None
+          self.lib.lammps_open.argtypes = [c_int, c_char_p*(narg+1), MPI_Comm, c_void_p]
         else:
           self.lib.lammps_open.argtypes = [c_int, c_char_p, MPI_Comm, c_void_p]
 
@@ -399,8 +400,9 @@ class lammps(object):
           for i in range(narg):
             if type(cmdargs[i]) is str:
               cmdargs[i] = cmdargs[i].encode()
-          cargs = (c_char_p*narg)(*cmdargs)
-          self.lib.lammps_open_no_mpi.argtypes = [c_int, c_char_p*narg, c_void_p]
+          cargs = (c_char_p*(narg+1))(*cmdargs)
+          cargs[narg] = None
+          self.lib.lammps_open_no_mpi.argtypes = [c_int, c_char_p*(narg+1), c_void_p]
           self.lmp = c_void_p(self.lib.lammps_open_no_mpi(narg,cargs,None))
         else:
           self.lib.lammps_open_no_mpi.argtypes = [c_int, c_char_p, c_void_p]

src/BOCS/fix_bocs.cpp

@@ -1024,7 +1024,10 @@ void FixBocs::final_integrate()
 
   if (pstat_flag) {
     if (pstyle == ISO) pressure->compute_scalar();
-    else pressure->compute_vector();
+    else {
+      temperature->compute_vector();
+      pressure->compute_vector();
+    }
     couple();
     pressure->addstep(update->ntimestep+1);
   }
@@ -1961,6 +1964,7 @@ void FixBocs::nhc_press_integrate()
   int ich,i,pdof;
   double expfac,factor_etap,kecurrent;
   double kt = boltz * t_target;
+  double lkt_press;
 
   // Update masses, to preserve initial freq, if flag set
 
@@ -2006,7 +2010,8 @@ void FixBocs::nhc_press_integrate()
     }
   }
 
-  double lkt_press = pdof * kt;
+  if (pstyle == ISO) lkt_press = kt;
+  else lkt_press = pdof * kt;
   etap_dotdot[0] = (kecurrent - lkt_press)/etap_mass[0];
 
   double ncfac = 1.0/nc_pchain;

src/DIFFRACTION/compute_xrd.cpp

@@ -261,7 +261,7 @@ void ComputeXRD::init()
   double ang = 0.0;
 
   double convf = 360 / MY_PI;
-  if (radflag ==1) convf = 1;
+  if (radflag == 1) convf = 2;
 
   int n = 0;
   for (int m = 0; m < mmax; m++) {

src/DPD-MESO/fix_mvv_dpd.cpp

@@ -63,13 +63,20 @@ int FixMvvDPD::setmask()
 void FixMvvDPD::init()
 {
   if (!atom->vest_flag)
-    error->all(FLERR,"Fix mvv/dpd requires atom attribute vest");
+    error->all(FLERR,"Fix mvv/dpd requires atom attribute vest e.g. from atom style mdpd");
+
+  if (!force->pair_match("^mdpd",0) && !force->pair_match("^dpd",0)) {
+    if (force->pair_match("^hybrid",0)) {
+      if (!(force->pair_match("^mdpd",0,1) || force->pair_match("^dpd",0),1)) {
+        error->all(FLERR, "Must use a dpd or mdpd pair style with fix mvv/dpd");
+      }
+    } else {
+      error->all(FLERR, "Must use a dpd or mdpd pair style with fix mvv/dpd");
+    }
+  }
 
   dtv = update->dt;
   dtf = 0.5 * update->dt * force->ftm2v;
-
-  if (!force->pair_match("^edpd",0) && !force->pair_match("^dpd",0))
-    error->all(FLERR, "Must use a dpd or edpd pair style with fix mvv/edpd");
 }
 
 /* ----------------------------------------------------------------------

src/DPD-MESO/fix_mvv_edpd.cpp

@@ -71,11 +71,20 @@ int FixMvvEDPD::setmask()
 
 void FixMvvEDPD::init()
 {
+  if (!atom->edpd_flag) error->all(FLERR,"Fix mvv/edpd requires atom style edpd");
+
+  if (!force->pair_match("^edpd",0)) {
+    if (force->pair_match("^hybrid",0)) {
+      if (!force->pair_match("^edpd",0,1)) {
+        error->all(FLERR, "Must use pair style edpd with fix mvv/edpd");
+      }
+    } else {
+      error->all(FLERR, "Must use pair style edpd with fix mvv/edpd");
+    }
+  }
+
   dtv = update->dt;
   dtf = 0.5 * update->dt * force->ftm2v;
-
-  if (!force->pair_match("^edpd",0))
-    error->all(FLERR, "Must use pair style edpd with fix mvv/edpd");
 }
 
 /* ----------------------------------------------------------------------

src/DPD-MESO/fix_mvv_tdpd.cpp

@@ -69,10 +69,20 @@ int FixMvvTDPD::setmask()
 
 void FixMvvTDPD::init()
 {
+  if (!atom->tdpd_flag) error->all(FLERR,"Fix mvv/tdpd requires atom style tdpd");
+
+  if (!force->pair_match("^tdpd",0)) {
+    if (force->pair_match("^hybrid",0)) {
+      if (!force->pair_match("^tdpd",0,1)) {
+        error->all(FLERR, "Must use pair style tdpd with fix mvv/tdpd");
+      }
+    } else {
+      error->all(FLERR, "Must use pair style tdpd with fix mvv/tdpd");
+    }
+  }
+
   dtv = update->dt;
   dtf = 0.5 * update->dt * force->ftm2v;
-  if (!force->pair_match("^tdpd",0))
-    error->all(FLERR, "Must use pair style tdpd with fix mvv/tdpd");
 }
 
 /* ----------------------------------------------------------------------

src/Depend.sh

@@ -64,6 +64,7 @@ fi
 
 if (test $1 = "COLLOID") then
   depend GPU
+  depend KOKKOS
   depend OPENMP
 fi
 

src/EFF/fix_langevin_eff.cpp

@@ -137,7 +137,7 @@ void FixLangevinEff::post_force_no_tally()
   dof = domain->dimension * particles;
   fix_dof = 0;
   for (int i = 0; i < modify->nfix; i++)
-    fix_dof += modify->fix[i]->dof(igroup);
+    fix_dof += (int)modify->fix[i]->dof(igroup);
 
   // extra_dof = domain->dimension
   dof -= domain->dimension + fix_dof;
@@ -306,7 +306,7 @@ void FixLangevinEff::post_force_tally()
   dof = domain->dimension * particles;
   fix_dof = 0;
   for (int i = 0; i < modify->nfix; i++)
-    fix_dof += modify->fix[i]->dof(igroup);
+    fix_dof += (int)modify->fix[i]->dof(igroup);
 
   // extra_dof = domain->dimension
   dof -= domain->dimension + fix_dof;

src/ELECTRODE/pppm_electrode.cpp

@@ -633,7 +633,9 @@ void PPPMElectrode::project_psi(double *vec, int sensor_grpbit)
   // project u_brick with weight matrix
   double **x = atom->x;
   int *mask = atom->mask;
-  double const scaleinv = 1.0 / (nx_pppm * ny_pppm * nz_pppm);
+  const bigint ngridtotal = (bigint) nx_pppm * ny_pppm * nz_pppm;
+  const double scaleinv = 1.0 / ngridtotal;
+
   for (int i = 0; i < atom->nlocal; i++) {
     if (!(mask[i] & sensor_grpbit)) continue;
     double v = 0.;
@@ -1362,7 +1364,7 @@ double PPPMElectrode::compute_qopt()
   // each proc calculates contributions from every Pth grid point
 
   bigint ngridtotal = (bigint) nx_pppm * ny_pppm * nz_pppm;
-  int nxy_pppm = nx_pppm * ny_pppm;
+  bigint nxy_pppm = (bigint) nx_pppm * ny_pppm;
 
   double qopt = 0.0;
 

src/EXTRA-DUMP/dump_yaml.cpp

@@ -24,6 +24,8 @@
 
 using namespace LAMMPS_NS;
 
+static constexpr char special_chars[] = "{}[],&:*#?|-<>=!%@\\";
+
 /* ---------------------------------------------------------------------- */
 DumpYAML::DumpYAML(class LAMMPS *_lmp, int narg, char **args) :
     DumpCustom(_lmp, narg, args), thermo(false)
@@ -67,7 +69,12 @@ void DumpYAML::write_header(bigint ndump)
       const auto &fields = th->get_fields();
 
       thermo_data += "thermo:\n  - keywords: [ ";
-      for (int i = 0; i < nfield; ++i) thermo_data += fmt::format("{}, ", keywords[i]);
+      for (int i = 0; i < nfield; ++i) {
+        if (keywords[i].find_first_of(special_chars) == std::string::npos)
+          thermo_data += fmt::format("{}, ", keywords[i]);
+        else
+          thermo_data += fmt::format("'{}', ", keywords[i]);
+      }
       thermo_data += "]\n  - data: [ ";
 
       for (int i = 0; i < nfield; ++i) {
@@ -107,7 +114,12 @@ void DumpYAML::write_header(bigint ndump)
     if (domain->triclinic) fmt::print(fp, "  - [ {}, {}, {} ]\n", boxxy, boxxz, boxyz);
 
     fmt::print(fp, "keywords: [ ");
-    for (const auto &item : utils::split_words(columns)) fmt::print(fp, "{}, ", item);
+    for (const auto &item : utils::split_words(columns)) {
+      if (item.find_first_of(special_chars) == std::string::npos)
+        fmt::print(fp, "{}, ", item);
+      else
+        fmt::print(fp, "'{}', ", item);
+    }
     fputs(" ]\ndata:\n", fp);
   } else    // reset so that the remainder of the output is not multi-proc
     filewriter = 0;

src/EXTRA-FIX/fix_viscosity.cpp

@@ -120,14 +120,14 @@ int FixViscosity::setmask()
 
 void FixViscosity::init()
 {
-  // warn if any fix ave/spatial comes after this fix
+  // warn if any fix ave/chunk comes after this fix
   // can cause glitch in averaging since ave will happen after swap
 
   int foundme = 0;
-  for (int i = 0; i < modify->nfix; i++) {
-    if (modify->fix[i] == this) foundme = 1;
-    if (foundme && strcmp(modify->fix[i]->style,"ave/spatial") == 0 && me == 0)
-      error->warning(FLERR,"Fix viscosity comes before fix ave/spatial");
+  for (const auto &ifix : modify->get_fix_list()) {
+    if (ifix == this) foundme = 1;
+    if (foundme && utils::strmatch(ifix->style,"^ave/chunk") && (me == 0))
+      error->warning(FLERR,"Fix viscosity comes before fix ave/chunk");
   }
 
   // set bounds of 2 slabs in pdim

src/EXTRA-MOLECULE/angle_cosine_periodic.cpp

@@ -117,7 +117,7 @@ void AngleCosinePeriodic::compute(int eflag, int vflag)
     tn = 1.0;
     tn_1 = 1.0;
     tn_2 = 0.0;
-    un = 1.0;
+    un = (m==1) ? 2.0 : 1.0;
     un_1 = 2.0;
     un_2 = 0.0;
 

src/FEP/fix_adapt_fep.cpp

@@ -40,17 +40,17 @@ using namespace LAMMPS_NS;
 using namespace FixConst;
 using namespace MathConst;
 
-enum{PAIR,KSPACE,ATOM};
-enum{DIAMETER,CHARGE};
+enum{PAIR, KSPACE, ATOM};
+enum{DIAMETER, CHARGE};
 
 /* ---------------------------------------------------------------------- */
 
 FixAdaptFEP::FixAdaptFEP(LAMMPS *lmp, int narg, char **arg) :
   Fix(lmp, narg, arg)
 {
-  if (narg < 5) error->all(FLERR,"Illegal fix adapt/fep command");
+  if (narg < 5) utils::missing_cmd_args(FLERR,"fix adapt/fep", error);
   nevery = utils::inumeric(FLERR,arg[3],false,lmp);
-  if (nevery < 0) error->all(FLERR,"Illegal fix adapt/fep command");
+  if (nevery < 0) error->all(FLERR,"Illegal fix adapt/fep every value {}", nevery);
 
   dynamic_group_allow = 1;
   create_attribute = 1;
@@ -62,21 +62,21 @@ FixAdaptFEP::FixAdaptFEP(LAMMPS *lmp, int narg, char **arg) :
   int iarg = 4;
   while (iarg < narg) {
     if (strcmp(arg[iarg],"pair") == 0) {
-      if (iarg+6 > narg) error->all(FLERR,"Illegal fix adapt/fep command");
+      if (iarg+6 > narg) utils::missing_cmd_args(FLERR,"fix adapt/fep pair", error);
       nadapt++;
       iarg += 6;
     } else if (strcmp(arg[iarg],"kspace") == 0) {
-      if (iarg+2 > narg) error->all(FLERR,"Illegal fix adapt/fep command");
+      if (iarg+2 > narg) utils::missing_cmd_args(FLERR,"fix adapt/fep kspace", error);
       nadapt++;
       iarg += 2;
     } else if (strcmp(arg[iarg],"atom") == 0) {
-      if (iarg+4 > narg) error->all(FLERR,"Illegal fix adapt/fep command");
+      if (iarg+4 > narg) utils::missing_cmd_args(FLERR,"fix adapt/fep atom", error);
       nadapt++;
       iarg += 4;
     } else break;
   }
 
-  if (nadapt == 0) error->all(FLERR,"Illegal fix adapt/fep command");
+  if (nadapt == 0) error->all(FLERR,"Nothing to adapt in fix adapt/fep command");
   adapt = new Adapt[nadapt];
 
   // parse keywords
@@ -136,11 +136,11 @@ FixAdaptFEP::FixAdaptFEP(LAMMPS *lmp, int narg, char **arg) :
 
   while (iarg < narg) {
     if (strcmp(arg[iarg],"reset") == 0) {
-      if (iarg+2 > narg) error->all(FLERR,"Illegal fix adapt/fep command");
+      if (iarg+2 > narg) utils::missing_cmd_args(FLERR,"fix adapt/fep reset", error);
       resetflag = utils::logical(FLERR,arg[iarg+1],false,lmp);
       iarg += 2;
     } else if (strcmp(arg[iarg],"scale") == 0) {
-      if (iarg+2 > narg) error->all(FLERR,"Illegal fix adapt/fep command");
+      if (iarg+2 > narg) utils::missing_cmd_args(FLERR,"fix adapt/fep scale", error);
       scaleflag = utils::logical(FLERR,arg[iarg+1],false,lmp);
       iarg += 2;
     } else if (strcmp(arg[iarg],"after") == 0) {
@@ -165,21 +165,21 @@ FixAdaptFEP::FixAdaptFEP(LAMMPS *lmp, int narg, char **arg) :
 FixAdaptFEP::~FixAdaptFEP()
 {
   for (int m = 0; m < nadapt; m++) {
-    delete [] adapt[m].var;
+    delete[] adapt[m].var;
     if (adapt[m].which == PAIR) {
-      delete [] adapt[m].pstyle;
-      delete [] adapt[m].pparam;
+      delete[] adapt[m].pstyle;
+      delete[] adapt[m].pparam;
       memory->destroy(adapt[m].array_orig);
     }
   }
-  delete [] adapt;
+  delete[] adapt;
 
   // check nfix in case all fixes have already been deleted
 
   if (id_fix_diam && modify->nfix) modify->delete_fix(id_fix_diam);
   if (id_fix_chg && modify->nfix) modify->delete_fix(id_fix_chg);
-  delete [] id_fix_diam;
-  delete [] id_fix_chg;
+  delete[] id_fix_diam;
+  delete[] id_fix_chg;
 }
 
 /* ---------------------------------------------------------------------- */
@@ -208,7 +208,7 @@ void FixAdaptFEP::post_constructor()
   id_fix_diam = nullptr;
   id_fix_chg = nullptr;
 
-  if (diam_flag) {
+  if (diam_flag && atom->radius_flag) {
     id_fix_diam = utils::strdup(id + std::string("_FIX_STORE_DIAM"));
     fix_diam = dynamic_cast<FixStoreAtom *>(
       modify->add_fix(fmt::format("{} {} STORE/ATOM 1 0 0 1", id_fix_diam,group->names[igroup])));
@@ -226,7 +226,7 @@ void FixAdaptFEP::post_constructor()
     }
   }
 
-  if (chgflag) {
+  if (chgflag && atom->q_flag) {
     id_fix_chg = utils::strdup(id + std::string("_FIX_STORE_CHG"));
     fix_chg = dynamic_cast<FixStoreAtom *>(
       modify->add_fix(fmt::format("{} {} STORE/ATOM 1 0 0 1",id_fix_chg,group->names[igroup])));
@@ -267,9 +267,9 @@ void FixAdaptFEP::init()
 
     ad->ivar = input->variable->find(ad->var);
     if (ad->ivar < 0)
-      error->all(FLERR,"Variable name for fix adapt/fep does not exist");
+      error->all(FLERR,"Variable name {} for fix adapt/fep does not exist", ad->var);
     if (!input->variable->equalstyle(ad->ivar))
-      error->all(FLERR,"Variable for fix adapt/fep is invalid style");
+      error->all(FLERR,"Variable {} for fix adapt/fep is invalid style", ad->var);
 
     if (ad->which == PAIR) {
       anypair = 1;
@@ -285,8 +285,9 @@ void FixAdaptFEP::init()
       if (ptr == nullptr)
         error->all(FLERR,"Fix adapt/fep pair style param not supported");
 
-      ad->pdim = 2;
-      if (ad->pdim == 0) ad->scalar = (double *) ptr;
+      if (ad->pdim != 2)
+        error->all(FLERR,"Pair style parameter {} is not compatible with fix adapt/fep", ad->pparam);
+
       if (ad->pdim == 2) ad->array = (double **) ptr;
 
       // if pair hybrid, test that ilo,ihi,jlo,jhi are valid for sub-style
@@ -434,6 +435,8 @@ void FixAdaptFEP::change_settings()
 
     } else if (ad->which == ATOM) {
 
+      if (scaleflag)
+        error->all(FLERR, "Keyword 'scale yes' is not supported with fix adapt/fep for 'atom'");
       // reset radius from diameter
       // also scale rmass to new value
 

src/GPU/pppm_gpu.cpp

@@ -405,7 +405,8 @@ void PPPMGPU::poisson_ik()
 
   // if requested, compute energy and virial contribution
 
-  double scaleinv = 1.0/(nx_pppm*ny_pppm*nz_pppm);
+  bigint ngridtotal = (bigint) nx_pppm * ny_pppm * nz_pppm;
+  double scaleinv = 1.0 / ngridtotal;
   double s2 = scaleinv*scaleinv;
 
   if (eflag_global || vflag_global) {

src/INTEL/fix_intel.cpp

@@ -551,6 +551,9 @@ void FixIntel::kspace_init_check()
 
   if (intel_pair == 0)
     error->all(FLERR,"Intel styles for kspace require intel pair style.");
+
+  if (utils::strmatch(update->integrate_style, "^verlet/split"))
+    error->all(FLERR,"Intel styles for kspace are not compatible with run_style verlet/split");
 }
 
 /* ---------------------------------------------------------------------- */

src/INTEL/pppm_electrode_intel.cpp

@@ -420,7 +420,9 @@ void PPPMElectrodeIntel::project_psi(IntelBuffers<flt_t, acc_t> *buffers, double
 #endif
   {
     int *mask = atom->mask;
-    const flt_t scaleinv = 1.0 / (nx_pppm * ny_pppm * nz_pppm);
+
+    const bigint ngridtotal = (bigint) nx_pppm * ny_pppm * nz_pppm;
+    const flt_t scaleinv = 1.0 / ngridtotal;
 
     const flt_t lo0 = boxlo[0];
     const flt_t lo1 = boxlo[1];

src/KOKKOS/atom_kokkos.cpp

@@ -44,6 +44,9 @@ AtomKokkos::AtomKokkos(LAMMPS *lmp) : Atom(lmp)
 
   h_tag_min = Kokkos::subview(h_tag_min_max,0);
   h_tag_max = Kokkos::subview(h_tag_min_max,1);
+
+  nprop_atom = 0;
+  fix_prop_atom = nullptr;
 }
 
 /* ---------------------------------------------------------------------- */
@@ -112,6 +115,7 @@ AtomKokkos::~AtomKokkos()
 
   memoryKK->destroy_kokkos(k_dvector, dvector);
   dvector = nullptr;
+  delete [] fix_prop_atom;
 }
 
 /* ---------------------------------------------------------------------- */
@@ -125,11 +129,37 @@ void AtomKokkos::init()
 
 /* ---------------------------------------------------------------------- */
 
+void AtomKokkos::update_property_atom()
+{
+  nprop_atom = 0;
+  std::vector<Fix *> prop_atom_fixes;
+  for (auto &ifix : modify->get_fix_by_style("^property/atom")) {
+    if (!ifix->kokkosable)
+      error->all(FLERR, "KOKKOS package requires a Kokkos-enabled version of fix property/atom");
+
+    ++nprop_atom;
+    prop_atom_fixes.push_back(ifix);
+  }
+
+  delete[] fix_prop_atom;
+  fix_prop_atom = new FixPropertyAtomKokkos *[nprop_atom];
+
+  int n = 0;
+  for (auto &ifix : prop_atom_fixes)
+    fix_prop_atom[n++] = dynamic_cast<FixPropertyAtomKokkos *>(ifix);
+}
+
+/* ---------------------------------------------------------------------- */
+
 void AtomKokkos::sync(const ExecutionSpace space, unsigned int mask)
 {
-  if (space == Device && lmp->kokkos->auto_sync) avecKK->modified(Host, mask);
+  if (space == Device && lmp->kokkos->auto_sync) {
+    avecKK->modified(Host, mask);
+    for (int n = 0; n < nprop_atom; n++) fix_prop_atom[n]->modified(Host, mask);
+  }
 
   avecKK->sync(space, mask);
+  for (int n = 0; n < nprop_atom; n++) fix_prop_atom[n]->sync(space, mask);
 }
 
 /* ---------------------------------------------------------------------- */
@@ -137,13 +167,20 @@ void AtomKokkos::sync(const ExecutionSpace space, unsigned int mask)
 void AtomKokkos::modified(const ExecutionSpace space, unsigned int mask)
 {
   avecKK->modified(space, mask);
+  for (int n = 0; n < nprop_atom; n++) fix_prop_atom[n]->modified(space, mask);
 
-  if (space == Device && lmp->kokkos->auto_sync) avecKK->sync(Host, mask);
+  if (space == Device && lmp->kokkos->auto_sync) {
+    avecKK->sync(Host, mask);
+    for (int n = 0; n < nprop_atom; n++) fix_prop_atom[n]->sync(Host, mask);
+  }
 }
 
+/* ---------------------------------------------------------------------- */
+
 void AtomKokkos::sync_overlapping_device(const ExecutionSpace space, unsigned int mask)
 {
   avecKK->sync_overlapping_device(space, mask);
+  for (int n = 0; n < nprop_atom; n++) fix_prop_atom[n]->sync_overlapping_device(space, mask);
 }
 /* ---------------------------------------------------------------------- */
 
@@ -375,7 +412,7 @@ AtomVec *AtomKokkos::new_avec(const std::string &style, int trysuffix, int &sfla
   int hybrid_substyle_flag = (avec != nullptr);
 
   AtomVec *avec = Atom::new_avec(style, trysuffix, sflag);
-  if (!avec->kokkosable) error->all(FLERR, "KOKKOS package requires a kokkos enabled atom_style");
+  if (!avec->kokkosable) error->all(FLERR, "KOKKOS package requires a Kokkos-enabled atom_style");
 
   if (!hybrid_substyle_flag)
     avecKK = dynamic_cast<AtomVecKokkos*>(avec);

src/KOKKOS/atom_kokkos.h

@@ -14,6 +14,7 @@
 
 #include "atom.h"               // IWYU pragma: export
 #include "kokkos_type.h"
+#include "fix_property_atom_kokkos.h"
 
 #include <Kokkos_Sort.hpp>
 
@@ -25,6 +26,8 @@ namespace LAMMPS_NS {
 class AtomKokkos : public Atom {
  public:
   bool sort_classic;
+  int nprop_atom;
+  FixPropertyAtomKokkos** fix_prop_atom;
 
   DAT::tdual_tagint_1d k_tag;
   DAT::tdual_int_1d k_type, k_mask;
@@ -144,6 +147,7 @@ class AtomKokkos : public Atom {
   }
 
   void init() override;
+  void update_property_atom();
   void allocate_type_arrays() override;
   void sync(const ExecutionSpace space, unsigned int mask);
   void modified(const ExecutionSpace space, unsigned int mask);

src/KOKKOS/atom_vec_dpd_kokkos.cpp

@@ -963,7 +963,6 @@ void AtomVecDPDKokkos::sync(ExecutionSpace space, unsigned int mask)
     if (mask & UCG_MASK) atomKK->k_uCG.sync<LMPDeviceType>();
     if (mask & UCGNEW_MASK) atomKK->k_uCGnew.sync<LMPDeviceType>();
     if (mask & DUCHEM_MASK) atomKK->k_duChem.sync<LMPDeviceType>();
-    if (mask & DVECTOR_MASK) atomKK->k_dvector.sync<LMPDeviceType>();
   } else {
     if (mask & X_MASK) atomKK->k_x.sync<LMPHostType>();
     if (mask & V_MASK) atomKK->k_v.sync<LMPHostType>();
@@ -980,7 +979,6 @@ void AtomVecDPDKokkos::sync(ExecutionSpace space, unsigned int mask)
     if (mask & UCG_MASK) atomKK->k_uCG.sync<LMPHostType>();
     if (mask & UCGNEW_MASK) atomKK->k_uCGnew.sync<LMPHostType>();
     if (mask & DUCHEM_MASK) atomKK->k_duChem.sync<LMPHostType>();
-    if (mask & DVECTOR_MASK) atomKK->k_dvector.sync<LMPHostType>();
   }
 }
 
@@ -1019,8 +1017,6 @@ void AtomVecDPDKokkos::sync_overlapping_device(ExecutionSpace space, unsigned in
       perform_async_copy<DAT::tdual_efloat_1d>(atomKK->k_uCGnew,space);
     if ((mask & DUCHEM_MASK) && atomKK->k_duChem.need_sync<LMPDeviceType>())
       perform_async_copy<DAT::tdual_efloat_1d>(atomKK->k_duChem,space);
-    if ((mask & DVECTOR_MASK) && atomKK->k_dvector.need_sync<LMPDeviceType>())
-      perform_async_copy<DAT::tdual_float_2d>(atomKK->k_dvector,space);
   } else {
     if ((mask & X_MASK) && atomKK->k_x.need_sync<LMPHostType>())
       perform_async_copy<DAT::tdual_x_array>(atomKK->k_x,space);
@@ -1052,8 +1048,6 @@ void AtomVecDPDKokkos::sync_overlapping_device(ExecutionSpace space, unsigned in
       perform_async_copy<DAT::tdual_efloat_1d>(atomKK->k_uCGnew,space);
     if ((mask & DUCHEM_MASK) && atomKK->k_duChem.need_sync<LMPHostType>())
       perform_async_copy<DAT::tdual_efloat_1d>(atomKK->k_duChem,space);
-    if ((mask & DVECTOR_MASK) && atomKK->k_dvector.need_sync<LMPHostType>())
-      perform_async_copy<DAT::tdual_float_2d>(atomKK->k_dvector,space);
   }
 }
 
@@ -1077,7 +1071,6 @@ void AtomVecDPDKokkos::modified(ExecutionSpace space, unsigned int mask)
     if (mask & UCG_MASK) atomKK->k_uCG.modify<LMPDeviceType>();
     if (mask & UCGNEW_MASK) atomKK->k_uCGnew.modify<LMPDeviceType>();
     if (mask & DUCHEM_MASK) atomKK->k_duChem.modify<LMPDeviceType>();
-    if (mask & DVECTOR_MASK) atomKK->k_dvector.modify<LMPDeviceType>();
   } else {
     if (mask & X_MASK) atomKK->k_x.modify<LMPHostType>();
     if (mask & V_MASK) atomKK->k_v.modify<LMPHostType>();
@@ -1094,6 +1087,5 @@ void AtomVecDPDKokkos::modified(ExecutionSpace space, unsigned int mask)
     if (mask & UCG_MASK) atomKK->k_uCG.modify<LMPHostType>();
     if (mask & UCGNEW_MASK) atomKK->k_uCGnew.modify<LMPHostType>();
     if (mask & DUCHEM_MASK) atomKK->k_duChem.modify<LMPHostType>();
-    if (mask & DVECTOR_MASK) atomKK->k_dvector.modify<LMPHostType>();
   }
 }

src/KOKKOS/atom_vec_kokkos.h

@@ -139,6 +139,8 @@ class AtomVecKokkos : virtual public AtomVec {
 
   DAT::tdual_int_1d k_count;
 
+ public:
+
   #ifdef LMP_KOKKOS_GPU
   template<class ViewType>
   Kokkos::View<typename ViewType::data_type,

src/KOKKOS/atom_vec_spin_kokkos.cpp

@@ -586,6 +586,38 @@ int AtomVecSpinKokkos::unpack_exchange_kokkos(DAT::tdual_xfloat_2d &k_buf, int n
   }
 }
 
+/* ----------------------------------------------------------------------
+   clear extra forces starting at atom N
+   nbytes = # of bytes to clear for a per-atom vector
+   include f b/c this is invoked from within SPIN pair styles
+------------------------------------------------------------------------- */
+
+void AtomVecSpinKokkos::force_clear(int n, size_t nbytes)
+{
+  int nzero = (double)nbytes/sizeof(double);
+
+  if (nzero) {
+    atomKK->k_fm.clear_sync_state(); // will be cleared below
+    atomKK->k_fm_long.clear_sync_state(); // will be cleared below
+
+    // local variables for lambda capture
+
+    auto l_fm = atomKK->k_fm.d_view;
+    auto l_fm_long = atomKK->k_fm_long.d_view;
+
+    Kokkos::parallel_for(nzero, LAMMPS_LAMBDA(int i) {
+      l_fm(i,0) = 0.0;
+      l_fm(i,1) = 0.0;
+      l_fm(i,2) = 0.0;
+      l_fm_long(i,0) = 0.0;
+      l_fm_long(i,1) = 0.0;
+      l_fm_long(i,2) = 0.0;
+    });
+
+    atomKK->modified(Device,FM_MASK|FML_MASK);
+  }
+}
+
 /* ---------------------------------------------------------------------- */
 
 void AtomVecSpinKokkos::sync(ExecutionSpace space, unsigned int mask)

src/KOKKOS/atom_vec_spin_kokkos.h

@@ -34,6 +34,7 @@ class AtomVecSpinKokkos : public AtomVecKokkos, public AtomVecSpin {
   AtomVecSpinKokkos(class LAMMPS *);
   void grow(int) override;
   void grow_pointers() override;
+  void force_clear(int, size_t) override;
   void sort_kokkos(Kokkos::BinSort<KeyViewType, BinOp> &Sorter) override;
   int pack_border_kokkos(int n, DAT::tdual_int_2d k_sendlist,
                          DAT::tdual_xfloat_2d buf,int iswap,

src/KOKKOS/comm_kokkos.cpp

@@ -931,6 +931,7 @@ void CommKokkos::exchange_device()
             if (nextrarecv) {
               kkbase->unpack_exchange_kokkos(
                 k_buf_recv,k_indices,nrecv/data_size,
+                nrecv1/data_size,nextrarecv1,
                 ExecutionSpaceFromDevice<DeviceType>::space);
               DeviceType().fence();
             }

src/KOKKOS/fix_dt_reset_kokkos.cpp

@@ -52,14 +52,6 @@ void FixDtResetKokkos<DeviceType>::init()
 {
   FixDtReset::init();
 
-  k_emax = Kokkos::DualView<double*, Kokkos::LayoutRight, DeviceType>("FixDtResetKokkos:gamma", 1);
-
-  k_emax.h_view(0) = emax;
-
-
-  k_emax.template modify<LMPHostType>();
-  k_emax.template sync<DeviceType>();
-
   if (utils::strmatch(update->integrate_style,"^respa"))
     error->all(FLERR,"Cannot (yet) use respa with Kokkos");
 }
@@ -93,43 +85,38 @@ void FixDtResetKokkos<DeviceType>::end_of_step()
 
   MPI_Allreduce(MPI_IN_PLACE, &dt, 1, MPI_DOUBLE, MPI_MIN, world);
 
-  atomKK->modified(execution_space, F_MASK);
+  if (minbound) dt = MAX(dt, tmin);
+  if (maxbound) dt = MIN(dt, tmax);
 
-   if (minbound) dt = MAX(dt, tmin);
-   if (maxbound) dt = MIN(dt, tmax);
+  // if timestep didn't change, just return
+  // else reset update->dt and other classes that depend on it
+  // rRESPA, pair style, fixes
 
-   // if timestep didn't change, just return
-   // else reset update->dt and other classes that depend on it
-   // rRESPA, pair style, fixes
+  if (dt == update->dt) return;
 
-   if (dt == update->dt) return;
+  laststep = update->ntimestep;
 
-   laststep = update->ntimestep;
-
-   // calls to other classes that need to know timestep size changed
-   // similar logic is in Input::timestep()
-
-   update->update_time();
-   update->dt = dt;
-   update->dt_default = 0;
-   if (force->pair) force->pair->reset_dt();
-   for (int i = 0; i < modify->nfix; i++) modify->fix[i]->reset_dt();
-   output->reset_dt();
+  // calls to other classes that need to know timestep size changed
+  // similar logic is in Input::timestep()
 
+  update->update_time();
+  update->dt = dt;
+  update->dt_default = 0;
+  if (force->pair) force->pair->reset_dt();
+  for (auto &ifix : modify->get_fix_list()) ifix->reset_dt();
+  output->reset_dt();
 }
 
 /* ---------------------------------------------------------------------- */
 
 template<class DeviceType>
 KOKKOS_INLINE_FUNCTION
-void FixDtResetKokkos<DeviceType>::operator()(TagFixDtResetMass, const int &i, double &k_dt) const {
+void FixDtResetKokkos<DeviceType>::operator()(TagFixDtResetMass, const int &i, double &dt_min) const {
 
-  double dtv, dtf, dte, dtsq;
+  double dt, dtv, dtf, dte, dtsq;
   double vsq, fsq, massinv;
   double delx, dely, delz, delr;
 
-  double emax = k_emax.d_view(0);
-
   if (mask[i] & groupbit) {
 
     massinv = 1.0 / mass[type[i]];
@@ -138,32 +125,31 @@ void FixDtResetKokkos<DeviceType>::operator()(TagFixDtResetMass, const int &i, d
     dtv = dtf = dte = BIG;
     if (vsq > 0.0) dtv = xmax / sqrt(vsq);
     if (fsq > 0.0) dtf = sqrt(2.0 * xmax / (ftm2v * sqrt(fsq) * massinv));
-    k_dt = MIN(dtv, dtf);
+    dt = MIN(dtv, dtf);
     if ((emax > 0.0) && (fsq * vsq > 0.0)) {
       dte = emax / sqrt(fsq * vsq) / sqrt(ftm2v * mvv2e);
-      k_dt = MIN(dt, dte);
+      dt = MIN(dt, dte);
     }
-    dtsq = k_dt * k_dt;
-    delx = k_dt * v(i,0) + 0.5 * dtsq * massinv * f(i,0) * ftm2v;
-    dely = k_dt * v(i,1) + 0.5 * dtsq * massinv * f(i,1) * ftm2v;
-    delz = k_dt * v(i,2) + 0.5 * dtsq * massinv * f(i,2) * ftm2v;
+    dtsq = dt * dt;
+    delx = dt * v(i,0) + 0.5 * dtsq * massinv * f(i,0) * ftm2v;
+    dely = dt * v(i,1) + 0.5 * dtsq * massinv * f(i,1) * ftm2v;
+    delz = dt * v(i,2) + 0.5 * dtsq * massinv * f(i,2) * ftm2v;
     delr = sqrt(delx * delx + dely * dely + delz * delz);
-    if (delr > xmax) k_dt *= xmax / delr;
+    if (delr > xmax) dt *= xmax / delr;
+    dt_min = MIN(dt_min,dt);
   }
- }
+}
 
 /* ---------------------------------------------------------------------- */
 
 template<class DeviceType>
 KOKKOS_INLINE_FUNCTION
-void FixDtResetKokkos<DeviceType>::operator()(TagFixDtResetRMass, const int &i, double &k_dt) const {
+void FixDtResetKokkos<DeviceType>::operator()(TagFixDtResetRMass, const int &i, double &dt_min) const {
 
-  double dtv, dtf, dte, dtsq;
+  double dt, dtv, dtf, dte, dtsq;
   double vsq, fsq, massinv;
   double delx, dely, delz, delr;
 
-  double emax = k_emax.d_view(0);
-
   if (mask[i] & groupbit) {
 
     massinv = 1.0 / rmass[i];
@@ -172,17 +158,18 @@ void FixDtResetKokkos<DeviceType>::operator()(TagFixDtResetRMass, const int &i,
     dtv = dtf = dte = BIG;
     if (vsq > 0.0) dtv = xmax / sqrt(vsq);
     if (fsq > 0.0) dtf = sqrt(2.0 * xmax / (ftm2v * sqrt(fsq) * massinv));
-    k_dt = MIN(dtv, dtf);
+    dt = MIN(dtv, dtf);
     if ((emax > 0.0) && (fsq * vsq > 0.0)) {
       dte = emax / sqrt(fsq * vsq) / sqrt(ftm2v * mvv2e);
-      k_dt = MIN(dt, dte);
+      dt = MIN(dt, dte);
     }
-    dtsq = k_dt * k_dt;
-    delx = k_dt * v(i,0) + 0.5 * dtsq * massinv * f(i,0) * ftm2v;
-    dely = k_dt * v(i,1) + 0.5 * dtsq * massinv * f(i,1) * ftm2v;
-    delz = k_dt * v(i,2) + 0.5 * dtsq * massinv * f(i,2) * ftm2v;
+    dtsq = dt * dt;
+    delx = dt * v(i,0) + 0.5 * dtsq * massinv * f(i,0) * ftm2v;
+    dely = dt * v(i,1) + 0.5 * dtsq * massinv * f(i,1) * ftm2v;
+    delz = dt * v(i,2) + 0.5 * dtsq * massinv * f(i,2) * ftm2v;
     delr = sqrt(delx * delx + dely * dely + delz * delz);
-    if (delr > xmax) k_dt *= xmax / delr;
+    if (delr > xmax) dt *= xmax / delr;
+    dt_min = MIN(dt_min,dt);
   }
 }
 

src/KOKKOS/fix_enforce2d_kokkos.cpp

@@ -17,15 +17,14 @@
 ------------------------------------------------------------------------- */
 
 #include "fix_enforce2d_kokkos.h"
+
 #include "atom_masks.h"
 #include "atom_kokkos.h"
 #include "comm.h"
 #include "error.h"
 
-
 using namespace LAMMPS_NS;
 
-
 template <class DeviceType>
 FixEnforce2DKokkos<DeviceType>::FixEnforce2DKokkos(LAMMPS *lmp, int narg, char **arg) :
   FixEnforce2D(lmp, narg, arg)
@@ -34,21 +33,16 @@ FixEnforce2DKokkos<DeviceType>::FixEnforce2DKokkos(LAMMPS *lmp, int narg, char *
   atomKK = (AtomKokkos *) atom;
   execution_space = ExecutionSpaceFromDevice<DeviceType>::space;
 
-  datamask_read   = V_MASK | F_MASK | OMEGA_MASK | MASK_MASK
-          | TORQUE_MASK | ANGMOM_MASK;
-
-  datamask_modify = V_MASK | F_MASK | OMEGA_MASK
-          | TORQUE_MASK | ANGMOM_MASK;
+  datamask_read   = V_MASK | F_MASK | OMEGA_MASK | MASK_MASK | TORQUE_MASK | ANGMOM_MASK;
+  datamask_modify = V_MASK | F_MASK | OMEGA_MASK | TORQUE_MASK | ANGMOM_MASK;
 }
 
-
 template <class DeviceType>
 void FixEnforce2DKokkos<DeviceType>::setup(int vflag)
 {
   post_force(vflag);
 }
 
-
 template <class DeviceType>
 void FixEnforce2DKokkos<DeviceType>::post_force(int /*vflag*/)
 {
@@ -66,7 +60,6 @@ void FixEnforce2DKokkos<DeviceType>::post_force(int /*vflag*/)
   if (atomKK->torque_flag)
     torque = atomKK->k_torque.view<DeviceType>();
 
-
   mask = atomKK->k_mask.view<DeviceType>();
 
   int nlocal = atomKK->nlocal;
@@ -125,13 +118,6 @@ void FixEnforce2DKokkos<DeviceType>::post_force(int /*vflag*/)
   copymode = 0;
 
   atomKK->modified(execution_space,datamask_modify);
-
-  for (int m = 0; m < nfixlist; m++) {
-    atomKK->sync(flist[m]->execution_space,flist[m]->datamask_read);
-    flist[m]->enforce2d();
-    atomKK->modified(flist[m]->execution_space,flist[m]->datamask_modify);
-  }
-
 }
 
 

src/KOKKOS/fix_neigh_history_kokkos.cpp

@@ -453,8 +453,12 @@ KOKKOS_INLINE_FUNCTION
 void FixNeighHistoryKokkos<DeviceType>::operator()(TagFixNeighHistoryUnpackExchange, const int &i) const
 {
   int index = d_indices(i);
+
   if (index > -1) {
     int m = (int) d_ubuf(d_buf(i)).i;
+    if (i >= nrecv1)
+      m = nextrarecv1 + (int) d_ubuf(d_buf(nextrarecv1 + i - nrecv1)).i;
+
     int n = (int) d_ubuf(d_buf(m++)).i;
     d_npartner(index) = n;
     for (int p = 0; p < n; p++) {
@@ -471,6 +475,7 @@ void FixNeighHistoryKokkos<DeviceType>::operator()(TagFixNeighHistoryUnpackExcha
 template<class DeviceType>
 void FixNeighHistoryKokkos<DeviceType>::unpack_exchange_kokkos(
   DAT::tdual_xfloat_2d &k_buf, DAT::tdual_int_1d &k_indices, int nrecv,
+  int nrecv1, int nextrarecv1,
   ExecutionSpace /*space*/)
 {
   d_buf = typename AT::t_xfloat_1d_um(
@@ -478,6 +483,9 @@ void FixNeighHistoryKokkos<DeviceType>::unpack_exchange_kokkos(
     k_buf.extent(0)*k_buf.extent(1));
   d_indices = k_indices.view<DeviceType>();
 
+  this->nrecv1 = nrecv1;
+  this->nextrarecv1 = nextrarecv1;
+
   d_npartner = k_npartner.template view<DeviceType>();
   d_partner = k_partner.template view<DeviceType>();
   d_valuepartner = k_valuepartner.template view<DeviceType>();

src/KOKKOS/fix_neigh_history_kokkos.h

@@ -72,12 +72,14 @@ class FixNeighHistoryKokkos : public FixNeighHistory, public KokkosBase {
 
   void unpack_exchange_kokkos(DAT::tdual_xfloat_2d &k_buf,
                               DAT::tdual_int_1d &indices,int nrecv,
+                              int nrecv1,int nrecv1extra,
                               ExecutionSpace space) override;
 
   typename DAT::tdual_int_2d k_firstflag;
   typename DAT::tdual_float_2d k_firstvalue;
 
  private:
+  int nrecv1,nextrarecv1;
   int nlocal,nsend,beyond_contact;
 
   typename AT::t_tagint_1d tag;

src/KOKKOS/fix_property_atom_kokkos.cpp

@@ -30,7 +30,46 @@ FixPropertyAtomKokkos::FixPropertyAtomKokkos(LAMMPS *lmp, int narg, char **arg)
   FixPropertyAtom(lmp, narg, arg)
 {
   atomKK = (AtomKokkos *) atom;
-  grow_arrays(atom->nmax);
+  kokkosable = 1;
+
+  dvector_flag = 0;
+  for (int nv = 0; nv < nvalue; nv++)
+    if (styles[nv] == DVEC) dvector_flag = 1;
+}
+
+/* ---------------------------------------------------------------------- */
+
+void FixPropertyAtomKokkos::post_constructor()
+{
+  atomKK->update_property_atom();
+
+  FixPropertyAtom::post_constructor();
+}
+
+/* ---------------------------------------------------------------------- */
+
+FixPropertyAtomKokkos::~FixPropertyAtomKokkos()
+{
+  // deallocate per-atom vectors in Atom class
+  // set ptrs to a null pointer, so they no longer exist for Atom class
+
+  for (int nv = 0; nv < nvalue; nv++) {
+    if (styles[nv] == MOLECULE) {
+      atom->molecule_flag = 0;
+      memoryKK->destroy_kokkos(atomKK->k_molecule,atom->molecule);
+      atom->molecule = nullptr;
+    } else if (styles[nv] == CHARGE) {
+      atom->q_flag = 0;
+      memoryKK->destroy_kokkos(atomKK->k_q,atom->q);
+      atom->q = nullptr;
+    } else if (styles[nv] == RMASS) {
+      atom->rmass_flag = 0;
+      memoryKK->destroy_kokkos(atomKK->k_rmass,atom->rmass);
+      atom->rmass = nullptr;
+    }
+  }
+
+  atomKK->update_property_atom();
 }
 
 /* ----------------------------------------------------------------------
@@ -44,17 +83,17 @@ void FixPropertyAtomKokkos::grow_arrays(int nmax)
 {
   for (int nv = 0; nv < nvalue; nv++) {
     if (styles[nv] == MOLECULE) {
-      memory->grow(atom->molecule,nmax,"atom:molecule");
-      size_t nbytes = (nmax-nmax_old) * sizeof(tagint);
-      memset(&atom->molecule[nmax_old],0,nbytes);
+      atomKK->sync(Device,MOLECULE_MASK);
+      memoryKK->grow_kokkos(atomKK->k_molecule,atom->molecule,nmax,"atom:molecule");
+      atomKK->modified(Device,MOLECULE_MASK);
     } else if (styles[nv] == CHARGE) {
-      memory->grow(atom->q,nmax,"atom:q");
-      size_t nbytes = (nmax-nmax_old) * sizeof(double);
-      memset(&atom->q[nmax_old],0,nbytes);
+      atomKK->sync(Device,Q_MASK);
+      memoryKK->grow_kokkos(atomKK->k_q,atom->q,nmax,"atom:q");
+      atomKK->modified(Device,Q_MASK);
     } else if (styles[nv] == RMASS) {
-      memory->grow(atom->rmass,nmax,"atom:rmass");
-      size_t nbytes = (nmax-nmax_old) * sizeof(double);
-      memset(&atom->rmass[nmax_old],0,nbytes);
+      atomKK->sync(Device,RMASS_MASK);
+      memoryKK->grow_kokkos(atomKK->k_rmass,atom->rmass,nmax,"atom:rmass");
+      atomKK->modified(Device,RMASS_MASK);
     } else if (styles[nv] == TEMPERATURE) {
       memory->grow(atom->temperature, nmax, "atom:temperature");
       size_t nbytes = (nmax - nmax_old) * sizeof(double);
@@ -69,7 +108,7 @@ void FixPropertyAtomKokkos::grow_arrays(int nmax)
       memset(&atom->ivector[index[nv]][nmax_old],0,nbytes);
     } else if (styles[nv] == DVEC) {
       atomKK->sync(Device,DVECTOR_MASK);
-      memoryKK->grow_kokkos(atomKK->k_dvector,atomKK->dvector,atomKK->k_dvector.extent(0),nmax,
+      memoryKK->grow_kokkos(atomKK->k_dvector,atom->dvector,atomKK->k_dvector.extent(0),nmax,
                           "atom:dvector");
       atomKK->modified(Device,DVECTOR_MASK);
     } else if (styles[nv] == IARRAY) {
@@ -84,3 +123,62 @@ void FixPropertyAtomKokkos::grow_arrays(int nmax)
   }
   nmax_old = nmax;
 }
+
+/* ---------------------------------------------------------------------- */
+
+void FixPropertyAtomKokkos::sync(ExecutionSpace space, unsigned int mask)
+{
+  if (space == Device) {
+    if (molecule_flag && (mask & MOLECULE_MASK)) atomKK->k_molecule.sync<LMPDeviceType>();
+    if (q_flag && (mask & Q_MASK)) atomKK->k_q.sync<LMPDeviceType>();
+    if (rmass_flag && (mask & RMASS_MASK)) {atomKK->k_rmass.sync<LMPDeviceType>();}
+    if (dvector_flag && (mask & DVECTOR_MASK)) atomKK->k_dvector.sync<LMPDeviceType>();
+  } else {
+    if (molecule_flag && (mask & MOLECULE_MASK)) atomKK->k_molecule.sync<LMPHostType>();
+    if (q_flag && (mask & Q_MASK)) atomKK->k_q.sync<LMPHostType>();
+    if (rmass_flag && (mask & RMASS_MASK)) atomKK->k_rmass.sync<LMPHostType>();
+    if (dvector_flag && (mask & DVECTOR_MASK)) atomKK->k_dvector.sync<LMPHostType>();
+  }
+}
+
+/* ---------------------------------------------------------------------- */
+
+void FixPropertyAtomKokkos::sync_overlapping_device(ExecutionSpace space, unsigned int mask)
+{
+  if (space == Device) {
+    if ((mask & MOLECULE_MASK) && atomKK->k_molecule.need_sync<LMPDeviceType>())
+      atomKK->avecKK->perform_async_copy<DAT::tdual_tagint_1d>(atomKK->k_molecule,space);
+    if ((mask & Q_MASK) && atomKK->k_q.need_sync<LMPDeviceType>())
+      atomKK->avecKK->perform_async_copy<DAT::tdual_float_1d>(atomKK->k_q,space);
+    if ((mask & RMASS_MASK) && atomKK->k_rmass.need_sync<LMPDeviceType>())
+      atomKK->avecKK->perform_async_copy<DAT::tdual_float_1d>(atomKK->k_rmass,space);
+    if ((mask & DVECTOR_MASK) && atomKK->k_dvector.need_sync<LMPDeviceType>())
+      atomKK->avecKK->perform_async_copy<DAT::tdual_float_2d>(atomKK->k_dvector,space);
+  } else {
+    if ((mask & MOLECULE_MASK) && atomKK->k_molecule.need_sync<LMPHostType>())
+      atomKK->avecKK->perform_async_copy<DAT::tdual_tagint_1d>(atomKK->k_molecule,space);
+    if ((mask & Q_MASK) && atomKK->k_q.need_sync<LMPHostType>())
+      atomKK->avecKK->perform_async_copy<DAT::tdual_float_1d>(atomKK->k_q,space);
+    if ((mask & RMASS_MASK) && atomKK->k_rmass.need_sync<LMPHostType>())
+      atomKK->avecKK->perform_async_copy<DAT::tdual_float_1d>(atomKK->k_rmass,space);
+    if ((mask & DVECTOR_MASK) && atomKK->k_dvector.need_sync<LMPHostType>())
+      atomKK->avecKK->perform_async_copy<DAT::tdual_float_2d>(atomKK->k_dvector,space);
+  }
+}
+
+/* ---------------------------------------------------------------------- */
+
+void FixPropertyAtomKokkos::modified(ExecutionSpace space, unsigned int mask)
+{
+  if (space == Device) {
+    if (molecule_flag && (mask & MOLECULE_MASK)) atomKK->k_molecule.modify<LMPDeviceType>();
+    if (q_flag && (mask & Q_MASK)) atomKK->k_q.modify<LMPDeviceType>();
+    if (rmass_flag && (mask & RMASS_MASK)) atomKK->k_rmass.modify<LMPDeviceType>();
+    if (dvector_flag && (mask & DVECTOR_MASK)) atomKK->k_dvector.modify<LMPDeviceType>();
+  } else {
+    if (molecule_flag && (mask & MOLECULE_MASK)) atomKK->k_molecule.modify<LMPHostType>();
+    if (q_flag && (mask & Q_MASK)) atomKK->k_q.modify<LMPHostType>();
+    if (rmass_flag && (mask & RMASS_MASK)) atomKK->k_rmass.modify<LMPHostType>();
+    if (dvector_flag && (mask & DVECTOR_MASK)) atomKK->k_dvector.modify<LMPHostType>();
+  }
+}

src/KOKKOS/fix_property_atom_kokkos.h

@@ -22,14 +22,23 @@ FixStyle(property/atom/kk,FixPropertyAtomKokkos);
 #define LMP_FIX_PROPERTY_ATOM_KOKKOS_H
 
 #include "fix_property_atom.h"
+#include "atom_vec_kokkos.h"
 
 namespace LAMMPS_NS {
 
 class FixPropertyAtomKokkos : public FixPropertyAtom {
  public:
   FixPropertyAtomKokkos(class LAMMPS *, int, char **);
-
+  void post_constructor() override;
+  ~FixPropertyAtomKokkos() override;
   void grow_arrays(int) override;
+
+  void sync(ExecutionSpace space, unsigned int mask);
+  void modified(ExecutionSpace space, unsigned int mask);
+  void sync_overlapping_device(ExecutionSpace space, unsigned int mask);
+
+ private:
+  int dvector_flag;
 };
 
 }

src/KOKKOS/fix_qeq_reaxff_kokkos.cpp

@@ -1416,6 +1416,7 @@ KOKKOS_INLINE_FUNCTION
 void FixQEqReaxFFKokkos<DeviceType>::operator()(TagQEqUnpackExchange, const int &i) const
 {
   int index = d_indices(i);
+
   if (index > -1) {
     for (int m = 0; m < nprev; m++) d_s_hist(index,m) = d_buf(i*nprev*2 + m);
     for (int m = 0; m < nprev; m++) d_t_hist(index,m) = d_buf(i*nprev*2 + nprev+m);
@@ -1427,6 +1428,7 @@ void FixQEqReaxFFKokkos<DeviceType>::operator()(TagQEqUnpackExchange, const int
 template <class DeviceType>
 void FixQEqReaxFFKokkos<DeviceType>::unpack_exchange_kokkos(
   DAT::tdual_xfloat_2d &k_buf, DAT::tdual_int_1d &k_indices, int nrecv,
+  int /*nrecv1*/, int /*nextrarecv1*/,
   ExecutionSpace /*space*/)
 {
   k_buf.sync<DeviceType>();

src/KOKKOS/fix_qeq_reaxff_kokkos.h

@@ -143,6 +143,7 @@ class FixQEqReaxFFKokkos : public FixQEqReaxFF, public KokkosBase {
 
   void unpack_exchange_kokkos(DAT::tdual_xfloat_2d &k_buf,
                               DAT::tdual_int_1d &indices,int nrecv,
+                              int nrecv1,int nextrarecv1,
                               ExecutionSpace space) override;
 
   struct params_qeq{

src/KOKKOS/fix_shake_kokkos.cpp

@@ -525,7 +525,7 @@ void FixShakeKokkos<DeviceType>::operator()(TagFixShakePostForce<NEIGHFLAG,EVFLA
 ------------------------------------------------------------------------- */
 
 template<class DeviceType>
-int FixShakeKokkos<DeviceType>::dof(int igroup)
+bigint FixShakeKokkos<DeviceType>::dof(int igroup)
 {
   d_mask = atomKK->k_mask.view<DeviceType>();
   d_tag = atomKK->k_tag.view<DeviceType>();
@@ -538,7 +538,7 @@ int FixShakeKokkos<DeviceType>::dof(int igroup)
   // count dof in a cluster if and only if
   // the central atom is in group and atom i is the central atom
 
-  int n = 0;
+  bigint n = 0;
   {
     // local variables for lambda capture
 
@@ -549,7 +549,7 @@ int FixShakeKokkos<DeviceType>::dof(int igroup)
     auto groupbit = group->bitmask[igroup];
 
     Kokkos::parallel_reduce(Kokkos::RangePolicy<DeviceType>(0,nlocal),
-     LAMMPS_LAMBDA(const int& i, int& n) {
+     LAMMPS_LAMBDA(const int& i, bigint& n) {
       if (!(mask[i] & groupbit)) return;
       if (d_shake_flag[i] == 0) return;
       if (d_shake_atom(i,0) != tag[i]) return;
@@ -560,8 +560,8 @@ int FixShakeKokkos<DeviceType>::dof(int igroup)
     },n);
   }
 
-  int nall;
-  MPI_Allreduce(&n,&nall,1,MPI_INT,MPI_SUM,world);
+  bigint nall;
+  MPI_Allreduce(&n,&nall,1,MPI_LMP_BIGINT,MPI_SUM,world);
   return nall;
 }
 
@@ -1581,8 +1581,8 @@ void FixShakeKokkos<DeviceType>::pack_exchange_item(const int &mysend, int &offs
     else offset++;
   } else {
 
-    d_buf[mysend] = nsend + offset;
     int m = nsend + offset;
+    d_buf[mysend] = m;
     d_buf[m++] = flag;
     if (flag == 1) {
       d_buf[m++] = d_shake_atom(i,0);
@@ -1703,6 +1703,8 @@ void FixShakeKokkos<DeviceType>::operator()(TagFixShakeUnpackExchange, const int
 
   if (index > -1) {
     int m = d_buf[i];
+    if (i >= nrecv1)
+      m = nextrarecv1 + d_buf[nextrarecv1 + i - nrecv1];
 
     int flag = d_shake_flag[index] = static_cast<int> (d_buf[m++]);
     if (flag == 1) {
@@ -1739,6 +1741,7 @@ void FixShakeKokkos<DeviceType>::operator()(TagFixShakeUnpackExchange, const int
 template <class DeviceType>
 void FixShakeKokkos<DeviceType>::unpack_exchange_kokkos(
   DAT::tdual_xfloat_2d &k_buf, DAT::tdual_int_1d &k_indices, int nrecv,
+  int nrecv1, int nextrarecv1,
   ExecutionSpace /*space*/)
 {
   k_buf.sync<DeviceType>();
@@ -1749,6 +1752,9 @@ void FixShakeKokkos<DeviceType>::unpack_exchange_kokkos(
     k_buf.extent(0)*k_buf.extent(1));
   d_indices = k_indices.view<DeviceType>();
 
+  this->nrecv1 = nrecv1;
+  this->nextrarecv1 = nextrarecv1;
+
   k_shake_flag.template sync<DeviceType>();
   k_shake_atom.template sync<DeviceType>();
   k_shake_type.template sync<DeviceType>();

src/KOKKOS/fix_shake_kokkos.h

@@ -44,8 +44,6 @@ struct TagFixShakeUnpackExchange{};
 template<class DeviceType>
 class FixShakeKokkos : public FixShake, public KokkosBase {
 
- //friend class FixEHEX;
-
  public:
   typedef DeviceType device_type;
   typedef EV_FLOAT value_type;
@@ -77,7 +75,7 @@ class FixShakeKokkos : public FixShake, public KokkosBase {
   void shake_end_of_step(int vflag) override;
   void correct_coordinates(int vflag) override;
 
-  int dof(int) override;
+  bigint dof(int) override;
 
   void unconstrained_update() override;
 
@@ -112,9 +110,12 @@ class FixShakeKokkos : public FixShake, public KokkosBase {
 
   void unpack_exchange_kokkos(DAT::tdual_xfloat_2d &k_buf,
                               DAT::tdual_int_1d &indices,int nrecv,
+                              int nrecv1,int nrecv1extra,
                               ExecutionSpace space) override;
 
  protected:
+  int nrecv1,nextrarecv1;
+
   typename AT::t_x_array d_x;
   typename AT::t_v_array d_v;
   typename AT::t_f_array d_f;
@@ -259,4 +260,3 @@ struct FixShakeKokkosPackExchangeFunctor {
 
 #endif
 #endif
-

src/KOKKOS/fix_wall_gran_kokkos.cpp

@@ -419,6 +419,7 @@ void FixWallGranKokkos<DeviceType>::operator()(TagFixWallGranUnpackExchange, con
 template<class DeviceType>
 void FixWallGranKokkos<DeviceType>::unpack_exchange_kokkos(
   DAT::tdual_xfloat_2d &k_buf, DAT::tdual_int_1d &k_indices, int nrecv,
+  int /*nrecv1*/, int /*nextrarecv1*/,
   ExecutionSpace /*space*/)
 {
   d_buf = typename ArrayTypes<DeviceType>::t_xfloat_1d_um(
@@ -430,7 +431,6 @@ void FixWallGranKokkos<DeviceType>::unpack_exchange_kokkos(
 
   copymode = 1;
 
-
   Kokkos::parallel_for(Kokkos::RangePolicy<DeviceType,TagFixWallGranUnpackExchange>(0,nrecv),*this);
 
   copymode = 0;

src/KOKKOS/fix_wall_gran_kokkos.h

@@ -62,12 +62,13 @@ class FixWallGranKokkos : public FixWallGranOld, public KokkosBase {
   void operator()(TagFixWallGranUnpackExchange, const int&) const;
 
   int pack_exchange_kokkos(const int &nsend,DAT::tdual_xfloat_2d &buf,
-			   DAT::tdual_int_1d k_sendlist,
-			   DAT::tdual_int_1d k_copylist,
-			   ExecutionSpace space) override;
+                           DAT::tdual_int_1d k_sendlist,
+                           DAT::tdual_int_1d k_copylist,
+                           ExecutionSpace space) override;
 
   void unpack_exchange_kokkos(DAT::tdual_xfloat_2d &k_buf,
                               DAT::tdual_int_1d &indices,int nrecv,
+                              int nrecv1,int nrecv1extra,
                               ExecutionSpace space) override;
 
  private:
@@ -91,6 +92,7 @@ class FixWallGranKokkos : public FixWallGranOld, public KokkosBase {
   typename AT::t_int_1d d_copylist;
   typename AT::t_int_1d d_indices;
 };
+
 }
 
 #endif

src/KOKKOS/grid3d_kokkos.cpp

@@ -640,7 +640,7 @@ void Grid3dKokkos<DeviceType>::forward_comm(int caller, void *ptr, int which, in
                             MPI_Datatype datatype)
 {
   if (caller == KSPACE) {
-    if (layout != Comm::LAYOUT_TILED)
+    if (comm->layout != Comm::LAYOUT_TILED)
     forward_comm_kspace_brick((KSpace *) ptr,which,nper,k_buf1,k_buf2,datatype);
   else
     forward_comm_kspace_tiled((KSpace *) ptr,which,nper,k_buf1,k_buf2,datatype);
@@ -780,7 +780,7 @@ void Grid3dKokkos<DeviceType>::reverse_comm(int caller, void *ptr, int which, in
                             MPI_Datatype datatype)
 {
   if (caller == KSPACE) {
-    if (layout != Comm::LAYOUT_TILED)
+    if (comm->layout != Comm::LAYOUT_TILED)
       reverse_comm_kspace_brick((KSpace *) ptr,which,nper,k_buf1,k_buf2,datatype);
     else
       reverse_comm_kspace_tiled((KSpace *) ptr,which,nper,k_buf1,k_buf2,datatype);

src/KOKKOS/kokkos.cpp

@@ -104,7 +104,9 @@ KokkosLMP::KokkosLMP(LAMMPS *lmp, int narg, char **arg) : Pointers(lmp)
 
   int me = 0;
   MPI_Comm_rank(world,&me);
-  if (me == 0) error->message(FLERR,"KOKKOS mode is enabled");
+  if (me == 0)
+    error->message(FLERR,"KOKKOS mode with Kokkos version {}.{}.{} is enabled",
+                   KOKKOS_VERSION / 10000, (KOKKOS_VERSION % 10000) / 100, KOKKOS_VERSION % 100);
 
   // process any command-line args that invoke Kokkos settings
 
@@ -143,6 +145,14 @@ KokkosLMP::KokkosLMP(LAMMPS *lmp, int narg, char **arg) : Pointers(lmp)
         if (device >= skip_gpu) device++;
         set_flag = 1;
       }
+      if ((str = getenv("FLUX_TASK_LOCAL_ID"))) {
+        if (ngpus > 0) {
+          int local_rank = atoi(str);
+          device = local_rank % ngpus;
+          if (device >= skip_gpu) device++;
+          set_flag = 1;
+        }
+      }
       if ((str = getenv("MPT_LRANK"))) {
         if (ngpus > 0) {
           int local_rank = atoi(str);

src/KOKKOS/kokkos_base.h

@@ -41,19 +41,18 @@ class KokkosBase {
                                            int, int *) {return 0;};
   virtual void unpack_forward_comm_fix_kokkos(int, int, DAT::tdual_xfloat_1d &) {}
 
-
-  // Region
-  virtual void match_all_kokkos(int, DAT::tdual_int_1d) {}
-
-  // Fix
   virtual int pack_exchange_kokkos(const int & /*nsend*/, DAT::tdual_xfloat_2d & /*k_buf*/,
                                    DAT::tdual_int_1d /*k_sendlist*/,
                                    DAT::tdual_int_1d /*k_copylist*/,
                                    ExecutionSpace /*space*/) { return 0; }
   virtual void unpack_exchange_kokkos(DAT::tdual_xfloat_2d & /*k_buf*/,
                                       DAT::tdual_int_1d & /*indices*/, int /*nrecv*/,
+                                      int /*nrecv1*/, int /*nextrarecv1*/,
                                       ExecutionSpace /*space*/) {}
 
+  // Region
+  virtual void match_all_kokkos(int, DAT::tdual_int_1d) {}
+
   using KeyViewType = DAT::t_x_array;
   using BinOp = BinOp3DLAMMPS<KeyViewType>;
   virtual void

src/KOKKOS/meam_force_kokkos.h

@@ -131,7 +131,7 @@ KOKKOS_INLINE_FUNCTION void MEAMKokkos<DeviceType>::operator()(TagMEAMForce<NEIG
   double drho3mdr1, drho3mdr2, drho3mds1, drho3mds2;
   double drho3mdrm1[3], drho3mdrm2[3];
 
-  // The f, etc. arrays are duplicated for OpenMP, atomic for CUDA, and neither for Serial
+  // The f, etc. arrays are duplicated for OpenMP, atomic for GPU, and neither for Serial
 
   auto v_f =
       ScatterViewHelper<NeedDup_v<NEIGHFLAG, DeviceType>, decltype(dup_f), decltype(ndup_f)>::get(
@@ -601,8 +601,31 @@ KOKKOS_INLINE_FUNCTION void MEAMKokkos<DeviceType>::operator()(TagMEAMForce<NEIG
             drho2mds2 = a2 * rhoa2mi * arg1j2m - 2.0 / 3.0 * d_arho2mb[j] * rhoa2mi;
             drho3mds1 = a3 * rhoa3mj * arg1i3m - a3a * rhoa3mj * arg3i3m;
             drho3mds2 = a3 * rhoa3mi * arg1j3m - a3a * rhoa3mi * arg3j3m;
+            drho1mds1 *= -1;
+            drho1mds2 *= -1;
             drho3mds1 *= -1;
             drho3mds2 *= -1;
+
+            t1i = 1.0;
+            t2i = 1.0;
+            t3i = 1.0;
+            t1j = 1.0;
+            t2j = 1.0;
+            t3j = 1.0;
+            dt1dr1 = 0.0;
+            dt1dr2 = 0.0;
+            dt2dr1 = 0.0;
+            dt2dr2 = 0.0;
+            dt3dr1 = 0.0;
+            dt3dr2 = 0.0;
+
+            // these formulae are simplifed by substituting t=1, dt=0 from above
+
+            drhods1 = d_dgamma1[i] * drho0ds1 + d_dgamma2[i]
+              * ((drho1ds1 - drho1mds1) + (drho2ds1 - drho2mds1) + (drho3ds1 - drho3mds1));
+            drhods2 = d_dgamma1[j] * drho0ds2 + d_dgamma2[j]
+              * ((drho1ds2 - drho1mds2) + (drho2ds2 - drho2mds2) + (drho3ds2 - drho3mds2));
+
           } else {
             drho1mds1 = 0.0;
             drho1mds2 = 0.0;
@@ -610,61 +633,49 @@ KOKKOS_INLINE_FUNCTION void MEAMKokkos<DeviceType>::operator()(TagMEAMForce<NEIG
             drho2mds2 = 0.0;
             drho3mds1 = 0.0;
             drho3mds2 = 0.0;
-          }
 
-          if (ialloy == 1) {
-            a1i = fdiv_zero_kk(rhoa0j, d_tsq_ave(i, 0));
-            a1j = fdiv_zero_kk(rhoa0i, d_tsq_ave(j, 0));
-            a2i = fdiv_zero_kk(rhoa0j, d_tsq_ave(i, 1));
-            a2j = fdiv_zero_kk(rhoa0i, d_tsq_ave(j, 1));
-            a3i = fdiv_zero_kk(rhoa0j, d_tsq_ave(i, 2));
-            a3j = fdiv_zero_kk(rhoa0i, d_tsq_ave(j, 2));
+            if (ialloy == 1) {
 
-            dt1ds1 = a1i * (t1mj - t1i * MathSpecialKokkos::square(t1mj));
-            dt1ds2 = a1j * (t1mi - t1j * MathSpecialKokkos::square(t1mi));
-            dt2ds1 = a2i * (t2mj - t2i * MathSpecialKokkos::square(t2mj));
-            dt2ds2 = a2j * (t2mi - t2j * MathSpecialKokkos::square(t2mi));
-            dt3ds1 = a3i * (t3mj - t3i * MathSpecialKokkos::square(t3mj));
-            dt3ds2 = a3j * (t3mi - t3j * MathSpecialKokkos::square(t3mi));
+              a1i = fdiv_zero_kk(rhoa0j, d_tsq_ave(i, 0));
+              a1j = fdiv_zero_kk(rhoa0i, d_tsq_ave(j, 0));
+              a2i = fdiv_zero_kk(rhoa0j, d_tsq_ave(i, 1));
+              a2j = fdiv_zero_kk(rhoa0i, d_tsq_ave(j, 1));
+              a3i = fdiv_zero_kk(rhoa0j, d_tsq_ave(i, 2));
+              a3j = fdiv_zero_kk(rhoa0i, d_tsq_ave(j, 2));
 
-          } else if (ialloy == 2) {
+              dt1ds1 = a1i * (t1mj - t1i * MathSpecialKokkos::square(t1mj));
+              dt1ds2 = a1j * (t1mi - t1j * MathSpecialKokkos::square(t1mi));
+              dt2ds1 = a2i * (t2mj - t2i * MathSpecialKokkos::square(t2mj));
+              dt2ds2 = a2j * (t2mi - t2j * MathSpecialKokkos::square(t2mi));
+              dt3ds1 = a3i * (t3mj - t3i * MathSpecialKokkos::square(t3mj));
+              dt3ds2 = a3j * (t3mi - t3j * MathSpecialKokkos::square(t3mi));
 
-            dt1ds1 = 0.0;
-            dt1ds2 = 0.0;
-            dt2ds1 = 0.0;
-            dt2ds2 = 0.0;
-            dt3ds1 = 0.0;
-            dt3ds2 = 0.0;
+            } else if (ialloy == 2) {
 
-          } else {
+              dt1ds1 = 0.0;
+              dt1ds2 = 0.0;
+              dt2ds1 = 0.0;
+              dt2ds2 = 0.0;
+              dt3ds1 = 0.0;
+              dt3ds2 = 0.0;
 
-            ai = 0.0;
-            if (!iszero_kk(d_rho0[i])) ai = rhoa0j / d_rho0[i];
-            aj = 0.0;
-            if (!iszero_kk(d_rho0[j])) aj = rhoa0i / d_rho0[j];
+            } else {
 
-            dt1ds1 = ai * (t1mj - t1i);
-            dt1ds2 = aj * (t1mi - t1j);
-            dt2ds1 = ai * (t2mj - t2i);
-            dt2ds2 = aj * (t2mi - t2j);
-            dt3ds1 = ai * (t3mj - t3i);
-            dt3ds2 = aj * (t3mi - t3j);
-          }
+              ai = 0.0;
+              if (!iszero_kk(d_rho0[i])) ai = rhoa0j / d_rho0[i];
+              aj = 0.0;
+              if (!iszero_kk(d_rho0[j])) aj = rhoa0i / d_rho0[j];
+
+              dt1ds1 = ai * (t1mj - t1i);
+              dt1ds2 = aj * (t1mi - t1j);
+              dt2ds1 = ai * (t2mj - t2i);
+              dt2ds2 = aj * (t2mi - t2j);
+              dt3ds1 = ai * (t3mj - t3i);
+              dt3ds2 = aj * (t3mi - t3j);
+            }
 
-          if (msmeamflag) {
             drhods1 = d_dgamma1[i] * drho0ds1 +
-              d_dgamma2[i] * (dt1ds1 * d_rho1[i] + t1i * (drho1ds1 - drho1mds1) +
-                              dt2ds1 * d_rho2[i] + t2i * (drho2ds1 - drho2mds1) +
-                              dt3ds1 * d_rho3[i] + t3i * (drho3ds1 - drho3mds1)) -
-              d_dgamma3[i] * (shpi[0] * dt1ds1 + shpi[1] * dt2ds1 + shpi[2] * dt3ds1);
-            drhods2 = d_dgamma1[j] * drho0ds2 +
-              d_dgamma2[j] * (dt1ds2 * d_rho1[j] + t1j * (drho1ds2 - drho1mds2) +
-                              dt2ds2 * d_rho2[j] + t2j * (drho2ds2 - drho2mds2) +
-                              dt3ds2 * d_rho3[j] + t3j * (drho3ds2 - drho3mds2)) -
-              d_dgamma3[j] * (shpj[0] * dt1ds2 + shpj[1] * dt2ds2 + shpj[2] * dt3ds2);
-          } else {
-            drhods1 = d_dgamma1[i] * drho0ds1 +
-                d_dgamma2[i] *
+              d_dgamma2[i] *
                     (dt1ds1 * d_rho1[i] + t1i * drho1ds1 + dt2ds1 * d_rho2[i] + t2i * drho2ds1 +
                      dt3ds1 * d_rho3[i] + t3i * drho3ds1) -
                 d_dgamma3[i] * (shpi[0] * dt1ds1 + shpi[1] * dt2ds1 + shpi[2] * dt3ds1);

src/KOKKOS/min_kokkos.cpp

@@ -59,6 +59,9 @@ void MinKokkos::init()
 {
   Min::init();
 
+  if (!fix_minimize->kokkosable)
+    error->all(FLERR,"KOKKOS package requires fix minimize/kk");
+
   fix_minimize_kk = (FixMinimizeKokkos*) fix_minimize;
 }
 
@@ -69,8 +72,7 @@ void MinKokkos::init()
 void MinKokkos::setup(int flag)
 {
   if (comm->me == 0 && screen) {
-    fmt::print(screen,"Setting up {} style minimization ...\n",
-               update->minimize_style);
+    fmt::print(screen,"Setting up {} style minimization ...\n", update->minimize_style);
     if (flag) {
       fmt::print(screen,"  Unit style    : {}\n", update->unit_style);
       fmt::print(screen,"  Current step  : {}\n", update->ntimestep);
@@ -89,14 +91,13 @@ void MinKokkos::setup(int flag)
     fextra = new double[nextra_global];
     if (comm->me == 0)
       error->warning(FLERR, "Energy due to {} extra global DOFs will"
-                     " be included in minimizer energies\n", nextra_global);
+                     " be included in minimizer energies\n",nextra_global);
   }
 
   // compute for potential energy
 
-  int id = modify->find_compute("thermo_pe");
-  if (id < 0) error->all(FLERR,"Minimization could not find thermo_pe compute");
-  pe_compute = modify->compute[id];
+  pe_compute = modify->get_compute_by_id("thermo_pe");
+  if (!pe_compute) error->all(FLERR,"Minimization could not find thermo_pe compute");
 
   // style-specific setup does two tasks
   // setup extra global dof vectors
@@ -534,6 +535,7 @@ double MinKokkos::energy_force(int resetflag)
     if (resetflag) fix_minimize_kk->reset_coords();
     reset_vectors();
   }
+
   return energy;
 }
 
@@ -572,7 +574,14 @@ void MinKokkos::force_clear()
         l_torque(i,2) = 0.0;
       }
     });
+
+    if (extraflag) {
+      size_t nbytes = sizeof(double) * atom->nlocal;
+      if (force->newton) nbytes += sizeof(double) * atom->nghost;
+      atom->avec->force_clear(0,nbytes);
+    }
   }
+
   atomKK->modified(Device,F_MASK|TORQUE_MASK);
 }
 

src/KOKKOS/mliap_data_kokkos.cpp

@@ -85,6 +85,7 @@ void MLIAPDataKokkos<DeviceType>::generate_neighdata(class NeighList *list_in, i
   // clear gradforce and elems arrays
 
   int nall = atom->nlocal + atom->nghost;
+  nlocal = atom->nlocal;
   ntotal = nall;
   if (gradgradflag > -1){
     auto d_gradforce = k_gradforce.template view<DeviceType>();

src/KOKKOS/mliap_data_kokkos.h

@@ -118,6 +118,7 @@ public:
     egradient(nullptr),
     ntotal(base.ntotal),
     nlistatoms(base.nlistatoms),
+    nlocal(base.nlocal),
     natomneigh(base.natomneigh),
     numneighs(base.numneighs),
     iatoms(base.k_iatoms.d_view.data()),
@@ -171,6 +172,7 @@ public:
   // Neighborlist stuff
   const int ntotal;
   const int nlistatoms;
+  const int nlocal;
   const int natomneigh;
   int *numneighs;
   int *iatoms;
@@ -191,7 +193,7 @@ public:
   int dev;
 
 #ifdef LMP_KOKKOS_GPU
-  MLIAPDataKokkosDevice(MLIAPDataKokkos<LMPHostType> &base) : ndescriptors(-1),nparams(-1),nelements(-1),ntotal(-1),nlistatoms(-1),natomneigh(-1),
+  MLIAPDataKokkosDevice(MLIAPDataKokkos<LMPHostType> &base) : ndescriptors(-1),nparams(-1),nelements(-1),ntotal(-1),nlistatoms(-1),nlocal(-1),natomneigh(-1),
       nneigh_max(-1),npairs(-1)
   {
     // It cannot get here, but needed for compilation

src/KOKKOS/mliap_model_python_couple_kokkos.pyx

@@ -25,6 +25,7 @@ cdef extern from "mliap_data_kokkos.h" namespace "LAMMPS_NS":
     cdef cppclass MLIAPDataKokkosDevice:
         # Array shapes
         int nlistatoms
+        int nlocal
         int ndescriptors
 
         # Input data
@@ -130,14 +131,14 @@ cdef create_array(device, void *pointer, shape,is_int):
                 return numpy.asarray(<int[:shape[0],:shape[1]]>pointer)
             else:
                 return numpy.asarray(<double[:shape[0],:shape[1]]>pointer)
-                
+
 
 cdef public void MLIAPPYKokkos_compute_gradients(MLIAPModelPythonKokkosDevice * c_model, MLIAPDataKokkosDevice * data) with gil:
 
     dev=data.dev
 
     torch.cuda.nvtx.range_push("set data fields")
-    model = retrieve(c_model) 
+    model = retrieve(c_model)
     n_d = data.ndescriptors
     n_a = data.nlistatoms
 
@@ -148,7 +149,7 @@ cdef public void MLIAPPYKokkos_compute_gradients(MLIAPModelPythonKokkosDevice *
     beta_cp = create_array(dev, data.betas, (n_a, n_d), False)
     desc_cp = create_array(dev, data.descriptors, (n_a, n_d), False)
     torch.cuda.nvtx.range_pop()
-    
+
     # Invoke python model on numpy arrays.
     torch.cuda.nvtx.range_push("call model")
     model(elem_cp,desc_cp,beta_cp,en_cp,dev==1)

src/KOKKOS/mliap_unified_couple_kokkos.pyx

@@ -59,6 +59,7 @@ cdef extern from "mliap_data_kokkos.h" namespace "LAMMPS_NS":
 
         int ntotal              # total number of owned and ghost atoms on this proc
         int nlistatoms          # current number of atoms in local atom lists
+        int nlocal
         int natomneigh          # current number of atoms and ghosts in atom neighbor arrays
         int * numneighs         # neighbors count for each atom
         int * iatoms            # index of each atom
@@ -133,14 +134,14 @@ cdef create_array(device, void *pointer, shape,is_int):
                 return np.asarray(<int[:shape[0],:shape[1]]>pointer)
             else:
                 return np.asarray(<double[:shape[0],:shape[1]]>pointer)
-            
+
 
 
 # Cython implementation of MLIAPData
 # Automatically converts between C arrays and numpy when needed
 cdef class MLIAPDataPy:
     cdef MLIAPDataKokkosDevice * data
-    
+
     def __cinit__(self):
         self.data = NULL
 
@@ -157,7 +158,7 @@ cdef class MLIAPDataPy:
             ptr = eij.data.ptr
         except:
             ptr = eij.data_ptr()
-        update_pair_energy(self.data, <double*>ptr)  
+        update_pair_energy(self.data, <double*>ptr)
     def update_pair_energy(self, eij):
         if self.data.dev==0:
             self.update_pair_energy_cpu(eij)
@@ -177,7 +178,7 @@ cdef class MLIAPDataPy:
             ptr = fij.data.ptr
         except:
             ptr = fij.data_ptr()
-        update_pair_forces(self.data, <double*>ptr)  
+        update_pair_forces(self.data, <double*>ptr)
     def update_pair_forces(self, fij):
         if self.data.dev==0:
             self.update_pair_forces_cpu(fij)
@@ -189,11 +190,11 @@ cdef class MLIAPDataPy:
             return None
         return create_array(self.data.dev, self.data.f, [self.ntotal, 3],False)
 
-    
+
     @property
     def size_gradforce(self):
         return self.data.size_gradforce
- 
+
     @write_only_property
     def gradforce(self, value):
         if self.data.gradforce is NULL:
@@ -202,7 +203,7 @@ cdef class MLIAPDataPy:
         cdef double[:, :] value_view = value
         gradforce_view[:] = value_view
         print("This code has not been tested or optimized for the GPU, if you are getting this warning optimize gradforce")
- 
+
     @write_only_property
     def betas(self, value):
         if self.data.betas is NULL:
@@ -280,7 +281,7 @@ cdef class MLIAPDataPy:
     @property
     def ntotal(self):
         return self.data.ntotal
-    
+
     @property
     def elems(self):
         if self.data.elems is NULL:
@@ -290,7 +291,11 @@ cdef class MLIAPDataPy:
     @property
     def nlistatoms(self):
         return self.data.nlistatoms
-    
+
+    @property
+    def nlocal(self):
+        return self.data.nlocal
+
     @property
     def natomneigh(self):
         return self.data.natomneigh
@@ -306,7 +311,7 @@ cdef class MLIAPDataPy:
         if self.data.iatoms is NULL:
             return None
         return create_array(self.data.dev, self.data.iatoms, [self.natomneigh],True)
-    
+
     @property
     def ielems(self):
         if self.data.ielems is NULL:
@@ -322,7 +327,7 @@ cdef class MLIAPDataPy:
         if self.data.pair_i is NULL:
             return None
         return create_array(self.data.dev, self.data.pair_i, [self.npairs],True)
-    
+
     @property
     def pair_j(self):
         return self.jatoms
@@ -332,7 +337,7 @@ cdef class MLIAPDataPy:
         if self.data.jatoms is NULL:
             return None
         return create_array(self.data.dev, self.data.jatoms, [self.npairs],True)
-    
+
     @property
     def jelems(self):
         if self.data.jelems is NULL:
@@ -383,13 +388,13 @@ cdef class MLIAPUnifiedInterfaceKokkos:
         self.model = NULL
         self.descriptor = NULL
         self.unified_impl = unified_impl
-    
+
     def compute_gradients(self, data):
         self.unified_impl.compute_gradients(data)
-    
+
     def compute_descriptors(self, data):
         self.unified_impl.compute_descriptors(data)
-    
+
     def compute_forces(self, data):
         self.unified_impl.compute_forces(data)
 
@@ -443,7 +448,7 @@ cdef public object mliap_unified_connect_kokkos(char *fname, MLIAPDummyModel * m
 
     if unified.element_types is None:
         raise ValueError("no element type set")
-    
+
     cdef int nelements = <int>len(unified.element_types)
     cdef char **elements = <char**>malloc(nelements * sizeof(char*))
 

src/KOKKOS/mliap_unified_kokkos.cpp

@@ -271,8 +271,8 @@ void LAMMPS_NS::update_pair_energy(MLIAPDataKokkosDevice *data, double *eij)
 {
   auto d_eatoms = data->eatoms;
   auto d_pair_i= data->pair_i;
-  const auto nlistatoms = data->nlistatoms;
-  Kokkos::parallel_for(nlistatoms, KOKKOS_LAMBDA(int ii){
+  const auto nlocal = data->nlocal;
+  Kokkos::parallel_for(nlocal, KOKKOS_LAMBDA(int ii){
     d_eatoms[ii] = 0;
   });
 
@@ -281,7 +281,7 @@ void LAMMPS_NS::update_pair_energy(MLIAPDataKokkosDevice *data, double *eij)
     double e = 0.5 * eij[ii];
 
     // must not count any contribution where i is not a local atom
-    if (i < nlistatoms) {
+    if (i < nlocal) {
       Kokkos::atomic_add(&d_eatoms[i], e);
       local_sum += e;
     }
@@ -294,7 +294,7 @@ void LAMMPS_NS::update_pair_energy(MLIAPDataKokkosDevice *data, double *eij)
 
 void LAMMPS_NS::update_pair_forces(MLIAPDataKokkosDevice *data, double *fij)
 {
-  const auto nlistatoms = data->nlistatoms;
+  const auto nlocal = data->nlocal;
   auto *f = data->f;
   auto pair_i = data->pair_i;
   auto j_atoms = data->jatoms;
@@ -315,7 +315,7 @@ void LAMMPS_NS::update_pair_forces(MLIAPDataKokkosDevice *data, double *fij)
     int i = pair_i[ii];
     int j = j_atoms[ii];
     // must not count any contribution where i is not a local atom
-    if (i < nlistatoms) {
+    if (i < nlocal) {
       Kokkos::atomic_add(&f[i*3+0], fij[ii3+0]);
       Kokkos::atomic_add(&f[i*3+1], fij[ii3+1]);
       Kokkos::atomic_add(&f[i*3+2], fij[ii3+2]);
@@ -378,12 +378,12 @@ void LAMMPS_NS::update_pair_forces(MLIAPDataKokkosDevice *data, double *fij)
 void LAMMPS_NS::update_atom_energy(MLIAPDataKokkosDevice *data, double *ei)
 {
   auto d_eatoms = data->eatoms;
-  const auto nlistatoms = data->nlistatoms;
+  const auto nlocal = data->nlocal;
 
-  Kokkos::parallel_reduce(nlistatoms, KOKKOS_LAMBDA(int i, double &local_sum){
+  Kokkos::parallel_reduce(nlocal, KOKKOS_LAMBDA(int i, double &local_sum){
     double e = ei[i];
     // must not count any contribution where i is not a local atom
-    if (i < nlistatoms) {
+    if (i < nlocal) {
       d_eatoms[i] = e;
       local_sum += e;
     }

src/KOKKOS/modify_kokkos.cpp

@@ -362,6 +362,17 @@ void ModifyKokkos::pre_reverse(int eflag, int vflag)
 
 void ModifyKokkos::post_force(int vflag)
 {
+  for (int i = 0; i < n_post_force_group; i++) {
+    atomKK->sync(fix[list_post_force_group[i]]->execution_space,
+                 fix[list_post_force_group[i]]->datamask_read);
+    int prev_auto_sync = lmp->kokkos->auto_sync;
+    if (!fix[list_post_force_group[i]]->kokkosable) lmp->kokkos->auto_sync = 1;
+    fix[list_post_force_group[i]]->post_force(vflag);
+    lmp->kokkos->auto_sync = prev_auto_sync;
+    atomKK->modified(fix[list_post_force_group[i]]->execution_space,
+                     fix[list_post_force_group[i]]->datamask_modify);
+  }
+
   for (int i = 0; i < n_post_force; i++) {
     atomKK->sync(fix[list_post_force[i]]->execution_space,
                  fix[list_post_force[i]]->datamask_read);

src/KOKKOS/neigh_bond_kokkos.cpp

@@ -112,9 +112,8 @@ void NeighBondKokkos<DeviceType>::init_topology_kk() {
   int i,m;
   int bond_off = 0;
   int angle_off = 0;
-  for (i = 0; i < modify->nfix; i++)
-    if ((strcmp(modify->fix[i]->style,"shake") == 0)
-        || (strcmp(modify->fix[i]->style,"rattle") == 0))
+  for (const auto &ifix : modify->get_fix_list())
+    if (utils::strmatch(ifix->style,"^shake") || utils::strmatch(ifix->style,"^rattle"))
       bond_off = angle_off = 1;
   if (force->bond && force->bond_match("quartic")) bond_off = 1;