Merge branch 'master' into plugin-loader

cmake/Modules/Packages/GPU.cmake

@@ -35,7 +35,7 @@ if(GPU_API STREQUAL "CUDA")
   if(NOT BIN2C)
     message(FATAL_ERROR "Could not find bin2c, use -DBIN2C=/path/to/bin2c to help cmake finding it.")
   endif()
-  option(CUDPP_OPT "Enable CUDPP_OPT" ON)
+  option(CUDPP_OPT "Enable GPU binning via CUDAPP (should be off for modern GPUs)" OFF)
   option(CUDA_MPS_SUPPORT "Enable tweaks to support CUDA Multi-process service (MPS)" OFF)
   if(CUDA_MPS_SUPPORT)
     if(CUDPP_OPT)

cmake/Modules/Packages/KOKKOS.cmake

@@ -51,10 +51,12 @@ if(DOWNLOAD_KOKKOS)
     INTERFACE_INCLUDE_DIRECTORIES "${INSTALL_DIR}/include"
     INTERFACE_LINK_LIBRARIES ${CMAKE_DL_LIBS})
   target_link_libraries(lammps PRIVATE LAMMPS::KOKKOS)
+  target_link_libraries(lmp PRIVATE LAMMPS::KOKKOS)
   add_dependencies(LAMMPS::KOKKOS kokkos_build)
 elseif(EXTERNAL_KOKKOS)
   find_package(Kokkos 3.3.01 REQUIRED CONFIG)
   target_link_libraries(lammps PRIVATE Kokkos::kokkos)
+  target_link_libraries(lmp PRIVATE Kokkos::kokkos)
 else()
   set(LAMMPS_LIB_KOKKOS_SRC_DIR ${LAMMPS_LIB_SOURCE_DIR}/kokkos)
   set(LAMMPS_LIB_KOKKOS_BIN_DIR ${LAMMPS_LIB_BINARY_DIR}/kokkos)
@@ -66,6 +68,7 @@ else()
                           ${LAMMPS_LIB_KOKKOS_BIN_DIR})
   target_include_directories(lammps PRIVATE ${Kokkos_INCLUDE_DIRS})
   target_link_libraries(lammps PRIVATE kokkos)
+  target_link_libraries(lmp PRIVATE kokkos)
 endif()
 target_compile_definitions(lammps PRIVATE -DLMP_KOKKOS)
 

cmake/presets/oneapi.cmake

@@ -0,0 +1,18 @@
+# preset that will enable the LLVM based Intel compilers with support for MPI and OpenMP (on Linux boxes)
+
+set(CMAKE_CXX_COMPILER "icpx" CACHE STRING "" FORCE)
+set(CMAKE_C_COMPILER "icx" CACHE STRING "" FORCE)
+set(CMAKE_Fortran_COMPILER "ifx" CACHE STRING "" FORCE)
+set(MPI_CXX "icpx" CACHE STRING "" FORCE)
+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" CACHE STRING "" FORCE)
+set(OpenMP_C_LIB_NAMES "omp" CACHE STRING "" FORCE)
+set(OpenMP_CXX "icpx" CACHE STRING "" FORCE)
+set(OpenMP_CXX_FLAGS "-qopenmp" CACHE STRING "" FORCE)
+set(OpenMP_CXX_LIB_NAMES "omp" CACHE STRING "" FORCE)
+set(OpenMP_Fortran_FLAGS "-qopenmp" CACHE STRING "" FORCE)
+set(OpenMP_omp_LIBRARY "libiomp5.so" CACHE PATH "" FORCE)
+

doc/src/Build_basics.rst

@@ -234,6 +234,8 @@ LAMMPS.
          cmake ../cmake -DCMAKE_C_COMPILER=gcc -DCMAKE_CXX_COMPILER=g++ -DCMAKE_Fortran_COMPILER=gfortran
          # Building with Intel Compilers:
          cmake ../cmake -DCMAKE_C_COMPILER=icc -DCMAKE_CXX_COMPILER=icpc -DCMAKE_Fortran_COMPILER=ifort
+         # Building with Intel oneAPI Compilers:
+         cmake ../cmake -DCMAKE_C_COMPILER=icx -DCMAKE_CXX_COMPILER=icpx -DCMAKE_Fortran_COMPILER=ifx
          # Building with LLVM/Clang Compilers:
          cmake ../cmake -DCMAKE_C_COMPILER=clang -DCMAKE_CXX_COMPILER=clang++ -DCMAKE_Fortran_COMPILER=flang
          # Building with PGI/Nvidia Compilers:
@@ -243,8 +245,10 @@ LAMMPS.
       provided that can be loaded with
       `-C ../cmake/presets/clang.cmake`.  Similarly,
       `-C ../cmake/presets/intel.cmake` should switch the compiler
-      toolchain to the Intel compilers and `-C ../cmake/presets/pgi.cmake`
+      toolchain to the legacy Intel compilers, `-C ../cmake/presets/oneapi.cmake`
-      should switch the compiler to the PGI compilers.
+      will switch to the LLVM based oneAPI Intel compilers,
+      and `-C ../cmake/presets/pgi.cmake`
+      will switch the compiler to the PGI compilers.
 
       In addition you can set ``CMAKE_TUNE_FLAGS`` to specifically add
       compiler flags to tune for optimal performance on given hosts. By

doc/src/Build_extras.rst

@@ -128,9 +128,9 @@ CMake build
                                 # default is 'gfx906' for HIP_PLATFORM=hcc and 'sm_50' for HIP_PLATFORM=nvcc
    -D HIP_USE_DEVICE_SORT=value # enables GPU sorting
                                 # value = yes (default) or no
-   -D CUDPP_OPT=value           # optimization setting for GPU_API=cuda
+   -D CUDPP_OPT=value           # use GPU binning on with CUDA (should be off for modern GPUs)
-                                # enables CUDA Performance Primitives Optimizations, must be "no" for CUDA_MPS_SUPPORT=yes
+                                # enables CUDA Performance Primitives, must be "no" for CUDA_MPS_SUPPORT=yes
-                                # value = yes (default) or no
+                                # value = yes or no (default)
    -D CUDA_MPS_SUPPORT=value    # enables some tweaks required to run with active nvidia-cuda-mps daemon
                                 # value = yes or no (default)
    -D USE_STATIC_OPENCL_LOADER=value  # downloads/includes OpenCL ICD loader library, no local OpenCL headers/libs needed
@@ -230,11 +230,12 @@ GPU architectures as supported by the CUDA toolkit in use. This is done
 through using the "--gencode " flag, which can be used multiple times and
 thus support all GPU architectures supported by your CUDA compiler.
 
-To include CUDA performance primitives set the Makefile variable
+To enable GPU binning via CUDA performance primitives set the Makefile variable
-``CUDPP_OPT = -DUSE_CUDPP -Icudpp_mini``.
+``CUDPP_OPT = -DUSE_CUDPP -Icudpp_mini``.  This should **not** be used with
+most modern GPUs.
 
 To support the CUDA multiprocessor server you can set the define
-``-DCUDA_PROXY``.  Please note that in this case you should **not** use
+``-DCUDA_PROXY``.  Please note that in this case you must **not** use
 the CUDA performance primitives and thus set the variable ``CUDPP_OPT``
 to empty.
 

doc/src/Howto_cmake.rst

@@ -411,10 +411,10 @@ interface (``ccmake`` or ``cmake-gui``).
 .. note::
 
    Using a preset to select a compiler package (``clang.cmake``,
-   ``gcc.cmake``, or ``intel.cmake``) are an exception to the option
+   ``gcc.cmake``, ``intel.cmake``, ``oneapi.cmake``, or ``pgi.cmake``)
-   of updating the configuration incrementally, as they will trigger
+   are an exception to the mechanism of updating the configuration incrementally,
-   a reset of cached internal CMake settings and thus reset them to
+   as they will trigger a reset of cached internal CMake settings and thus
-   their default values.
+   reset settings to their default values.
 
 Compilation and build targets
 -----------------------------

doc/src/Tools.rst

@@ -93,6 +93,7 @@ Miscellaneous tools
    * :ref:`i-pi <ipi>`
    * :ref:`kate <kate>`
    * :ref:`LAMMPS shell <lammps_shell>`
+   * :ref:`LAMMPS magic patterns for file(1) <magic>`
    * :ref:`singularity <singularity_tool>`
    * :ref:`SWIG interface <swig>`
    * :ref:`vim <vim>`
@@ -641,6 +642,39 @@ This tool was written by Ara Kooser at Sandia (askoose at sandia.gov).
 
 ----------
 
+.. _magic:
+
+Magic patterns for the "file" command
+-------------------------------------
+
+.. versionadded:: 10Mar2021
+
+The file ``magic`` contains patterns that are used by the
+`file program <https://en.wikipedia.org/wiki/File_(command)>`_
+available on most Unix-like operating systems which enables it
+to detect various LAMMPS files and print some useful information
+about them.  To enable these patterns, append or copy the contents
+of the file to either the file ``.magic`` in your home directory
+or (as administrator) to ``/etc/magic`` (for a system-wide
+installation).  Afterwards the ``file`` command should be able to
+detect most LAMMPS restarts, dump, data and log files. Examples:
+
+.. code-block:: bash
+
+   $ file *.*
+   dihedral-quadratic.restart:   LAMMPS binary restart file (rev 2), Version 10 Mar 2021, Little Endian
+   mol-pair-wf_cut.restart:      LAMMPS binary restart file (rev 2), Version 24 Dec 2020, Little Endian
+   atom.bin:                     LAMMPS atom style binary dump (rev 2), Little Endian, First time step: 445570
+   custom.bin:                   LAMMPS custom style binary dump (rev 2), Little Endian, First time step: 100
+   bn1.lammpstrj:                LAMMPS text mode dump, First time step: 5000
+   data.fourmol:                 LAMMPS data file written by LAMMPS
+   pnc.data:                     LAMMPS data file written by msi2lmp
+   data.spce:                    LAMMPS data file written by TopoTools
+   B.data:                       LAMMPS data file written by OVITO
+   log.lammps:                   LAMMPS log file written by version 10 Feb 2021
+
+----------
+
 .. _matlab:
 
 matlab tool

doc/utils/sphinx-config/false_positives.txt

@@ -160,6 +160,7 @@ atwt
 augt
 AuO
 automagically
+Auvergne
 Avalos
 avalue
 aveforce
@@ -1740,6 +1741,7 @@ lwsock
 lx
 ly
 Lybrand
+lyon
 Lyulin
 lz
 Maaravi
@@ -3256,6 +3258,7 @@ Tz
 Tzou
 ub
 Uberuaga
+uca
 uChem
 uCond
 uef

examples/ELASTIC_T/displace.mod

@@ -100,22 +100,12 @@ run ${nrun}
 
 # Obtain new stress tensor
  
-variable tmp equal pe
+variable pxx1 equal f_avp[1]
-variable e1 equal ${tmp}
+variable pyy1 equal f_avp[2]
-variable tmp equal press
+variable pzz1 equal f_avp[3]
-variable p1 equal ${tmp}
+variable pxy1 equal f_avp[4]
-variable tmp equal pxx
+variable pxz1 equal f_avp[5]
-variable pxx1 equal ${tmp}
+variable pyz1 equal f_avp[6]
-variable tmp equal pyy
-variable pyy1 equal ${tmp}
-variable tmp equal pzz
-variable pzz1 equal ${tmp}
-variable tmp equal pxy
-variable pxy1 equal ${tmp}
-variable tmp equal pxz
-variable pxz1 equal ${tmp}
-variable tmp equal pyz
-variable pyz1 equal ${tmp}
 
 # Compute elastic constant from pressure tensor
 

lib/gpu/Makefile.cuda

@@ -29,7 +29,8 @@ LIB_DIR = ./
 AR = ar
 BSH = /bin/sh
 
-CUDPP_OPT = -DUSE_CUDPP -Icudpp_mini
+# GPU binning not recommended for modern GPUs
+CUDPP_OPT = # -DUSE_CUDPP -Icudpp_mini
 CUDA_MPS  =
 
 # device code compiler and settings

lib/gpu/Makefile.linux

@@ -70,7 +70,8 @@ LIB_DIR = ./
 AR = ar
 BSH = /bin/sh
 
-CUDPP_OPT = -DUSE_CUDPP -Icudpp_mini
+# GPU binning not recommended with modern GPUs
+CUDPP_OPT = #-DUSE_CUDPP -Icudpp_mini
 
 include Nvidia.makefile
 

lib/gpu/Makefile.linux_multi

@@ -77,7 +77,8 @@ LIB_DIR = ./
 AR = ar
 BSH = /bin/sh
 
-CUDPP_OPT = -DUSE_CUDPP -Icudpp_mini
+# GPU binning not recommended with modern GPUs
+CUDPP_OPT = #-DUSE_CUDPP -Icudpp_mini
 
 include Nvidia.makefile_multi
 

lib/gpu/Makefile.serial

@@ -66,7 +66,8 @@ LIB_DIR = ./
 AR = ar
 BSH = /bin/sh
 
-CUDPP_OPT = -DUSE_CUDPP -Icudpp_mini
+# GPU binning not recommended for most modern GPUs
+CUDPP_OPT = #-DUSE_CUDPP -Icudpp_mini
 
 include Nvidia.makefile
 

lib/gpu/Makefile.turing

@@ -16,6 +16,7 @@ LIB_DIR = .
 AR = ar
 BSH = /bin/sh
 
+# GPU binning not recommended with most modern GPUs
 CUDPP_OPT = #-DUSE_CUDPP -Icudpp_mini
 
 include Nvidia.makefile

lib/gpu/geryon/nvd_device.h

@@ -24,6 +24,10 @@
 #ifndef NVD_DEVICE
 #define NVD_DEVICE
 
+// workaround after GPU package Feb2021 update
+// todo: make new neighbor code work with CUDA
+#define LAL_USE_OLD_NEIGHBOR
+
 #include <string>
 #include <vector>
 #include <iostream>

python/examples/pylammps/README

@@ -1,28 +0,0 @@
-# Compile LAMMPS as shared library
-
-git clone https://github.com/lammps/lammps.git
-cd lammps/src
-python Make.py -m mpi -png -s ffmpeg exceptions -a file
-
-make -j 4 mode=shlib auto
-cd ../..
-
-# Install Python package
-
-virtualenv testing
-source testing/bin/activate
-
-(testing) cd lammps/python
-(testing) python install.py
-(testing) pip install jupyter matplotlib mpi4py
-
-(testing) cd ../../examples
-
-# Launch jupter and work inside browser
-
-(testing) jupyter notebook
-
-# Use Ctrl+c to stop jupyter
-
-# finally exit the virtualenv
-(testing) deactivate

python/examples/pylammps/README.md

