Merge branch 'develop' into quick-regression

.github/workflows/check-vla.yml

@@ -0,0 +1,89 @@
+# GitHub action to build LAMMPS on Linux with gcc and -Werror=vla
+name: "Check for Variable Length Arrays"
+
+on:
+  push:
+    branches:
+      - develop
+  pull_request:
+    branches:
+      - develop
+
+  workflow_dispatch:
+
+jobs:
+  build:
+    name: Build with -Werror=vla
+    if: ${{ github.repository == 'lammps/lammps' }}
+    runs-on: ubuntu-latest
+    env:
+      CCACHE_DIR: ${{ github.workspace }}/.ccache
+
+    steps:
+    - name: Checkout repository
+      uses: actions/checkout@v4
+      with:
+        fetch-depth: 2
+
+    - name: Install extra packages
+      run: |
+        sudo apt-get install -y ccache \
+                                libeigen3-dev \
+                                libgsl-dev \
+                                libcurl4-openssl-dev \
+                                mold \
+                                mpi-default-bin \
+                                mpi-default-dev \
+                                ninja-build \
+                                python3-dev
+
+    - name: Create Build Environment
+      run: mkdir build
+
+    - name: Set up ccache
+      uses: actions/cache@v4
+      with:
+        path: ${{ env.CCACHE_DIR }}
+        key: linux-vla-ccache-${{ github.sha }}
+        restore-keys: linux-vla-ccache-
+
+    - name: Building LAMMPS via CMake
+      shell: bash
+      run: |
+        ccache -z
+        python3 -m venv linuxenv
+        source linuxenv/bin/activate
+        python3 -m pip install numpy
+        python3 -m pip install pyyaml
+        cmake -S cmake -B build \
+              -C cmake/presets/most.cmake \
+              -D CMAKE_CXX_COMPILER=g++ \
+              -D CMAKE_C_COMPILER=gcc \
+              -D CMAKE_CXX_COMPILER_LAUNCHER=ccache \
+              -D CMAKE_C_COMPILER_LAUNCHER=ccache \
+              -D CMAKE_BUILD_TYPE=Debug \
+              -D CMAKE_CXX_FLAGS_DEBUG="-Og -g -Werror=vla" \
+              -D DOWNLOAD_POTENTIALS=off \
+              -D BUILD_MPI=on \
+              -D BUILD_SHARED_LIBS=off \
+              -D BUILD_TOOLS=off \
+              -D ENABLE_TESTING=off \
+              -D MLIAP_ENABLE_ACE=on \
+              -D MLIAP_ENABLE_PYTHON=off \
+              -D PKG_AWPMD=on \
+              -D PKG_GPU=on \
+              -D GPU_API=opencl \
+              -D PKG_LATBOLTZ=on \
+              -D PKG_MDI=on \
+              -D PKG_MANIFOLD=on \
+              -D PKG_ML-PACE=on \
+              -D PKG_ML-RANN=off \
+              -D PKG_MOLFILE=on \
+              -D PKG_RHEO=on \
+              -D PKG_PTM=on \
+              -D PKG_PYTHON=on \
+              -D PKG_QTB=on \
+              -D PKG_SMTBQ=on \
+              -G Ninja
+        cmake --build build
+        ccache -s

.github/workflows/compile-msvc.yml

@@ -1,5 +1,5 @@
-# GitHub action to build LAMMPS on Windows with Visual C++ and run unit tests
+# GitHub action to test LAMMPS on Windows with Visual C++
-name: "Native Windows Unit Tests"
+name: "Windows Unit Tests"
 
 on:
   push:
@@ -16,6 +16,8 @@ jobs:
     name: Windows Compilation Test
     if: ${{ github.repository == 'lammps/lammps' }}
     runs-on: windows-latest
+    env:
+      CCACHE_DIR: ${{ github.workspace }}/.ccache
 
     steps:
     - name: Checkout repository
@@ -23,36 +25,41 @@ jobs:
       with:
         fetch-depth: 2
 
+    - name: Enable MSVC++
+      uses: lammps/setup-msvc-dev@v3
+      with:
+        arch: x64
+
+    - name: Install Ccache
+      run: |
+        choco install ccache ninja -y
+
+    - name: Set up ccache
+      uses: actions/cache@v4
+      with:
+        path: ${{ env.CCACHE_DIR }}
+        key: win-unit-ccache-${{ github.sha }}
+        restore-keys: win-unit-ccache-
+
     - name: Select Python version
       uses: actions/setup-python@v5
       with:
         python-version: '3.11'
 
     - name: Building LAMMPS via CMake
-      shell: bash
       run: |
+        ccache -z
         python3 -m pip install numpy
         python3 -m pip install pyyaml
-        nuget install MSMPIsdk
+        cmake -C cmake\presets\windows.cmake -D CMAKE_CXX_COMPILER=cl -D CMAKE_CXX_COMPILER_LAUNCHER=ccache -D CMAKE_C_COMPILER=cl -D CMAKE_C_COMPILER_LAUNCHER=ccache -D CMAKE_Fortran_COMPILER="" -D DOWNLOAD_POTENTIALS=off -D PKG_PYTHON=on -D WITH_PNG=off -D WITH_JPEG=off -S cmake -B build -D BUILD_SHARED_LIBS=on -D ENABLE_TESTING=on -D CMAKE_BUILD_TYPE=Release -G Ninja
-        nuget install MSMPIDIST
+        cmake --build build
-        cmake -C cmake/presets/windows.cmake \
+        ccache -s
-              -D DOWNLOAD_POTENTIALS=off \
-              -D PKG_PYTHON=on \
-              -D WITH_PNG=off \
-              -D WITH_JPEG=off \
-              -S cmake -B build \
-              -D BUILD_SHARED_LIBS=on \
-              -D LAMMPS_EXCEPTIONS=on \
-              -D ENABLE_TESTING=on
-        cmake --build build --config Release --parallel 2
 
     - name: Run LAMMPS executable
-      shell: bash
       run: |
-        ./build/Release/lmp.exe -h
+        build\lmp.exe -h
-        ./build/Release/lmp.exe -in bench/in.lj
+        build\lmp.exe -in bench\in.lj
 
     - name: Run Unit Tests
       working-directory: build
-      shell: bash
+      run: ctest -V -E FixTimestep:python_move_nve
-      run: ctest -V -C Release -E FixTimestep:python_move_nve

.github/workflows/unittest-linux.yml

@@ -0,0 +1,82 @@
+# GitHub action to build LAMMPS on Linux and run standard unit tests
+name: "Unittest for Linux /w LAMMPS_BIGBIG"
+
+on:
+  push:
+    branches:
+      - develop
+  pull_request:
+    branches:
+      - develop
+
+  workflow_dispatch:
+
+jobs:
+  build:
+    name: Linux Unit Test
+    if: ${{ github.repository == 'lammps/lammps' }}
+    runs-on: ubuntu-latest
+    env:
+      CCACHE_DIR: ${{ github.workspace }}/.ccache
+
+    steps:
+    - name: Checkout repository
+      uses: actions/checkout@v4
+      with:
+        fetch-depth: 2
+
+    - name: Install extra packages
+      run: |
+        sudo apt-get install -y ccache \
+                                libeigen3-dev \
+                                libgsl-dev \
+                                libcurl4-openssl-dev \
+                                mold \
+                                ninja-build \
+                                python3-dev
+
+    - name: Create Build Environment
+      run: mkdir build
+
+    - name: Set up ccache
+      uses: actions/cache@v4
+      with:
+        path: ${{ env.CCACHE_DIR }}
+        key: linux-unit-ccache-${{ github.sha }}
+        restore-keys: linux-unit-ccache-
+
+    - name: Building LAMMPS via CMake
+      shell: bash
+      run: |
+        ccache -z
+        python3 -m venv linuxenv
+        source linuxenv/bin/activate
+        python3 -m pip install numpy
+        python3 -m pip install pyyaml
+        cmake -S cmake -B build \
+              -C cmake/presets/gcc.cmake \
+              -C cmake/presets/most.cmake \
+              -D CMAKE_CXX_COMPILER_LAUNCHER=ccache \
+              -D CMAKE_C_COMPILER_LAUNCHER=ccache \
+              -D BUILD_SHARED_LIBS=on \
+              -D LAMMPS_SIZES=bigbig \
+              -D DOWNLOAD_POTENTIALS=off \
+              -D ENABLE_TESTING=on \
+              -D MLIAP_ENABLE_ACE=on \
+              -D MLIAP_ENABLE_PYTHON=off \
+              -D PKG_MANIFOLD=on \
+              -D PKG_ML-PACE=on \
+              -D PKG_ML-RANN=on \
+              -D PKG_RHEO=on \
+              -D PKG_PTM=on \
+              -D PKG_PYTHON=on \
+              -D PKG_QTB=on \
+              -D PKG_SMTBQ=on \
+              -G Ninja
+        cmake --build build
+        ccache -s
+
+    - name: Run Tests
+      working-directory: build
+      shell: bash
+      run: ctest -V

cmake/CMakeLists.txt

@@ -474,13 +474,13 @@ if(BUILD_OMP)
   if(CMAKE_VERSION VERSION_LESS 3.28)
     get_filename_component(_exe "${CMAKE_CXX_COMPILER}" NAME)
     if((CMAKE_CXX_COMPILER_ID STREQUAL "Clang") AND (_exe STREQUAL "crayCC"))
-      set(CMAKE_SHARED_LINKER_FLAGS_${BTYPE} "${CMAKE_SHARED_LINKER_FLAGS_${BTYPE} -fopenmp")
+      set(CMAKE_SHARED_LINKER_FLAGS_${BTYPE} "${CMAKE_SHARED_LINKER_FLAGS_${BTYPE}} -fopenmp")
-      set(CMAKE_STATIC_LINKER_FLAGS_${BTYPE} "${CMAKE_STATIC_LINKER_FLAGS_${BTYPE} -fopenmp")
+      set(CMAKE_STATIC_LINKER_FLAGS_${BTYPE} "${CMAKE_STATIC_LINKER_FLAGS_${BTYPE}} -fopenmp")
     endif()
   else()
     if(CMAKE_CXX_COMPILER_ID STREQUAL "CrayClang")
-      set(CMAKE_SHARED_LINKER_FLAGS_${BTYPE} "${CMAKE_SHARED_LINKER_FLAGS_${BTYPE} -fopenmp")
+      set(CMAKE_SHARED_LINKER_FLAGS_${BTYPE} "${CMAKE_SHARED_LINKER_FLAGS_${BTYPE}} -fopenmp")
-      set(CMAKE_STATIC_LINKER_FLAGS_${BTYPE} "${CMAKE_STATIC_LINKER_FLAGS_${BTYPE} -fopenmp")
+      set(CMAKE_STATIC_LINKER_FLAGS_${BTYPE} "${CMAKE_STATIC_LINKER_FLAGS_${BTYPE}} -fopenmp")
     endif()
   endif()
 endif()