@@ -0,0 +1,89 @@
+# PyLammps and Jupyter Notebooks
+
+This folder contains examples showcasing the usage of the PyLammps Python
+interface and Jupyter notebooks. To use this you will need LAMMPS compiled as
+a shared library and the LAMMPS Python package installed.
+
+An extensive guide on how to achieve this is documented in the [LAMMPS manual](https://lammps.sandia.gov/doc/Python_install.html). There is also a [PyLammps tutorial](https://lammps.sandia.gov/doc/Howto_pylammps.html).
+
+The following will show one way of creating a Python virtual environment
+which has both LAMMPS and its Python package installed:
+
+1. Clone the LAMMPS source code
+
+   ```shell
+   $ git clone -b stable https://github.com/lammps/lammps.git
+   $ cd lammps
+   ```
+
+2. Create a build folder
+
+   ```shell
+   $ mkdir build
+   $ cd build
+   ```
+
+3. Create a virtual environment for Python
+
+   ```shell
+   $ python3 -m venv myenv
+   ```
+
+4. Extend `LD_LIBRARY_PATH` (Unix/Linux) or `DYLD_LIBRARY_PATH` (MacOS)
+
+   On Unix/Linux:
+   ```shell
+   $ echo 'export LD_LIBRARY_PATH=$VIRTUAL_ENV/lib:$LD_LIBRARY_PATH' >> myenv/bin/activate
+   ```
+
+   On MacOS:
+   ```shell
+   echo 'export DYLD_LIBRARY_PATH=$VIRTUAL_ENV/lib:$DYLD_LIBRARY_PATH' >> myenv/bin/activate
+   ```
+
+5. Activate the virtual environment
+
+   ```shell
+   $ source myenv/bin/activate
+   (myenv)$
+   ```
+
+6. Configure LAMMPS compilation (CMake)
+
+   ```shell
+   (myenv)$ cmake -C ../cmake/presets/minimal.cmake \
+                  -D BUILD_SHARED_LIBS=on \
+                  -D LAMMPS_EXCEPTIONS=on -D PKG_PYTHON=on \
+                  -D CMAKE_INSTALL_PREFIX=$VIRTUAL_ENV \
+                  ../cmake
+   ```
+
+7. Compile LAMMPS
+
+   ```shell
+   (myenv)$ cmake --build .
+   ```
+
+8. Install LAMMPS and Python package into virtual environment
+
+   ```shell
+   (myenv)$ cmake --install .
+   ```
+
+9. Install other Python packages into virtual environment
+
+   ```shell
+   (myenv)$ pip install jupyter matplotlib mpi4py
+   ```
+
+10. Navigate to pylammps examples folder
+
+    ```shell
+    (myenv)$ cd ../python/examples/pylammmps
+    ```
+
+11. Launch Jupyter and work inside browser
+
+    ```shell
+    (myenv)$ jupyter notebook
+    ```

python/examples/pylammps/dihedrals/dihedral.ipynb

@@ -19,9 +19,7 @@
   {
    "cell_type": "code",
    "execution_count": null,
-   "metadata": {
+   "metadata": {},
-    "collapsed": true
-   },
    "outputs": [],
    "source": [
     "import matplotlib.pyplot as plt"
@@ -30,9 +28,7 @@
   {
    "cell_type": "code",
    "execution_count": null,
-   "metadata": {
+   "metadata": {},
-    "collapsed": true
-   },
    "outputs": [],
    "source": [
     "from lammps import IPyLammps"
@@ -194,9 +190,7 @@
   {
    "cell_type": "code",
    "execution_count": null,
-   "metadata": {
+   "metadata": {},
-    "collapsed": true
-   },
    "outputs": [],
    "source": [
     "K = 80.0\n",
@@ -232,9 +226,7 @@
   {
    "cell_type": "code",
    "execution_count": null,
-   "metadata": {
+   "metadata": {},
-    "collapsed": true
-   },
    "outputs": [],
    "source": []
   }
@@ -255,7 +247,7 @@
    "name": "python",
    "nbconvert_exporter": "python",
    "pygments_lexer": "ipython3",
-   "version": "3.6.2"
+   "version": "3.9.2"
   }
  },
  "nbformat": 4,

python/examples/pylammps/interface_usage.ipynb

@@ -4,71 +4,29 @@
    "cell_type": "markdown",
    "metadata": {},
    "source": [
-    "# Using LAMMPS with iPython and Jupyter"
+    "# Example 2: Using the PyLammps interface"
    ]
   },
   {
    "cell_type": "markdown",
    "metadata": {},
    "source": [
-    "LAMMPS can be run interactively using iPython easily. This tutorial shows how to set this up."
+    "## Prerequisites\n",
+    "\n",
+    "Before running this example, make sure your Python environment can find the LAMMPS shared library (`liblammps.so`) and the LAMMPS Python package is installed. If you followed the [README](README.md) in this folder, this should already be the case. You can also find more information about how to compile LAMMPS and install the LAMMPS Python package in the [LAMMPS manual](https://lammps.sandia.gov/doc/Python_install.html). There is also a dedicated [PyLammps HowTo](https://lammps.sandia.gov/doc/Howto_pylammps.html)."
    ]
   },
   {
    "cell_type": "markdown",
    "metadata": {},
    "source": [
-    "## Installation"
+    "## Setup system"
-   ]
-  },
-  {
-   "cell_type": "markdown",
-   "metadata": {},
-   "source": [
-    "1. Download the latest version of LAMMPS into a folder (we will calls this `$LAMMPS_DIR` from now on)\n",
-    "2. Compile LAMMPS as a shared library and enable exceptions and PNG support\n",
-    "   ```bash\n",
-    "   cd $LAMMPS_DIR/src\n",
-    "   make mpi mode=shlib LMP_INC=\"-DLAMMPS_PNG -DLAMMPS_EXCEPTIONS\" JPG_LIB=\"-lpng\"\n",
-    "   ```\n",
-    "\n",
-    "3. Create a python virtualenv\n",
-    "   ```bash\n",
-    "   virtualenv testing\n",
-    "   source testing/bin/activate\n",
-    "   ```\n",
-    "\n",
-    "4. Inside the virtualenv install the lammps package\n",
-    "   ```\n",
-    "   (testing) cd $LAMMPS_DIR/python\n",
-    "   (testing) python install.py\n",
-    "   (testing) cd   # move to your working directory\n",
-    "   ```\n",
-    "\n",
-    "5. Install jupyter and ipython in the virtualenv\n",
-    "   ```bash\n",
-    "   (testing) pip install ipython jupyter\n",
-    "   ```\n",
-    "\n",
-    "6. Run jupyter notebook\n",
-    "   ```bash\n",
-    "   (testing) jupyter notebook\n",
-    "   ```"
-   ]
-  },
-  {
-   "cell_type": "markdown",
-   "metadata": {},
-   "source": [
-    "## Example"
    ]
   },
   {
    "cell_type": "code",
    "execution_count": null,
-   "metadata": {
+   "metadata": {},
-    "collapsed": true
-   },
    "outputs": [],
    "source": [
     "from lammps import IPyLammps"
@@ -90,7 +48,6 @@
    "outputs": [],
    "source": [
     "# 3d Lennard-Jones melt\n",
-    "\n",
     "L.units(\"lj\")\n",
     "L.atom_style(\"atomic\")\n",
     "L.atom_modify(\"map array\")\n",
@@ -116,6 +73,13 @@
     "L.run(10)"
    ]
   },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "## Visualize the initial state"
+   ]
+  },
   {
    "cell_type": "code",
    "execution_count": null,
@@ -205,9 +169,7 @@
   {
    "cell_type": "code",
    "execution_count": null,
-   "metadata": {
+   "metadata": {},
-    "collapsed": true
-   },
    "outputs": [],
    "source": [
     "L.variable(\"a index 2\")"
@@ -225,9 +187,7 @@
   {
    "cell_type": "code",
    "execution_count": null,
-   "metadata": {
+   "metadata": {},
-    "collapsed": true
-   },
    "outputs": [],
    "source": [
     "L.variable(\"t equal temp\")"
@@ -283,9 +243,7 @@
   {
    "cell_type": "code",
    "execution_count": null,
-   "metadata": {
+   "metadata": {},
-    "collapsed": true
-   },
    "outputs": [],
    "source": [
     "L.variable(\"b index a b c\")"
@@ -321,9 +279,16 @@
   {
    "cell_type": "code",
    "execution_count": null,
-   "metadata": {
+   "metadata": {},
-    "collapsed": true
+   "outputs": [],
-   },
+   "source": [
+    "L.lmp.command('variable i loop 10')"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
    "outputs": [],
    "source": [
     "L.variable(\"i loop 10\")"
@@ -379,7 +344,7 @@
    "metadata": {},
    "outputs": [],
    "source": [
-    "[x for x in dir(L.atoms[0]) if not x.startswith('__')]"
+    "dir(L.atoms[0])"
    ]
   },
   {
@@ -479,6 +444,15 @@
     "L.runs[0].thermo"
    ]
   },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "dir(L.runs[0].thermo)"
+   ]
+  },
   {
    "cell_type": "markdown",
    "metadata": {},
@@ -486,14 +460,48 @@
     "## Saving session to as LAMMPS input file"
    ]
   },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "PyLammps can keep track of all LAMMPS commands that are executed. This allows you to prototype a script and then later on save it as a regular input script:"
+   ]
+  },
   {
    "cell_type": "code",
    "execution_count": null,
-   "metadata": {
+   "metadata": {},
-    "collapsed": true
-   },
    "outputs": [],
    "source": [
+    "L = IPyLammps()\n",
+    "\n",
+    "# enable command history\n",
+    "L.enable_cmd_history = True\n",
+    "\n",
+    "# 3d Lennard-Jones melt\n",
+    "L.units(\"lj\")\n",
+    "L.atom_style(\"atomic\")\n",
+    "L.atom_modify(\"map array\")\n",
+    "\n",
+    "L.lattice(\"fcc\", 0.8442)\n",
+    "L.region(\"box block\", 0, 4, 0, 4, 0, 4)\n",
+    "L.create_box(1, \"box\")\n",
+    "L.create_atoms(1, \"box\")\n",
+    "L.mass(1, 1.0)\n",
+    "\n",
+    "L.velocity(\"all create\", 1.44, 87287, \"loop geom\")\n",
+    "\n",
+    "L.pair_style(\"lj/cut\", 2.5)\n",
+    "L.pair_coeff(1, 1, 1.0, 1.0, 2.5)\n",
+    "\n",
+    "L.neighbor(0.3, \"bin\")\n",
+    "L.neigh_modify(\"delay 0 every 20 check no\")\n",
+    "\n",
+    "L.fix(\"1 all nve\")\n",
+    "\n",
+    "L.run(10)\n",
+    "\n",
+    "# write LAMMPS input script with all commands executed so far (including implicit ones)\n",
     "L.write_script(\"in.output\")"
    ]
   },
@@ -503,7 +511,7 @@
    "metadata": {},
    "outputs": [],
    "source": [
-    "dir(L.runs[0].thermo)"
+    "!cat in.output"
    ]
   },
   {
@@ -530,7 +538,7 @@
    "name": "python",
    "nbconvert_exporter": "python",
    "pygments_lexer": "ipython3",
-   "version": "3.6.2"
+   "version": "3.9.2"
   }
  },
  "nbformat": 4,

python/examples/pylammps/interface_usage_bonds.ipynb

@@ -4,72 +4,29 @@
    "cell_type": "markdown",
    "metadata": {},
    "source": [
-    "# Using LAMMPS with iPython and Jupyter"
+    "# Example 3: 2D circle of particles inside of box with LJ walls"
    ]
   },
   {
    "cell_type": "markdown",
    "metadata": {},
    "source": [
-    "LAMMPS can be run interactively using iPython easily. This tutorial shows how to set this up."
+    "## Prerequisites\n",
+    "\n",
+    "Before running this example, make sure your Python environment can find the LAMMPS shared library (`liblammps.so`) and the LAMMPS Python package is installed. If you followed the [README](README.md) in this folder, this should already be the case. You can also find more information about how to compile LAMMPS and install the LAMMPS Python package in the [LAMMPS manual](https://lammps.sandia.gov/doc/Python_install.html). There is also a dedicated [PyLammps HowTo](https://lammps.sandia.gov/doc/Howto_pylammps.html)."
    ]
   },
   {
    "cell_type": "markdown",
    "metadata": {},
    "source": [
-    "## Installation"
+    "## Setup system"
-   ]
-  },
-  {
-   "cell_type": "markdown",
-   "metadata": {},
-   "source": [
-    "1. Download the latest version of LAMMPS into a folder (we will calls this `$LAMMPS_DIR` from now on)\n",
-    "2. Compile LAMMPS as a shared library and enable exceptions and PNG support\n",
-    "   ```bash\n",
-    "   cd $LAMMPS_DIR/src\n",
-    "   make yes-molecule\n",
-    "   make mpi mode=shlib LMP_INC=\"-DLAMMPS_PNG -DLAMMPS_EXCEPTIONS\" JPG_LIB=\"-lpng\"\n",
-    "   ```\n",
-    "\n",
-    "3. Create a python virtualenv\n",
-    "   ```bash\n",
-    "   virtualenv testing\n",
-    "   source testing/bin/activate\n",
-    "   ```\n",
-    "\n",
-    "4. Inside the virtualenv install the lammps package\n",
-    "   ```\n",
-    "   (testing) cd $LAMMPS_DIR/python\n",
-    "   (testing) python install.py\n",
-    "   (testing) cd   # move to your working directory\n",
-    "   ```\n",
-    "\n",
-    "5. Install jupyter and ipython in the virtualenv\n",
-    "   ```bash\n",
-    "   (testing) pip install ipython jupyter\n",
-    "   ```\n",
-    "\n",
-    "6. Run jupyter notebook\n",
-    "   ```bash\n",
-    "   (testing) jupyter notebook\n",
-    "   ```"
-   ]
-  },
-  {
-   "cell_type": "markdown",
-   "metadata": {},
-   "source": [
-    "## Example"
    ]
   },
   {
    "cell_type": "code",
    "execution_count": null,
-   "metadata": {
+   "metadata": {},
-    "collapsed": true
-   },
    "outputs": [],
    "source": [
     "from lammps import IPyLammps"
@@ -135,6 +92,13 @@
     "L.fix(3, \"all wall/lj93 ylo 0.0 1 1 2.5 yhi\", y, \"1 1 2.5\")"
    ]
   },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "## Visualize initial state"
+   ]
+  },
   {
    "cell_type": "code",
    "execution_count": null,
@@ -144,6 +108,13 @@
     "L.image(zoom=1.8)"
    ]
   },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "## Run simulation and visualize new state"
+   ]
+  },
   {
    "cell_type": "code",
    "execution_count": null,
@@ -254,9 +225,7 @@
   {
    "cell_type": "code",
    "execution_count": null,
-   "metadata": {
+   "metadata": {},
-    "collapsed": true
-   },
    "outputs": [],
    "source": [
     "L.variable(\"a index 2\")"
@@ -274,9 +243,7 @@
   {
    "cell_type": "code",
    "execution_count": null,
-   "metadata": {
+   "metadata": {},
-    "collapsed": true
-   },
    "outputs": [],
    "source": [
     "L.variable(\"t equal temp\")"
@@ -332,9 +299,7 @@
   {
    "cell_type": "code",
    "execution_count": null,
-   "metadata": {
+   "metadata": {},
-    "collapsed": true
-   },
    "outputs": [],
    "source": [
     "L.variable(\"b index a b c\")"
@@ -370,9 +335,7 @@
   {
    "cell_type": "code",
    "execution_count": null,
-   "metadata": {
+   "metadata": {},
-    "collapsed": true
-   },
    "outputs": [],
    "source": [
     "L.variable(\"i loop 10\")"
@@ -428,7 +391,7 @@
    "metadata": {},
    "outputs": [],
    "source": [
-    "[x for x in dir(L.atoms[0]) if not x.startswith('__')]"
+    "dir(L.atoms[0])"
    ]
   },
   {
@@ -479,9 +442,7 @@
   {
    "cell_type": "code",
    "execution_count": null,
-   "metadata": {
+   "metadata": {},
-    "collapsed": true
-   },
    "outputs": [],
    "source": []
   }
@@ -502,7 +463,7 @@
    "name": "python",
    "nbconvert_exporter": "python",
    "pygments_lexer": "ipython3",
-   "version": "3.6.2"
+   "version": "3.9.2"
   }
  },
  "nbformat": 4,

python/examples/pylammps/simple.ipynb

@@ -4,74 +4,32 @@
    "cell_type": "markdown",
    "metadata": {},
    "source": [
-    "# Using LAMMPS with iPython and Jupyter"
+    "# Example 1: Using LAMMPS with PyLammps"
    ]
   },
   {
    "cell_type": "markdown",
    "metadata": {},
    "source": [
-    "LAMMPS can be run interactively using iPython easily. This tutorial shows how to set this up."
+    "The LAMMPS Python package provides multiple interfaces. The `PyLammps` interface is a high-level abstration of the low-level `lammps` interface. `IPyLammps` further extends this interface with functions that are useful for Jupyter notebooks to enable embedding generated graphics and videos."
    ]
   },
   {
    "cell_type": "markdown",
    "metadata": {},
    "source": [
-    "## Installation"
+    "## Prerequisites\n",
+    "\n",
+    "Before Running this example, make sure your Python environment can find the LAMMPS shared library (`liblammps.so`) and the LAMMPS Python package is installed. If you followed the [README](README.md) in this folder, this should already be the case. You can also find more information about how to compile LAMMPS and install the LAMMPS Python package in the [LAMMPS manual](https://lammps.sandia.gov/doc/Python_install.html). There is also a dedicated [PyLammps HowTo](https://lammps.sandia.gov/doc/Howto_pylammps.html)."
    ]
   },
   {
    "cell_type": "markdown",
    "metadata": {},
    "source": [
-    "1. Download the latest version of LAMMPS into a folder (we will calls this `$LAMMPS_DIR` from now on)\n",
+    "## Creating a new simulation\n",
-    "2. Compile LAMMPS as a shared library and enable exceptions and PNG support\n",
-    "   ```bash\n",
-    "   cd $LAMMPS_DIR/src\n",
-    "   make mpi mode=shlib LMP_INC=\"-DLAMMPS_PNG -DLAMMPS_EXCEPTIONS\" JPG_LIB=\"-lpng\"\n",
-    "   ```\n",
     "\n",
-    "3. Create a python virtualenv\n",
+    "Once the LAMMPS shared library and the LAMMPS Python package are installed, you can create a new LAMMMPS instance in your Python interpreter as follows:"
-    "   ```bash\n",
-    "   virtualenv testing\n",
-    "   source testing/bin/activate\n",
-    "   ```\n",
-    "\n",
-    "4. Inside the virtualenv install the lammps package\n",
-    "   ```\n",
-    "   (testing) cd $LAMMPS_DIR/python\n",
-    "   (testing) python install.py\n",
-    "   (testing) cd   # move to your working directory\n",
-    "   ```\n",
-    "\n",
-    "5. Install jupyter and ipython in the virtualenv\n",
-    "   ```bash\n",
-    "   (testing) pip install ipython jupyter\n",
-    "   ```\n",
-    "\n",
-    "6. Run jupyter notebook\n",
-    "   ```bash\n",
-    "   (testing) jupyter notebook\n",
-    "   ```"
-   ]
-  },
-  {
-   "cell_type": "markdown",
-   "metadata": {},
-   "source": [
-    "## Example"
-   ]
-  },
-  {
-   "cell_type": "code",
-   "execution_count": null,
-   "metadata": {
-    "collapsed": true
-   },
-   "outputs": [],
-   "source": [
-    "from lammps import IPyLammps"
    ]
   },
   {
@@ -80,22 +38,53 @@
    "metadata": {},
    "outputs": [],
    "source": [
+    "from lammps import IPyLammps\n",
     "L = IPyLammps()"
    ]
   },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "With `PyLammps`/`IPyLammps` you can write LAMMPS simulations similar to the input script language. Take the following LAMMPS input script:\n",
+    "\n",
+    "```bash\n",
+    "# 3d Lennard-Jones melt\n",
+    "\n",
+    "units        lj\n",
+    "atom_style   atomic\n",
+    "\n",
+    "lattice      fcc 0.8442\n",
+    "region       box block 0 4 0 4 0 4\n",
+    "create_box   1 box\n",
+    "create_atoms 1 box\n",
+    "mass         1 1.0\n",
+    "\n",
+    "velocity     all create 1.44 87287 loop geom\n",
+    "\n",
+    "pair_style   lj/cut 2.5\n",
+    "pair_coeff   1 1 1.0 1.0 2.5\n",
+    "\n",
+    "neighbor     0.3 bin\n",
+    "neigh_modify delay 0 every 20 check no\n",
+    "\n",
+    "fix          1 all nve\n",
+    "\n",
+    "thermo       50\n",
+    "```\n",
+    "The equivalent can be written with `PyLammps`/`IPyLammps`:"
+   ]
+  },
   {
    "cell_type": "code",
    "execution_count": null,
    "metadata": {},
    "outputs": [],
    "source": [
-    "import math\n",
-    "\n",
     "# 3d Lennard-Jones melt\n",
     "\n",
     "L.units(\"lj\")\n",
     "L.atom_style(\"atomic\")\n",
-    "L.atom_modify(\"map array\")\n",
     "\n",
     "L.lattice(\"fcc\", 0.8442)\n",
     "L.region(\"box\", \"block\", 0, 4, 0, 4, 0, 4)\n",
@@ -111,13 +100,18 @@
     "L.neighbor(0.3, \"bin\")\n",
     "L.neigh_modify(\"delay\", 0, \"every\", 20, \"check no\")\n",
     "\n",
-    "L.fix(\"1 all nve\")\n",
+    "L.fix(\"1\", \"all\", \"nve\")\n",
     "\n",
-    "L.variable(\"fx atom fx\")\n",
+    "L.thermo(50)"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "## Visualizing the initial state\n",
     "\n",
-    "L.info(\"all\")\n",
+    "`IPyLammps` allows you to visualize the current simulation state with the [image](https://lammps.sandia.gov/doc/Python_module.html#lammps.IPyLammps.image) command. Here we use it to create an image of the initial state of the system."
-    "\n",
-    "L.run(10)"
    ]
   },
   {
@@ -128,6 +122,170 @@
    "source": [
     "L.image(zoom=1.0)"
    ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "## Running simulations\n",
+    "\n",
+    "Use the `run` command to start the simulation. In Jupyter the return value of the last command will be displayed. The `run` command will return the output of the simulation."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "L.run(150)"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "You can suppress it by adding a semicolon `;`."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "L.run(100);"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "Visualizing the system will now show us how the atoms have moved."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "L.image(zoom=1.0)"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "## Post-processing thermo output\n",
+    "\n",
+    "Independent of whether or not you suppress or show the output of the `run` command, `PyLammps` will record the output. Each `run` command creates a new entry in the `L.runs` list. So far our PyLammps instance `L` executed two `run` commands:"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "len(L.runs)"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "Each entry contains information about the simulation run, including the thermo output for the printed out time steps.\n",
+    "\n",
+    "```bash\n",
+    "# thermo output of a LAMMPS simulation run\n",
+    "Step Temp E_pair E_mol TotEng Press\n",
+    "       0         1.44   -6.7733681            0   -4.6218056   -5.0244179\n",
+    "      50   0.70303849   -5.6796164            0    -4.629178   0.50453907\n",
+    "     100   0.72628044   -5.7150774            0   -4.6299123   0.29765862\n",
+    "     150   0.78441711    -5.805142            0   -4.6331125 -0.086709661\n",
+    "```\n",
+    "\n",
+    "`PyLammps` already parses this information and makes it available as dictionaries and arrays."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "L.runs[0]"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "L.runs[1]"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "For example, the first run was 150 time steps, with printing out a line every 50 steps. You can access the list of time steps using `{entry}.thermo.Step`."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "L.runs[0].thermo.Step"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "The corresponding values of each thermo quantity are also accessed this way:"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "L.runs[0].thermo.TotEng"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "Together you can use this information to run post-processing on these values or even plot it using `matplotlib`:"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "%matplotlib inline\n",
+    "import matplotlib.pyplot as plt\n",
+    "\n",
+    "plt.xlabel('time step')\n",
+    "plt.ylabel('Total Energy')\n",
+    "plt.plot(L.runs[0].thermo.Step, L.runs[0].thermo.TotEng)"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": []
   }
  ],
  "metadata": {
@@ -146,7 +304,7 @@
    "name": "python",
    "nbconvert_exporter": "python",
    "pygments_lexer": "ipython3",
-   "version": "3.6.2"
+   "version": "3.9.2"
   }
  },
  "nbformat": 4,

python/lammps/pylammps.py

@@ -132,6 +132,9 @@ class Atom(object):
     self._pylmp = pylammps_instance
     self.index = index
 
+  def __dir__(self):
+    return [k for k in super().__dir__() if not k.startswith('_')]
+
   @property
   def id(self):
     """
@@ -298,6 +301,9 @@ class variable_set:
     def __str__(self):
         return "{}({})".format(self._name, ','.join(["{}={}".format(k, self.__dict__[k]) for k in self.__dict__.keys() if not k.startswith('_')]))
 
+    def __dir__(self):
+        return [k for k in self.__dict__.keys() if not k.startswith('_')]
+
     def __repr__(self):
         return self.__str__()
 

src/GPU/pair_beck_gpu.cpp

@@ -158,7 +158,7 @@ void PairBeckGPU::init_style()
   double cell_size = sqrt(maxcut) + neighbor->skin;
 
   int maxspecial=0;
-  if (atom->molecular)
+  if (atom->molecular != Atom::ATOMIC)
     maxspecial=atom->maxspecial;
   int mnf = 5e-2 * neighbor->oneatom;
   int success = beck_gpu_init(atom->ntypes+1, cutsq, aa, alpha, beta,

src/GPU/pair_born_coul_long_cs_gpu.cpp

@@ -194,7 +194,7 @@ void PairBornCoulLongCSGPU::init_style()
   g_ewald = force->kspace->g_ewald;
 
   int maxspecial=0;
-  if (atom->molecular)
+  if (atom->molecular != Atom::ATOMIC)
     maxspecial=atom->maxspecial;
   int mnf = 5e-2 * neighbor->oneatom;
   int success = bornclcs_gpu_init(atom->ntypes+1, cutsq,  rhoinv,

src/GPU/pair_born_coul_long_gpu.cpp

@@ -193,7 +193,7 @@ void PairBornCoulLongGPU::init_style()
   if (ncoultablebits) init_tables(cut_coul,cut_respa);
 
   int maxspecial=0;
-  if (atom->molecular)
+  if (atom->molecular != Atom::ATOMIC)
     maxspecial=atom->maxspecial;
   int mnf = 5e-2 * neighbor->oneatom;
   int success = borncl_gpu_init(atom->ntypes+1, cutsq,  rhoinv,

src/GPU/pair_born_coul_wolf_cs_gpu.cpp

@@ -179,7 +179,7 @@ void PairBornCoulWolfCSGPU::init_style()
     cut_coul;
 
   int maxspecial=0;
-  if (atom->molecular)
+  if (atom->molecular != Atom::ATOMIC)
     maxspecial=atom->maxspecial;
   int mnf = 5e-2 * neighbor->oneatom;
   int success = borncwcs_gpu_init(atom->ntypes+1, cutsq, rhoinv,

src/GPU/pair_born_coul_wolf_gpu.cpp

@@ -177,7 +177,7 @@ void PairBornCoulWolfGPU::init_style()
     cut_coul;
 
   int maxspecial=0;
-  if (atom->molecular)
+  if (atom->molecular != Atom::ATOMIC)
     maxspecial=atom->maxspecial;
   int mnf = 5e-2 * neighbor->oneatom;
   int success = borncw_gpu_init(atom->ntypes+1, cutsq, rhoinv,

src/GPU/pair_born_gpu.cpp

@@ -161,7 +161,7 @@ void PairBornGPU::init_style()
   double cell_size = sqrt(maxcut) + neighbor->skin;
 
   int maxspecial=0;
-  if (atom->molecular)
+  if (atom->molecular != Atom::ATOMIC)
     maxspecial=atom->maxspecial;
   int mnf = 5e-2 * neighbor->oneatom;
   int success = born_gpu_init(atom->ntypes+1, cutsq, rhoinv,

src/GPU/pair_buck_coul_cut_gpu.cpp

@@ -165,7 +165,7 @@ void PairBuckCoulCutGPU::init_style()
   double cell_size = sqrt(maxcut) + neighbor->skin;
 
   int maxspecial=0;
-  if (atom->molecular)
+  if (atom->molecular != Atom::ATOMIC)
     maxspecial=atom->maxspecial;
   int mnf = 5e-2 * neighbor->oneatom;
   int success = buckc_gpu_init(atom->ntypes+1, cutsq, rhoinv, buck1, buck2,

src/GPU/pair_buck_coul_long_gpu.cpp

@@ -189,7 +189,7 @@ void PairBuckCoulLongGPU::init_style()
   if (ncoultablebits) init_tables(cut_coul,cut_respa);
 
   int maxspecial=0;
-  if (atom->molecular)
+  if (atom->molecular != Atom::ATOMIC)
     maxspecial=atom->maxspecial;
   int mnf = 5e-2 * neighbor->oneatom;
   int success = buckcl_gpu_init(atom->ntypes+1, cutsq,  rhoinv, buck1, buck2,

src/GPU/pair_buck_gpu.cpp

@@ -159,7 +159,7 @@ void PairBuckGPU::init_style()
   double cell_size = sqrt(maxcut) + neighbor->skin;
 
   int maxspecial=0;
-  if (atom->molecular)
+  if (atom->molecular != Atom::ATOMIC)
     maxspecial=atom->maxspecial;
   int mnf = 5e-2 * neighbor->oneatom;
   int success = buck_gpu_init(atom->ntypes+1, cutsq, rhoinv, buck1, buck2,

src/GPU/pair_colloid_gpu.cpp

@@ -169,7 +169,7 @@ void PairColloidGPU::init_style()
     }
   }
   int maxspecial=0;
-  if (atom->molecular)
+  if (atom->molecular != Atom::ATOMIC)
     maxspecial=atom->maxspecial;
   int mnf = 5e-2 * neighbor->oneatom;
   int success = colloid_gpu_init(atom->ntypes+1, cutsq, lj1, lj2, lj3, lj4,

src/GPU/pair_coul_cut_gpu.cpp

@@ -164,7 +164,7 @@ void PairCoulCutGPU::init_style()
   double cell_size = sqrt(maxcut) + neighbor->skin;
 
   int maxspecial=0;
-  if (atom->molecular)
+  if (atom->molecular != Atom::ATOMIC)
     maxspecial=atom->maxspecial;
   int mnf = 5e-2 * neighbor->oneatom;
   int success = coul_gpu_init(atom->ntypes+1, scale, cutsq,

src/GPU/pair_coul_debye_gpu.cpp

@@ -165,7 +165,7 @@ void PairCoulDebyeGPU::init_style()
   double cell_size = sqrt(maxcut) + neighbor->skin;
 
   int maxspecial=0;
-  if (atom->molecular)
+  if (atom->molecular != Atom::ATOMIC)
     maxspecial=atom->maxspecial;
   int mnf = 5e-2 * neighbor->oneatom;
   int success = cdebye_gpu_init(atom->ntypes+1, scale, cutsq,

src/GPU/pair_coul_dsf_gpu.cpp

@@ -182,7 +182,7 @@ void PairCoulDSFGPU::init_style()
   e_shift = erfcc/cut_coul - f_shift*cut_coul;
 
   int maxspecial=0;
-  if (atom->molecular)
+  if (atom->molecular != Atom::ATOMIC)
     maxspecial=atom->maxspecial;
   int mnf = 5e-2 * neighbor->oneatom;
   int success = cdsf_gpu_init(atom->ntypes+1, atom->nlocal,

src/GPU/pair_coul_long_cs_gpu.cpp

@@ -184,7 +184,7 @@ void PairCoulLongCSGPU::init_style()
   if (ncoultablebits) init_tables(cut_coul,cut_respa);
 
   int maxspecial=0;
-  if (atom->molecular)
+  if (atom->molecular != Atom::ATOMIC)
     maxspecial=atom->maxspecial;
   int mnf = 5e-2 * neighbor->oneatom;
   int success = clcs_gpu_init(atom->ntypes+1, scale,

src/GPU/pair_coul_long_gpu.cpp

@@ -179,7 +179,7 @@ void PairCoulLongGPU::init_style()
   if (ncoultablebits) init_tables(cut_coul,cut_respa);
 
   int maxspecial=0;
-  if (atom->molecular)
+  if (atom->molecular != Atom::ATOMIC)
     maxspecial=atom->maxspecial;
   int mnf = 5e-2 * neighbor->oneatom;
   int success = cl_gpu_init(atom->ntypes+1, scale,

src/GPU/pair_dpd_gpu.cpp

@@ -306,7 +306,7 @@ void PairDPDGPU::init_style()
   double cell_size = sqrt(maxcut) + neighbor->skin;
 
   int maxspecial=0;
-  if (atom->molecular)
+  if (atom->molecular != Atom::ATOMIC)
     maxspecial=atom->maxspecial;
   int mnf = 5e-2 * neighbor->oneatom;
   int success = dpd_gpu_init(atom->ntypes+1, cutsq, a0, gamma, sigma,

src/GPU/pair_dpd_tstat_gpu.cpp

@@ -325,7 +325,7 @@ void PairDPDTstatGPU::init_style()
   double cell_size = sqrt(maxcut) + neighbor->skin;
 
   int maxspecial=0;
-  if (atom->molecular)
+  if (atom->molecular != Atom::ATOMIC)
     maxspecial=atom->maxspecial;
   int mnf = 5e-2 * neighbor->oneatom;
   int success = dpd_tstat_gpu_init(atom->ntypes+1, cutsq, a0, gamma, sigma,

src/GPU/pair_eam_alloy_gpu.cpp

@@ -181,7 +181,7 @@ void PairEAMAlloyGPU::init_style()
   double cell_size = sqrt(maxcut) + neighbor->skin;
 
   int maxspecial=0;
-  if (atom->molecular)
+  if (atom->molecular != Atom::ATOMIC)
     maxspecial=atom->maxspecial;
   int fp_size;
   int mnf = 5e-2 * neighbor->oneatom;

src/GPU/pair_eam_fs_gpu.cpp

@@ -180,7 +180,7 @@ void PairEAMFSGPU::init_style()
   double cell_size = sqrt(maxcut) + neighbor->skin;
 
   int maxspecial=0;
-  if (atom->molecular)
+  if (atom->molecular != Atom::ATOMIC)
     maxspecial=atom->maxspecial;
   int fp_size;
   int mnf = 5e-2 * neighbor->oneatom;

src/GPU/pair_eam_gpu.cpp

@@ -182,7 +182,7 @@ void PairEAMGPU::init_style()
   double cell_size = sqrt(maxcut) + neighbor->skin;
 
   int maxspecial=0;
-  if (atom->molecular)
+  if (atom->molecular != Atom::ATOMIC)
     maxspecial=atom->maxspecial;
   int fp_size;
   int mnf = 5e-2 * neighbor->oneatom;

src/GPU/pair_gauss_gpu.cpp

@@ -157,7 +157,7 @@ void PairGaussGPU::init_style()
   double cell_size = sqrt(maxcut) + neighbor->skin;
 
   int maxspecial=0;
-  if (atom->molecular)
+  if (atom->molecular != Atom::ATOMIC)
     maxspecial=atom->maxspecial;
   int mnf = 5e-2 * neighbor->oneatom;
   int success = gauss_gpu_init(atom->ntypes+1, cutsq, a, b,

src/GPU/pair_gayberne_gpu.cpp

@@ -205,7 +205,7 @@ void PairGayBerneGPU::init_style()
   double cell_size = sqrt(maxcut) + neighbor->skin;
 
   int maxspecial=0;
-  if (atom->molecular)
+  if (atom->molecular != Atom::ATOMIC)
     maxspecial=atom->maxspecial;
   int mnf = 5e-2 * neighbor->oneatom;
   int success = gb_gpu_init(atom->ntypes+1, gamma, upsilon, mu,

src/GPU/pair_lj96_cut_gpu.cpp

@@ -158,7 +158,7 @@ void PairLJ96CutGPU::init_style()
   double cell_size = sqrt(maxcut) + neighbor->skin;
 
   int maxspecial=0;
-  if (atom->molecular)
+  if (atom->molecular != Atom::ATOMIC)
     maxspecial=atom->maxspecial;
   int mnf = 5e-2 * neighbor->oneatom;
   int success = lj96_gpu_init(atom->ntypes+1, cutsq, lj1, lj2, lj3, lj4,

src/GPU/pair_lj_charmm_coul_charmm_gpu.cpp

@@ -165,7 +165,7 @@ void PairLJCharmmCoulCharmmGPU::init_style()
   double cell_size = sqrt(cut_bothsq) + neighbor->skin;
 
   int maxspecial=0;
-  if (atom->molecular)
+  if (atom->molecular != Atom::ATOMIC)
     maxspecial=atom->maxspecial;
 
   bool arithmetic = true;

src/GPU/pair_lj_charmm_coul_long_gpu.cpp

@@ -191,7 +191,7 @@ void PairLJCharmmCoulLongGPU::init_style()
   if (ncoultablebits) init_tables(cut_coul,cut_respa);
 
   int maxspecial=0;
-  if (atom->molecular)
+  if (atom->molecular != Atom::ATOMIC)
     maxspecial=atom->maxspecial;
 
   bool arithmetic = true;

src/GPU/pair_lj_class2_coul_long_gpu.cpp

@@ -186,7 +186,7 @@ void PairLJClass2CoulLongGPU::init_style()
   if (ncoultablebits) init_tables(cut_coul,cut_respa);
 
   int maxspecial=0;
-  if (atom->molecular)
+  if (atom->molecular != Atom::ATOMIC)
     maxspecial=atom->maxspecial;
   int mnf = 5e-2 * neighbor->oneatom;
   int success = c2cl_gpu_init(atom->ntypes+1, cutsq, lj1, lj2, lj3, lj4,

src/GPU/pair_lj_class2_gpu.cpp

@@ -155,7 +155,7 @@ void PairLJClass2GPU::init_style()
   double cell_size = sqrt(maxcut) + neighbor->skin;
 
   int maxspecial=0;
-  if (atom->molecular)
+  if (atom->molecular != Atom::ATOMIC)
     maxspecial=atom->maxspecial;
   int mnf = 5e-2 * neighbor->oneatom;
   int success = lj96_gpu_init(atom->ntypes+1, cutsq, lj1, lj2, lj3, lj4,

src/GPU/pair_lj_cubic_gpu.cpp

@@ -163,7 +163,7 @@ void PairLJCubicGPU::init_style()
   double cell_size = sqrt(maxcut) + neighbor->skin;
 
   int maxspecial=0;
-  if (atom->molecular)
+  if (atom->molecular != Atom::ATOMIC)
     maxspecial=atom->maxspecial;
   int mnf = 5e-2 * neighbor->oneatom;
   int success = ljcb_gpu_init(atom->ntypes+1, cutsq, cut_inner_sq,

src/GPU/pair_lj_cut_coul_cut_gpu.cpp

@@ -166,7 +166,7 @@ void PairLJCutCoulCutGPU::init_style()
   double cell_size = sqrt(maxcut) + neighbor->skin;
 
   int maxspecial=0;
-  if (atom->molecular)
+  if (atom->molecular != Atom::ATOMIC)
     maxspecial=atom->maxspecial;
   int mnf = 5e-2 * neighbor->oneatom;
   int success = ljc_gpu_init(atom->ntypes+1, cutsq, lj1, lj2, lj3, lj4,

src/GPU/pair_lj_cut_coul_debye_gpu.cpp

@@ -168,7 +168,7 @@ void PairLJCutCoulDebyeGPU::init_style()
   double cell_size = sqrt(maxcut) + neighbor->skin;
 
   int maxspecial=0;
-  if (atom->molecular)
+  if (atom->molecular != Atom::ATOMIC)
     maxspecial=atom->maxspecial;
   int mnf = 5e-2 * neighbor->oneatom;
   int success = ljcd_gpu_init(atom->ntypes+1, cutsq, lj1, lj2, lj3, lj4,

src/GPU/pair_lj_cut_coul_dsf_gpu.cpp

@@ -183,7 +183,7 @@ void PairLJCutCoulDSFGPU::init_style()
   e_shift = erfcc/cut_coul - f_shift*cut_coul;
 
   int maxspecial=0;
-  if (atom->molecular)
+  if (atom->molecular != Atom::ATOMIC)
     maxspecial=atom->maxspecial;
   int mnf = 5e-2 * neighbor->oneatom;
   int success = ljd_gpu_init(atom->ntypes+1, cutsq, lj1, lj2, lj3, lj4,

src/GPU/pair_lj_cut_coul_long_gpu.cpp

@@ -191,7 +191,7 @@ void PairLJCutCoulLongGPU::init_style()
   if (ncoultablebits) init_tables(cut_coul,cut_respa);
 
   int maxspecial=0;
-  if (atom->molecular)
+  if (atom->molecular != Atom::ATOMIC)
     maxspecial=atom->maxspecial;
   int mnf = 5e-2 * neighbor->oneatom;
   int success = ljcl_gpu_init(atom->ntypes+1, cutsq, lj1, lj2, lj3, lj4,

src/GPU/pair_lj_cut_coul_msm_gpu.cpp

@@ -177,7 +177,7 @@ void PairLJCutCoulMSMGPU::init_style()
   if (ncoultablebits) init_tables(cut_coul,cut_respa);
 
   int maxspecial=0;
-  if (atom->molecular)
+  if (atom->molecular != Atom::ATOMIC)
     maxspecial=atom->maxspecial;
   int mnf = 5e-2 * neighbor->oneatom;
   int success = ljcm_gpu_init(atom->ntypes+1, cutsq, lj1, lj2, lj3, lj4,

src/GPU/pair_lj_cut_dipole_cut_gpu.cpp

@@ -171,7 +171,7 @@ void PairLJCutDipoleCutGPU::init_style()
   double cell_size = sqrt(maxcut) + neighbor->skin;
 
   int maxspecial=0;
-  if (atom->molecular)
+  if (atom->molecular != Atom::ATOMIC)
     maxspecial=atom->maxspecial;
   int mnf = 5e-2 * neighbor->oneatom;
   int success = dpl_gpu_init(atom->ntypes+1, cutsq, lj1, lj2, lj3, lj4,

src/GPU/pair_lj_cut_dipole_long_gpu.cpp

@@ -195,7 +195,7 @@ void PairLJCutDipoleLongGPU::init_style()
   if (ncoultablebits) init_tables(cut_coul,nullptr);
 
   int maxspecial=0;
-  if (atom->molecular)
+  if (atom->molecular != Atom::ATOMIC)
     maxspecial=atom->maxspecial;
   int mnf = 5e-2 * neighbor->oneatom;
   int success = dplj_gpu_init(atom->ntypes+1, cutsq, lj1, lj2, lj3, lj4,

src/GPU/pair_lj_cut_gpu.cpp

@@ -162,7 +162,7 @@ void PairLJCutGPU::init_style()
   double cell_size = sqrt(maxcut) + neighbor->skin;
 
   int maxspecial=0;
-  if (atom->molecular)
+  if (atom->molecular != Atom::ATOMIC)
     maxspecial=atom->maxspecial;
   int mnf = 5e-2 * neighbor->oneatom;
   int success = ljl_gpu_init(atom->ntypes+1, cutsq, lj1, lj2, lj3, lj4,

src/GPU/pair_lj_cut_tip4p_long_gpu.cpp

@@ -210,7 +210,7 @@ void PairLJCutTIP4PLongGPU::init_style()
   if (ncoultablebits) init_tables(cut_coul,cut_respa);
 
   int maxspecial=0;
-  if (atom->molecular)
+  if (atom->molecular != Atom::ATOMIC)
     maxspecial=atom->maxspecial;
 
   // set alpha parameter

src/GPU/pair_lj_expand_coul_long_gpu.cpp

@@ -191,7 +191,7 @@ void PairLJExpandCoulLongGPU::init_style()
   if (ncoultablebits) init_tables(cut_coul,cut_respa);
 
   int maxspecial=0;
-  if (atom->molecular)
+  if (atom->molecular != Atom::ATOMIC)
     maxspecial=atom->maxspecial;
   int mnf = 5e-2 * neighbor->oneatom;
   int success = ljecl_gpu_init(atom->ntypes+1, cutsq, lj1, lj2, lj3, lj4,

src/GPU/pair_lj_expand_gpu.cpp

@@ -159,7 +159,7 @@ void PairLJExpandGPU::init_style()
   double cell_size = sqrt(maxcut) + neighbor->skin;
 
   int maxspecial=0;
-  if (atom->molecular)
+  if (atom->molecular != Atom::ATOMIC)
     maxspecial=atom->maxspecial;
   int mnf = 5e-2 * neighbor->oneatom;
   int success = lje_gpu_init(atom->ntypes+1, cutsq, lj1, lj2, lj3, lj4,

src/GPU/pair_lj_gromacs_gpu.cpp

@@ -162,7 +162,7 @@ void PairLJGromacsGPU::init_style()
   double cell_size = sqrt(maxcut) + neighbor->skin;
 
   int maxspecial=0;
-  if (atom->molecular)
+  if (atom->molecular != Atom::ATOMIC)
     maxspecial=atom->maxspecial;
 
   int mnf = 5e-2 * neighbor->oneatom;

src/GPU/pair_lj_sdk_coul_long_gpu.cpp

@@ -195,7 +195,7 @@ void PairLJSDKCoulLongGPU::init_style()
   if (ncoultablebits) init_tables(cut_coul,nullptr);
 
   int maxspecial=0;
-  if (atom->molecular)
+  if (atom->molecular != Atom::ATOMIC)
     maxspecial=atom->maxspecial;
   int mnf = 5e-2 * neighbor->oneatom;
   int success = sdkl_gpu_init(atom->ntypes+1, cutsq, lj_type, lj1, lj2, lj3,

src/GPU/pair_lj_sdk_gpu.cpp

@@ -164,7 +164,7 @@ void PairLJSDKGPU::init_style()
   double cell_size = sqrt(maxcut) + neighbor->skin;
 
   int maxspecial=0;
-  if (atom->molecular)
+  if (atom->molecular != Atom::ATOMIC)
     maxspecial=atom->maxspecial;
   int mnf = 5e-2 * neighbor->oneatom;
   int success = sdk_gpu_init(atom->ntypes+1,cutsq,lj_type,lj1,lj2,lj3,lj4,

src/GPU/pair_lj_sf_dipole_sf_gpu.cpp

@@ -170,7 +170,7 @@ void PairLJSFDipoleSFGPU::init_style()
   double cell_size = sqrt(maxcut) + neighbor->skin;
 
   int maxspecial=0;
-  if (atom->molecular)
+  if (atom->molecular != Atom::ATOMIC)
     maxspecial=atom->maxspecial;
   int mnf = 5e-2 * neighbor->oneatom;
   int success = dplsf_gpu_init(atom->ntypes+1, cutsq, lj1, lj2, lj3, lj4,

src/GPU/pair_mie_cut_gpu.cpp

@@ -159,7 +159,7 @@ void PairMIECutGPU::init_style()
   double cell_size = sqrt(maxcut) + neighbor->skin;
 
   int maxspecial=0;
-  if (atom->molecular)
+  if (atom->molecular != Atom::ATOMIC)
     maxspecial=atom->maxspecial;
   int mnf = 5e-2 * neighbor->oneatom;
   int success = mie_gpu_init(atom->ntypes+1, cutsq, mie1, mie2, mie3, mie4,

src/GPU/pair_morse_gpu.cpp

@@ -155,7 +155,7 @@ void PairMorseGPU::init_style()
   double cell_size = sqrt(maxcut) + neighbor->skin;
 
   int maxspecial=0;
-  if (atom->molecular)
+  if (atom->molecular != Atom::ATOMIC)
     maxspecial=atom->maxspecial;
   int mnf = 5e-2 * neighbor->oneatom;
   int success = mor_gpu_init(atom->ntypes+1, cutsq, morse1, r0, alpha, d0,

src/GPU/pair_resquared_gpu.cpp

@@ -203,7 +203,7 @@ void PairRESquaredGPU::init_style()
   double cell_size = sqrt(maxcut) + neighbor->skin;
 
   int maxspecial=0;
-  if (atom->molecular)
+  if (atom->molecular != Atom::ATOMIC)
     maxspecial=atom->maxspecial;
   int mnf = 5e-2 * neighbor->oneatom;
   int success = re_gpu_init(atom->ntypes+1, shape1, well, cutsq, sigma,

src/GPU/pair_soft_gpu.cpp

@@ -160,7 +160,7 @@ void PairSoftGPU::init_style()
   double cell_size = sqrt(maxcut) + neighbor->skin;
 
   int maxspecial=0;
-  if (atom->molecular)
+  if (atom->molecular != Atom::ATOMIC)
     maxspecial=atom->maxspecial;
   int mnf = 5e-2 * neighbor->oneatom;
   int success = soft_gpu_init(atom->ntypes+1, cutsq, prefactor, cut,

src/GPU/pair_table_gpu.cpp

@@ -229,7 +229,7 @@ void PairTableGPU::init_style()
   }
 
   int maxspecial=0;
-  if (atom->molecular)
+  if (atom->molecular != Atom::ATOMIC)
     maxspecial=atom->maxspecial;
   int mnf = 5e-2 * neighbor->oneatom;
   int success = table_gpu_init(atom->ntypes+1, cutsq, table_coeffs, table_data,

src/GPU/pair_ufm_gpu.cpp

@@ -163,7 +163,7 @@ void PairUFMGPU::init_style()
   double cell_size = sqrt(maxcut) + neighbor->skin;
 
   int maxspecial=0;
-  if (atom->molecular)
+  if (atom->molecular != Atom::ATOMIC)
     maxspecial=atom->maxspecial;
   int mnf = 5e-2 * neighbor->oneatom;
   int success = ufml_gpu_init(atom->ntypes+1, cutsq, uf1, uf2, uf3,

src/GPU/pair_yukawa_colloid_gpu.cpp

@@ -168,7 +168,7 @@ void PairYukawaColloidGPU::init_style()
   double cell_size = sqrt(maxcut) + neighbor->skin;
 
   int maxspecial=0;
-  if (atom->molecular)
+  if (atom->molecular != Atom::ATOMIC)
     maxspecial=atom->maxspecial;
   int mnf = 5e-2 * neighbor->oneatom;
   int success = ykcolloid_gpu_init(atom->ntypes+1, cutsq, a,

src/GPU/pair_yukawa_gpu.cpp

@@ -157,7 +157,7 @@ void PairYukawaGPU::init_style()
   double cell_size = sqrt(maxcut) + neighbor->skin;
 
   int maxspecial=0;
-  if (atom->molecular)
+  if (atom->molecular != Atom::ATOMIC)
     maxspecial=atom->maxspecial;
   int mnf = 5e-2 * neighbor->oneatom;
   int success = yukawa_gpu_init(atom->ntypes+1, cutsq, kappa, a,

src/GPU/pair_zbl_gpu.cpp

@@ -163,7 +163,7 @@ void PairZBLGPU::init_style()
   cut_globalsq = cut_global * cut_global;
 
   int maxspecial=0;
-  if (atom->molecular)
+  if (atom->molecular != Atom::ATOMIC)
     maxspecial=atom->maxspecial;
   int mnf = 5e-2 * neighbor->oneatom;
   int success = zbl_gpu_init(atom->ntypes+1, cutsq, sw1, sw2, sw3, sw4,

src/KOKKOS/min_kokkos.cpp

@@ -462,7 +462,7 @@ double MinKokkos::energy_force(int resetflag)
     timer->stamp(Timer::PAIR);
   }
 
-  if (atom->molecular) {
+  if (atom->molecular != Atom::ATOMIC) {
     if (force->bond) {
       atomKK->sync(force->bond->execution_space,force->bond->datamask_read);
       force->bond->compute(eflag,vflag);

src/KOKKOS/npair_kokkos.cpp

@@ -189,7 +189,7 @@ void NPairKokkos<DeviceType,HALF_NEIGH,GHOST,TRI,SIZE>::build(NeighList *list_)
   k_bins.sync<DeviceType>();
   k_atom2bin.sync<DeviceType>();
 
-  if (atom->molecular) {
+  if (atom->molecular != Atom::ATOMIC) {
     if (exclude)
       atomKK->sync(Device,X_MASK|RADIUS_MASK|TYPE_MASK|MASK_MASK|MOLECULE_MASK|TAG_MASK|SPECIAL_MASK);
     else
@@ -545,14 +545,14 @@ __device__ __forceinline__ int __syncthreads_count(int predicate) {
 
 #ifdef LMP_KOKKOS_GPU
 template<class DeviceType> template<int HalfNeigh,int Newton,int Tri>
-__device__ inline
+#if defined(KOKKOS_ENABLE_CUDA) || defined(KOKKOS_ENABLE_HIP)
-void NeighborKokkosExecute<DeviceType>::build_ItemCuda(typename Kokkos::TeamPolicy<DeviceType>::member_type dev) const
+__device__
-{
-#ifdef KOKKOS_ENABLE_HIP
-  HIP_DYNAMIC_SHARED(X_FLOAT, sharedmem);
-#else
-  extern __shared__ X_FLOAT sharedmem[];
 #endif
+inline
+void NeighborKokkosExecute<DeviceType>::build_ItemGPU(typename Kokkos::TeamPolicy<DeviceType>::member_type dev,
+                                                      size_t sharedsize) const
+{
+  auto* sharedmem = static_cast<X_FLOAT *>(dev.team_shmem().get_shmem(sharedsize));
   /* loop over atoms in i's bin,
   */
   const int atoms_per_bin = c_bins.extent(1);
@@ -657,7 +657,7 @@ void NeighborKokkosExecute<DeviceType>::build_ItemCuda(typename Kokkos::TeamPoli
 
     }
   }
-  __syncthreads();
+  dev.team_barrier();
 
   const typename ArrayTypes<DeviceType>::t_int_1d_const_um stencil
     = d_stencil;
@@ -678,7 +678,7 @@ void NeighborKokkosExecute<DeviceType>::build_ItemCuda(typename Kokkos::TeamPoli
 
     other_id[MY_II] = j;
 
-    __syncthreads();
+    dev.team_barrier();
 
     if (i >= 0 && i < nlocal) {
       #pragma unroll 8
@@ -735,7 +735,7 @@ void NeighborKokkosExecute<DeviceType>::build_ItemCuda(typename Kokkos::TeamPoli
 
       }
     }
-    __syncthreads();
+    dev.team_barrier();
   }
 
   if (i >= 0 && i < nlocal) {
@@ -988,14 +988,14 @@ void NeighborKokkosExecute<DeviceType>::
 
 #ifdef LMP_KOKKOS_GPU
 template<class DeviceType> template<int HalfNeigh,int Newton,int Tri>
-__device__ inline
+#if defined(KOKKOS_ENABLE_CUDA) || defined(KOKKOS_ENABLE_HIP)
-void NeighborKokkosExecute<DeviceType>::build_ItemSizeCuda(typename Kokkos::TeamPolicy<DeviceType>::member_type dev) const
+__device__
-{
-#ifdef KOKKOS_ENABLE_HIP
-  HIP_DYNAMIC_SHARED(X_FLOAT, sharedmem);
-#else
-  extern __shared__ X_FLOAT sharedmem[];
 #endif
+inline
+void NeighborKokkosExecute<DeviceType>::build_ItemSizeGPU(typename Kokkos::TeamPolicy<DeviceType>::member_type dev,
+                                                          size_t sharedsize) const
+{
+  auto* sharedmem = static_cast<X_FLOAT *>(dev.team_shmem().get_shmem(sharedsize));
   /* loop over atoms in i's bin,
    */
   const int atoms_per_bin = c_bins.extent(1);
@@ -1086,7 +1086,7 @@ void NeighborKokkosExecute<DeviceType>::build_ItemSizeCuda(typename Kokkos::Team
         }
       }
     }
-    __syncthreads();
+    dev.team_barrier();
 
     const typename ArrayTypes<DeviceType>::t_int_1d_const_um stencil
       = d_stencil;
@@ -1108,7 +1108,7 @@ void NeighborKokkosExecute<DeviceType>::build_ItemSizeCuda(typename Kokkos::Team
 
       other_id[MY_II] = j;
 
-      __syncthreads();
+      dev.team_barrier();
 
       if (i >= 0 && i < nlocal) {
         #pragma unroll 8
@@ -1147,7 +1147,7 @@ void NeighborKokkosExecute<DeviceType>::build_ItemSizeCuda(typename Kokkos::Team
           }
         }
       }
-      __syncthreads();
+      dev.team_barrier();
     }
 
     if (i >= 0 && i < nlocal) {

src/KOKKOS/npair_kokkos.h

@@ -312,12 +312,20 @@ class NeighborKokkosExecute
 
 #ifdef LMP_KOKKOS_GPU
   template<int HalfNeigh, int Newton, int Tri>
-  __device__ inline
+#if defined(KOKKOS_ENABLE_CUDA) || defined(KOKKOS_ENABLE_HIP)
-  void build_ItemCuda(typename Kokkos::TeamPolicy<DeviceType>::member_type dev) const;
+  __device__
+#endif
+  inline
+  void build_ItemGPU(typename Kokkos::TeamPolicy<DeviceType>::member_type dev,
+                     size_t sharedsize) const;
 
   template<int HalfNeigh, int Newton, int Tri>
-  __device__ inline
+#if defined(KOKKOS_ENABLE_CUDA) || defined(KOKKOS_ENABLE_HIP)
-  void build_ItemSizeCuda(typename Kokkos::TeamPolicy<DeviceType>::member_type dev) const;
+  __device__
+#endif
+  inline
+  void build_ItemSizeGPU(typename Kokkos::TeamPolicy<DeviceType>::member_type dev,
+                         size_t sharedsize) const;
 #endif
 
   KOKKOS_INLINE_FUNCTION
@@ -388,10 +396,12 @@ struct NPairKokkosBuildFunctor {
     c.template build_Item<HALF_NEIGH,GHOST_NEWTON,TRI>(i);
   }
 #ifdef LMP_KOKKOS_GPU
-  __device__ inline
+#if defined(KOKKOS_ENABLE_CUDA) || defined(KOKKOS_ENABLE_HIP)
-
+  __device__
+#endif
+  inline
   void operator() (typename Kokkos::TeamPolicy<DeviceType>::member_type dev) const {
-    c.template build_ItemCuda<HALF_NEIGH,GHOST_NEWTON,TRI>(dev);
+    c.template build_ItemGPU<HALF_NEIGH,GHOST_NEWTON,TRI>(dev, sharedsize);
   }
   size_t team_shmem_size(const int team_size) const { (void) team_size; return sharedsize; }
 #endif
@@ -446,9 +456,12 @@ struct NPairKokkosBuildFunctorSize {
   }
 
 #ifdef LMP_KOKKOS_GPU
-  __device__ inline
+  #if defined(KOKKOS_ENABLE_CUDA) || defined(KOKKOS_ENABLE_HIP)
+  __device__
+#endif
+  inline
   void operator() (typename Kokkos::TeamPolicy<DeviceType>::member_type dev) const {
-    c.template build_ItemSizeCuda<HALF_NEIGH,GHOST_NEWTON,TRI>(dev);
+    c.template build_ItemSizeGPU<HALF_NEIGH,GHOST_NEWTON,TRI>(dev, sharedsize);
   }
   size_t team_shmem_size(const int team_size) const { (void) team_size; return sharedsize; }
 #endif

src/MC/fix_atom_swap.cpp

@@ -504,7 +504,7 @@ double FixAtomSwap::energy_full()
 
   if (force->pair) force->pair->compute(eflag,vflag);
 
-  if (atom->molecular) {
+  if (atom->molecular != Atom::ATOMIC) {
     if (force->bond) force->bond->compute(eflag,vflag);
     if (force->angle) force->angle->compute(eflag,vflag);
     if (force->dihedral) force->dihedral->compute(eflag,vflag);

src/MC/fix_bond_break.cpp

@@ -80,7 +80,7 @@ FixBondBreak::FixBondBreak(LAMMPS *lmp, int narg, char **arg) :
 
   // error check
 
-  if (atom->molecular != 1)
+  if (atom->molecular != Atom::MOLECULAR)
     error->all(FLERR,"Cannot use fix bond/break with non-molecular systems");
 
   // initialize Marsaglia RNG with processor-unique seed

src/MC/fix_bond_create.cpp

@@ -148,7 +148,7 @@ FixBondCreate::FixBondCreate(LAMMPS *lmp, int narg, char **arg) :
 
   // error check
 
-  if (atom->molecular != 1)
+  if (atom->molecular != Atom::MOLECULAR)
     error->all(FLERR,"Cannot use fix bond/create with non-molecular systems");
   if (iatomtype == jatomtype &&
       ((imaxbond != jmaxbond) || (inewtype != jnewtype)))

src/MC/fix_bond_swap.cpp

@@ -80,7 +80,7 @@ FixBondSwap::FixBondSwap(LAMMPS *lmp, int narg, char **arg) :
 
   // error check
 
-  if (atom->molecular != 1)
+  if (atom->molecular != Atom::MOLECULAR)
     error->all(FLERR,"Cannot use fix bond/swap with non-molecular systems");
 
   // create a new compute temp style

src/MC/fix_gcmc.cpp

@@ -2318,7 +2318,7 @@ double FixGCMC::energy_full()
 
   if (force->pair) force->pair->compute(eflag,vflag);
 
-  if (atom->molecular) {
+  if (atom->molecular != Atom::ATOMIC) {
     if (force->bond) force->bond->compute(eflag,vflag);
     if (force->angle) force->angle->compute(eflag,vflag);
     if (force->dihedral) force->dihedral->compute(eflag,vflag);

src/MC/fix_widom.cpp

@@ -1045,7 +1045,7 @@ double FixWidom::energy_full()
 
   if (force->pair) force->pair->compute(eflag,vflag);
 
-  if (atom->molecular) {
+  if (atom->molecular != Atom::ATOMIC) {
     if (force->bond) force->bond->compute(eflag,vflag);
     if (force->angle) force->angle->compute(eflag,vflag);
     if (force->dihedral) force->dihedral->compute(eflag,vflag);

src/MLIAP/mliap_model_nn.cpp

@@ -378,7 +378,7 @@ void MLIAPModelNN::compute_gradients(MLIAPData* data)
    egradient is derivative of energy w.r.t. parameters
    ---------------------------------------------------------------------- */
 
-void MLIAPModelNN::compute_gradgrads(class MLIAPData* data)
+void MLIAPModelNN::compute_gradgrads(class MLIAPData * /*data*/)
 {
   error->all(FLERR,"compute_gradgrads not implemented");
 }
@@ -388,7 +388,7 @@ void MLIAPModelNN::compute_gradgrads(class MLIAPData* data)
    egradient is derivative of energy w.r.t. parameters
    ---------------------------------------------------------------------- */
 
-void MLIAPModelNN::compute_force_gradients(class MLIAPData* data)
+void MLIAPModelNN::compute_force_gradients(class MLIAPData * /*data*/)
 {
   error->all(FLERR,"compute_force_gradients not implemented");
 }
@@ -397,7 +397,7 @@ void MLIAPModelNN::compute_force_gradients(class MLIAPData* data)
    count the number of non-zero entries in gamma matrix
    ---------------------------------------------------------------------- */
 
-int MLIAPModelNN::get_gamma_nnz(class MLIAPData* data)
+int MLIAPModelNN::get_gamma_nnz(class MLIAPData * /*data*/)
 {
   // todo: get_gamma_nnz
   return 0;

src/REPLICA/verlet_split.cpp

@@ -367,7 +367,7 @@ void VerletSplit::run(int n)
         timer->stamp(Timer::PAIR);
       }
 
-      if (atom->molecular) {
+      if (atom->molecular != Atom::ATOMIC) {
         if (force->bond) force->bond->compute(eflag,vflag);
         if (force->angle) force->angle->compute(eflag,vflag);
         if (force->dihedral) force->dihedral->compute(eflag,vflag);

src/USER-INTEL/npair_full_bin_intel.cpp

@@ -74,7 +74,7 @@ fbi(NeighList *list, IntelBuffers<flt_t,acc_t> *buffers) {
                      _fix->nbor_pack_width());
 
   int need_ic = 0;
-  if (atom->molecular)
+  if (atom->molecular != Atom::ATOMIC)
     dminimum_image_check(need_ic, neighbor->cutneighmax, neighbor->cutneighmax,
                          neighbor->cutneighmax);
 

src/USER-INTEL/npair_half_bin_newton_intel.cpp

@@ -73,7 +73,7 @@ hbni(NeighList *list, IntelBuffers<flt_t,acc_t> *buffers) {
   buffers->grow_list(list, atom->nlocal, comm->nthreads, 0, off_end);
 
   int need_ic = 0;
-  if (atom->molecular)
+  if (atom->molecular != Atom::ATOMIC)
     dminimum_image_check(need_ic, neighbor->cutneighmax, neighbor->cutneighmax,
                          neighbor->cutneighmax);
 

src/USER-INTEL/npair_half_bin_newton_tri_intel.cpp

@@ -73,7 +73,7 @@ hbnti(NeighList *list, IntelBuffers<flt_t,acc_t> *buffers) {
   buffers->grow_list(list, atom->nlocal, comm->nthreads, 0, off_end);
 
   int need_ic = 0;
-  if (atom->molecular)
+  if (atom->molecular != Atom::ATOMIC)
     dminimum_image_check(need_ic, neighbor->cutneighmax, neighbor->cutneighmax,
                          neighbor->cutneighmax);
 

src/USER-INTEL/verlet_lrt_intel.cpp

@@ -316,7 +316,7 @@ void VerletLRTIntel::run(int n)
       timer->stamp(Timer::PAIR);
     }
 
-    if (atom->molecular) {
+    if (atom->molecular != Atom::ATOMIC) {
       if (force->bond) force->bond->compute(eflag,vflag);
       if (force->angle) force->angle->compute(eflag,vflag);
       if (force->dihedral) force->dihedral->compute(eflag,vflag);

src/USER-MEAMC/meam.h

@@ -17,7 +17,7 @@
 #include "math_const.h"         // IWYU pragma: export
 
 #include <cmath>
-#include <cstring>
+#include <string>
 
 #define maxelt 5
 
@@ -242,29 +242,29 @@ public:
   // only use single-element lattices if single=true
   // return false on failure
   // return true and set lat on success
-  static bool str_to_lat(const char* str, bool single, lattice_t& lat)
+  static bool str_to_lat(const std::string & str, bool single, lattice_t& lat)
   {
-    if (strcmp(str,"fcc") == 0) lat = FCC;
+    if (str == "fcc") lat = FCC;
-    else if (strcmp(str,"bcc") == 0) lat = BCC;
+    else if (str == "bcc") lat = BCC;
-    else if (strcmp(str,"hcp") == 0) lat = HCP;
+    else if (str == "hcp") lat = HCP;
-    else if (strcmp(str,"dim") == 0) lat = DIM;
+    else if (str == "dim") lat = DIM;
-    else if (strcmp(str,"dia") == 0) lat = DIA;
+    else if (str == "dia") lat = DIA;
-    else if (strcmp(str,"dia3") == 0) lat = DIA3;
+    else if (str == "dia3") lat = DIA3;
-    else if (strcmp(str,"lin") == 0) lat = LIN;
+    else if (str == "lin") lat = LIN;
-    else if (strcmp(str,"zig") == 0) lat = ZIG;
+    else if (str == "zig") lat = ZIG;
-    else if (strcmp(str,"tri") == 0) lat = TRI;
+    else if (str == "tri") lat = TRI;
     else {
       if (single)
         return false;
 
-      if (strcmp(str,"b1")  == 0) lat = B1;
+      if (str == "b1") lat = B1;
-      else if (strcmp(str,"c11") == 0) lat = C11;
+      else if (str == "c11") lat = C11;
-      else if (strcmp(str,"l12") == 0) lat = L12;
+      else if (str == "l12") lat = L12;
-      else if (strcmp(str,"b2")  == 0) lat = B2;
+      else if (str == "b2") lat = B2;
-      else if (strcmp(str,"ch4")  == 0) lat = CH4;
+      else if (str == "ch4") lat = CH4;
-      else if (strcmp(str,"lin")  == 0) lat =LIN;
+      else if (str == "lin") lat =LIN;
-      else if (strcmp(str,"zig")  == 0) lat = ZIG;
+      else if (str == "zig") lat = ZIG;
-      else if (strcmp(str,"tri")  == 0) lat = TRI;
+      else if (str == "tri") lat = TRI;
       else return false;
     }
     return true;

src/USER-MEAMC/pair_meamc.cpp

@@ -17,7 +17,7 @@
 
 #include "pair_meamc.h"
 
-#include <cstring>
+#include <memory>
 
 #include "meam.h"
 #include "atom.h"
@@ -28,6 +28,8 @@
 #include "neigh_request.h"
 #include "memory.h"
 #include "error.h"
+#include "potential_file_reader.h"
+#include "tokenizer.h"
 
 using namespace LAMMPS_NS;
 
@@ -54,8 +56,6 @@ PairMEAMC::PairMEAMC(LAMMPS *lmp) : Pair(lmp)
   allocated = 0;
 
   nelements = 0;
-  elements = nullptr;
-  mass = nullptr;
   meam_inst = new MEAM(memory);
   scale = nullptr;
 
@@ -74,10 +74,6 @@ PairMEAMC::~PairMEAMC()
 {
   delete meam_inst;
 
-  for (int i = 0; i < nelements; i++) delete [] elements[i];
-  delete [] elements;
-  delete [] mass;
-
   if (allocated) {
     memory->destroy(setflag);
     memory->destroy(cutsq);
@@ -192,7 +188,7 @@ void PairMEAMC::allocate()
    global settings
 ------------------------------------------------------------------------- */
 
-void PairMEAMC::settings(int narg, char **/*arg*/)
+void PairMEAMC::settings(int narg, char ** /*arg*/)
 {
   if (narg != 0) error->all(FLERR,"Illegal pair_style command");
 }
@@ -247,9 +243,8 @@ void PairMEAMC::coeff(int narg, char **arg)
   // elements = list of unique element names
 
   if (nelements) {
-    for (int i = 0; i < nelements; i++) delete [] elements[i];
+    elements.clear();
-    delete [] elements;
+    mass.clear();
-    delete [] mass;
   }
 
   nelements = paridx - 3;
@@ -258,13 +253,10 @@ void PairMEAMC::coeff(int narg, char **arg)
     error->all(FLERR,fmt::format("Too many elements extracted from MEAM "
                                  "library (current limit: {}). Increase "
                                  "'maxelt' in meam.h and recompile.", maxelt));
-  elements = new char*[nelements];
-  mass = new double[nelements];
 
   for (int i = 0; i < nelements; i++) {
-    n = strlen(arg[i+3]) + 1;
+    elements.push_back(arg[i+3]);
-    elements[i] = new char[n];
+    mass.push_back(0.0);
-    strcpy(elements[i],arg[i+3]);
   }
 
   // read MEAM library and parameter files
@@ -281,7 +273,7 @@ void PairMEAMC::coeff(int narg, char **arg)
     m = i - (4+nelements) + 1;
     int j;
     for (j = 0; j < nelements; j++)
-      if (strcmp(arg[i],elements[j]) == 0) break;
+      if (elements[j] == arg[i]) break;
     if (j < nelements) map[m] = j;
     else if (strcmp(arg[i],"NULL") == 0) map[m] = -1;
     else error->all(FLERR,"Incorrect args for pair coefficients");
@@ -358,272 +350,214 @@ double PairMEAMC::init_one(int i, int j)
 void PairMEAMC::read_files(const std::string &globalfile,
                            const std::string &userfile)
 {
+  read_global_meamc_file(globalfile);
+  read_user_meamc_file(userfile);
+}
+
+/* ---------------------------------------------------------------------- */
+
+void PairMEAMC::read_global_meamc_file(const std::string &globalfile)
+{
+  // allocate parameter arrays
+  std::vector<lattice_t> lat(nelements);
+  std::vector<int> ielement(nelements);
+  std::vector<int> ibar(nelements);
+  std::vector<double> z(nelements);
+  std::vector<double> atwt(nelements);
+  std::vector<double> alpha(nelements);
+  std::vector<double> b0(nelements);
+  std::vector<double> b1(nelements);
+  std::vector<double> b2(nelements);
+  std::vector<double> b3(nelements);
+  std::vector<double> alat(nelements);
+  std::vector<double> esub(nelements);
+  std::vector<double> asub(nelements);
+  std::vector<double> t0(nelements);
+  std::vector<double> t1(nelements);
+  std::vector<double> t2(nelements);
+  std::vector<double> t3(nelements);
+  std::vector<double> rozero(nelements);
+  std::vector<bool> found(nelements, false);
+
   // open global meamf file on proc 0
 
-  FILE *fp;
   if (comm->me == 0) {
-    fp = utils::open_potential(globalfile,lmp,nullptr);
+    PotentialFileReader reader(lmp, globalfile, "MEAM", " library");
-    if (fp == nullptr)
+    char * line;
-      error->one(FLERR,fmt::format("Cannot open MEAM potential file {}",
-                                   globalfile));
-  }
 
-  // allocate parameter arrays
+    const int params_per_line = 19;
-
-  int params_per_line = 19;
-
-  lattice_t *lat = new lattice_t[nelements];
-  int *ielement = new int[nelements];
-  int *ibar = new int[nelements];
-  double *z = new double[nelements];
-  double *atwt = new double[nelements];
-  double *alpha = new double[nelements];
-  double *b0 = new double[nelements];
-  double *b1 = new double[nelements];
-  double *b2 = new double[nelements];
-  double *b3 = new double[nelements];
-  double *alat = new double[nelements];
-  double *esub = new double[nelements];
-  double *asub = new double[nelements];
-  double *t0 = new double[nelements];
-  double *t1 = new double[nelements];
-  double *t2 = new double[nelements];
-  double *t3 = new double[nelements];
-  double *rozero = new double[nelements];
-
-  bool *found = new bool[nelements];
-  for (int i = 0; i < nelements; i++) found[i] = false;
-
-  // read each set of params from global MEAM file
-  // one set of params can span multiple lines
-  // store params if element name is in element list
-  // if element name appears multiple times, only store 1st entry
-
-  if (comm->me == 0) {
     int nset = 0;
-    char **words = new char*[params_per_line+1];
+
-    char line[MAXLINE];
+    while ((line = reader.next_line(params_per_line))) {
-    while (1) {
+      try {
-      char *ptr;
+        ValueTokenizer values(line, "' \t\n\r\f");
-      ptr = fgets(line,MAXLINE,fp);
+
-      if (ptr == nullptr) {
+        // read each set of params from global MEAM file
-        fclose(fp);
+        // one set of params can span multiple lines
-        break;
+        // store params if element name is in element list
+        // if element name appears multiple times, only store 1st entry
+        std::string element = values.next_string();
+
+        // skip if element name isn't in element list
+
+        int index;
+        for (index = 0; index < nelements; index++)
+          if (elements[index] == element) break;
+        if (index == nelements) continue;
+
+        // skip if element already appeared (technically error in library file, but always ignored)
+
+        if (found[index]) continue;
+        found[index] = true;
+
+        // map lat string to an integer
+        std::string lattice_type = values.next_string();
+
+        if (!MEAM::str_to_lat(lattice_type.c_str(), true, lat[index]))
+          error->one(FLERR,fmt::format("Unrecognized lattice type in MEAM "
+                                       "library file: {}", lattice_type));
+
+        // store parameters
+
+        z[index] = values.next_double();
+        ielement[index] = values.next_int();
+        atwt[index] = values.next_double();
+        alpha[index] = values.next_double();
+        b0[index] = values.next_double();
+        b1[index] = values.next_double();
+        b2[index] = values.next_double();
+        b3[index] = values.next_double();
+        alat[index] = values.next_double();
+        esub[index] = values.next_double();
+        asub[index] = values.next_double();
+        t0[index] = values.next_double();
+        t1[index] = values.next_double();
+        t2[index] = values.next_double();
+        t3[index] = values.next_double();
+        rozero[index] = values.next_double();
+        ibar[index] = values.next_int();
+
+        if (!isone(t0[index]))
+          error->one(FLERR,"Unsupported parameter in MEAM library file: t0!=1");
+
+        // z given is ignored: if this is mismatched, we definitely won't do what the user said -> fatal error
+        if (z[index] != MEAM::get_Zij(lat[index]))
+          error->one(FLERR,"Mismatched parameter in MEAM library file: z!=lat");
+
+        nset++;
+      } catch (TokenizerException &e) {
+        error->one(FLERR, e.what());
       }
-
-      // strip comment, skip line if blank
-
-      if ((ptr = strchr(line,'#'))) *ptr = '\0';
-      int nwords = utils::count_words(line);
-      if (nwords == 0) continue;
-
-      // concatenate additional lines until have params_per_line words
-
-      while (nwords < params_per_line) {
-        int n = strlen(line);
-        ptr = fgets(&line[n],MAXLINE-n,fp);
-        if (ptr == nullptr) {
-          fclose(fp);
-          break;
-        }
-        if ((ptr = strchr(line,'#'))) *ptr = '\0';
-        nwords = utils::count_words(line);
-      }
-
-      if (nwords != params_per_line)
-        error->one(FLERR,"Incorrect format in MEAM library file");
-
-      // words = ptrs to all words in line
-      // strip single and double quotes from words
-
-      nwords = 0;
-      words[nwords++] = strtok(line,"' \t\n\r\f");
-      while ((words[nwords++] = strtok(nullptr,"' \t\n\r\f"))) continue;
-
-      // skip if element name isn't in element list
-
-      int index;
-      for (index = 0; index < nelements; index++)
-        if (strcmp(words[0],elements[index]) == 0) break;
-      if (index == nelements) continue;
-
-      // skip if element already appeared (technically error in library file, but always ignored)
-
-      if (found[index] == true) continue;
-      found[index] = true;
-
-      // map lat string to an integer
-
-      if (!MEAM::str_to_lat(words[1], true, lat[index]))
-        error->one(FLERR,fmt::format("Unrecognized lattice type in MEAM "
-                                     "library file: {}", words[1]));
-
-      // store parameters
-
-      z[index] = atof(words[2]);
-      ielement[index] = atoi(words[3]);
-      atwt[index] = atof(words[4]);
-      alpha[index] = atof(words[5]);
-      b0[index] = atof(words[6]);
-      b1[index] = atof(words[7]);
-      b2[index] = atof(words[8]);
-      b3[index] = atof(words[9]);
-      alat[index] = atof(words[10]);
-      esub[index] = atof(words[11]);
-      asub[index] = atof(words[12]);
-      t0[index] = atof(words[13]);
-      t1[index] = atof(words[14]);
-      t2[index] = atof(words[15]);
-      t3[index] = atof(words[16]);
-      rozero[index] = atof(words[17]);
-      ibar[index] = atoi(words[18]);
-
-      if (!isone(t0[index]))
-        error->one(FLERR,"Unsupported parameter in MEAM library file: t0!=1");
-
-      // z given is ignored: if this is mismatched, we definitely won't do what the user said -> fatal error
-      if (z[index] != MEAM::get_Zij(lat[index]))
-        error->one(FLERR,"Mismatched parameter in MEAM library file: z!=lat");
-
-      nset++;
     }
 
     // error if didn't find all elements in file
 
     if (nset != nelements) {
-      char str[128] = "Did not find all elements in MEAM library file, missing:";
+      std::string msg = "Did not find all elements in MEAM library file, missing:";
       for (int i = 0; i < nelements; i++)
         if (!found[i]) {
-          strcat(str," ");
+          msg += " ";
-          strcat(str,elements[i]);
+          msg += elements[i];
         }
-      error->one(FLERR,str);
+      error->one(FLERR,msg);
     }
-
-    delete [] words;
   }
 
   // distribute complete parameter sets
-  MPI_Bcast(lat, nelements, MPI_INT, 0, world);
+  MPI_Bcast(lat.data(), nelements, MPI_INT, 0, world);
-  MPI_Bcast(ielement, nelements, MPI_INT, 0, world);
+  MPI_Bcast(ielement.data(), nelements, MPI_INT, 0, world);
-  MPI_Bcast(ibar, nelements, MPI_INT, 0, world);
+  MPI_Bcast(ibar.data(), nelements, MPI_INT, 0, world);
-  MPI_Bcast(z, nelements, MPI_DOUBLE, 0, world);
+  MPI_Bcast(z.data(), nelements, MPI_DOUBLE, 0, world);
-  MPI_Bcast(atwt, nelements, MPI_DOUBLE, 0, world);
+  MPI_Bcast(atwt.data(), nelements, MPI_DOUBLE, 0, world);
-  MPI_Bcast(alpha, nelements, MPI_DOUBLE, 0, world);
+  MPI_Bcast(alpha.data(), nelements, MPI_DOUBLE, 0, world);
-  MPI_Bcast(b0, nelements, MPI_DOUBLE, 0, world);
+  MPI_Bcast(b0.data(), nelements, MPI_DOUBLE, 0, world);
-  MPI_Bcast(b1, nelements, MPI_DOUBLE, 0, world);
+  MPI_Bcast(b1.data(), nelements, MPI_DOUBLE, 0, world);
-  MPI_Bcast(b2, nelements, MPI_DOUBLE, 0, world);
+  MPI_Bcast(b2.data(), nelements, MPI_DOUBLE, 0, world);
-  MPI_Bcast(b3, nelements, MPI_DOUBLE, 0, world);
+  MPI_Bcast(b3.data(), nelements, MPI_DOUBLE, 0, world);
-  MPI_Bcast(alat, nelements, MPI_DOUBLE, 0, world);
+  MPI_Bcast(alat.data(), nelements, MPI_DOUBLE, 0, world);
-  MPI_Bcast(esub, nelements, MPI_DOUBLE, 0, world);
+  MPI_Bcast(esub.data(), nelements, MPI_DOUBLE, 0, world);
-  MPI_Bcast(asub, nelements, MPI_DOUBLE, 0, world);
+  MPI_Bcast(asub.data(), nelements, MPI_DOUBLE, 0, world);
-  MPI_Bcast(t0, nelements, MPI_DOUBLE, 0, world);
+  MPI_Bcast(t0.data(), nelements, MPI_DOUBLE, 0, world);
-  MPI_Bcast(t1, nelements, MPI_DOUBLE, 0, world);
+  MPI_Bcast(t1.data(), nelements, MPI_DOUBLE, 0, world);
-  MPI_Bcast(t2, nelements, MPI_DOUBLE, 0, world);
+  MPI_Bcast(t2.data(), nelements, MPI_DOUBLE, 0, world);
-  MPI_Bcast(t3, nelements, MPI_DOUBLE, 0, world);
+  MPI_Bcast(t3.data(), nelements, MPI_DOUBLE, 0, world);
-  MPI_Bcast(rozero, nelements, MPI_DOUBLE, 0, world);
+  MPI_Bcast(rozero.data(), nelements, MPI_DOUBLE, 0, world);
 
   // pass element parameters to MEAM package
 
-  meam_inst->meam_setup_global(nelements,lat,ielement,atwt,alpha,b0,b1,b2,b3,
+  meam_inst->meam_setup_global(nelements, lat.data(), ielement.data(), atwt.data(),
-                       alat,esub,asub,t0,t1,t2,t3,rozero,ibar);
+                               alpha.data(), b0.data(), b1.data(), b2.data(), b3.data(),
+                               alat.data(), esub.data(), asub.data(), t0.data(), t1.data(),
+                               t2.data(), t3.data(), rozero.data(), ibar.data());
 
   // set element masses
 
   for (int i = 0; i < nelements; i++) mass[i] = atwt[i];
+}
 
-  // clean-up memory
+/* ---------------------------------------------------------------------- */
-
-  delete [] lat;
-  delete [] ielement;
-  delete [] ibar;
-  delete [] z;
-  delete [] atwt;
-  delete [] alpha;
-  delete [] b0;
-  delete [] b1;
-  delete [] b2;
-  delete [] b3;
-  delete [] alat;
-  delete [] esub;
-  delete [] asub;
-  delete [] t0;
-  delete [] t1;
-  delete [] t2;
-  delete [] t3;
-  delete [] rozero;
-  delete [] found;
 
+void PairMEAMC::read_user_meamc_file(const std::string &userfile)
+{
   // done if user param file is "NULL"
 
   if (userfile == "NULL") return;
 
   // open user param file on proc 0
 
+  std::shared_ptr<PotentialFileReader> reader;
+
   if (comm->me == 0) {
-    fp = utils::open_potential(userfile,lmp,nullptr);
+    reader = std::make_shared<PotentialFileReader>(lmp, userfile, "MEAM");
-    if (fp == nullptr)
-      error->one(FLERR,fmt::format("Cannot open MEAM potential file {}",
-                                   userfile));
   }
 
   // read settings
   // pass them one at a time to MEAM package
   // match strings to list of corresponding ints
+  char * line = nullptr;
+  char buffer[MAXLINE];
 
-  int maxparams = 6;
-  char **params = new char*[maxparams];
   while (1) {
     int which;
     int nindex, index[3];
     double value;
-    char line[MAXLINE];
     int nline;
-    char *ptr;
     if (comm->me == 0) {
-      ptr = fgets(line,MAXLINE,fp);
+      line = reader->next_line();
-      if (ptr == nullptr) {
+      if (line == nullptr) {
-        fclose(fp);
         nline = -1;
       } else nline = strlen(line) + 1;
+    } else {
+      line = buffer;
     }
+
     MPI_Bcast(&nline,1,MPI_INT,0,world);
     if (nline<0) break;
     MPI_Bcast(line,nline,MPI_CHAR,0,world);
 
-    // strip comment, skip line if blank
+    ValueTokenizer values(line, "=(), '\t\n\r\f");
-
+    int nparams = values.count();
-    if ((ptr = strchr(line,'#'))) *ptr = '\0';
+    std::string keyword = values.next_string();
-    if (utils::count_words(line) == 0) continue;
-
-    // params = ptrs to all fields in line
-
-    int nparams = 0;
-    params[nparams++] = strtok(line,"=(), '\t\n\r\f");
-    while (nparams < maxparams &&
-           (params[nparams++] = strtok(nullptr,"=(), '\t\n\r\f")))
-      continue;
-    nparams--;
 
     for (which = 0; which < nkeywords; which++)
-      if (strcmp(params[0],keywords[which]) == 0) break;
+      if (keyword == keywords[which]) break;
     if (which == nkeywords)
       error->all(FLERR,fmt::format("Keyword {} in MEAM parameter file not "
-                                   "recognized", params[0]));
+                                   "recognized", keyword));
 
     nindex = nparams - 2;
-    for (int i = 0; i < nindex; i++) index[i] = atoi(params[i+1]) - 1;
+    for (int i = 0; i < nindex; i++) index[i] = values.next_int() - 1;
 
     // map lattce_meam value to an integer
-
     if (which == 4) {
+      std::string lattice_type = values.next_string();
       lattice_t latt;
-      if (!MEAM::str_to_lat(params[nparams-1], false, latt))
+      if (!MEAM::str_to_lat(lattice_type, false, latt))
         error->all(FLERR, fmt::format("Unrecognized lattice type in MEAM "
-                                      "parameter file: {}", params[nparams-1]));
+                                      "parameter file: {}", lattice_type));
       value = latt;
     }
-    else value = atof(params[nparams-1]);
+    else value = values.next_double();
 
     // pass single setting to MEAM package
 
@@ -634,13 +568,10 @@ void PairMEAMC::read_files(const std::string &globalfile,
               "is out of range (please report a bug)",
               "expected more indices",
               "has out of range element index"};
-      char str[256];
       if ((errorflag < 0) || (errorflag > 3)) errorflag = 0;
-      snprintf(str,256,"Error in MEAM parameter file: keyword %s %s",params[0],descr[errorflag]);
+      error->all(FLERR, fmt::format("Error in MEAM parameter file: keyword {} {}", keyword, descr[errorflag]));
-      error->all(FLERR,str);
     }
   }
-  delete [] params;
 }
 
 /* ---------------------------------------------------------------------- */

src/USER-MEAMC/pair_meamc.h

@@ -22,7 +22,8 @@ PairStyle(meam,PairMEAMC)
 #define LMP_PAIR_MEAMC_H
 
 #include "pair.h"
-
+#include <vector>
+#include <string>
 
 namespace LAMMPS_NS {
 
@@ -48,14 +49,16 @@ class PairMEAMC : public Pair {
   class MEAM *meam_inst;
   double cutmax;                // max cutoff for all elements
   int nelements;                // # of unique elements
-  char **elements;              // names of unique elements
+  std::vector<std::string> elements; // names of unique elements
-  double *mass;                 // mass of each element
+  std::vector<double> mass;          // mass of each element
 
   int *map;                     // mapping from atom types (1-indexed) to elements (1-indexed)
   double **scale;               // scaling factor for adapt
 
   void allocate();
   void read_files(const std::string &, const std::string &);
+  void read_global_meamc_file(const std::string &);
+  void read_user_meamc_file(const std::string &);
   void neigh_strip(int, int *, int *, int **);
 };
 

src/USER-OMP/thr_data.cpp

@@ -139,20 +139,6 @@ void ThrData::init_adp(int nall, double *rho, double **mu, double **lambda)
 
 /* ---------------------------------------------------------------------- */
 
-void ThrData::init_cdeam(int nall, double *rho, double *rhoB, double *D_values)
-{
-  init_eam(nall, rho);
-
-  if (nall >= 0 && rhoB && D_values) {
-    _rhoB = rhoB + _tid*nall;
-    _D_values = D_values + _tid*nall;
-    memset(_rhoB, 0, nall*sizeof(double));
-    memset(_D_values, 0, nall*sizeof(double));
-  }
-}
-
-/* ---------------------------------------------------------------------- */
-
 void ThrData::init_eim(int nall, double *rho, double *fp)
 {
   init_eam(nall, rho);

src/USER-OMP/thr_data.h

@@ -52,7 +52,6 @@ class ThrData {
 
   // setup and erase per atom arrays
   void init_adp(int, double *, double **, double **); // ADP (+ EAM)
-  void init_cdeam(int, double *, double *, double *); // CDEAM (+ EAM)
   void init_eam(int, double *);                       // EAM
   void init_eim(int, double *, double *);             // EIM (+ EAM)
 

src/arg_info.h

@@ -67,14 +67,14 @@ namespace LAMMPS_NS {
     /*! get dimension of reference
      *
      * This will return either 0, 1, 2 depending on whether the
-     * reference has no, one or two "[<number>]" postfixes.
+     * reference has no, one or two "[{number}]" postfixes.
      *
      * \return  integer with the dimensionality of the reference */
     int  get_dim()  const { return dim; }
 
     /*! get index of first dimension
      *
-     * This will return the <number> in the first "[<number>]"
+     * This will return the number in the first "[{number}]"
      * postfix or 0 if there is no postfix.
      *
      * \return  integer with index or the postfix or 0 */
@@ -82,7 +82,7 @@ namespace LAMMPS_NS {
 
     /*! get index of second dimension
      *
-     * This will return the <number> in the second "[<number>]"
+     * This will return the number in the second "[{number}]"
      * postfix or -1 if there is no second postfix.
      *
      * \return  integer with index of the postfix or -1 */

src/atom.cpp

@@ -341,11 +341,8 @@ void Atom::settings(Atom *old)
   map_style = old->map_style;
   sortfreq = old->sortfreq;
   userbinsize = old->userbinsize;
-  if (old->firstgroupname) {
+  if (old->firstgroupname)
-    int n = strlen(old->firstgroupname) + 1;
+    firstgroupname = utils::strdup(old->firstgroupname);
-    firstgroupname = new char[n];
-    strcpy(firstgroupname,old->firstgroupname);
-  }
 }
 
 /* ----------------------------------------------------------------------
@@ -532,9 +529,7 @@ void Atom::add_peratom(const char *name, void *address,
       memory->srealloc(peratom,maxperatom*sizeof(PerAtom),"atom:peratom");
   }
 
-  int n = strlen(name) + 1;
+  peratom[nperatom].name = utils::strdup(name);
-  peratom[nperatom].name = new char[n];
-  strcpy(peratom[nperatom].name,name);
   peratom[nperatom].address = address;
   peratom[nperatom].datatype = datatype;
   peratom[nperatom].cols = cols;
@@ -582,9 +577,7 @@ void Atom::add_peratom_vary(const char *name, void *address,
       memory->srealloc(peratom,maxperatom*sizeof(PerAtom),"atom:peratom");
   }
 
-  int n = strlen(name) + 1;
+  peratom[nperatom].name = utils::strdup(name);
-  peratom[nperatom].name = new char[n];
-  strcpy(peratom[nperatom].name,name);
   peratom[nperatom].address = address;
   peratom[nperatom].datatype = datatype;
   peratom[nperatom].cols = -1;
@@ -670,7 +663,7 @@ void Atom::create_avec(const std::string &style, int narg, char **arg, int trysu
   // map style will be reset to array vs hash to by map_init()
 
   molecular = avec->molecular;
-  if (molecular && tag_enable == 0)
+  if ((molecular != Atom::ATOMIC) && (tag_enable == 0))
     error->all(FLERR,"Atom IDs must be used for molecular systems");
   if (molecular != Atom::ATOMIC) map_style = MAP_YES;
 }
@@ -815,9 +808,7 @@ void Atom::modify_params(int narg, char **arg)
         delete [] firstgroupname;
         firstgroupname = nullptr;
       } else {
-        int n = strlen(arg[iarg+1]) + 1;
+        firstgroupname = utils::strdup(arg[iarg+1]);
-        firstgroupname = new char[n];
-        strcpy(firstgroupname,arg[iarg+1]);
         sortfreq = 0;
       }
       iarg += 2;
@@ -1114,13 +1105,13 @@ void Atom::data_atoms(int n, char *buf, tagint id_offset, tagint mol_offset,
   for (int i = 0; i < n; i++) {
     next = strchr(buf,'\n');
 
-    values[0] = strtok(buf," \t\n\r\f");
+    for (m = 0; m < nwords; m++) {
-    if (values[0] == nullptr)
+      buf += strspn(buf," \t\n\r\f");
-      error->all(FLERR,"Incorrect atom format in data file");
+      buf[strcspn(buf," \t\n\r\f")] = '\0';
-    for (m = 1; m < nwords; m++) {
+      if (strlen(buf) == 0)
-      values[m] = strtok(nullptr," \t\n\r\f");
-      if (values[m] == nullptr)
         error->all(FLERR,"Incorrect atom format in data file");
+      values[m] = buf;
+      buf += strlen(buf)+1;
     }
 
     int imx = 0, imy = 0, imz = 0;
@@ -1217,9 +1208,12 @@ void Atom::data_vels(int n, char *buf, tagint id_offset)
   for (int i = 0; i < n; i++) {
     next = strchr(buf,'\n');
 
-    values[0] = strtok(buf," \t\n\r\f");
+    for (j = 0; j < nwords; j++) {
-    for (j = 1; j < nwords; j++)
+      buf += strspn(buf," \t\n\r\f");
-      values[j] = strtok(nullptr," \t\n\r\f");
+      buf[strcspn(buf," \t\n\r\f")] = '\0';
+      values[j] = buf;
+      buf += strlen(buf)+1;
+    }
 
     tagdata = ATOTAGINT(values[0]) + id_offset;
     if (tagdata <= 0 || tagdata > map_tag_max)
@@ -1569,9 +1563,12 @@ void Atom::data_bonus(int n, char *buf, AtomVec *avec_bonus, tagint id_offset)
   for (int i = 0; i < n; i++) {
     next = strchr(buf,'\n');
 
-    values[0] = strtok(buf," \t\n\r\f");
+    for (j = 0; j < nwords; j++) {
-    for (j = 1; j < nwords; j++)
+      buf += strspn(buf," \t\n\r\f");
-      values[j] = strtok(nullptr," \t\n\r\f");
+      buf[strcspn(buf," \t\n\r\f")] = '\0';
+      values[j] = buf;
+      buf += strlen(buf)+1;
+    }
 
     tagdata = ATOTAGINT(values[0]) + id_offset;
     if (tagdata <= 0 || tagdata > map_tag_max)
@@ -1611,8 +1608,10 @@ void Atom::data_bodies(int n, char *buf, AtomVec *avec_body, tagint id_offset)
   // else skip values
 
   for (int i = 0; i < n; i++) {
-    if (i == 0) tagdata = ATOTAGINT(strtok(buf," \t\n\r\f")) + id_offset;
+    buf += strspn(buf," \t\n\r\f");
-    else tagdata = ATOTAGINT(strtok(nullptr," \t\n\r\f")) + id_offset;
+    buf[strcspn(buf," \t\n\r\f")] = '\0';
+    tagdata = utils::tnumeric(FLERR,buf,false,lmp) + id_offset;
+    buf += strlen(buf)+1;
 
     if (tagdata <= 0 || tagdata > map_tag_max)
       error->one(FLERR,"Invalid atom ID in Bodies section of data file");
@@ -1622,8 +1621,15 @@ void Atom::data_bodies(int n, char *buf, AtomVec *avec_body, tagint id_offset)
     else
       error->one(FLERR,"Duplicate atom ID in Bodies section of data file");
 
-    ninteger = utils::inumeric(FLERR,strtok(nullptr," \t\n\r\f"),false,lmp);
+    buf += strspn(buf," \t\n\r\f");
-    ndouble = utils::inumeric(FLERR,strtok(nullptr," \t\n\r\f"),false,lmp);
+    buf[strcspn(buf," \t\n\r\f")] = '\0';
+    ninteger = utils::inumeric(FLERR,buf,false,lmp);
+    buf += strlen(buf)+1;
+
+    buf += strspn(buf," \t\n\r\f");
+    buf[strcspn(buf," \t\n\r\f")] = '\0';
+    ndouble = utils::inumeric(FLERR,buf,false,lmp);
+    buf += strlen(buf)+1;
 
     if ((m = map(tagdata)) >= 0) {
       if (ninteger > maxint) {
@@ -1637,17 +1643,29 @@ void Atom::data_bodies(int n, char *buf, AtomVec *avec_body, tagint id_offset)
         dvalues = new double[maxdouble];
       }
 
-      for (j = 0; j < ninteger; j++)
+      for (j = 0; j < ninteger; j++) {
-        ivalues[j] = utils::inumeric(FLERR,strtok(nullptr," \t\n\r\f"),false,lmp);
+        buf += strspn(buf," \t\n\r\f");
-      for (j = 0; j < ndouble; j++)
+        buf[strcspn(buf," \t\n\r\f")] = '\0';
-        dvalues[j] = utils::numeric(FLERR,strtok(nullptr," \t\n\r\f"),false,lmp);
+        ivalues[j] = utils::inumeric(FLERR,buf,false,lmp);
+        buf += strlen(buf)+1;
+      }
+
+      for (j = 0; j < ndouble; j++) {
+        buf += strspn(buf," \t\n\r\f");
+        buf[strcspn(buf," \t\n\r\f")] = '\0';
+        dvalues[j] = utils::numeric(FLERR,buf,false,lmp);
+        buf += strlen(buf)+1;
+      }
 
       avec_body->data_body(m,ninteger,ndouble,ivalues,dvalues);
 
     } else {
       nvalues = ninteger + ndouble;    // number of values to skip
-      for (j = 0; j < nvalues; j++)
+      for (j = 0; j < nvalues; j++) {
-        strtok(nullptr," \t\n\r\f");
+        buf += strspn(buf," \t\n\r\f");
+        buf[strcspn(buf," \t\n\r\f")] = '\0';
+        buf += strlen(buf)+1;
+      }
     }
   }
 
@@ -2400,9 +2418,7 @@ int Atom::add_custom(const char *name, int flag)
     nivector++;
     iname = (char **) memory->srealloc(iname,nivector*sizeof(char *),
                                        "atom:iname");
-    int n = strlen(name) + 1;
+    iname[index] = utils::strdup(name);
-    iname[index] = new char[n];
-    strcpy(iname[index],name);
     ivector = (int **) memory->srealloc(ivector,nivector*sizeof(int *),
                                         "atom:ivector");
     memory->create(ivector[index],nmax,"atom:ivector");
@@ -2411,9 +2427,7 @@ int Atom::add_custom(const char *name, int flag)
     ndvector++;
     dname = (char **) memory->srealloc(dname,ndvector*sizeof(char *),
                                        "atom:dname");
-    int n = strlen(name) + 1;
+    dname[index] = utils::strdup(name);
-    dname[index] = new char[n];
-    strcpy(dname[index],name);
     dvector = (double **) memory->srealloc(dvector,ndvector*sizeof(double *),
                                            "atom:dvector");
     memory->create(dvector[index],nmax,"atom:dvector");

src/atom_vec.cpp

@@ -154,11 +154,8 @@ void AtomVec::store_args(int narg, char **arg)
   nargcopy = narg;
   if (nargcopy) argcopy = new char*[nargcopy];
   else argcopy = nullptr;
-  for (int i = 0; i < nargcopy; i++) {
+  for (int i = 0; i < nargcopy; i++)
-    int n = strlen(arg[i]) + 1;
+    argcopy[i] = utils::strdup(arg[i]);
-    argcopy[i] = new char[n];
-    strcpy(argcopy[i],arg[i]);
-  }
 }
 
 /* ----------------------------------------------------------------------

src/atom_vec_hybrid.cpp

@@ -109,8 +109,7 @@ void AtomVecHybrid::process_args(int narg, char **arg)
       if (strcmp(arg[iarg],keywords[i]) == 0)
         error->all(FLERR,"Atom style hybrid cannot use same atom style twice");
     styles[nstyles] = atom->new_avec(arg[iarg],1,dummy);
-    keywords[nstyles] = new char[strlen(arg[iarg])+1];
+    keywords[nstyles] = utils::strdup(arg[iarg]);
-    strcpy(keywords[nstyles],arg[iarg]);
     jarg = iarg + 1;
     while (jarg < narg && !known_style(arg[jarg])) jarg++;
     styles[nstyles]->process_args(jarg-iarg-1,&arg[iarg+1]);
@@ -601,13 +600,10 @@ void AtomVecHybrid::build_styles()
 
   allstyles = new char*[nallstyles];
 
-  int n;
   nallstyles = 0;
 #define ATOM_CLASS
-#define AtomStyle(key,Class)                \
+#define AtomStyle(key,Class)                   \
-  n = strlen(#key) + 1;                     \
+  allstyles[nallstyles] = utils::strdup(#key); \
-  allstyles[nallstyles] = new char[n];      \
-  strcpy(allstyles[nallstyles],#key);       \
   nallstyles++;
 #include "style_atom.h"  // IWYU pragma: keep
 #undef AtomStyle

src/body.cpp

@@ -20,9 +20,7 @@ using namespace LAMMPS_NS;
 
 Body::Body(LAMMPS *lmp, int /*narg*/, char **arg) : Pointers(lmp)
 {
-  int n = strlen(arg[0]) + 1;
+  style = utils::strdup(arg[0]);
-  style = new char[n];
-  strcpy(style,arg[0]);
   icp = nullptr;
   dcp = nullptr;
 }