@@ -973,6 +973,9 @@ message(STATUS "<<< Compilers and Flags: >>>
       C++ Standard:  ${CMAKE_CXX_STANDARD}
       C++ Flags:    ${CMAKE_CXX_FLAGS} ${CMAKE_CXX_FLAGS_${BTYPE}}
       Defines:       ${DEFINES}")
+if(CMAKE_CXX_COMPILER_LAUNCHER)
+  message(STATUS "      Launcher:      ${CMAKE_CXX_COMPILER_LAUNCHER}")
+endif()
 get_target_property(OPTIONS lammps COMPILE_OPTIONS)
 if(OPTIONS)
   message("      Options:       ${OPTIONS}")
@@ -991,6 +994,9 @@ if(_index GREATER -1)
       Type:          ${CMAKE_C_COMPILER_ID}
       Version:       ${CMAKE_C_COMPILER_VERSION}
       C Flags:      ${CMAKE_C_FLAGS} ${CMAKE_C_FLAGS_${BTYPE}}")
+  if(CMAKE_C_COMPILER_LAUNCHER)
+    message(STATUS "      Launcher:      ${CMAKE_C_COMPILER_LAUNCHER}")
+  endif()
 endif()
 message(STATUS "<<< Linker flags: >>>")
 message(STATUS "Executable name:  ${LAMMPS_BINARY}")

doc/src/Build_development.rst

@@ -630,11 +630,35 @@ The following target are available for both, GNU make and CMake:
 GitHub command line interface
 -----------------------------
 
-GitHub is developing a `tool for the command line
+GitHub has developed a `command line tool <https://cli.github.com>`_
-<https://cli.github.com>`_ that interacts with the GitHub website via a
+to interact with the GitHub website via a command called ``gh``.
-command called ``gh``.  This can be extremely convenient when working
+This is extremely convenient when working with a Git repository hosted
-with a Git repository hosted on GitHub (like LAMMPS).  It is thus highly
+on GitHub (like LAMMPS).  It is thus highly recommended to install it
-recommended to install it when doing LAMMPS development.
+when doing LAMMPS development.  To use ``gh`` you must be within a git
+checkout of a repository and you must obtain an authentication token
+to connect your checkout with a GitHub user.  This is done with the
+command: ``gh auth login`` where you then have to follow the prompts.
+Here are some examples:
 
-The capabilities of the ``gh`` command is continually expanding, so
+.. list-table::
-please see the documentation at https://cli.github.com/manual/
+   :header-rows: 1
+   :widths: 34 66
+
+   * - Command
+     - Description
+   * - ``gh pr list``
+     - List currently open pull requests
+   * - ``gh pr checks 404``
+     - Shows the status of all checks for pull request #404
+   * - ``gh pr view 404``
+     - Shows the description and recent comments for pull request #404
+   * - ``gh co 404``
+     - Check out the branch from pull request #404; set up for pushing changes
+   * - ``gh issue list``
+     - List currently open issues
+   * - ``gh issue view 430 --comments``
+     - Shows the description and all comments for issue #430
+
+The capabilities of the ``gh`` command are continually expanding, so
+for more details please see the documentation at https://cli.github.com/manual/
+or use ``gh --help`` or ``gh <command> --help`` for embedded help.

doc/src/Library_atoms.rst

@@ -4,6 +4,7 @@ Per-atom properties
 This section documents the following functions:
 
 - :cpp:func:`lammps_extract_atom_datatype`
+- :cpp:func:`lammps_extract_atom_size`
 - :cpp:func:`lammps_extract_atom`
 
 -----------------------
@@ -13,6 +14,11 @@ This section documents the following functions:
 
 -----------------------
 
+.. doxygenfunction:: lammps_extract_atom_size
+   :project: progguide
+
+-----------------------
+
 .. doxygenfunction:: lammps_extract_atom
    :project: progguide
 

doc/src/pair_hbond_dreiding.rst

@@ -18,28 +18,27 @@ Syntax
 
 .. code-block:: LAMMPS
 
-   pair_style style N inner_distance_cutoff outer_distance_cutoff angle_cutof
+   pair_style style N inner_distance_cutoff outer_distance_cutoff angle_cutoff
 
 * style = *hbond/dreiding/lj* or *hbond/dreiding/morse*
-* n = cosine angle periodicity
+* N = power of cosine of angle theta (integer)
 * inner_distance_cutoff = global inner cutoff for Donor-Acceptor interactions (distance units)
 * outer_distance_cutoff = global cutoff for Donor-Acceptor interactions (distance units)
-* angle_cutoff = global angle cutoff for Acceptor-Hydrogen-Donor
+* angle_cutoff = global angle cutoff for Acceptor-Hydrogen-Donor interactions (degrees)
-* interactions (degrees)
 
 Examples
 """"""""
 
 .. code-block:: LAMMPS
 
-   pair_style hybrid/overlay lj/cut 10.0 hbond/dreiding/lj 4 9.0 11.0 90
+   pair_style hybrid/overlay lj/cut 10.0 hbond/dreiding/lj 4 9.0 11.0 90.0
    pair_coeff 1 2 hbond/dreiding/lj 3 i 9.5 2.75 4 9.0 11.0 90.0
 
-   pair_style hybrid/overlay lj/cut 10.0 hbond/dreiding/morse 2 9.0 11.0 90
+   pair_style hybrid/overlay lj/cut 10.0 hbond/dreiding/morse 2 9.0 11.0 90.0
-   pair_coeff 1 2 hbond/dreiding/morse 3 i 3.88 1.7241379 2.9 2 9 11 90
+   pair_coeff 1 2 hbond/dreiding/morse 3 i 3.88 1.7241379 2.9 2 9.0 11.0 90.0
 
    labelmap atom 1 C 2 O 3 H
-   pair_coeff C O hbond/dreiding/morse H i 3.88 1.7241379 2.9 2 9 11 90
+   pair_coeff C O hbond/dreiding/morse H i 3.88 1.7241379 2.9 2 9.0 11.0 90.0
 
 Description
 """""""""""
@@ -65,7 +64,8 @@ force field, given by:
 
 where :math:`r_{\rm in}` is the inner spline distance cutoff,
 :math:`r_{\rm out}` is the outer distance cutoff, :math:`\theta_c` is
-the angle cutoff, and *n* is the cosine periodicity.
+the angle cutoff, and :math:`n` is the power of the cosine of the angle
+:math:`\theta`.
 
 Here, *r* is the radial distance between the donor (D) and acceptor
 (A) atoms and :math:`\theta` is the bond angle between the acceptor, the
@@ -217,7 +217,8 @@ These pair styles do not support the :doc:`pair_modify <pair_modify>`
 tail option for adding long-range tail corrections to energy and
 pressure.
 
-These pair styles do not write their information to :doc:`binary restart files <restart>`, so pair_style and pair_coeff commands need to be
+These pair styles do not write their information to :doc:`binary restart
+files <restart>`, so pair_style and pair_coeff commands need to be
 re-specified in an input script that reads a restart file.
 
 These pair styles can only be used via the *pair* keyword of the

examples/COUPLE/plugin/liblammpsplugin.c

@@ -105,6 +105,7 @@ liblammpsplugin_t *liblammpsplugin_load(const char *lib)
   ADDSYM(map_atom);
 
   ADDSYM(extract_atom_datatype);
+  ADDSYM(extract_atom_size);
   ADDSYM(extract_atom);
 
   ADDSYM(extract_compute);

examples/COUPLE/plugin/liblammpsplugin.h

@@ -151,6 +151,7 @@ struct _liblammpsplugin {
   int (*map_atom)(void *, const void *);
 
   int (*extract_atom_datatype)(void *, const char *);
+  int (*extract_atom_size)(void *, const char *, int);
   void *(*extract_atom)(void *, const char *);
 
   void *(*extract_compute)(void *, const char *, int, int);

fortran/lammps.f90

@@ -542,6 +542,14 @@ MODULE LIBLAMMPS
       INTEGER(c_int) :: lammps_extract_atom_datatype
     END FUNCTION lammps_extract_atom_datatype
 
+    FUNCTION lammps_extract_atom_size(handle, name, dtype) BIND(C)
+      IMPORT :: c_ptr, c_int
+      IMPLICIT NONE
+      TYPE(c_ptr), INTENT(IN), VALUE :: handle, name
+      INTEGER(c_int), INTENT(IN), VALUE :: dtype
+      INTEGER(c_int) :: lammps_extract_atom_size
+    END FUNCTION lammps_extract_atom_size
+
     FUNCTION lammps_extract_atom(handle, name) BIND(C)
       IMPORT :: c_ptr
       IMPLICIT NONE
@@ -1461,43 +1469,35 @@ CONTAINS
     ntypes = lmp_extract_setting(self, 'ntypes')
     Cname = f2c_string(name)
     datatype = lammps_extract_atom_datatype(self%handle, Cname)
+    nrows = lammps_extract_atom_size(self%handle, Cname, LMP_SIZE_ROWS)
+    ncols = lammps_extract_atom_size(self%handle, Cname, LMP_SIZE_COLS)
     Cptr = lammps_extract_atom(self%handle, Cname)
     CALL lammps_free(Cname)
 
-    SELECT CASE (name)
-      CASE ('mass')
-        ncols = ntypes + 1
-        nrows = 1
-      CASE ('x','v','f','mu','omega','torque','angmom')
-        ncols = nmax
-        nrows = 3
-      CASE DEFAULT
-        ncols = nmax
-        nrows = 1
-    END SELECT
-
     peratom_data%lammps_instance => self
     SELECT CASE (datatype)
       CASE (LAMMPS_INT)
         peratom_data%datatype = DATA_INT_1D
-        CALL C_F_POINTER(Cptr, peratom_data%i32_vec, [ncols])
+        CALL C_F_POINTER(Cptr, peratom_data%i32_vec, [nrows])
       CASE (LAMMPS_INT64)
         peratom_data%datatype = DATA_INT64_1D
-        CALL C_F_POINTER(Cptr, peratom_data%i64_vec, [ncols])
+        CALL C_F_POINTER(Cptr, peratom_data%i64_vec, [nrows])
       CASE (LAMMPS_DOUBLE)
         peratom_data%datatype = DATA_DOUBLE_1D
+        ! The mass array is allocated from 0, but only used from 1. We also want to use it from 1.
         IF (name == 'mass') THEN
-          CALL C_F_POINTER(Cptr, dummy, [ncols])
+          CALL C_F_POINTER(Cptr, dummy, [nrows])
           peratom_data%r64_vec(0:) => dummy
         ELSE
-          CALL C_F_POINTER(Cptr, peratom_data%r64_vec, [ncols])
+          CALL C_F_POINTER(Cptr, peratom_data%r64_vec, [nrows])
         END IF
       CASE (LAMMPS_DOUBLE_2D)
         peratom_data%datatype = DATA_DOUBLE_2D
         ! First, we dereference the void** pointer to point to the void*
-        CALL C_F_POINTER(Cptr, Catomptr, [ncols])
+        CALL C_F_POINTER(Cptr, Catomptr, [nrows])
         ! Catomptr(1) now points to the first element of the array
-        CALL C_F_POINTER(Catomptr(1), peratom_data%r64_mat, [nrows,ncols])
+        ! rows and columns are swapped in Fortran
+        CALL C_F_POINTER(Catomptr(1), peratom_data%r64_mat, [ncols,nrows])
       CASE (-1)
         CALL lmp_error(self, LMP_ERROR_ALL + LMP_ERROR_WORLD, &
           'per-atom property ' // name // ' not found in extract_setting')

python/lammps/core.py

@@ -318,6 +318,8 @@ class lammps(object):
     self.lib.lammps_extract_atom.argtypes = [c_void_p, c_char_p]
     self.lib.lammps_extract_atom_datatype.argtypes = [c_void_p, c_char_p]
     self.lib.lammps_extract_atom_datatype.restype = c_int
+    self.lib.lammps_extract_atom_size.argtypes = [c_void_p, c_char_p, c_int]
+    self.lib.lammps_extract_atom_size.restype = c_int
 
     self.lib.lammps_extract_fix.argtypes = [c_void_p, c_char_p, c_int, c_int, c_int, c_int]
 
@@ -1070,31 +1072,59 @@ class lammps(object):
     else: return None
     return self.lib.lammps_extract_atom_datatype(self.lmp, newname)
 
+  # -------------------------------------------------------------------------
+  # extract per-atom info datatype
+
+  def extract_atom_size(self, name, dtype):
+    """Retrieve per-atom property dimensions from LAMMPS
+
+    This is a wrapper around the :cpp:func:`lammps_extract_atom_size`
+    function of the C-library interface. Its documentation includes a
+    list of the supported keywords.
+    This function returns ``None`` if the keyword is not
+    recognized. Otherwise it will return an integer value with the size
+    of the per-atom vector or array.  If *name* corresponds to a per-atom
+    array, the *dtype* keyword must be either LMP_SIZE_ROWS or LMP_SIZE_COLS
+    from the :ref:`type <py_type_constants>` constants defined in the
+    :py:mod:`lammps` module.  The return value is the requested size.
+    If *name* corresponds to a per-atom vector the *dtype* keyword is ignored.
+
+    :param name: name of the property
+    :type name:  string
+    :param type: either LMP_SIZE_ROWS or LMP_SIZE_COLS for arrays, otherwise ignored
+    :type type:  int
+    :return: data type of per-atom property (see :ref:`py_datatype_constants`)
+    :rtype: int
+    """
+    if name: newname = name.encode()
+    else: return None
+    return self.lib.lammps_extract_atom_size(self.lmp, newname, dtype)
+
   # -------------------------------------------------------------------------
   # extract per-atom info
 
   def extract_atom(self, name, dtype=LAMMPS_AUTODETECT):
     """Retrieve per-atom properties from LAMMPS
 
-    This is a wrapper around the :cpp:func:`lammps_extract_atom`
+    This is a wrapper around the :cpp:func:`lammps_extract_atom` function of the
-    function of the C-library interface. Its documentation includes a
+    C-library interface. Its documentation includes a list of the supported
-    list of the supported keywords and their data types.
+    keywords and their data types.  Since Python needs to know the data type to
-    Since Python needs to know the data type to be able to interpret
+    be able to interpret the result, by default, this function will try to
-    the result, by default, this function will try to auto-detect the data type
+    auto-detect the data type by asking the library. You can also force a
-    by asking the library. You can also force a specific data type by setting ``dtype``
+    specific data type by setting ``dtype`` to one of the :ref:`data type
-    to one of the :ref:`data type <py_datatype_constants>` constants defined in the
+    <py_datatype_constants>` constants defined in the :py:mod:`lammps` module.
-    :py:mod:`lammps` module.
+    This function returns ``None`` if either the keyword is not recognized, or
-    This function returns ``None`` if either the keyword is not
+    an invalid data type constant is used.
-    recognized, or an invalid data type constant is used.
 
     .. note::
 
-       While the returned arrays of per-atom data are dimensioned
+       While the returned vectors or arrays of per-atom data are dimensioned for
-       for the range [0:nmax] - as is the underlying storage -
+       the range [0:nmax] - as is the underlying storage - the data is usually
-       the data is usually only valid for the range of [0:nlocal],
+       only valid for the range of [0:nlocal], unless the property of interest
-       unless the property of interest is also updated for ghost
+       is also updated for ghost atoms.  In some cases, this depends on a LAMMPS
-       atoms.  In some cases, this depends on a LAMMPS setting, see
+       setting, see for example :doc:`comm_modify vel yes <comm_modify>`.
-       for example :doc:`comm_modify vel yes <comm_modify>`.
+       The actual size can be determined by calling
+       py:meth:`extract_atom_size() <lammps.lammps.extract_atom_size>`.
 
     :param name: name of the property
     :type name:  string
@@ -1105,6 +1135,7 @@ class lammps(object):
             ctypes.POINTER(ctypes.c_int64), ctypes.POINTER(ctypes.POINTER(ctypes.c_int64)),
             ctypes.POINTER(ctypes.c_double), ctypes.POINTER(ctypes.POINTER(ctypes.c_double)),
             or NoneType
+
     """
     if dtype == LAMMPS_AUTODETECT:
       dtype = self.extract_atom_datatype(name)
@@ -2522,3 +2553,7 @@ class lammps(object):
     newcomputeid = computeid.encode()
     idx = self.lib.lammps_find_compute_neighlist(self.lmp, newcomputeid, reqid)
     return idx
+
+# Local Variables:
+# fill-column: 80
+# End:

python/lammps/numpy_wrapper.py

@@ -54,7 +54,8 @@ class numpy_wrapper:
 
   # -------------------------------------------------------------------------
 
-  def extract_atom(self, name, dtype=LAMMPS_AUTODETECT, nelem=LAMMPS_AUTODETECT, dim=LAMMPS_AUTODETECT):
+  def extract_atom(self, name, dtype=LAMMPS_AUTODETECT, nelem=LAMMPS_AUTODETECT,
+                   dim=LAMMPS_AUTODETECT):
     """Retrieve per-atom properties from LAMMPS as NumPy arrays
 
     This is a wrapper around the :py:meth:`lammps.extract_atom()` method.
@@ -63,16 +64,16 @@ class numpy_wrapper:
 
     .. note::
 
-       The returned arrays of per-atom data are by default dimensioned
+       The returned vectors or arrays of per-atom data are dimensioned
-       for the range [0:nlocal] since that data is *always* valid.  The
+       according to the return value of :py:meth:`lammps.extract_atom_size()`.
-       underlying storage for the data, however, is typically allocated
+       Except for the "mass" property, the underlying storage will always be
-       for the range of [0:nmax]. Whether there is valid data in the range
+       dimensioned for the range [0:nmax].  The actual usable data may be
-       [nlocal:nlocal+nghost] depends on whether the property of interest
+       only in the range [0:nlocal] or [0:nlocal][0:dim].  Whether there is
-       is also updated for ghost atoms.  This is not often the case.  In
+       valid data in the range [nlocal:nlocal+nghost] or [nlocal:local+nghost][0:dim]
-       some cases, it depends on a LAMMPS setting, see for example
+       depends on whether the property of interest is also updated for ghost atoms.
-       :doc:`comm_modify vel yes <comm_modify>`.  By using the optional
+       Also the value of *dim* depends on the value of *name*.  By using the optional
-       *nelem* parameter the size of the returned NumPy can be overridden.
+       *nelem* and *dim* parameters the dimensions of the returned NumPy array can
-       There is no check whether the number of elements chosen is valid.
+       be overridden.  There is no check whether the number of elements chosen is valid.
 
     :param name: name of the property
     :type name:  string
@@ -89,21 +90,10 @@ class numpy_wrapper:
       dtype = self.lmp.extract_atom_datatype(name)
 
     if nelem == LAMMPS_AUTODETECT:
-      if name == "mass":
+      nelem = self.lmp.extract_atom_size(name, LMP_SIZE_ROWS)
-        nelem = self.lmp.extract_global("ntypes") + 1
-      else:
-        nelem = self.lmp.extract_global("nlocal")
     if dim == LAMMPS_AUTODETECT:
       if dtype in (LAMMPS_INT_2D, LAMMPS_DOUBLE_2D, LAMMPS_INT64_2D):
-        # TODO add other fields
+        dim = self.lmp.extract_atom_size(name, LMP_SIZE_COLS)
-        if name in ("x", "v", "f", "x0","omega", "angmom", "torque", "csforce", "vforce", "vest"):
-          dim = 3
-        elif name == "smd_data_9":
-          dim = 9
-        elif name == "smd_stress":
-          dim = 6
-        else:
-          dim = 2
       else:
         dim = 1
 
@@ -119,37 +109,6 @@ class numpy_wrapper:
 
   # -------------------------------------------------------------------------
 
-  def extract_atom_iarray(self, name, nelem, dim=1):
-    warnings.warn("deprecated, use extract_atom instead", DeprecationWarning)
-
-    if name in ['id', 'molecule']:
-      c_int_type = self.lmp.c_tagint
-    elif name in ['image']:
-      c_int_type = self.lmp.c_imageint
-    else:
-      c_int_type = c_int
-
-    if dim == 1:
-      raw_ptr = self.lmp.extract_atom(name, LAMMPS_INT)
-    else:
-      raw_ptr = self.lmp.extract_atom(name, LAMMPS_INT_2D)
-
-    return self.iarray(c_int_type, raw_ptr, nelem, dim)
-
-  # -------------------------------------------------------------------------
-
-  def extract_atom_darray(self, name, nelem, dim=1):
-    warnings.warn("deprecated, use extract_atom instead", DeprecationWarning)
-
-    if dim == 1:
-      raw_ptr = self.lmp.extract_atom(name, LAMMPS_DOUBLE)
-    else:
-      raw_ptr = self.lmp.extract_atom(name, LAMMPS_DOUBLE_2D)
-
-    return self.darray(raw_ptr, nelem, dim)
-
-  # -------------------------------------------------------------------------
-
   def extract_compute(self, cid, cstyle, ctype):
     """Retrieve data from a LAMMPS compute
 

src/CG-DNA/atom_vec_oxdna.cpp

@@ -37,6 +37,7 @@ AtomVecOxdna::AtomVecOxdna(LAMMPS *lmp) : AtomVec(lmp)
   fields_grow = {"id5p"};
   fields_copy = {"id5p"};
   fields_border = {"id5p"};
+  fields_border_vel = {"id5p"};
   fields_exchange = {"id5p"};
   fields_restart = {"id5p"};
   fields_data_atom = {"id", "type", "x"};

src/DIFFRACTION/compute_saed.cpp

@@ -401,7 +401,7 @@ void ComputeSAED::compute_vector()
 
   // Setting up OMP
 #if defined(_OPENMP)
-  if (me == 0 && echo) utils::logmesg(lmp," using {}OMP threads\n",comm->nthreads);
+  if (me == 0 && echo) utils::logmesg(lmp," using {} OMP thread(s)\n",comm->nthreads);
 #endif
 
   if (me == 0 && echo) utils::logmesg(lmp,"\n");
@@ -478,7 +478,7 @@ void ComputeSAED::compute_vector()
         }
       }
     } // End of pragma omp for region
-    delete [] f;
+    delete[] f;
   }
 
   auto scratch = new double[2*nRows];
@@ -499,10 +499,10 @@ void ComputeSAED::compute_vector()
     utils::logmesg(lmp," 100% \nTime elapsed during compute_saed = {:.2f} sec "
                    "using {:.2f} Mbytes/processor\n-----\n", t2-t0, bytes/1024.0/1024.0);
 
-  delete [] xlocal;
+  delete[] xlocal;
-  delete [] typelocal;
+  delete[] typelocal;
-  delete [] scratch;
+  delete[] scratch;
-  delete [] Fvec;
+  delete[] Fvec;
 }
 
 /* ----------------------------------------------------------------------

src/DIFFRACTION/compute_xrd.cpp

@@ -332,7 +332,7 @@ void ComputeXRD::compute_array()
 
 // Setting up OMP
 #if defined(_OPENMP)
-  if ((me == 0) && echo) utils::logmesg(lmp," using {} OMP threads\n",comm->nthreads);
+  if ((me == 0) && echo) utils::logmesg(lmp," using {} OMP thread(s)\n",comm->nthreads);
 #endif
 
   if ((me == 0) && echo) {
@@ -482,7 +482,7 @@ void ComputeXRD::compute_array()
         }
       } // End of pragma omp for region
     } // End of if LP=1 check
-    delete [] f;
+    delete[] f;
   } // End of pragma omp parallel region
 
   auto scratch = new double[2*size_array_rows];
@@ -503,10 +503,10 @@ void ComputeXRD::compute_array()
     utils::logmesg(lmp," 100% \nTime elapsed during compute_xrd = {:.2f} sec "
                    "using {:.2f} Mbytes/processor\n-----\n", t2-t0, bytes/1024.0/1024.0);
 
-  delete [] scratch;
+  delete[] scratch;
-  delete [] Fvec;
+  delete[] Fvec;
-  delete [] xlocal;
+  delete[] xlocal;
-  delete [] typelocal;
+  delete[] typelocal;
 }
 
 /* ----------------------------------------------------------------------

src/DIFFRACTION/fix_saed_vtk.cpp

@@ -114,6 +114,7 @@ FixSAEDVTK::FixSAEDVTK(LAMMPS *lmp, int narg, char **arg) :
   memory->create(vector_total,nrows,"saed/vtk:vector_total");
 
   vector_flag = 1;
+  extvector = 0;
   size_vector = nrows;
 
   if (nOutput == 0) {
@@ -248,8 +249,8 @@ FixSAEDVTK::FixSAEDVTK(LAMMPS *lmp, int narg, char **arg) :
 
 FixSAEDVTK::~FixSAEDVTK()
 {
-  delete [] filename;
+  delete[] filename;
-  delete [] ids;
+  delete[] ids;
   memory->destroy(vector);
   memory->destroy(vector_total);
   if (fp && comm->me == 0) fclose(fp);

src/ML-QUIP/pair_quip.cpp

@@ -286,7 +286,7 @@ void PairQUIP::coeff(int narg, char **arg)
   // and returns the necessary size of quip_potential. This behavior
   // is invoked by setting n_potential_quip to 0.
   n_quip_potential = 0;
-  quip_potential = new int[0];
+  quip_potential = new int[1];
   quip_lammps_potential_initialise(quip_potential, &n_quip_potential, &cutoff, quip_file,
                                    &n_quip_file, quip_string, &n_quip_string);
   delete[] quip_potential;

src/MOLECULE/pair_hbond_dreiding_lj.cpp

@@ -396,14 +396,14 @@ void PairHbondDreidingLJ::init_style()
   //   and computing forces on A,H which may be on different procs
 
   if (atom->molecular == Atom::ATOMIC)
-    error->all(FLERR,"Pair style hbond/dreiding requires molecular system");
+    error->all(FLERR,"Pair style hbond/dreiding/lj requires molecular system");
   if (atom->tag_enable == 0)
-    error->all(FLERR,"Pair style hbond/dreiding requires atom IDs");
+    error->all(FLERR,"Pair style hbond/dreiding/lj requires atom IDs");
   if (atom->map_style == Atom::MAP_NONE)
-    error->all(FLERR,"Pair style hbond/dreiding requires an atom map, "
+    error->all(FLERR,"Pair style hbond/dreiding/lj requires an atom map, "
                "see atom_modify");
   if (force->newton_pair == 0)
-    error->all(FLERR,"Pair style hbond/dreiding requires newton pair on");
+    error->all(FLERR,"Pair style hbond/dreiding/lj requires newton pair on");
 
   // set donor[M]/acceptor[M] if any atom of type M is a donor/acceptor
 
@@ -419,7 +419,7 @@ void PairHbondDreidingLJ::init_style()
           acceptor[j] = 1;
         }
 
-  if (!anyflag) error->all(FLERR,"No pair hbond/dreiding coefficients set");
+  if (!anyflag) error->all(FLERR,"No pair hbond/dreiding/lj coefficients set");
 
   // set additional param values
   // offset is for LJ only, angle term is not included

src/MOLECULE/pair_hbond_dreiding_morse.cpp

@@ -323,14 +323,14 @@ void PairHbondDreidingMorse::init_style()
   //   and computing forces on A,H which may be on different procs
 
   if (atom->molecular == Atom::ATOMIC)
-    error->all(FLERR,"Pair style hbond/dreiding requires molecular system");
+    error->all(FLERR,"Pair style hbond/dreiding/morse requires molecular system");
   if (atom->tag_enable == 0)
-    error->all(FLERR,"Pair style hbond/dreiding requires atom IDs");
+    error->all(FLERR,"Pair style hbond/dreiding/morse requires atom IDs");
   if (atom->map_style == Atom::MAP_NONE)
-    error->all(FLERR,"Pair style hbond/dreiding requires an atom map, "
+    error->all(FLERR,"Pair style hbond/dreiding/morse requires an atom map, "
                "see atom_modify");
   if (force->newton_pair == 0)
-    error->all(FLERR,"Pair style hbond/dreiding requires newton pair on");
+    error->all(FLERR,"Pair style hbond/dreiding/morse requires newton pair on");
 
   // set donor[M]/acceptor[M] if any atom of type M is a donor/acceptor
 
@@ -346,7 +346,7 @@ void PairHbondDreidingMorse::init_style()
           acceptor[j] = 1;
         }
 
-  if (!anyflag) error->all(FLERR,"No pair hbond/dreiding coefficients set");
+  if (!anyflag) error->all(FLERR,"No pair hbond/dreiding/morse coefficients set");
 
   // set additional param values
   // offset is for Morse only, angle term is not included

src/atom.cpp

@@ -1911,7 +1911,11 @@ void Atom::allocate_type_arrays()
   if (avec->mass_type == AtomVec::PER_TYPE) {
     mass = new double[ntypes+1];
     mass_setflag = new int[ntypes+1];
-    for (int itype = 1; itype <= ntypes; itype++) mass_setflag[itype] = 0;
+    // start loop from 0 to avoid uninitialized access when operating on the whole array
+    for (int itype = 0; itype <= ntypes; itype++) {
+      mass_setflag[itype] = 0;
+      mass[itype] = 0.0;
+    }
   }
 }
 
@@ -2739,20 +2743,21 @@ Classes rarely need to check on ghost communication and so `find_custom`
 is typically preferred to this function. See :doc:`pair amoeba <pair_amoeba>`
 for an example where checking ghost communication is necessary.
 \endverbatim
- * \param name Name of the property (w/o a "i_" or "d_" or "i2_" or "d2_" prefix)
+ * \param name    Name of the property (w/o a "i_" or "d_" or "i2_" or "d2_" prefix)
- * \param &flag Returns data type of property: 0 for int, 1 for double
+ * \param &flag   Returns data type of property: 0 for int, 1 for double
- * \param &cols Returns number of values: 0 for a single value, 1 or more for a vector of values
+ * \param &cols   Returns number of values: 0 for a single value, 1 or more for a vector of values
- * \param &ghost Returns whether property is communicated to ghost atoms: 0 for no, 1 for yes
+ * \param &ghost  Returns whether property is communicated to ghost atoms: 0 for no, 1 for yes
  * \return index of property in the respective list of properties
  */
 int Atom::find_custom_ghost(const char *name, int &flag, int &cols, int &ghost)
 {
   int i = find_custom(name, flag, cols);
+  ghost = 0;
   if (i == -1) return i;
   if ((flag == 0) && (cols == 0)) ghost = ivghost[i];
   else if ((flag == 1) && (cols == 0)) ghost = dvghost[i];
-  else if ((flag == 0) && (cols == 1)) ghost = iaghost[i];
+  else if ((flag == 0) && (cols > 0)) ghost = iaghost[i];
-  else if ((flag == 1) && (cols == 1)) ghost = daghost[i];
+  else if ((flag == 1) && (cols > 0)) ghost = daghost[i];
   return i;
 }
 
@@ -2999,11 +3004,13 @@ length of the data area, and a short description.
      - N double values defined by fix property/atom array name
 
 *See also*
-   :cpp:func:`lammps_extract_atom`
+   :cpp:func:`lammps_extract_atom`,
+   :cpp:func:`lammps_extract_atom_datatype`,
+   :cpp:func:`lammps_extract_atom_size`
 
 \endverbatim
  *
- * \sa extract_datatype
+ * \sa extract_datatype, extract_size
  *
  * \param  name  string with the keyword of the desired property.
                  Typically the name of the pointer variable returned
@@ -3142,7 +3149,7 @@ void *Atom::extract(const char *name)
 
 \endverbatim
  *
- * \sa extract
+ * \sa extract extract_size
  *
  * \param  name  string with the keyword of the desired property.
  * \return       data type constant for desired property or -1 */
@@ -3177,10 +3184,14 @@ int Atom::extract_datatype(const char *name)
   if (strcmp(name,"temperature") == 0) return LAMMPS_DOUBLE;
   if (strcmp(name,"heatflow") == 0) return LAMMPS_DOUBLE;
 
+  // PERI package (and in part MACHDYN)
+
   if (strcmp(name,"vfrac") == 0) return LAMMPS_DOUBLE;
   if (strcmp(name,"s0") == 0) return LAMMPS_DOUBLE;
   if (strcmp(name,"x0") == 0) return LAMMPS_DOUBLE_2D;
 
+  // AWPMD package (and in part EFF and ELECTRODE)
+
   if (strcmp(name,"espin") == 0) return LAMMPS_INT;
   if (strcmp(name,"spin") == 0) return LAMMPS_INT;   // backwards compatibility
   if (strcmp(name,"eradius") == 0) return LAMMPS_DOUBLE;
@@ -3248,14 +3259,255 @@ int Atom::extract_datatype(const char *name)
     if (!array) index = find_custom(&name[2],flag,cols);
     else index = find_custom(&name[3],flag,cols);
 
+    // consistency checks
     if (index < 0) return -1;
     if (which != flag) return -1;
     if ((!array && cols) || (array && !cols)) return -1;
 
-    if (which == 0) return LAMMPS_INT;
+    if (!which && !array) return LAMMPS_INT;
-    else return LAMMPS_DOUBLE;
+    if (which && !array) return LAMMPS_DOUBLE;
+    if (!which && array) return LAMMPS_INT_2D;
+    if (which && array) return LAMMPS_DOUBLE_2D;
   }
+  return -1;
+}
 
+/** Provide vector or array size info of internal data of the Atom class
+ *
+\verbatim embed:rst
+
+.. versionadded:: TBD
+
+\endverbatim
+ *
+ * \sa extract extract_datatype
+ *
+ * \param  name  string with the keyword of the desired property.
+ * \param  type  either LMP_SIZE_ROWS or LMP_SIZE_COLS for per-atom array or ignored
+ * \return       size of the vector or size of the array for the requested dimension or -1 */
+
+int Atom::extract_size(const char *name, int type)
+{
+  // --------------------------------------------------------------------
+  // 6th customization section: customize by adding new variable name
+
+  const auto datatype = extract_datatype(name);
+  const auto nall = nlocal + nghost;
+  const auto ghost_vel = comm->ghost_velocity;
+
+  if ((datatype == LAMMPS_DOUBLE_2D) || (datatype == LAMMPS_INT_2D)) {
+    if (type == LMP_SIZE_ROWS) {
+      if (strcmp(name,"x") == 0) return nall;
+      if (strcmp(name,"v") == 0) {
+        if (ghost_vel) return nall;
+        else return nlocal;
+      }
+      if (strcmp(name,"f") == 0) return nall;
+      if (strcmp(name,"mu") == 0) return nall;
+      if (strcmp(name,"omega") == 0) {
+        if (ghost_vel) return nall;
+        else return nlocal;
+      }
+      if (strcmp(name,"angmom") == 0) {
+        if (ghost_vel) return nall;
+        else return nlocal;
+      }
+      if (strcmp(name,"torque") == 0) return nlocal;
+      if (strcmp(name,"quat") == 0) {
+        if (ghost_vel) return nall;
+        else return nlocal;
+      }
+
+      // PERI package
+
+      if (strcmp(name,"x0") == 0) return nall;
+
+      // SPIN package
+
+      if (strcmp(name,"sp") == 0) return nall;
+      if (strcmp(name,"fm") == 0) return nlocal;
+      if (strcmp(name,"fm_long") == 0) return nlocal;
+
+      // AWPMD package
+
+      if (strcmp(name,"cs") == 0) {
+        if (ghost_vel) return nall;
+        else return nlocal;
+      }
+      if (strcmp(name,"csforce") == 0) return nlocal;
+      if (strcmp(name,"vforce") == 0) return nlocal;
+
+      // SPH package
+
+      if (strcmp(name,"vest") == 0) return nall;
+
+      // MACHDYN package
+
+      if (strcmp(name, "smd_data_9") == 0) return LAMMPS_DOUBLE_2D;
+      if (strcmp(name, "smd_stress") == 0) return LAMMPS_DOUBLE_2D;
+
+    } else if (type == LMP_SIZE_COLS) {
+
+      if (strcmp(name,"x") == 0) return 3;
+      if (strcmp(name,"v") == 0) return 3;
+      if (strcmp(name,"f") == 0) return 3;
+      if (strcmp(name,"mu") == 0) return 4;
+      if (strcmp(name,"omega") == 0) return 3;
+      if (strcmp(name,"angmom") == 0) return 3;
+      if (strcmp(name,"torque") == 0) return 3;
+      if (strcmp(name,"quat") == 0) return 4;
+
+      // PERI package
+
+      if (strcmp(name,"x0") == 0) return 3;
+
+      // SPIN package
+
+      if (strcmp(name,"sp") == 0) return 4;
+      if (strcmp(name,"fm") == 0) return 3;
+      if (strcmp(name,"fm_long") == 0) return 3;
+
+      // AWPMD package
+
+      if (strcmp(name,"cs") == 0) return 2;
+      if (strcmp(name,"csforce") == 0) return 2;
+      if (strcmp(name,"vforce") == 0) return 3;
+
+      // SPH package
+
+      if (strcmp(name,"vest") == 0) return 3;
+
+      // MACHDYN package
+
+      if (strcmp(name, "smd_data_9") == 0) return 9;
+      if (strcmp(name, "smd_stress") == 0) return 6;
+    }
+
+    // custom arrays
+
+    if (utils::strmatch(name,"^[id]2_")) {
+      int which = 0;
+      if (name[0] == 'd') which = 1;
+
+      int index,flag,cols,ghost;
+      index = find_custom_ghost(&name[3],flag,cols,ghost);
+
+      // consistency checks
+      if (index < 0) return -1;
+      if (which != flag) return -1;
+      if (!cols) return -1;
+
+      if (type == LMP_SIZE_ROWS) {
+        if (ghost) return nall;
+        else return nlocal;
+      } else if (type == LMP_SIZE_COLS) {
+        return cols;
+      }
+    }
+  } else {
+
+    if (strcmp(name,"mass") == 0) return ntypes + 1;
+
+    if (strcmp(name,"id") == 0) return nall;
+    if (strcmp(name,"type") == 0) return nall;
+    if (strcmp(name,"mask") == 0) return nall;
+    if (strcmp(name,"image") == 0) return nlocal;
+    if (strcmp(name,"molecule") == 0) return nall;
+    if (strcmp(name,"q") == 0) return nall;
+    if (strcmp(name,"radius") == 0) return nall;
+    if (strcmp(name,"rmass") == 0) return nall;
+
+    // ASPHERE package
+
+    if (strcmp(name,"ellipsoid") == 0) return nlocal;
+
+    // BODY package
+
+    if (strcmp(name,"line") == 0) return nlocal;
+    if (strcmp(name,"tri") == 0) return nlocal;
+    if (strcmp(name,"body") == 0) return nlocal;
+
+    // PERI package (and in part MACHDYN)
+
+    if (strcmp(name,"vfrac") == 0) return nall;
+    if (strcmp(name,"s0") == 0) return nall;
+
+    // AWPMD package (and in part EFF and ELECTRODE)
+
+    if (strcmp(name,"espin") == 0) return nall;
+    if (strcmp(name,"spin") == 0) return nall;   // backwards compatibility
+    if (strcmp(name,"eradius") == 0) return nall;
+    if (strcmp(name,"ervel") == 0) return nlocal;
+    if (strcmp(name,"erforce") == 0) return nlocal;
+    if (strcmp(name,"ervelforce") == 0) return nlocal;
+    if (strcmp(name,"etag") == 0) return nall;
+
+    // CG-DNA package
+
+    if (strcmp(name,"id5p") == 0) return nall;
+
+    // RHEO package
+
+    if (strcmp(name,"temperature") == 0) return nlocal;
+    if (strcmp(name,"heatflow") == 0) return nlocal;
+    if (strcmp(name,"rheo_status") == 0) return nall;
+    if (strcmp(name,"conductivity") == 0) return nlocal;
+    if (strcmp(name,"pressure") == 0) return nlocal;
+    if (strcmp(name,"viscosity") == 0) return nlocal;
+
+    // SPH package
+
+    if (strcmp(name,"rho") == 0) return nall;
+    if (strcmp(name,"drho") == 0) return nlocal;
+    if (strcmp(name,"esph") == 0) return nall;
+    if (strcmp(name,"desph") == 0) return nlocal;
+    if (strcmp(name,"cv") == 0) return nall;
+
+    // MACHDYN package
+
+    if (strcmp(name, "contact_radius") == 0) return nall;
+    if (strcmp(name, "eff_plastic_strain") == 0) return nlocal;
+    if (strcmp(name, "eff_plastic_strain_rate") == 0) return nlocal;
+    if (strcmp(name, "damage") == 0) return nlocal;
+
+    // DPD-REACT package
+
+    if (strcmp(name,"dpdTheta") == 0) return nall;
+
+    // DPD-MESO package
+
+    if (strcmp(name,"edpd_temp") == 0) return nall;
+
+    // DIELECTRIC package
+
+    if (strcmp(name,"area") == 0) return nall;
+    if (strcmp(name,"ed") == 0) return nall;
+    if (strcmp(name,"em") == 0) return nall;
+    if (strcmp(name,"epsilon") == 0) return nall;
+    if (strcmp(name,"curvature") == 0) return nall;
+    if (strcmp(name,"q_unscaled") == 0) return nall;
+
+    // end of customization section
+    // --------------------------------------------------------------------
+
+    // custom vectors
+
+    if (utils::strmatch(name,"^[id]_")) {
+      int which = 0;
+      if (name[0] == 'd') which = 1;
+
+      int index,flag,cols,ghost;
+      index = find_custom_ghost(&name[2],flag,cols,ghost);
+
+      // consistency checks
+      if (index < 0) return -1;
+      if (which != flag) return -1;
+      if (cols) return -1;
+
+      if (ghost) return nall;
+      else return nlocal;
+    }
+  }
   return -1;
 }
 

src/atom.h

@@ -378,6 +378,7 @@ class Atom : protected Pointers {
 
   void *extract(const char *);
   int extract_datatype(const char *);
+  int extract_size(const char *, int);
 
   inline int *get_map_array() { return map_array; };
   inline int get_map_size() { return map_tag_max + 1; };

src/error.h

@@ -29,14 +29,14 @@ class Error : protected Pointers {
 
   [[noreturn]] void all(const std::string &, int, const std::string &);
   template <typename... Args>
-  void all(const std::string &file, int line, const std::string &format, Args &&...args)
+  [[noreturn]] void all(const std::string &file, int line, const std::string &format, Args &&...args)
   {
     _all(file, line, format, fmt::make_format_args(args...));
   }
 
   [[noreturn]] void one(const std::string &, int, const std::string &);
   template <typename... Args>
-  void one(const std::string &file, int line, const std::string &format, Args &&...args)
+  [[noreturn]] void one(const std::string &file, int line, const std::string &format, Args &&...args)
   {
     _one(file, line, format, fmt::make_format_args(args...));
   }

src/library.cpp

@@ -2087,10 +2087,13 @@ int lammps_map_atom(void *handle, const void *id)
 
 .. versionadded:: 18Sep2020
 
-This function returns an integer that encodes the data type of the per-atom
+This function returns an integer that encodes the data type of the
-property with the specified name. See :cpp:enum:`_LMP_DATATYPE_CONST` for valid
+per-atom property with the specified name. See
-values. Callers of :cpp:func:`lammps_extract_atom` can use this information
+:cpp:enum:`_LMP_DATATYPE_CONST` for valid values. Callers of
-to then decide how to cast the ``void *`` pointer and access the data.
+:cpp:func:`lammps_extract_atom` can use this information to decide how
+to cast the ``void *`` pointer and access the data.  In addition,
+:cpp:func:`lammps_extract_atom_size` can be used to get information
+about the vector or array dimensions.
 
 \endverbatim
  *
@@ -2108,18 +2111,53 @@ int lammps_extract_atom_datatype(void *handle, const char *name)
 
 /* ---------------------------------------------------------------------- */
 
+/** Get dimension info of a LAMMPS per-atom property
+ *
+\verbatim embed:rst
+
+.. versionadded:: TBD
+
+This function returns an integer with the size of the per-atom
+property with the specified name.  This allows to accurately determine
+the size of the per-atom data vectors or arrays.  For per-atom arrays,
+the *type* argument is required to return either the number of rows or the
+number of columns.  It is ignored for per-atom vectors.
+
+Callers of :cpp:func:`lammps_extract_atom` can use this information in
+combination with the result from :cpp:func:`lammps_extract_atom_datatype`
+to decide how to cast the ``void *`` pointer and access the data.
+
+\endverbatim
+ *
+ * \param  handle  pointer to a previously created LAMMPS instance
+ * \param  name    string with the name of the extracted property
+ * \param  type    either LMP_SIZE_ROWS or LMP_SIZE_COLS if *name* refers
+                   to a per-atom array otherwise ignored
+ * \return         integer with the size of the vector or array dimension or -1
+ * */
+
+int lammps_extract_atom_size(void *handle, const char *name, int type)
+{
+  auto lmp = (LAMMPS *) handle;
+  return lmp->atom->extract_size(name, type);
+}
+
+/* ---------------------------------------------------------------------- */
+
 /** Get pointer to a LAMMPS per-atom property.
  *
 \verbatim embed:rst
 
-This function returns a pointer to the location of per-atom properties
+This function returns a pointer to the location of per-atom properties (and
-(and per-atom-type properties in the case of the 'mass' keyword).
+per-atom-type properties in the case of the 'mass' keyword).  Per-atom data is
-Per-atom data is distributed across sub-domains and thus MPI ranks.  The
+distributed across sub-domains and thus MPI ranks.  The returned pointer is cast
-returned pointer is cast to ``void *`` and needs to be cast to a pointer
+to ``void *`` and needs to be cast to a pointer of data type that the entity
-of data type that the entity represents.
+represents.  You can use the functions :cpp:func:`lammps_extract_atom_datatype`
+and :cpp:func:`lammps_extract_atom_size` to determine data type, dimensions and
+sizes of the storage pointed to by the returned pointer.
 
-A table with supported keywords is included in the documentation
+A table with supported keywords is included in the documentation of the
-of the :cpp:func:`Atom::extract() <LAMMPS_NS::Atom::extract>` function.
+:cpp:func:`Atom::extract() <LAMMPS_NS::Atom::extract>` function.
 
 .. warning::
 
@@ -7027,5 +7065,5 @@ int lammps_python_api_version() {
 }
 
 // Local Variables:
-// fill-column: 72
+// fill-column: 80
 // End:

src/library.h

@@ -172,6 +172,7 @@ int lammps_map_atom(void *handle, const void *id);
  * ---------------------------------------------------------------------- */
 
 int lammps_extract_atom_datatype(void *handle, const char *name);
+int lammps_extract_atom_size(void *handle, const char *name, int type);
 void *lammps_extract_atom(void *handle, const char *name);
 
 /* ----------------------------------------------------------------------

src/version.h

@@ -1 +1,2 @@
 #define LAMMPS_VERSION "29 Aug 2024"
+#define LAMMPS_UPDATE "Development"

tools/lammps-gui/CMakeLists.txt

@@ -1,6 +1,6 @@
 cmake_minimum_required(VERSION 3.16)
 
-project(lammps-gui VERSION 1.6.10 LANGUAGES CXX)
+project(lammps-gui VERSION 1.6.11 LANGUAGES CXX)
 
 set(CMAKE_AUTOUIC ON)
 set(CMAKE_AUTOMOC ON)
@@ -86,9 +86,7 @@ else()
       add_compile_options(/wd4244)
       add_compile_options(/wd4267)
       add_compile_options(/wd4250)
-      if(LAMMPS_EXCEPTIONS)
+      add_compile_options(/EHsc)
-        add_compile_options(/EHsc)
-      endif()
     endif()
     add_compile_definitions(_CRT_SECURE_NO_WARNINGS)
   endif()

tools/lammps-gui/lammps-gui.appdata.xml

@@ -54,6 +54,10 @@
   </provides>
 
   <releases>
+    <release version="1.6.11" timestamp="1725080055">
+      <description>
+      </description>
+    </release>
     <release version="1.6.10" timestamp="1724585189">
       <description>
         Resolve plugin mode issues.

tools/swig/lammps.i

@@ -130,6 +130,7 @@ extern void  *lammps_extract_pair(void *handle, const char *name);
 extern int    lammps_map_atom(void *handle, const void *id);
 
 extern int    lammps_extract_atom_datatype(void *handle, const char *name);
+extern int    lammps_extract_atom_size(void *handle, const char *name, int type);
 extern void  *lammps_extract_atom(void *handle, const char *name);
 
 extern void  *lammps_extract_compute(void *handle, const char *id, int, int);
@@ -319,6 +320,7 @@ extern void  *lammps_extract_pair(void *handle, const char *name);
 extern int    lammps_map_atom(void *handle, const void *id);
 
 extern int    lammps_extract_atom_datatype(void *handle, const char *name);
+extern int    lammps_extract_atom_size(void *handle, const char *name, int type);
 extern void  *lammps_extract_atom(void *handle, const char *name);
 
 extern void  *lammps_extract_compute(void *handle, const char *id, int, int);

unittest/c-library/test_library_properties.cpp

@@ -49,6 +49,7 @@ protected:
         if (verbose) std::cout << output;
         EXPECT_THAT(output, StartsWith("LAMMPS ("));
     }
+
     void TearDown() override
     {
         ::testing::internal::CaptureStdout();
@@ -470,9 +471,9 @@ TEST_F(LibraryProperties, global)
     EXPECT_EQ(lammps_extract_global_datatype(lmp, "xlattice"), LAMMPS_DOUBLE);
     EXPECT_EQ(lammps_extract_global_datatype(lmp, "ylattice"), LAMMPS_DOUBLE);
     EXPECT_EQ(lammps_extract_global_datatype(lmp, "zlattice"), LAMMPS_DOUBLE);
-    auto *xlattice   = (double *)lammps_extract_global(lmp, "xlattice");
+    auto *xlattice = (double *)lammps_extract_global(lmp, "xlattice");
-    auto *ylattice   = (double *)lammps_extract_global(lmp, "ylattice");
+    auto *ylattice = (double *)lammps_extract_global(lmp, "ylattice");
-    auto *zlattice   = (double *)lammps_extract_global(lmp, "zlattice");
+    auto *zlattice = (double *)lammps_extract_global(lmp, "zlattice");
     EXPECT_NE(xlattice, nullptr);
     EXPECT_NE(ylattice, nullptr);
     EXPECT_NE(zlattice, nullptr);
@@ -484,9 +485,9 @@ TEST_F(LibraryProperties, global)
     lammps_command(lmp, "units real");
     lammps_command(lmp, "lattice fcc 2.0");
     if (!verbose) ::testing::internal::GetCapturedStdout();
-    xlattice   = (double *)lammps_extract_global(lmp, "xlattice");
+    xlattice = (double *)lammps_extract_global(lmp, "xlattice");
-    ylattice   = (double *)lammps_extract_global(lmp, "ylattice");
+    ylattice = (double *)lammps_extract_global(lmp, "ylattice");
-    zlattice   = (double *)lammps_extract_global(lmp, "zlattice");
+    zlattice = (double *)lammps_extract_global(lmp, "zlattice");
     EXPECT_NE(xlattice, nullptr);
     EXPECT_NE(ylattice, nullptr);
     EXPECT_NE(zlattice, nullptr);
@@ -694,11 +695,10 @@ TEST_F(LibraryProperties, has_error)
 class AtomProperties : public ::testing::Test {
 protected:
     void *lmp;
+    int ntypes, nlocal, nall;
 
-    AtomProperties() = default;
+    AtomProperties()           = default;
-    ;
     ~AtomProperties() override = default;
-    ;
 
     void SetUp() override
     {
@@ -713,11 +713,30 @@ protected:
         if (verbose) std::cout << output;
         EXPECT_THAT(output, StartsWith("LAMMPS ("));
         ::testing::internal::CaptureStdout();
+        lammps_command(lmp, "fix props all property/atom i_one i2_two 2 d_three d2_four 2");
+        lammps_command(lmp, "fix rmass all property/atom mol q rmass ghost yes");
         lammps_command(lmp, "region box block 0 2 0 2 0 2");
         lammps_command(lmp, "create_box 1 box");
         lammps_command(lmp, "mass 1 3.0");
         lammps_command(lmp, "create_atoms 1 single 1.0 1.0 1.5");
         lammps_command(lmp, "create_atoms 1 single 0.2 0.1 0.1");
+        lammps_command(lmp, "set group all mass 2.0");
+        lammps_command(lmp, "set atom 1 charge -1");
+        lammps_command(lmp, "set atom 2 charge  1");
+        lammps_command(lmp, "set atom 1 mol 2");
+        lammps_command(lmp, "set atom 2 mol 1");
+        lammps_command(lmp, "set atom 1 i_one -3");
+        lammps_command(lmp, "set atom 2 i_one  3");
+        lammps_command(lmp, "set atom 1 d_three -1.3");
+        lammps_command(lmp, "set atom 2 d_three  3.5");
+        lammps_command(lmp, "set atom 1 i_two[1] -3");
+        lammps_command(lmp, "set atom 2 i_two[2]  3");
+        lammps_command(lmp, "set atom * d_four[1] -1.3");
+        lammps_command(lmp, "set atom * d_four[2]  3.5");
+        ntypes = lammps_extract_setting(lmp, "ntypes");
+        nlocal = lammps_extract_setting(lmp, "nlocal");
+        nall   = lammps_extract_setting(lmp, "nall");
+
         output = ::testing::internal::GetCapturedStdout();
         if (verbose) std::cout << output;
     }
@@ -740,14 +759,42 @@ TEST_F(AtomProperties, invalid)
 TEST_F(AtomProperties, mass)
 {
     EXPECT_EQ(lammps_extract_atom_datatype(lmp, "mass"), LAMMPS_DOUBLE);
+    EXPECT_EQ(lammps_extract_atom_size(lmp, "mass", 0), ntypes + 1);
     auto *mass = (double *)lammps_extract_atom(lmp, "mass");
     ASSERT_NE(mass, nullptr);
     ASSERT_DOUBLE_EQ(mass[1], 3.0);
+    EXPECT_EQ(lammps_extract_atom_datatype(lmp, "rmass"), LAMMPS_DOUBLE);
+    EXPECT_EQ(lammps_extract_atom_size(lmp, "rmass", 0), nall);
+    mass = (double *)lammps_extract_atom(lmp, "rmass");
+    ASSERT_NE(mass, nullptr);
+    ASSERT_DOUBLE_EQ(mass[0], 2.0);
+    ASSERT_DOUBLE_EQ(mass[1], 2.0);
+}
+
+TEST_F(AtomProperties, charge)
+{
+    EXPECT_EQ(lammps_extract_atom_datatype(lmp, "q"), LAMMPS_DOUBLE);
+    EXPECT_EQ(lammps_extract_atom_size(lmp, "rmass", 0), nall);
+    auto *charge = (double *)lammps_extract_atom(lmp, "q");
+    ASSERT_NE(charge, nullptr);
+    ASSERT_DOUBLE_EQ(charge[0], -1.0);
+    ASSERT_DOUBLE_EQ(charge[1], 1.0);
+}
+
+TEST_F(AtomProperties, molecule)
+{
+    EXPECT_EQ(lammps_extract_atom_datatype(lmp, "molecule"), LAMMPS_TAGINT);
+    EXPECT_EQ(lammps_extract_atom_size(lmp, "molecule", 0), nall);
+    auto *molecule = (tagint *)lammps_extract_atom(lmp, "molecule");
+    ASSERT_NE(molecule, nullptr);
+    ASSERT_EQ(molecule[0], 2);
+    ASSERT_EQ(molecule[1], 1);
 }
 
 TEST_F(AtomProperties, id)
 {
     EXPECT_EQ(lammps_extract_atom_datatype(lmp, "id"), LAMMPS_TAGINT);
+    EXPECT_EQ(lammps_extract_atom_size(lmp, "id", 0), nall);
     auto *id = (tagint *)lammps_extract_atom(lmp, "id");
     ASSERT_NE(id, nullptr);
     ASSERT_EQ(id[0], 1);
@@ -757,6 +804,7 @@ TEST_F(AtomProperties, id)
 TEST_F(AtomProperties, type)
 {
     EXPECT_EQ(lammps_extract_atom_datatype(lmp, "type"), LAMMPS_INT);
+    EXPECT_EQ(lammps_extract_atom_size(lmp, "type", 0), nall);
     int *type = (int *)lammps_extract_atom(lmp, "type");
     ASSERT_NE(type, nullptr);
     ASSERT_EQ(type[0], 1);
@@ -766,6 +814,8 @@ TEST_F(AtomProperties, type)
 TEST_F(AtomProperties, position)
 {
     EXPECT_EQ(lammps_extract_atom_datatype(lmp, "x"), LAMMPS_DOUBLE_2D);
+    EXPECT_EQ(lammps_extract_atom_size(lmp, "x", LMP_SIZE_ROWS), nall);
+    EXPECT_EQ(lammps_extract_atom_size(lmp, "x", LMP_SIZE_COLS), 3);
     auto **x = (double **)lammps_extract_atom(lmp, "x");
     ASSERT_NE(x, nullptr);
     EXPECT_DOUBLE_EQ(x[0][0], 1.0);
@@ -776,6 +826,41 @@ TEST_F(AtomProperties, position)
     EXPECT_DOUBLE_EQ(x[1][2], 0.1);
 }
 
+TEST_F(AtomProperties, custom)
+{
+    EXPECT_EQ(lammps_extract_atom_datatype(lmp, "i_one"), LAMMPS_INT);
+    EXPECT_EQ(lammps_extract_atom_size(lmp, "i_one", 0), nlocal);
+    auto *one = (int *)lammps_extract_atom(lmp, "i_one");
+    ASSERT_NE(one, nullptr);
+    EXPECT_EQ(lammps_extract_atom_datatype(lmp, "i2_two"), LAMMPS_INT_2D);
+    EXPECT_EQ(lammps_extract_atom_size(lmp, "i2_two", LMP_SIZE_ROWS), nlocal);
+    EXPECT_EQ(lammps_extract_atom_size(lmp, "i2_two", LMP_SIZE_COLS), 2);
+    auto **two = (int **)lammps_extract_atom(lmp, "i2_two");
+    ASSERT_NE(two, nullptr);
+    EXPECT_EQ(lammps_extract_atom_datatype(lmp, "d_three"), LAMMPS_DOUBLE);
+    EXPECT_EQ(lammps_extract_atom_size(lmp, "d_three", 0), nlocal);
+    auto *three = (double *)lammps_extract_atom(lmp, "d_three");
+    ASSERT_NE(three, nullptr);
+    EXPECT_EQ(lammps_extract_atom_datatype(lmp, "d2_four"), LAMMPS_DOUBLE_2D);
+    EXPECT_EQ(lammps_extract_atom_size(lmp, "d2_four", LMP_SIZE_ROWS), nlocal);
+    EXPECT_EQ(lammps_extract_atom_size(lmp, "d2_four", LMP_SIZE_COLS), 2);
+    auto **four = (double **)lammps_extract_atom(lmp, "d2_four");
+    ASSERT_NE(four, nullptr);
+
+    EXPECT_EQ(one[0], -3);
+    EXPECT_EQ(one[1], 3);
+    EXPECT_EQ(two[0][0], -3);
+    EXPECT_EQ(two[0][1], 0);
+    EXPECT_EQ(two[1][0], 0);
+    EXPECT_EQ(two[1][1], 3);
+    EXPECT_DOUBLE_EQ(three[0], -1.3);
+    EXPECT_DOUBLE_EQ(three[1], 3.5);
+    EXPECT_DOUBLE_EQ(four[0][0], -1.3);
+    EXPECT_DOUBLE_EQ(four[0][1], 3.5);
+    EXPECT_DOUBLE_EQ(four[1][0], -1.3);
+    EXPECT_DOUBLE_EQ(four[1][1], 3.5);
+}
+
 TEST(SystemSettings, kokkos)
 {
     if (!lammps_config_has_package("KOKKOS")) GTEST_SKIP();

unittest/force-styles/check_tests.py

@@ -88,10 +88,9 @@ for header in headers:
                 style = m[1]
                 if upper.match(style):
                     continue
-                if style in ['reax/c', 'reax/c/omp', 'reax/c/kk',
+                if style in ['lj/sdk', 'lj/sdk/coul/long', 'lj/sdk/coul/msm', 'sdk', 'lj/sdk/gpu',
-                        'reax/c/kk/device', 'reax/c/kk/host',
+                             'lj/sdk/coul/long/gpu', 'lj/sdk/omp', 'lj/sdk/coul/long/omp', 'sdk/omp',
-                        'reax/c/species', 'reax/c/bonds',
+                             'lj/sdk/coul/msm/omp', 'lj/sdk/kk', 'lj/sdk/coul/long/kk', 'sdk/kk']:
-                        'reax/c/species/kk', 'reax/c/bonds/kk', 'meam/c']:
                     continue
 
                 # detect, process, and flag suffix styles:
@@ -176,11 +175,12 @@ def check_tests(name,styles,yaml,search,skip=()):
 
 counter = 0
 counter += check_tests('pair',pair,'*-pair-*.yaml',
-                       '.*pair_style:\\s*((\\S+).*)?',skip=('meam','lj/sf'))
+                       '.*pair_style:\\s*((\\S+).*)?',
+                       skip=('lj/sf','lj/sdk', 'lj/sdk/coul/long', 'lj/sdk/coul/msm'))
 counter += check_tests('bond',bond,'bond-*.yaml',
                        '.*bond_style:\\s*((\\S+).*)?')
 counter += check_tests('angle',angle,'angle-*.yaml',
-                       '.*angle_style:\\s*((\\S+).*)?')
+                       '.*angle_style:\\s*((\\S+).*)?', skip=('sdk'))
 counter += check_tests('dihedral',dihedral,'dihedral-*.yaml',
                        '.*dihedral_style:\\s*((\\S+).*)?')
 counter += check_tests('improper',improper,'improper-*.yaml',

unittest/fortran/keepstuff.f90

@@ -4,9 +4,9 @@ MODULE keepstuff
   TYPE(LAMMPS), SAVE :: lmp
   INTEGER, SAVE :: mycomm
   CHARACTER(LEN=40), DIMENSION(3), PARAMETER :: demo_input = &
-      [ CHARACTER(LEN=40) ::                                &
+      [ CHARACTER(LEN=40) ::                                 &
-      'region       box block 0 $x 0 2 0 2',                &
+      'region       box block 0 $x 0 2 0 2',                 &
-      'create_box 1 box',                                   &
+      'create_box 1 box',                                    &
       'create_atoms 1 single 1.0 1.0 ${zpos}' ]
   CHARACTER(LEN=40), DIMENSION(3), PARAMETER :: big_input = &
       [ CHARACTER(LEN=40) ::                                &
@@ -14,15 +14,26 @@ MODULE keepstuff
       'create_box 1 box',                                   &
       'create_atoms 1 single 1.0 1.0 ${zpos}' ]
   CHARACTER(LEN=40), DIMENSION(2), PARAMETER :: cont_input = &
-      [ CHARACTER(LEN=40) ::                                &
+      [ CHARACTER(LEN=40) ::                                 &
-      'create_atoms 1 single &',                            &
+      'create_atoms 1 single &',                             &
       ' 0.2 0.1 0.1' ]
+  CHARACTER(LEN=60), DIMENSION(18), PARAMETER :: prop_input =         &
+      [ CHARACTER(LEN=60) ::  'fix 1 all nve', 'mass 1 3.0',          &
+      'fix 2 all property/atom mol q rmass ghost yes',                &
+      'fix 3 all property/atom i_one i2_two 2 d_three d2_four 2',     &
+      'set group all mass 2.0', 'set atom 1 charge -1',               &
+      'set atom 2 charge  1', 'set atom 1 mol 2', 'set atom 2 mol 1', &
+      'set atom 1 i_one -3', 'set atom 2 i_one  3',                   &
+      'set atom 1 d_three -1.3', 'set atom 2 d_three  3.5',           &
+      'set atom 1 i_two[1] -3', 'set atom 2 i_two[2]  3',             &
+      'set atom * d_four[1] -1.3', 'set atom * d_four[2]  3.5',       &
+      'run 0 post no' ]
   CHARACTER(LEN=40), DIMENSION(1), PARAMETER :: more_input = &
       [ CHARACTER(LEN=40) :: 'create_atoms 1 single 0.5 0.5 0.5' ]
   CHARACTER(LEN=40), DIMENSION(3), PARAMETER :: pair_input = &
-      [ CHARACTER(LEN=40) ::                                &
+      [ CHARACTER(LEN=40) ::                                 &
-      'pair_style lj/cut 2.5',                              &
+      'pair_style lj/cut 2.5',                               &
-      'pair_coeff 1 1 1.0 1.0',                             &
+      'pair_coeff 1 1 1.0 1.0',                              &
       'mass 1 2.0' ]
 
   INTERFACE
@@ -63,4 +74,3 @@ CONTAINS
   END FUNCTION f2c_string
 
 END MODULE keepstuff
-

unittest/fortran/test_fortran_extract_atom.f90

@@ -24,12 +24,13 @@ END SUBROUTINE f_lammps_close
 
 SUBROUTINE f_lammps_setup_extract_atom() BIND(C)
    USE LIBLAMMPS
-   USE keepstuff, ONLY : lmp, big_input, cont_input, pair_input
+   USE keepstuff, ONLY : lmp, big_input, cont_input, pair_input, prop_input
    IMPLICIT NONE
 
    CALL lmp%commands_list(big_input)
    CALL lmp%commands_list(cont_input)
    CALL lmp%commands_list(pair_input)
+   CALL lmp%commands_list(prop_input)
 END SUBROUTINE f_lammps_setup_extract_atom
 
 FUNCTION f_lammps_extract_atom_mass() BIND(C)
@@ -44,6 +45,19 @@ FUNCTION f_lammps_extract_atom_mass() BIND(C)
    f_lammps_extract_atom_mass = mass(1)
 END FUNCTION f_lammps_extract_atom_mass
 
+FUNCTION f_lammps_extract_atom_mass_size() BIND(C)
+   USE, INTRINSIC :: ISO_C_BINDING, ONLY : c_double, c_int
+   USE LIBLAMMPS
+   USE keepstuff, ONLY : lmp
+   IMPLICIT NONE
+   INTEGER(c_int) :: f_lammps_extract_atom_mass_size, ntypes
+   REAL(c_double), DIMENSION(:), POINTER :: mass => NULL()
+
+   ntypes = lmp%extract_setting('ntypes')
+   mass = lmp%extract_atom('mass')
+   f_lammps_extract_atom_mass_size = SIZE(mass)
+END FUNCTION f_lammps_extract_atom_mass_size
+
 FUNCTION f_lammps_extract_atom_tag_int(i) BIND(C)
    USE, INTRINSIC :: ISO_C_BINDING, ONLY : c_double, c_int
    USE LIBLAMMPS
@@ -83,6 +97,18 @@ FUNCTION f_lammps_extract_atom_type(i) BIND(C)
    f_lammps_extract_atom_type = atype(i)
 END FUNCTION f_lammps_extract_atom_type
 
+FUNCTION f_lammps_extract_atom_type_size() BIND(C)
+   USE, INTRINSIC :: ISO_C_BINDING, ONLY : c_int
+   USE LIBLAMMPS
+   USE keepstuff, ONLY : lmp
+   IMPLICIT NONE
+   INTEGER(c_int) :: f_lammps_extract_atom_type_size
+   INTEGER(c_int), DIMENSION(:), POINTER :: atype => NULL()
+
+   atype = lmp%extract_atom('type')
+   f_lammps_extract_atom_type_size = size(atype)
+END FUNCTION f_lammps_extract_atom_type_size
+
 FUNCTION f_lammps_extract_atom_mask(i) BIND(C)
    USE, INTRINSIC :: ISO_C_BINDING, ONLY : c_int
    USE LIBLAMMPS
@@ -109,6 +135,19 @@ SUBROUTINE f_lammps_extract_atom_x(i, x) BIND(C)
    x = xptr(:,i)
 END SUBROUTINE f_lammps_extract_atom_x
 
+FUNCTION f_lammps_extract_atom_x_size(i) BIND(C)
+   USE, INTRINSIC :: ISO_C_BINDING, ONLY : c_double, c_int
+   USE LIBLAMMPS
+   USE keepstuff, ONLY : lmp
+   IMPLICIT NONE
+   INTEGER(c_int), INTENT(IN), VALUE :: i
+   INTEGER(c_int) :: f_lammps_extract_atom_x_size
+   REAL(c_double), DIMENSION(:,:), POINTER :: xptr => NULL()
+
+   xptr = lmp%extract_atom('x')
+   f_lammps_extract_atom_x_size = SIZE(xptr, i)
+END FUNCTION f_lammps_extract_atom_x_size
+
 SUBROUTINE f_lammps_extract_atom_v(i, v) BIND(C)
    USE, INTRINSIC :: ISO_C_BINDING, ONLY : c_double, c_int
    USE LIBLAMMPS
@@ -121,3 +160,16 @@ SUBROUTINE f_lammps_extract_atom_v(i, v) BIND(C)
    vptr = lmp%extract_atom('v')
    v = vptr(:,i)
 END SUBROUTINE f_lammps_extract_atom_v
+
+FUNCTION f_lammps_extract_atom_v_size(i) BIND(C)
+   USE, INTRINSIC :: ISO_C_BINDING, ONLY : c_double, c_int
+   USE LIBLAMMPS
+   USE keepstuff, ONLY : lmp
+   IMPLICIT NONE
+   INTEGER(c_int), INTENT(IN), VALUE :: i
+   INTEGER(c_int) :: f_lammps_extract_atom_v_size
+   REAL(c_double), DIMENSION(:,:), POINTER :: xptr => NULL()
+
+   xptr = lmp%extract_atom('v')
+   f_lammps_extract_atom_v_size = SIZE(xptr, i)
+END FUNCTION f_lammps_extract_atom_v_size

unittest/fortran/wrap_extract_atom.cpp

@@ -1,6 +1,7 @@
 // unit tests for extracting Atom class data from a LAMMPS instance through the
 // Fortran wrapper
 
+#include "atom.h"
 #include "lammps.h"
 #include "library.h"
 #include <cstdint>
@@ -16,12 +17,16 @@ void *f_lammps_with_args();
 void f_lammps_close();
 void f_lammps_setup_extract_atom();
 double f_lammps_extract_atom_mass();
+int f_lammps_extract_atom_mass_size();
 int f_lammps_extract_atom_tag_int(int);
 int64_t f_lammps_extract_atom_tag_int64(int64_t);
 int f_lammps_extract_atom_type(int);
+int f_lammps_extract_atom_type_size();
 int f_lammps_extract_atom_mask(int);
 void f_lammps_extract_atom_x(int, double *);
+int f_lammps_extract_atom_x_size(int);
 void f_lammps_extract_atom_v(int, double *);
+int f_lammps_extract_atom_v_size(int);
 }
 
 class LAMMPS_extract_atom : public ::testing::Test {
@@ -50,7 +55,9 @@ protected:
 TEST_F(LAMMPS_extract_atom, mass)
 {
     f_lammps_setup_extract_atom();
-    EXPECT_DOUBLE_EQ(f_lammps_extract_atom_mass(), 2.0);
+    int ntypes = lmp->atom->ntypes;
+    EXPECT_DOUBLE_EQ(f_lammps_extract_atom_mass(), 3.0);
+    EXPECT_EQ(f_lammps_extract_atom_mass_size(), ntypes + 1);
 };
 
 TEST_F(LAMMPS_extract_atom, tag)
@@ -68,8 +75,10 @@ TEST_F(LAMMPS_extract_atom, tag)
 TEST_F(LAMMPS_extract_atom, type)
 {
     f_lammps_setup_extract_atom();
+    int nall = lmp->atom->nlocal + lmp->atom->nghost;
     EXPECT_EQ(f_lammps_extract_atom_type(1), 1);
     EXPECT_EQ(f_lammps_extract_atom_type(2), 1);
+    EXPECT_EQ(f_lammps_extract_atom_type_size(), nall);
 };
 
 TEST_F(LAMMPS_extract_atom, mask)
@@ -86,6 +95,7 @@ TEST_F(LAMMPS_extract_atom, mask)
 TEST_F(LAMMPS_extract_atom, x)
 {
     f_lammps_setup_extract_atom();
+    int nall = lmp->atom->nlocal + lmp->atom->nghost;
     double x1[3];
     double x2[3];
     f_lammps_extract_atom_x(1, x1);
@@ -96,11 +106,15 @@ TEST_F(LAMMPS_extract_atom, x)
     EXPECT_DOUBLE_EQ(x2[0], 0.2);
     EXPECT_DOUBLE_EQ(x2[1], 0.1);
     EXPECT_DOUBLE_EQ(x2[2], 0.1);
+    // in Fortran row and column are swapped
+    EXPECT_EQ(f_lammps_extract_atom_x_size(1), 3);
+    EXPECT_EQ(f_lammps_extract_atom_x_size(2), nall);
 }
 
 TEST_F(LAMMPS_extract_atom, v)
 {
     f_lammps_setup_extract_atom();
+    int nall = lmp->atom->nlocal + lmp->atom->nghost;
     double v1[3];
     double v2[3];
     f_lammps_extract_atom_v(1, v1);
@@ -117,4 +131,13 @@ TEST_F(LAMMPS_extract_atom, v)
     EXPECT_DOUBLE_EQ(v1[0], 1.0);
     EXPECT_DOUBLE_EQ(v1[1], 2.0);
     EXPECT_DOUBLE_EQ(v1[2], 3.0);
+    // in Fortran row and column are swapped!
+    EXPECT_EQ(f_lammps_extract_atom_v_size(1), 3);
+    EXPECT_EQ(f_lammps_extract_atom_v_size(2), lmp->atom->nlocal);
+    lammps_command(lmp, "comm_modify vel yes");
+    lammps_command(lmp, "run 0 post no");
+    EXPECT_EQ(f_lammps_extract_atom_v_size(1), 3);
+    EXPECT_EQ(f_lammps_extract_atom_v_size(2), nall);
 }
+
+// TODO: write tests for custom properties

unittest/python/python-numpy.py

@@ -155,67 +155,104 @@ class PythonNumpy(unittest.TestCase):
         self.assertEqual(values[1,0], 1.5)
         self.assertEqual(values[1,3], 1.5)
 
-    def testExtractAtomDeprecated(self):
-        self.lmp.command("units lj")
-        self.lmp.command("atom_style atomic")
-        self.lmp.command("atom_modify map array")
-        self.lmp.command("region box block 0 2 0 2 0 2")
-        self.lmp.command("create_box 1 box")
-
-        x = [
-          1.0, 1.0, 1.0,
-          1.0, 1.0, 1.5
-        ]
-
-        types = [1, 1]
-
-        self.assertEqual(self.lmp.create_atoms(2, id=None, type=types, x=x), 2)
-        nlocal = self.lmp.extract_global("nlocal", LAMMPS_INT)
-        self.assertEqual(nlocal, 2)
-
-        ident = self.lmp.numpy.extract_atom_iarray("id", nlocal, dim=1)
-        self.assertEqual(len(ident), 2)
-
-        ntypes = self.lmp.extract_global("ntypes", LAMMPS_INT)
-        self.assertEqual(ntypes, 1)
-
-        x = self.lmp.numpy.extract_atom_darray("x", nlocal, dim=3)
-        v = self.lmp.numpy.extract_atom_darray("v", nlocal, dim=3)
-        self.assertEqual(len(x), 2)
-        self.assertTrue((x[0] == (1.0, 1.0, 1.0)).all())
-        self.assertTrue((x[1] == (1.0, 1.0, 1.5)).all())
-        self.assertEqual(len(v), 2)
-
     def testExtractAtom(self):
         self.lmp.command("units lj")
         self.lmp.command("atom_style atomic")
         self.lmp.command("atom_modify map array")
         self.lmp.command("region box block 0 2 0 2 0 2")
-        self.lmp.command("create_box 1 box")
+        self.lmp.command("create_box 2 box")
 
-        x = [
+        x = [ 1.0, 1.0, 1.0,  1.0, 1.0, 1.5,  1.5, 1.0, 1.0 ]
-          1.0, 1.0, 1.0,
+        types = [1, 2, 1]
-          1.0, 1.0, 1.5
+        ids = [1, 2, 3]
-        ]
+        self.assertEqual(self.lmp.create_atoms(3, id=ids, type=types, x=x), 3)
+        self.lmp.command("mass * 2.0")
+        self.lmp.command("pair_style zero 1.1")
+        self.lmp.command("pair_coeff * *")
+        self.lmp.command("fix props all property/atom i_one i2_two 2 d_three d2_four 2");
+        self.lmp.command("fix rmass all property/atom mol q rmass ghost yes");
+        self.lmp.command("fix 1 all nve")
+        self.lmp.command("run 0 post no")
+        ntypes = self.lmp.extract_setting("ntypes");
+        nlocal = self.lmp.extract_setting("nlocal");
+        nall   = self.lmp.extract_setting("nall");
+        self.assertEqual(nlocal, 3)
+        self.assertEqual(ntypes, 2)
+        self.assertEqual(nall, 63)
 
-        types = [1, 1]
+        self.lmp.command("set atom 1 charge -1");
+        self.lmp.command("set atom 2 charge  1");
+        self.lmp.command("set atom 3 charge  0");
+        self.lmp.command("set atom * mol 2");
+        self.lmp.command("set atom 2 mol 1");
+        self.lmp.command("set atom 1 i_one -3");
+        self.lmp.command("set atom 2 i_one  3");
+        self.lmp.command("set atom 2 d_three -1.3");
+        self.lmp.command("set atom 3 d_three  3.5");
+        self.lmp.command("set atom 1 i_two[1] -3");
+        self.lmp.command("set atom 2 i_two[2]  3");
+        self.lmp.command("set atom * d_four[1] -1.3");
+        self.lmp.command("set atom * d_four[2]  3.5");
+        self.lmp.command("run 0 post no")
 
-        self.assertEqual(self.lmp.create_atoms(2, id=None, type=types, x=x), 2)
+        mass = self.lmp.numpy.extract_atom("mass")
-        nlocal = self.lmp.extract_global("nlocal")
+        self.assertEqual(len(mass), ntypes + 1)
-        self.assertEqual(nlocal, 2)
+        self.assertTrue((mass == (0.0, 2.0, 2.0)).all())
+
+        rmass = self.lmp.numpy.extract_atom("rmass")
+        self.assertEqual(len(rmass), nall)
+        self.assertTrue((rmass[0:3] == (0.0, 0.0, 0.0)).all())
+
+        charge = self.lmp.numpy.extract_atom("q")
+        self.assertEqual(len(charge), nall)
+        self.assertTrue((charge[0:3] == (-1.0, 1.0, 0.0)).all())
+
+        molecule = self.lmp.numpy.extract_atom("molecule")
+        self.assertEqual(len(molecule), nall)
+        self.assertTrue((molecule[0:3] == (2, 1, 2)).all())
 
         ident = self.lmp.numpy.extract_atom("id")
-        self.assertEqual(len(ident), 2)
+        self.assertEqual(len(ident), nall)
+        self.assertTrue((ident[0:3] == (1, 2, 3)).all())
 
-        ntypes = self.lmp.extract_global("ntypes")
+        atype = self.lmp.numpy.extract_atom("type")
-        self.assertEqual(ntypes, 1)
+        self.assertEqual(len(atype), nall)
+        self.assertTrue((atype[0:3] == (1, 2, 1)).all())
 
         x = self.lmp.numpy.extract_atom("x")
         v = self.lmp.numpy.extract_atom("v")
-        self.assertEqual(len(x), 2)
+        self.assertEqual(len(x), nall)
+        self.assertEqual(len(x[0]), 3)
         self.assertTrue((x[0] == (1.0, 1.0, 1.0)).all())
         self.assertTrue((x[1] == (1.0, 1.0, 1.5)).all())
-        self.assertEqual(len(v), 2)
+        self.assertTrue((x[2] == (1.5, 1.0, 1.0)).all())
+        self.assertEqual(len(v), nlocal)
+        self.assertEqual(len(v[0]), 3)
+
+        self.lmp.command("comm_modify vel yes");
+        self.lmp.command("run 0 post no")
+
+        v = self.lmp.numpy.extract_atom("v")
+        self.assertEqual(len(v), nall)
+
+        one = self.lmp.numpy.extract_atom("i_one")
+        two = self.lmp.numpy.extract_atom("i2_two")
+        three = self.lmp.numpy.extract_atom("d_three")
+        four = self.lmp.numpy.extract_atom("d2_four")
+        self.assertEqual(len(one), nlocal)
+        self.assertTrue((one == (-3, 3, 0)).all())
+        self.assertEqual(len(two), nlocal)
+        self.assertEqual(len(two[0]), 2)
+        self.assertTrue((two[0] == (-3, 0)).all())
+        self.assertTrue((two[1] == (0, 3)).all())
+        self.assertTrue((two[2] == (0, 0)).all())
+        self.assertEqual(len(three), nlocal)
+        self.assertTrue((three == (0.0, -1.3, 3.5)).all())
+        self.assertEqual(len(four), nlocal)
+        self.assertEqual(len(four[0]), 2)
+        self.assertTrue((four[0] == (-1.3, 3.5)).all())
+        self.assertTrue((four[1] == (-1.3, 3.5)).all())
+        self.assertTrue((four[2] == (-1.3, 3.5)).all())
 
     @unittest.skipIf(not has_full,"Gather bonds test")
     def testGatherBond_newton_on(self):