diff --git a/.github/dependabot.yml b/.github/dependabot.yml new file mode 100644 index 0000000000..5ace4600a1 --- /dev/null +++ b/.github/dependabot.yml @@ -0,0 +1,6 @@ +version: 2 +updates: + - package-ecosystem: "github-actions" + directory: "/" + schedule: + interval: "weekly" diff --git a/.github/workflows/codeql-analysis.yml b/.github/workflows/codeql-analysis.yml index 7ea31fcf14..07f7564727 100644 --- a/.github/workflows/codeql-analysis.yml +++ b/.github/workflows/codeql-analysis.yml @@ -13,6 +13,11 @@ jobs: if: ${{ github.repository == 'lammps/lammps' }} runs-on: ubuntu-latest + permissions: + security-events: write + actions: read + contents: read + strategy: fail-fast: false matrix: @@ -20,17 +25,17 @@ jobs: steps: - name: Checkout repository - uses: actions/checkout@v2 + uses: actions/checkout@v3 with: fetch-depth: 2 - name: Setup Python - uses: actions/setup-python@v2 + uses: actions/setup-python@v4 with: python-version: '3.x' - name: Initialize CodeQL - uses: github/codeql-action/init@v1 + uses: github/codeql-action/init@v2 with: languages: ${{ matrix.language }} config-file: ./.github/codeql/${{ matrix.language }}.yml @@ -48,4 +53,4 @@ jobs: cmake --build . --parallel 2 - name: Perform CodeQL Analysis - uses: github/codeql-action/analyze@v1 + uses: github/codeql-action/analyze@v2 diff --git a/.github/workflows/compile-msvc.yml b/.github/workflows/compile-msvc.yml index 15fcf1099d..7747be7b46 100644 --- a/.github/workflows/compile-msvc.yml +++ b/.github/workflows/compile-msvc.yml @@ -15,12 +15,12 @@ jobs: steps: - name: Checkout repository - uses: actions/checkout@v2 + uses: actions/checkout@v3 with: fetch-depth: 2 - name: Select Python version - uses: actions/setup-python@v2 + uses: actions/setup-python@v4 with: python-version: '3.10' diff --git a/.github/workflows/unittest-macos.yml b/.github/workflows/unittest-macos.yml index 2d903af646..e6e5ccfdc8 100644 --- a/.github/workflows/unittest-macos.yml +++ b/.github/workflows/unittest-macos.yml @@ -17,7 +17,7 @@ jobs: steps: - name: Checkout repository - uses: actions/checkout@v2 + uses: actions/checkout@v3 with: fetch-depth: 2 @@ -28,7 +28,7 @@ jobs: run: mkdir build - name: Set up ccache - uses: actions/cache@v2 + uses: actions/cache@v3 with: path: ${{ env.CCACHE_DIR }} key: macos-ccache-${{ github.sha }} diff --git a/cmake/CMakeLists.txt b/cmake/CMakeLists.txt index 2566497c0e..b0b8bfd363 100644 --- a/cmake/CMakeLists.txt +++ b/cmake/CMakeLists.txt @@ -135,11 +135,13 @@ set(CMAKE_CXX_EXTENSIONS OFF CACHE BOOL "Use compiler extensions") # ugly hacks for MSVC which by default always reports an old C++ standard in the __cplusplus macro # and prints lots of pointless warnings about "unsafe" functions if(MSVC) - add_compile_options(/Zc:__cplusplus) - add_compile_options(/wd4244) - add_compile_options(/wd4267) - if(LAMMPS_EXCEPTIONS) - add_compile_options(/EHsc) + if(CMAKE_CXX_COMPILER_ID STREQUAL "MSVC") + add_compile_options(/Zc:__cplusplus) + add_compile_options(/wd4244) + add_compile_options(/wd4267) + if(LAMMPS_EXCEPTIONS) + add_compile_options(/EHsc) + endif() endif() add_compile_definitions(_CRT_SECURE_NO_WARNINGS) endif() @@ -784,14 +786,16 @@ if(BUILD_SHARED_LIBS) find_package(Python COMPONENTS Interpreter) endif() if(BUILD_IS_MULTI_CONFIG) - set(LIBLAMMPS_SHARED_BINARY ${CMAKE_BINARY_DIR}/$/liblammps${LAMMPS_MACHINE}${CMAKE_SHARED_LIBRARY_SUFFIX}) + set(MY_BUILD_DIR ${CMAKE_BINARY_DIR}/$) else() - set(LIBLAMMPS_SHARED_BINARY ${CMAKE_BINARY_DIR}/liblammps${LAMMPS_MACHINE}${CMAKE_SHARED_LIBRARY_SUFFIX}) + set(MY_BUILD_DIR ${CMAKE_BINARY_DIR}) endif() + set(LIBLAMMPS_SHARED_BINARY ${MY_BUILD_DIR}/liblammps${LAMMPS_MACHINE}${CMAKE_SHARED_LIBRARY_SUFFIX}) if(Python_EXECUTABLE) add_custom_target( install-python ${CMAKE_COMMAND} -E remove_directory build - COMMAND ${Python_EXECUTABLE} ${LAMMPS_PYTHON_DIR}/install.py -p ${LAMMPS_PYTHON_DIR}/lammps -l ${LIBLAMMPS_SHARED_BINARY} + COMMAND ${Python_EXECUTABLE} ${LAMMPS_PYTHON_DIR}/install.py -p ${LAMMPS_PYTHON_DIR}/lammps + -l ${LIBLAMMPS_SHARED_BINARY} -w ${MY_BUILD_DIR} COMMENT "Installing LAMMPS Python module") else() add_custom_target( diff --git a/cmake/Modules/LAMMPSInterfacePlugin.cmake b/cmake/Modules/LAMMPSInterfacePlugin.cmake new file mode 100644 index 0000000000..ec7a739785 --- /dev/null +++ b/cmake/Modules/LAMMPSInterfacePlugin.cmake @@ -0,0 +1,197 @@ +# CMake script code to define LAMMPS settings required for building LAMMPS plugins + +# enforce out-of-source build +if(${CMAKE_SOURCE_DIR} STREQUAL ${CMAKE_BINARY_DIR}) + message(FATAL_ERROR "In-source builds are not allowed. You must create and use a build directory. " + "Please remove CMakeCache.txt and CMakeFiles first.") +endif() + +# global LAMMPS/plugin build settings +set(LAMMPS_SOURCE_DIR "" CACHE PATH "Location of LAMMPS sources folder") +if(NOT LAMMPS_SOURCE_DIR) + message(FATAL_ERROR "Must set LAMMPS_SOURCE_DIR") +endif() + +# by default, install into $HOME/.local (not /usr/local), +# so that no root access (and sudo) is needed +if(CMAKE_INSTALL_PREFIX_INITIALIZED_TO_DEFAULT) + set(CMAKE_INSTALL_PREFIX "$ENV{HOME}/.local" CACHE PATH "Default install path" FORCE) +endif() + +# ugly hacks for MSVC which by default always reports an old C++ standard in the __cplusplus macro +# and prints lots of pointless warnings about "unsafe" functions +if(MSVC) + if(CMAKE_CXX_COMPILER_ID STREQUAL "MSVC") + add_compile_options(/Zc:__cplusplus) + add_compile_options(/wd4244) + add_compile_options(/wd4267) + if(LAMMPS_EXCEPTIONS) + add_compile_options(/EHsc) + endif() + endif() + add_compile_definitions(_CRT_SECURE_NO_WARNINGS) +endif() + +# C++11 is required +set(CMAKE_CXX_STANDARD 11) +set(CMAKE_CXX_STANDARD_REQUIRED ON) + +# Need -restrict with Intel compilers +if(CMAKE_CXX_COMPILER_ID STREQUAL "Intel") + set(CMAKE_CXX_FLAGS "${CMAKE_CXX_FLAGS} -restrict") +endif() +set(CMAKE_POSITION_INDEPENDENT_CODE TRUE) + +####### +# helper functions from LAMMPSUtils.cmake +function(validate_option name values) + string(TOLOWER ${${name}} needle_lower) + string(TOUPPER ${${name}} needle_upper) + list(FIND ${values} ${needle_lower} IDX_LOWER) + list(FIND ${values} ${needle_upper} IDX_UPPER) + if(${IDX_LOWER} LESS 0 AND ${IDX_UPPER} LESS 0) + list_to_bulletpoints(POSSIBLE_VALUE_LIST ${${values}}) + message(FATAL_ERROR "\n########################################################################\n" + "Invalid value '${${name}}' for option ${name}\n" + "\n" + "Possible values are:\n" + "${POSSIBLE_VALUE_LIST}" + "########################################################################") + endif() +endfunction(validate_option) + +# helper function for getting the most recently modified file or folder from a glob pattern +function(get_newest_file path variable) + file(GLOB _dirs ${path}) + set(_besttime 2000-01-01T00:00:00) + set(_bestfile "") + foreach(_dir ${_dirs}) + file(TIMESTAMP ${_dir} _newtime) + if(_newtime IS_NEWER_THAN _besttime) + set(_bestfile ${_dir}) + set(_besttime ${_newtime}) + endif() + endforeach() + if(_bestfile STREQUAL "") + message(FATAL_ERROR "Could not find valid path at: ${path}") + endif() + set(${variable} ${_bestfile} PARENT_SCOPE) +endfunction() + +################################################################################# +# LAMMPS C++ interface. We only need the header related parts except on windows. +add_library(lammps INTERFACE) +target_include_directories(lammps INTERFACE ${LAMMPS_SOURCE_DIR}) +if((CMAKE_SYSTEM_NAME STREQUAL "Windows") AND CMAKE_CROSSCOMPILING) + target_link_libraries(lammps INTERFACE ${CMAKE_BINARY_DIR}/../liblammps.dll.a) +endif() + +################################################################################ +# MPI configuration +if(NOT CMAKE_CROSSCOMPILING) + find_package(MPI QUIET) + option(BUILD_MPI "Build MPI version" ${MPI_FOUND}) +else() + option(BUILD_MPI "Build MPI version" OFF) +endif() + +if(BUILD_MPI) + # do not include the (obsolete) MPI C++ bindings which makes + # for leaner object files and avoids namespace conflicts + set(MPI_CXX_SKIP_MPICXX TRUE) + # We use a non-standard procedure to cross-compile with MPI on Windows + if((CMAKE_SYSTEM_NAME STREQUAL "Windows") AND CMAKE_CROSSCOMPILING) + # Download and configure custom MPICH files for Windows + message(STATUS "Downloading and configuring MPICH-1.4.1 for Windows") + set(MPICH2_WIN64_DEVEL_URL "${LAMMPS_THIRDPARTY_URL}/mpich2-win64-devel.tar.gz" CACHE STRING "URL for MPICH2 (win64) tarball") + set(MPICH2_WIN32_DEVEL_URL "${LAMMPS_THIRDPARTY_URL}/mpich2-win32-devel.tar.gz" CACHE STRING "URL for MPICH2 (win32) tarball") + set(MPICH2_WIN64_DEVEL_MD5 "4939fdb59d13182fd5dd65211e469f14" CACHE STRING "MD5 checksum of MPICH2 (win64) tarball") + set(MPICH2_WIN32_DEVEL_MD5 "a61d153500dce44e21b755ee7257e031" CACHE STRING "MD5 checksum of MPICH2 (win32) tarball") + mark_as_advanced(MPICH2_WIN64_DEVEL_URL) + mark_as_advanced(MPICH2_WIN32_DEVEL_URL) + mark_as_advanced(MPICH2_WIN64_DEVEL_MD5) + mark_as_advanced(MPICH2_WIN32_DEVEL_MD5) + + include(ExternalProject) + if(CMAKE_SYSTEM_PROCESSOR STREQUAL "x86_64") + ExternalProject_Add(mpi4win_build + URL ${MPICH2_WIN64_DEVEL_URL} + URL_MD5 ${MPICH2_WIN64_DEVEL_MD5} + CONFIGURE_COMMAND "" BUILD_COMMAND "" INSTALL_COMMAND "" + BUILD_BYPRODUCTS /lib/libmpi.a) + else() + ExternalProject_Add(mpi4win_build + URL ${MPICH2_WIN32_DEVEL_URL} + URL_MD5 ${MPICH2_WIN32_DEVEL_MD5} + CONFIGURE_COMMAND "" BUILD_COMMAND "" INSTALL_COMMAND "" + BUILD_BYPRODUCTS /lib/libmpi.a) + endif() + + ExternalProject_get_property(mpi4win_build SOURCE_DIR) + file(MAKE_DIRECTORY "${SOURCE_DIR}/include") + add_library(MPI::MPI_CXX UNKNOWN IMPORTED) + set_target_properties(MPI::MPI_CXX PROPERTIES + IMPORTED_LOCATION "${SOURCE_DIR}/lib/libmpi.a" + INTERFACE_INCLUDE_DIRECTORIES "${SOURCE_DIR}/include" + INTERFACE_COMPILE_DEFINITIONS "MPICH_SKIP_MPICXX") + add_dependencies(MPI::MPI_CXX mpi4win_build) + + # set variables for status reporting at the end of CMake run + set(MPI_CXX_INCLUDE_PATH "${SOURCE_DIR}/include") + set(MPI_CXX_COMPILE_DEFINITIONS "MPICH_SKIP_MPICXX") + set(MPI_CXX_LIBRARIES "${SOURCE_DIR}/lib/libmpi.a") + else() + find_package(MPI REQUIRED) + option(LAMMPS_LONGLONG_TO_LONG "Workaround if your system or MPI version does not recognize 'long long' data types" OFF) + if(LAMMPS_LONGLONG_TO_LONG) + target_compile_definitions(lammps INTERFACE -DLAMMPS_LONGLONG_TO_LONG) + endif() + endif() + target_link_libraries(lammps INTERFACE MPI::MPI_CXX) +else() + add_library(mpi_stubs INTERFACE) + target_include_directories(mpi_stubs INTERFACE $) + target_link_libraries(lammps INTERFACE mpi_stubs) +endif() + +################################################################################ +# detect if we may enable OpenMP support by default +set(BUILD_OMP_DEFAULT OFF) +find_package(OpenMP QUIET) +if(OpenMP_FOUND) + check_include_file_cxx(omp.h HAVE_OMP_H_INCLUDE) + if(HAVE_OMP_H_INCLUDE) + set(BUILD_OMP_DEFAULT ON) + endif() +endif() + +option(BUILD_OMP "Build with OpenMP support" ${BUILD_OMP_DEFAULT}) + +if(BUILD_OMP) + find_package(OpenMP REQUIRED) + check_include_file_cxx(omp.h HAVE_OMP_H_INCLUDE) + if(NOT HAVE_OMP_H_INCLUDE) + message(FATAL_ERROR "Cannot find the 'omp.h' header file required for full OpenMP support") + endif() + + if (((CMAKE_CXX_COMPILER_ID STREQUAL "GNU") AND (CMAKE_CXX_COMPILER_VERSION VERSION_GREATER_EQUAL 9.0)) OR + (CMAKE_CXX_COMPILER_ID STREQUAL "PGI") OR + ((CMAKE_CXX_COMPILER_ID STREQUAL "Clang") AND (CMAKE_CXX_COMPILER_VERSION VERSION_GREATER_EQUAL 10.0)) OR + ((CMAKE_CXX_COMPILER_ID STREQUAL "Intel") AND (CMAKE_CXX_COMPILER_VERSION VERSION_GREATER_EQUAL 19.0))) + # GCC 9.x and later plus Clang 10.x and later implement strict OpenMP 4.0 semantics for consts. + # Intel 18.0 was tested to support both, so we switch to OpenMP 4+ from 19.x onward to be safe. + target_compile_definitions(lammps INTERFACE -DLAMMPS_OMP_COMPAT=4) + else() + target_compile_definitions(lammps INTERFACE -DLAMMPS_OMP_COMPAT=3) + endif() + target_link_libraries(lammps INTERFACE OpenMP::OpenMP_CXX) +endif() + +################ +# integer size selection +set(LAMMPS_SIZES "smallbig" CACHE STRING "LAMMPS integer sizes (smallsmall: all 32-bit, smallbig: 64-bit #atoms #timesteps, bigbig: also 64-bit imageint, 64-bit atom ids)") +set(LAMMPS_SIZES_VALUES smallbig bigbig smallsmall) +set_property(CACHE LAMMPS_SIZES PROPERTY STRINGS ${LAMMPS_SIZES_VALUES}) +validate_option(LAMMPS_SIZES LAMMPS_SIZES_VALUES) +string(TOUPPER ${LAMMPS_SIZES} LAMMPS_SIZES) +target_compile_definitions(lammps INTERFACE -DLAMMPS_${LAMMPS_SIZES}) diff --git a/cmake/Modules/Packages/MDI.cmake b/cmake/Modules/Packages/MDI.cmake index 88859f0285..90188f91a8 100644 --- a/cmake/Modules/Packages/MDI.cmake +++ b/cmake/Modules/Packages/MDI.cmake @@ -8,8 +8,8 @@ option(DOWNLOAD_MDI "Download and compile the MDI library instead of using an al if(DOWNLOAD_MDI) message(STATUS "MDI download requested - we will build our own") - set(MDI_URL "https://github.com/MolSSI-MDI/MDI_Library/archive/v1.3.2.tar.gz" CACHE STRING "URL for MDI tarball") - set(MDI_MD5 "836f5da400d8cff0f0e4435640f9454f" CACHE STRING "MD5 checksum for MDI tarball") + set(MDI_URL "https://github.com/MolSSI-MDI/MDI_Library/archive/v1.4.1.tar.gz" CACHE STRING "URL for MDI tarball") + set(MDI_MD5 "f9505fccd4c79301a619f6452dad4ad9" CACHE STRING "MD5 checksum for MDI tarball") mark_as_advanced(MDI_URL) mark_as_advanced(MDI_MD5) enable_language(C) diff --git a/cmake/Modules/Packages/ML-PACE.cmake b/cmake/Modules/Packages/ML-PACE.cmake index c647c8873a..c553809ff1 100644 --- a/cmake/Modules/Packages/ML-PACE.cmake +++ b/cmake/Modules/Packages/ML-PACE.cmake @@ -32,5 +32,6 @@ target_include_directories(pace PUBLIC ${PACE_EVALUATOR_INCLUDE_DIR} ${YAML_CPP_ target_link_libraries(pace PRIVATE yaml-cpp-pace) - -target_link_libraries(lammps PRIVATE pace) +if(CMAKE_PROJECT_NAME STREQUAL "lammps") + target_link_libraries(lammps PRIVATE pace) +endif() diff --git a/cmake/presets/mingw-cross.cmake b/cmake/presets/mingw-cross.cmake index 81db71729a..13cfee9018 100644 --- a/cmake/presets/mingw-cross.cmake +++ b/cmake/presets/mingw-cross.cmake @@ -46,8 +46,8 @@ set(WIN_PACKAGES MISC ML-HDNNP ML-IAP - ML-SNAP ML-RANN + ML-SNAP MOFFF MOLECULE MOLFILE @@ -56,6 +56,7 @@ set(WIN_PACKAGES ORIENT PERI PHONON + PLUGIN POEMS PTM QEQ diff --git a/doc/Makefile b/doc/Makefile index 07b201b07e..8841ae4825 100644 --- a/doc/Makefile +++ b/doc/Makefile @@ -13,7 +13,7 @@ VENV = $(BUILDDIR)/docenv ANCHORCHECK = $(VENV)/bin/rst_anchor_check SPHINXCONFIG = $(BUILDDIR)/utils/sphinx-config MATHJAX = $(SPHINXCONFIG)/_static/mathjax -MATHJAXTAG = 3.2.1 +MATHJAXTAG = 3.2.2 PYTHON = $(word 3,$(shell type python3)) DOXYGEN = $(word 3,$(shell type doxygen)) diff --git a/doc/lammps.1 b/doc/lammps.1 index cfeb23aadb..5f1c25867e 100644 --- a/doc/lammps.1 +++ b/doc/lammps.1 @@ -1,7 +1,7 @@ -.TH LAMMPS "1" "2 June 2022" "2022-6-2" +.TH LAMMPS "1" "23 June 2022" "2022-6-23" .SH NAME .B LAMMPS -\- Molecular Dynamics Simulator. Version 2 June 2022 +\- Molecular Dynamics Simulator. Version 23 June 2022 .SH SYNOPSIS .B lmp diff --git a/doc/src/Commands_compute.rst b/doc/src/Commands_compute.rst index 61c5e83eda..bb0c3e9c2a 100644 --- a/doc/src/Commands_compute.rst +++ b/doc/src/Commands_compute.rst @@ -138,6 +138,8 @@ KOKKOS, o = OPENMP, t = OPT. * :doc:`smd/vol ` * :doc:`snap ` * :doc:`sna/atom ` + * :doc:`sna/grid ` + * :doc:`sna/grid/local ` * :doc:`snad/atom ` * :doc:`snav/atom ` * :doc:`sph/e/atom ` diff --git a/doc/src/Developer_plugins.rst b/doc/src/Developer_plugins.rst index 9bf52801a7..36fdd010b3 100644 --- a/doc/src/Developer_plugins.rst +++ b/doc/src/Developer_plugins.rst @@ -276,10 +276,27 @@ Compilation of the plugin can be managed via both, CMake or traditional GNU makefiles. Some examples that can be used as a template are in the ``examples/plugins`` folder. The CMake script code has some small adjustments to allow building the plugins for running unit tests with -them. Another example that converts the KIM package into a plugin can be -found in the ``examples/kim/plugin`` folder. No changes to the sources -of the KIM package themselves are needed; only the plugin interface and -loader code needs to be added. This example only supports building with -CMake, but is probably a more typical example. To compile you need to -run CMake with -DLAMMPS_SOURCE_DIR=. Other +them. + +Another example that converts the KIM package into a plugin can be found +in the ``examples/kim/plugin`` folder. No changes to the sources of the +KIM package themselves are needed; only the plugin interface and loader +code needs to be added. This example only supports building with CMake, +but is probably a more typical example. To compile you need to run CMake +with -DLAMMPS_SOURCE_DIR=. Other configuration setting are identical to those for compiling LAMMPS. + +A second example for a plugin from a package is in the +``examples/PACKAGES/pace/plugin`` folder that will create a plugin from +the ML-PACE package. In this case the bulk of the code is in a static +external library that is being downloaded and compiled first and then +combined with the pair style wrapper and the plugin loader. This +example also contains a NSIS script that can be used to create an +Installer package for Windows (the mutual licensing terms of the +external library and LAMMPS conflict when distributing binaries, so the +ML-PACE package cannot be linked statically, but the LAMMPS headers +required to build the plugin are also available under a less restrictive +license). This will automatically set the required environment variable +and launching a (compatible) LAMMPS binary will load and register the +plugin and the ML-PACE package can then be used as it was linked into +LAMMPS. diff --git a/doc/src/Errors_warnings.rst b/doc/src/Errors_warnings.rst index 2d588a1b77..ab06ac523c 100644 --- a/doc/src/Errors_warnings.rst +++ b/doc/src/Errors_warnings.rst @@ -470,6 +470,12 @@ This will most likely cause errors in kinetic fluctuations. *More than one compute sna/atom* Self-explanatory. +*More than one compute sna/grid* + Self-explanatory. + +*More than one compute sna/grid/local* + Self-explanatory. + *More than one compute snad/atom* Self-explanatory. diff --git a/doc/src/Howto_mdi.rst b/doc/src/Howto_mdi.rst index 578f98e61e..14bef1a546 100644 --- a/doc/src/Howto_mdi.rst +++ b/doc/src/Howto_mdi.rst @@ -72,7 +72,7 @@ either a stand-alone code or a plugin library. As explained on the `fix mdi/qm ` command doc page, it can be used to perform *ab initio* MD simulations or energy minimizations, -or to evalute the quantum energy and forces for a series of +or to evaluate the quantum energy and forces for a series of independent systems. The examples/mdi directory has example input scripts for all of these use cases. diff --git a/doc/src/Howto_structured_data.rst b/doc/src/Howto_structured_data.rst index 4ea6c28086..18a5dfd775 100644 --- a/doc/src/Howto_structured_data.rst +++ b/doc/src/Howto_structured_data.rst @@ -184,7 +184,7 @@ frame. .. code-block:: python - import re, yaml + import yaml import pandas as pd try: @@ -193,7 +193,7 @@ frame. from yaml import SafeLoader as Loader with open("ave.yaml") as f: - ave = yaml.load(docs, Loader=Loader) + ave = yaml.load(f, Loader=Loader) keys = ave['keywords'] df = {} diff --git a/doc/src/Packages_details.rst b/doc/src/Packages_details.rst index aaac06a4a8..59922bd939 100644 --- a/doc/src/Packages_details.rst +++ b/doc/src/Packages_details.rst @@ -657,7 +657,7 @@ advection-diffusion-reaction systems. The equations of motion of these DPD extensions are integrated through a modified velocity-Verlet (MVV) algorithm. -**Author:** Zhen Li (Division of Applied Mathematics, Brown University) +**Author:** Zhen Li (Department of Mechanical Engineering, Clemson University) **Supporting info:** @@ -1479,7 +1479,7 @@ the :doc:`Build extras ` page. * lib/mdi/README * :doc:`Howto MDI ` * :doc:`mdi ` -* :doc:`fix mdi/aimd ` +* :doc:`fix mdi/qm ` * examples/PACKAGES/mdi ---------- @@ -1801,6 +1801,8 @@ computes which analyze attributes of the potential. * src/ML-SNAP: filenames -> commands * :doc:`pair_style snap ` * :doc:`compute sna/atom ` +* :doc:`compute sna/grid ` +* :doc:`compute sna/grid/local ` * :doc:`compute snad/atom ` * :doc:`compute snav/atom ` * examples/snap diff --git a/doc/src/Speed_measure.rst b/doc/src/Speed_measure.rst index 1daf49f4c4..888e8d9790 100644 --- a/doc/src/Speed_measure.rst +++ b/doc/src/Speed_measure.rst @@ -42,5 +42,4 @@ inaccurate relative timing data, because processors have to wait when communication occurs for other processors to catch up. Thus the reported times for "Communication" or "Other" may be higher than they really are, due to load-imbalance. If this is an issue, you can -uncomment the MPI_Barrier() lines in src/timer.cpp, and re-compile -LAMMPS, to obtain synchronized timings. +use the :doc:`timer sync ` command to obtain synchronized timings. diff --git a/doc/src/Tools.rst b/doc/src/Tools.rst index ef403daa84..b57c91ffee 100644 --- a/doc/src/Tools.rst +++ b/doc/src/Tools.rst @@ -95,7 +95,7 @@ Miscellaneous tools * :ref:`LAMMPS shell ` * :ref:`LAMMPS magic patterns for file(1) ` * :ref:`Offline build tool ` - * :ref:`singularity ` + * :ref:`singularity/apptainer ` * :ref:`SWIG interface ` * :ref:`vim ` @@ -1007,14 +1007,15 @@ Ivanov, at University of Iceland (ali5 at hi.is). .. _singularity_tool: -singularity tool ----------------------------------------- +singularity/apptainer tool +-------------------------- -The singularity sub-directory contains container definitions files -that can be used to build container images for building and testing -LAMMPS on specific OS variants using the `Singularity `_ -container software. Contributions for additional variants are welcome. -For more details please see the README.md file in that folder. +The singularity sub-directory contains container definitions files that +can be used to build container images for building and testing LAMMPS on +specific OS variants using the `Apptainer `_ or +`Singularity `_ container software. Contributions for +additional variants are welcome. For more details please see the +README.md file in that folder. ---------- diff --git a/doc/src/compute.rst b/doc/src/compute.rst index 6fdedbbb95..508c440e78 100644 --- a/doc/src/compute.rst +++ b/doc/src/compute.rst @@ -284,6 +284,8 @@ The individual style names on the :doc:`Commands compute ` pag * :doc:`smd/vol ` - per-particle volumes and their sum in Smooth Mach Dynamics * :doc:`snap ` - gradients of SNAP energy and forces w.r.t. linear coefficients and related quantities for fitting SNAP potentials * :doc:`sna/atom ` - bispectrum components for each atom +* :doc:`sna/grid ` - global array of bispectrum components on a regular grid +* :doc:`sna/grid/local ` - local array of bispectrum components on a regular grid * :doc:`snad/atom ` - derivative of bispectrum components for each atom * :doc:`snav/atom ` - virial contribution from bispectrum components for each atom * :doc:`sph/e/atom ` - per-atom internal energy of Smooth-Particle Hydrodynamics atoms diff --git a/doc/src/compute_ave_sphere_atom.rst b/doc/src/compute_ave_sphere_atom.rst index db04682865..48dbf377cb 100644 --- a/doc/src/compute_ave_sphere_atom.rst +++ b/doc/src/compute_ave_sphere_atom.rst @@ -35,16 +35,24 @@ Examples Description """"""""""" -Define a computation that calculates the local density and temperature -for each atom and neighbors inside a spherical cutoff. +Define a computation that calculates the local mass density and +temperature for each atom based on its neighbors inside a spherical +cutoff. If an atom has M neighbors, then its local mass density is +calculated as the sum of its mass and its M neighbor masses, divided +by the volume of the cutoff sphere (or circle in 2d). The local +temperature of the atom is calculated as the temperature of the +collection of M+1 atoms, after subtracting the center-of-mass velocity +of the M+1 atoms from each of the M+1 atom's velocities. This is +effectively the thermal velocity of the neighborhood of the central +atom, similar to :doc:`compute temp/com `. -The optional keyword *cutoff* defines the distance cutoff -used when searching for neighbors. The default value is the cutoff -specified by the pair style. If no pair style is defined, then a cutoff -must be defined using this keyword. If the specified cutoff is larger than -that of the pair_style plus neighbor skin (or no pair style is defined), -the *comm_modify cutoff* option must also be set to match that of the -*cutoff* keyword. +The optional keyword *cutoff* defines the distance cutoff used when +searching for neighbors. The default value is the cutoff specified by +the pair style. If no pair style is defined, then a cutoff must be +defined using this keyword. If the specified cutoff is larger than +that of the pair_style plus neighbor skin (or no pair style is +defined), the *comm_modify cutoff* option must also be set to match +that of the *cutoff* keyword. The neighbor list needed to compute this quantity is constructed each time the calculation is performed (i.e. each time a snapshot of atoms @@ -55,16 +63,16 @@ too frequently. If you have a bonded system, then the settings of :doc:`special_bonds ` command can remove pairwise - interactions between atoms in the same bond, angle, or dihedral. This - is the default setting for the :doc:`special_bonds ` - command, and means those pairwise interactions do not appear in the - neighbor list. Because this fix uses the neighbor list, it also means - those pairs will not be included in the order parameter. This - difficulty can be circumvented by writing a dump file, and using the - :doc:`rerun ` command to compute the order parameter for - snapshots in the dump file. The rerun script can use a - :doc:`special_bonds ` command that includes all pairs in - the neighbor list. + interactions between atoms in the same bond, angle, or dihedral. + This is the default setting for the :doc:`special_bonds + ` command, and means those pairwise interactions do + not appear in the neighbor list. Because this compute uses the + neighbor list, it also means those pairs will not be included in + the order parameter. This difficulty can be circumvented by + writing a dump file, and using the :doc:`rerun ` command to + compute the order parameter for snapshots in the dump file. The + rerun script can use a :doc:`special_bonds ` command + that includes all pairs in the neighbor list. ---------- @@ -77,17 +85,20 @@ too frequently. Output info """"""""""" -This compute calculates a per-atom array with two columns: density and temperature. +This compute calculates a per-atom array with two columns: mass +density in density :doc:`units ` and temperature in temperature +:doc:`units `. These values can be accessed by any command that uses per-atom values -from a compute as input. See the :doc:`Howto output ` doc -page for an overview of LAMMPS output options. +from a compute as input. See the :doc:`Howto output ` +doc page for an overview of LAMMPS output options. Restrictions """""""""""" -This compute is part of the EXTRA-COMPUTE package. It is only enabled if -LAMMPS was built with that package. See the :doc:`Build package ` page for more info. +This compute is part of the EXTRA-COMPUTE package. It is only enabled +if LAMMPS was built with that package. See the :doc:`Build package +` page for more info. Related commands """""""""""""""" @@ -97,5 +108,5 @@ Related commands Default """"""" -The option defaults are *cutoff* = pair style cutoff +The option defaults are *cutoff* = pair style cutoff. diff --git a/doc/src/compute_sna_atom.rst b/doc/src/compute_sna_atom.rst index 54a6df02a2..deb717d84d 100644 --- a/doc/src/compute_sna_atom.rst +++ b/doc/src/compute_sna_atom.rst @@ -2,6 +2,8 @@ .. index:: compute snad/atom .. index:: compute snav/atom .. index:: compute snap +.. index:: compute sna/grid +.. index:: compute sna/grid/local compute sna/atom command ======================== @@ -15,6 +17,12 @@ compute snav/atom command compute snap command ==================== +compute sna/grid command +======================== + +compute sna/grid/local command +============================== + Syntax """""" @@ -24,6 +32,9 @@ Syntax compute ID group-ID snad/atom rcutfac rfac0 twojmax R_1 R_2 ... w_1 w_2 ... keyword values ... compute ID group-ID snav/atom rcutfac rfac0 twojmax R_1 R_2 ... w_1 w_2 ... keyword values ... compute ID group-ID snap rcutfac rfac0 twojmax R_1 R_2 ... w_1 w_2 ... keyword values ... + compute ID group-ID snap rcutfac rfac0 twojmax R_1 R_2 ... w_1 w_2 ... keyword values ... + compute ID group-ID sna/grid nx ny nz rcutfac rfac0 twojmax R_1 R_2 ... w_1 w_2 ... keyword values ... + compute ID group-ID sna/grid/local nx ny nz rcutfac rfac0 twojmax R_1 R_2 ... w_1 w_2 ... keyword values ... * ID, group-ID are documented in :doc:`compute ` command * sna/atom = style name of this compute command @@ -32,6 +43,7 @@ Syntax * twojmax = band limit for bispectrum components (non-negative integer) * R_1, R_2,... = list of cutoff radii, one for each type (distance units) * w_1, w_2,... = list of neighbor weights, one for each type +* nx, ny, nz = number of grid points in x, y, and z directions (positive integer) * zero or more keyword/value pairs may be appended * keyword = *rmin0* or *switchflag* or *bzeroflag* or *quadraticflag* or *chem* or *bnormflag* or *wselfallflag* or *bikflag* or *switchinnerflag* or *sinner* or *dinner* @@ -78,6 +90,7 @@ Examples compute snap all snap 1.4 0.95 6 2.0 1.0 compute snap all snap 1.0 0.99363 6 3.81 3.83 1.0 0.93 chem 2 0 1 compute snap all snap 1.0 0.99363 6 3.81 3.83 1.0 0.93 switchinnerflag 1 sinner 1.35 1.6 dinner 0.25 0.3 + compute bgrid all sna/grid/local 200 200 200 1.4 0.95 6 2.0 1.0 Description """"""""""" @@ -212,6 +225,46 @@ command: See section below on output for a detailed explanation of the data layout in the global array. +The compute *sna/grid* and *sna/grid/local* commands calculate +bispectrum components for a regular grid of points. +These are calculated from the local density of nearby atoms *i'* +around each grid point, as if there was a central atom *i* +at the grid point. This is useful for characterizing fine-scale +structure in a configuration of atoms, and it is used +in the `MALA package `_ +to build machine-learning surrogates for finite-temperature Kohn-Sham +density functional theory (:ref:`Ellis et al. `) +Neighbor atoms not in the group do not contribute to the +bispectrum components of the grid points. The distance cutoff :math:`R_{ii'}` +assumes that *i* has the same type as the neighbor atom *i'*. + +Compute *sna/grid* calculates a global array containing bispectrum +components for a regular grid of points. +The grid is aligned with the current box dimensions, with the +first point at the box origin, and forming a regular 3d array with +*nx*, *ny*, and *nz* points in the x, y, and z directions. For triclinic +boxes, the array is congruent with the periodic lattice vectors +a, b, and c. The array contains one row for each of the +:math:`nx \times ny \times nz` grid points, looping over the index for *ix* fastest, +then *iy*, and *iz* slowest. Each row of the array contains the *x*, *y*, +and *z* coordinates of the grid point, followed by the bispectrum +components. See section below on output for a detailed explanation of the data +layout in the global array. + +Compute *sna/grid/local* calculates bispectrum components of a regular +grid of points similarly to compute *sna/grid* described above. +However, because the array is local, it contains only rows for grid points +that are local to the processor sub-domain. The global grid +of :math:`nx \times ny \times nz` points is still laid out in space the same as for *sna/grid*, +but grid points are strictly partitioned, so that every grid point appears in +one and only one local array. The array contains one row for each of the +local grid points, looping over the global index *ix* fastest, +then *iy*, and *iz* slowest. Each row of the array contains +the global indexes *ix*, *iy*, and *iz* first, followed by the *x*, *y*, +and *z* coordinates of the grid point, followed by the bispectrum +components. See section below on output for a detailed explanation of the data +layout in the global array. + The value of all bispectrum components will be zero for atoms not in the group. Neighbor atoms not in the group do not contribute to the bispectrum of atoms in the group. @@ -414,6 +467,21 @@ number of columns in the global array generated by *snap* are 31, and 931, respectively, while the number of rows is 1+3\*\ *N*\ +6, where *N* is the total number of atoms. +Compute *sna/grid* evaluates a global array. +The array contains one row for each of the +:math:`nx \times ny \times nz` grid points, looping over the index for *ix* fastest, +then *iy*, and *iz* slowest. Each row of the array contains the *x*, *y*, +and *z* coordinates of the grid point, followed by the bispectrum +components. + +Compute *sna/grid/local* evaluates a local array. +The array contains one row for each of the +local grid points, looping over the global index *ix* fastest, +then *iy*, and *iz* slowest. Each row of the array contains +the global indexes *ix*, *iy*, and *iz* first, followed by the *x*, *y*, +and *z* coordinates of the grid point, followed by the bispectrum +components. + If the *quadratic* keyword value is set to 1, then additional columns are generated, corresponding to the products of all distinct pairs of bispectrum components. If the number of bispectrum components is *K*, @@ -464,8 +532,7 @@ The optional keyword defaults are *rmin0* = 0, .. _Thompson20141: -**(Thompson)** Thompson, Swiler, Trott, Foiles, Tucker, under review, preprint -available at `arXiv:1409.3880 `_ +**(Thompson)** Thompson, Swiler, Trott, Foiles, Tucker, J Comp Phys, 285, 316, (2015). .. _Bartok20101: @@ -486,4 +553,8 @@ of Angular Momentum, World Scientific, Singapore (1987). .. _Cusentino2020: -**(Cusentino)** Cusentino, Wood, and Thompson, J Phys Chem A, xxx, xxxxx, (2020) +**(Cusentino)** Cusentino, Wood, Thompson, J Phys Chem A, 124, 5456, (2020) + +.. _Ellis2021: + +**(Ellis)** Ellis, Fiedler, Popoola, Modine, Stephens, Thompson, Cangi, Rajamanickam, Phys Rev B, 104, 035120, (2021) diff --git a/doc/src/fix_mdi_qm.rst b/doc/src/fix_mdi_qm.rst index a00f6e2edc..ca47db5e7a 100644 --- a/doc/src/fix_mdi_qm.rst +++ b/doc/src/fix_mdi_qm.rst @@ -57,9 +57,9 @@ explained below. These are example use cases for this fix, discussed further below: -* perform an ab intitio MD (AIMD) simulation with quantum forces -* perform an energy minimziation with quantum forces -* perform a nudged elatic band (NEB) calculation with quantum forces +* perform an ab initio MD (AIMD) simulation with quantum forces +* perform an energy minimization with quantum forces +* perform a nudged elastic band (NEB) calculation with quantum forces * perform a QM calculation for a series of independent systems which LAMMPS reads or generates The code coupling performed by this command is done via the `MDI @@ -71,7 +71,7 @@ information about how LAMMPS can operate as either an MDI driver or engine. The examples/mdi directory contains input scripts using this fix in -the various use cases dicussed below. In each case, two instances of +the various use cases discussed below. In each case, two instances of LAMMPS are used, once as an MDI driver, once as an MDI engine (surrogate for a QM code). The examples/mdi/README file explains how to launch two codes so that they communicate via the MDI library using @@ -141,7 +141,7 @@ directory. See its README file for more details. (1) To run an ab initio MD (AIMD) dynamics simulation, or an energy minimization with QM forces, or a multi-replica NEB calculation, use *add yes* and *every 1* (the defaults). This is so that every time -LAMMPS needs energy and forces, the QM code will be invoked. +LAMMPS needs energy and forces, the QM code will be invoked. Both LAMMPS and the QM code should define the same system (simulation box, atoms and their types) in their respective input scripts. Note @@ -157,7 +157,7 @@ atom coordinates. The engine code computes quantum forces on each atom and the total energy of the system and returns them to LAMMPS. Note that if the AIMD simulation is an NPT or NPH model, or the energy -minimization includesf :doc:`fix box relax ` to +minimization includes :doc:`fix box relax ` to equilibrate the box size/shape, then LAMMPS computes a pressure. This means the *virial* keyword should be set to *yes* so that the QM contribution to the pressure can be included. @@ -172,7 +172,7 @@ atom coordinates and QM forces to a file. Likewise the QM energy and virial could be output with the :doc:`thermo_style custom ` command. -(3) To do a QM evaulation of energy and forces for a series of *N* +(3) To do a QM evaluation of energy and forces for a series of *N* independent systems (simulation box and atoms), use *add no* and *every 1*. Write a LAMMPS input script which loops over the *N* systems. See the :doc:`Howto multiple ` doc page for @@ -180,7 +180,7 @@ details on looping and removing old systems. The series of systems could be initialized by reading them from data files with :doc:`read_data ` commands. Or, for example, by using the :doc:`lattice ` , :doc:`create_atoms `, -:doc:`delete_atoms `, and/or :doc:`displace_atoms +:doc:`delete_atoms `, and/or :doc:`displace_atoms random ` commands to generate a series of different systems. At the end of the loop perform :doc:`run 0 ` and :doc:`write_dump ` commands to invoke the QM code and diff --git a/doc/src/kspace_style.rst b/doc/src/kspace_style.rst index 59efbdc5bd..30f3e550c5 100644 --- a/doc/src/kspace_style.rst +++ b/doc/src/kspace_style.rst @@ -129,8 +129,8 @@ Examples kspace_style pppm 1.0e-4 kspace_style pppm/cg 1.0e-5 1.0e-6 - kspace style msm 1.0e-4 - kspace style scafacos fmm 1.0e-4 + kspace_style msm 1.0e-4 + kspace_style scafacos fmm 1.0e-4 kspace_style none Used in input scripts: diff --git a/doc/src/mdi.rst b/doc/src/mdi.rst index 2591b3ea16..b5ee0098bd 100644 --- a/doc/src/mdi.rst +++ b/doc/src/mdi.rst @@ -326,15 +326,15 @@ able to initiate and terminate the connection to the engine code. The only current MDI driver command in LAMMPS is the :doc:`fix mdi/qm ` command. If it is only used once in an input script then it can initiate and terminate the connection. But if it is being -issuedd multiple times, e.g. in a loop that issues a :doc:`clear -` command, then it cannot initiate/terminate the connection +issued multiple times, e.g. in a loop that issues a :doc:`clear +` command, then it cannot initiate or terminate the connection multiple times. Instead, the *mdi connect* and *mdi exit* commands -should be used outside the loop to intiate/terminate the connection. +should be used outside the loop to initiate or terminate the connection. See the examples/mdi/in.series.driver script for an example of how this is done. The LOOP in that script is reading a series of data file configurations and passing them to an MDI engine (e.g. quantum -code) for enery and force evaluation. A *clear* command inside the +code) for energy and force evaluation. A *clear* command inside the loop wipes out the current system so a new one can be defined. This operation also destroys all fixes. So the :doc:`fix mdi/qm ` command is issued once per loop iteration. Note that it diff --git a/doc/src/pair_e3b.rst b/doc/src/pair_e3b.rst index b75fb8450c..8ae6d10b82 100644 --- a/doc/src/pair_e3b.rst +++ b/doc/src/pair_e3b.rst @@ -50,6 +50,12 @@ Examples pair_style hybrid/overlay e3b 1 lj/cut/tip4p/long 1 2 1 1 0.15 8.5 pair_coeff * * e3b preset 2011 +Used in example input script: + +.. parsed-literal:: + + examples/PACKAGES/e3b/in.e3b-tip4p2005 + Description """"""""""" @@ -68,21 +74,27 @@ The *e3b* style computes an \"explicit three-body\" (E3B) potential for water :r 0 & r>R_f\\ \end{cases} -This potential was developed as a water model that includes the three-body cooperativity of hydrogen bonding explicitly. -To use it in this way, it must be applied in conjunction with a conventional two-body water model, through *pair_style hybrid/overlay*. -The three body interactions are split into three types: A, B, and C. -Type A corresponds to anti-cooperative double hydrogen bond donor interactions. -Type B corresponds to the cooperative interaction of molecules that both donate and accept a hydrogen bond. -Type C corresponds to anti-cooperative double hydrogen bond acceptor interactions. -The three-body interactions are smoothly cutoff by the switching function s(r) between Rs and Rc3. -The two-body interactions are designed to correct for the effective many-body interactions implicitly included in the conventional two-body potential. -The two-body interactions are cut off sharply at Rc2, because K3 is typically significantly smaller than K2. -See :ref:`(Kumar 2008) ` for more details. +This potential was developed as a water model that includes the +three-body cooperativity of hydrogen bonding explicitly. To use it in +this way, it must be applied in conjunction with a conventional two-body +water model, through pair style :doc:`hybrid/overlay `. The +three body interactions are split into three types: A, B, and C. Type A +corresponds to anti-cooperative double hydrogen bond donor interactions. +Type B corresponds to the cooperative interaction of molecules that both +donate and accept a hydrogen bond. Type C corresponds to +anti-cooperative double hydrogen bond acceptor interactions. The +three-body interactions are smoothly cutoff by the switching function +s(r) between Rs and Rc3. The two-body interactions are designed to +correct for the effective many-body interactions implicitly included in +the conventional two-body potential. The two-body interactions are cut +off sharply at Rc2, because K3 is typically significantly smaller than +K2. See :ref:`(Kumar 2008) ` for more details. -Only a single *pair_coeff* command is used with the *e3b* style. -The first two arguments must be \* \*. -The oxygen atom type for the pair style is passed as the only argument to the *pair_style* command, not in the *pair_coeff* command. -The hydrogen atom type is inferred by the ordering of the atoms. +Only a single :doc:`pair_coeff ` command is used with the +*e3b* style and the first two arguments must be \* \*. The oxygen atom +type for the pair style is passed as the only argument to the +*pair_style* command, not in the *pair_coeff* command. The hydrogen +atom type is inferred from the ordering of the atoms. .. note:: @@ -90,26 +102,41 @@ The hydrogen atom type is inferred by the ordering of the atoms. Each water molecule must have consecutive IDs with the oxygen first. This pair style does not test that this criteria is met. -The *pair_coeff* command must have at least one keyword/value pair, as described above. -The *preset* keyword sets the potential parameters to the values used in :ref:`(Tainter 2011) ` or :ref:`(Tainter 2015) `. -To use the water models defined in those references, the *e3b* style should always be used in conjunction with an *lj/cut/tip4p/long* style through *pair_style hybrid/overlay*, as demonstrated in the second example above. -The *preset 2011* option should be used with the :doc:`TIP4P water model `. -The *preset 2015* option should be used with the :doc:`TIP4P/2005 water model `. -If the *preset* keyword is used, no other keyword is needed. -Changes to the preset parameters can be made by specifying the *preset* keyword followed by the specific parameter to change, like *Ea*\ . -Note that the other keywords must come after *preset* in the pair_style command. -The *e3b* style can also be used to implement any three-body potential of the same form by specifying all the keywords except *neigh*\ : *Ea*, *Eb*, *Ec*, *E2*, *K3*, *K2*, *Rc3*, *Rc2*, *Rs*, and *bondL*\ . -The keyword *bondL* specifies the intramolecular OH bond length of the water model being used. -This is needed to include H atoms that are within the cutoff even when the attached oxygen atom is not. +The *pair_coeff* command must have at least one keyword/value pair, as +described above. The *preset* keyword sets the potential parameters to +the values used in :ref:`(Tainter 2011) ` or +:ref:`(Tainter 2015) `. To use the water models defined in +those references, the *e3b* style should always be used in conjunction +with an *lj/cut/tip4p/long* style through *pair_style hybrid/overlay*, +as demonstrated in the second example above. The *preset 2011* option +should be used with the :doc:`TIP4P water model `. The +*preset 2015* option should be used with the :doc:`TIP4P/2005 water +model `. If the *preset* keyword is used, no other keyword +is needed. Changes to the preset parameters can be made by specifying +the *preset* keyword followed by the specific parameter to change, like +*Ea*\ . Note that the other keywords must come after *preset* in the +pair_style command. The *e3b* style can also be used to implement any +three-body potential of the same form by specifying all the keywords +except *neigh*\ : *Ea*, *Eb*, *Ec*, *E2*, *K3*, *K2*, *Rc3*, *Rc2*, +*Rs*, and *bondL*\ . The keyword *bondL* specifies the intramolecular +OH bond length of the water model being used. This is needed to include +H atoms that are within the cutoff even when the attached oxygen atom is +not. -This pair style allocates arrays sized according to the number of pairwise interactions within Rc3. -To do this it needs an estimate for the number of water molecules within Rc3 of an oxygen atom. -This estimate defaults to 10 and can be changed using the *neigh* keyword, which takes an integer as an argument. -If the neigh setting is too small, the simulation will fail with the error "neigh is too small". -If the neigh setting is too large, the pair style will use more memory than necessary. +This pair style allocates arrays sized according to the number of +pairwise interactions within Rc3. To do this it needs an estimate for +the number of water molecules within Rc3 of an oxygen atom. This +estimate defaults to 10 and can be changed using the *neigh* keyword, +which takes an integer as an argument. If the neigh setting is too +small, the simulation will fail with the error "neigh is too small". If +the neigh setting is too large, the pair style will use more memory than +necessary. -This pair style tallies a breakdown of the total E3B potential energy into sub-categories, which can be accessed via the :doc:`compute pair ` command as a vector of values of length 4. -The 4 values correspond to the terms in the first equation above: the E2 term, the Ea term, the Eb term, and the Ec term. +This pair style tallies a breakdown of the total E3B potential energy +into sub-categories, which can be accessed via the :doc:`compute pair +` command as a vector of values of length 4. The 4 values +correspond to the terms in the first equation above: the E2 term, the Ea +term, the Eb term, and the Ec term. See the examples/PACKAGES/e3b directory for a complete example script. diff --git a/doc/src/pair_sw_angle_table.rst b/doc/src/pair_sw_angle_table.rst index 6431917d67..ff88711d7a 100644 --- a/doc/src/pair_sw_angle_table.rst +++ b/doc/src/pair_sw_angle_table.rst @@ -23,9 +23,9 @@ Examples Used in example input script: - .. parsed-literal:: +.. parsed-literal:: - examples/PACKAGES/manybody_table/in.spce_sw + examples/PACKAGES/manybody_table/in.spce_sw Description diff --git a/doc/src/pair_threebody_table.rst b/doc/src/pair_threebody_table.rst index 68661270c9..19c90feccd 100644 --- a/doc/src/pair_threebody_table.rst +++ b/doc/src/pair_threebody_table.rst @@ -27,10 +27,10 @@ Examples Used in example input scripts: - .. parsed-literal:: +.. parsed-literal:: - examples/PACKAGES/manybody_table/in.spce - examples/PACKAGES/manybody_table/in.spce2 + examples/PACKAGES/manybody_table/in.spce + examples/PACKAGES/manybody_table/in.spce2 Description """"""""""" diff --git a/doc/src/read_data.rst b/doc/src/read_data.rst index 68e1efc496..853d86785e 100644 --- a/doc/src/read_data.rst +++ b/doc/src/read_data.rst @@ -56,7 +56,7 @@ Examples read_data ../run7/data.polymer.gz read_data data.protein fix mycmap crossterm CMAP read_data data.water add append offset 3 1 1 1 1 shift 0.0 0.0 50.0 - read_data data.water add merge 1 group solvent + read_data data.water add merge group solvent Description """"""""""" @@ -622,6 +622,8 @@ of analysis. - atom-ID molecule-ID atom-type x y z * - charge - atom-ID atom-type q x y z + * - dielectric + - atom-ID atom-type q x y z normx normy normz area ed em epsilon curvature * - dipole - atom-ID atom-type q x y z mux muy muz * - dpd diff --git a/doc/src/read_restart.rst b/doc/src/read_restart.rst index 4d68167560..f3065a2a5c 100644 --- a/doc/src/read_restart.rst +++ b/doc/src/read_restart.rst @@ -11,7 +11,6 @@ Syntax read_restart file flag * file = name of binary restart file to read in -* flag = remap (optional) Examples """""""" @@ -19,44 +18,40 @@ Examples .. code-block:: LAMMPS read_restart save.10000 - read_restart save.10000 remap read_restart restart.* read_restart restart.*.mpiio - read_restart poly.*.% remap Description """"""""""" Read in a previously saved system configuration from a restart file. This allows continuation of a previous run. Details about what -information is stored (and not stored) in a restart file is given -below. Basically this operation will re-create the simulation box -with all its atoms and their attributes as well as some related global -settings, at the point in time it was written to the restart file by a -previous simulation. The simulation box will be partitioned into a -regular 3d grid of rectangular bricks, one per processor, based on the -number of processors in the current simulation and the settings of the +information is stored (and not stored) in a restart file is given below. +Basically this operation will re-create the simulation box with all its +atoms and their attributes as well as some related global settings, at +the point in time it was written to the restart file by a previous +simulation. The simulation box will be partitioned into a regular 3d +grid of rectangular bricks, one per processor, based on the number of +processors in the current simulation and the settings of the :doc:`processors ` command. The partitioning can later be -changed by the :doc:`balance ` or :doc:`fix balance ` commands. - -.. note:: - - Normally, restart files are written by the - :doc:`restart ` or :doc:`write_restart ` commands - so that all atoms in the restart file are inside the simulation box. - If this is not the case, the read_restart command will print an error - that atoms were "lost" when the file is read. This error should be - reported to the LAMMPS developers so the invalid writing of the - restart file can be fixed. If you still wish to use the restart file, - the optional *remap* flag can be appended to the read_restart command. - This should avoid the error, by explicitly remapping each atom back - into the simulation box, updating image flags for the atom - appropriately. +changed by the :doc:`balance ` or :doc:`fix balance +` commands. Restart files are saved in binary format to enable exact restarts, meaning that the trajectories of a restarted run will precisely match those produced by the original run had it continued on. +The binary restart file format was not designed with backward, forward, +or cross-platform compatibility in mind, so the files are only expected +to be read correctly by the same LAMMPS executable on the same platform. +Changes to the architecture, compilation settings, or LAMMPS version can +render a restart file unreadable or it may read the data incorrectly. +If you want a more portable format, you can use the data file format as +created by the :doc:`write_data ` command. Binary restart +files can also be converted into a data file from the command line by +the LAMMPS executable that wrote them using the :ref:`-restart2data +` command line flag. + Several things can prevent exact restarts due to round-off effects, in which case the trajectories in the 2 runs will slowly diverge. These include running on a different number of processors or changing @@ -65,7 +60,8 @@ certain settings such as those set by the :doc:`newton ` or these cases. Certain fixes will not restart exactly, though they should provide -statistically similar results. These include :doc:`fix shake ` and :doc:`fix langevin `. +statistically similar results. These include :doc:`fix shake +` and :doc:`fix langevin `. Certain pair styles will not restart exactly, though they should provide statistically similar results. This is because the forces @@ -81,18 +77,19 @@ produced the restart file, it could be a LAMMPS bug, so consider :doc:`reporting it ` if you think the behavior is a bug. Because restart files are binary, they may not be portable to other -machines. In this case, you can use the :doc:`-restart command-line switch ` to convert a restart file to a data file. +machines. In this case, you can use the :doc:`-restart command-line +switch ` to convert a restart file to a data file. -Similar to how restart files are written (see the -:doc:`write_restart ` and :doc:`restart ` -commands), the restart filename can contain two wild-card characters. -If a "\*" appears in the filename, the directory is searched for all -filenames that match the pattern where "\*" is replaced with a timestep -value. The file with the largest timestep value is read in. Thus, -this effectively means, read the latest restart file. It's useful if -you want your script to continue a run from where it left off. See -the :doc:`run ` command and its "upto" option for how to specify -the run command so it does not need to be changed either. +Similar to how restart files are written (see the :doc:`write_restart +` and :doc:`restart ` commands), the restart +filename can contain two wild-card characters. If a "\*" appears in the +filename, the directory is searched for all filenames that match the +pattern where "\*" is replaced with a timestep value. The file with the +largest timestep value is read in. Thus, this effectively means, read +the latest restart file. It's useful if you want your script to +continue a run from where it left off. See the :doc:`run ` command +and its "upto" option for how to specify the run command so it does not +need to be changed either. If a "%" character appears in the restart filename, LAMMPS expects a set of multiple files to exist. The :doc:`restart ` and @@ -222,17 +219,17 @@ its calculations in a consistent manner. .. note:: - There are a handful of commands which can be used before or - between runs which may require a system initialization. Examples - include the "balance", "displace_atoms", "delete_atoms", "set" (some - options), and "velocity" (some options) commands. This is because - they can migrate atoms to new processors. Thus they will also discard - unused "state" information from fixes. You will know the discard has + There are a handful of commands which can be used before or between + runs which may require a system initialization. Examples include the + "balance", "displace_atoms", "delete_atoms", "set" (some options), + and "velocity" (some options) commands. This is because they can + migrate atoms to new processors. Thus they will also discard unused + "state" information from fixes. You will know the discard has occurred because a list of discarded fixes will be printed to the screen and log file, as explained above. This means that if you wish to retain that info in a restarted run, you must re-specify the - relevant fixes and computes (which create fixes) before those commands - are used. + relevant fixes and computes (which create fixes) before those + commands are used. Some pair styles, like the :doc:`granular pair styles `, also use a fix to store "state" information that persists from timestep to @@ -245,18 +242,19 @@ LAMMPS allows bond interactions (angle, etc) to be turned off or deleted in various ways, which can affect how their info is stored in a restart file. -If bonds (angles, etc) have been turned off by the :doc:`fix shake ` or :doc:`delete_bonds ` command, -their info will be written to a restart file as if they are turned on. -This means they will need to be turned off again in a new run after -the restart file is read. +If bonds (angles, etc) have been turned off by the :doc:`fix shake +` or :doc:`delete_bonds ` command, their info +will be written to a restart file as if they are turned on. This means +they will need to be turned off again in a new run after the restart +file is read. Bonds that are broken (e.g. by a bond-breaking potential) are written to the restart file as broken bonds with a type of 0. Thus these bonds will still be broken when the restart file is read. -Bonds that have been broken by the :doc:`fix bond/break ` command have disappeared from the -system. No information about these bonds is written to the restart -file. +Bonds that have been broken by the :doc:`fix bond/break +` command have disappeared from the system. No +information about these bonds is written to the restart file. ---------- diff --git a/doc/utils/sphinx-config/false_positives.txt b/doc/utils/sphinx-config/false_positives.txt index 836dc4efa8..77a1da6643 100644 --- a/doc/utils/sphinx-config/false_positives.txt +++ b/doc/utils/sphinx-config/false_positives.txt @@ -120,6 +120,7 @@ Antonelli api Apoorva Appl +apptainer Apu arallel arccos @@ -372,6 +373,7 @@ Caltech Camilloni Camiloni Campana +Cangi Cao Capolungo Caro @@ -2368,7 +2370,11 @@ nopreliminary Nord norder Nordlund +noremap normals +normx +normy +normz Noskov noslip noticable @@ -2676,6 +2682,7 @@ polyhedra Polym polymorphism popen +Popoola Popov popstore Poresag @@ -2819,6 +2826,7 @@ radians Rafferty rahman Rahman +Rajamanickam Ralf Raman ramped @@ -3671,13 +3679,19 @@ vx Vx vxcm vxmu +vxx +vxy +vxz vy Vy vycm +vyy +vyz vz Vz vzcm vzi +vzz Waals Wadley wallstyle diff --git a/examples/COUPLE/fortran2/LAMMPS-wrapper.cpp b/examples/COUPLE/fortran2/LAMMPS-wrapper.cpp index 8b6a257755..eb6f421606 100644 --- a/examples/COUPLE/fortran2/LAMMPS-wrapper.cpp +++ b/examples/COUPLE/fortran2/LAMMPS-wrapper.cpp @@ -12,8 +12,8 @@ ------------------------------------------------------------------------- */ /* ------------------------------------------------------------------------ - Contributing author: Karl D. Hammond - University of Tennessee, Knoxville (USA), 2012 + Contributing author: Karl D. Hammond + University of Missouri (USA), 2012 ------------------------------------------------------------------------- */ /* This is set of "wrapper" functions to assist LAMMPS.F90, which itself @@ -23,6 +23,7 @@ #include #include "LAMMPS-wrapper.h" +#define LAMMPS_LIB_MPI 1 #include #include #include @@ -56,181 +57,40 @@ void lammps_error_all (void *ptr, const char *file, int line, const char *str) int lammps_extract_compute_vectorsize (void *ptr, char *id, int style) { - class LAMMPS *lmp = (class LAMMPS *) ptr; - int icompute = lmp->modify->find_compute(id); - if ( icompute < 0 ) return 0; - class Compute *compute = lmp->modify->compute[icompute]; - - if ( style == 0 ) - { - if ( !compute->vector_flag ) - return 0; - else - return compute->size_vector; - } - else if ( style == 1 ) - { - return lammps_get_natoms (ptr); - } - else if ( style == 2 ) - { - if ( !compute->local_flag ) - return 0; - else - return compute->size_local_rows; - } - else - return 0; + int *size; + size = (int *) lammps_extract_compute(ptr, id, style, LMP_SIZE_VECTOR); + if (size) return *size; + return 0; } void lammps_extract_compute_arraysize (void *ptr, char *id, int style, int *nrows, int *ncols) { - class LAMMPS *lmp = (class LAMMPS *) ptr; - int icompute = lmp->modify->find_compute(id); - if ( icompute < 0 ) - { - *nrows = 0; - *ncols = 0; - } - class Compute *compute = lmp->modify->compute[icompute]; - - if ( style == 0 ) - { - if ( !compute->array_flag ) - { - *nrows = 0; - *ncols = 0; - } - else - { - *nrows = compute->size_array_rows; - *ncols = compute->size_array_cols; - } - } - else if ( style == 1 ) - { - if ( !compute->peratom_flag ) - { - *nrows = 0; - *ncols = 0; - } - else - { - *nrows = lammps_get_natoms (ptr); - *ncols = compute->size_peratom_cols; - } - } - else if ( style == 2 ) - { - if ( !compute->local_flag ) - { - *nrows = 0; - *ncols = 0; - } - else - { - *nrows = compute->size_local_rows; - *ncols = compute->size_local_cols; - } - } - else - { - *nrows = 0; - *ncols = 0; - } - + int *tmp; + tmp = (int *) lammps_extract_compute(ptr, id, style, LMP_SIZE_ROWS); + if (tmp) *nrows = *tmp; + tmp = (int *) lammps_extract_compute(ptr, id, style, LMP_SIZE_COLS); + if (tmp) *ncols = *tmp; return; } int lammps_extract_fix_vectorsize (void *ptr, char *id, int style) { - class LAMMPS *lmp = (class LAMMPS *) ptr; - int ifix = lmp->modify->find_fix(id); - if ( ifix < 0 ) return 0; - class Fix *fix = lmp->modify->fix[ifix]; - - if ( style == 0 ) - { - if ( !fix->vector_flag ) - return 0; - else - return fix->size_vector; - } - else if ( style == 1 ) - { - return lammps_get_natoms (ptr); - } - else if ( style == 2 ) - { - if ( !fix->local_flag ) - return 0; - else - return fix->size_local_rows; - } - else - return 0; + int *size; + size = (int *) lammps_extract_fix(ptr, id, style, LMP_SIZE_VECTOR, 0, 0); + if (size) return *size; + return 0; } void lammps_extract_fix_arraysize (void *ptr, char *id, int style, int *nrows, int *ncols) { - class LAMMPS *lmp = (class LAMMPS *) ptr; - int ifix = lmp->modify->find_fix(id); - if ( ifix < 0 ) - { - *nrows = 0; - *ncols = 0; - } - class Fix *fix = lmp->modify->fix[ifix]; - - if ( style == 0 ) - { - if ( !fix->array_flag ) - { - *nrows = 0; - *ncols = 0; - } - else - { - *nrows = fix->size_array_rows; - *ncols = fix->size_array_cols; - } - } - else if ( style == 1 ) - { - if ( !fix->peratom_flag ) - { - *nrows = 0; - *ncols = 0; - } - else - { - *nrows = lammps_get_natoms (ptr); - *ncols = fix->size_peratom_cols; - } - } - else if ( style == 2 ) - { - if ( !fix->local_flag ) - { - *nrows = 0; - *ncols = 0; - } - else - { - *nrows = fix->size_local_rows; - *ncols = fix->size_local_cols; - } - } - else - { - *nrows = 0; - *ncols = 0; - } - + int *tmp; + tmp = (int *) lammps_extract_fix(ptr, id, style, LMP_SIZE_ROWS, 0, 0); + if (tmp) *nrows = *tmp; + tmp = (int *) lammps_extract_fix(ptr, id, style, LMP_SIZE_COLS, 0, 0); + if (tmp) *ncols = *tmp; return; - } /* vim: set ts=3 sts=3 expandtab: */ diff --git a/examples/COUPLE/fortran2/LAMMPS-wrapper.h b/examples/COUPLE/fortran2/LAMMPS-wrapper.h index 68e03ae05a..479af91738 100644 --- a/examples/COUPLE/fortran2/LAMMPS-wrapper.h +++ b/examples/COUPLE/fortran2/LAMMPS-wrapper.h @@ -12,8 +12,8 @@ ------------------------------------------------------------------------- */ /* ------------------------------------------------------------------------ - Contributing author: Karl D. Hammond - University of Tennessee, Knoxville (USA), 2012 + Contributing author: Karl D. Hammond + University of Missouri (USA), 2012 ------------------------------------------------------------------------- */ /* This is set of "wrapper" functions to assist LAMMPS.F90, which itself diff --git a/examples/COUPLE/fortran2/LAMMPS.F90 b/examples/COUPLE/fortran2/LAMMPS.F90 index 8b74a38721..251b56f48f 100644 --- a/examples/COUPLE/fortran2/LAMMPS.F90 +++ b/examples/COUPLE/fortran2/LAMMPS.F90 @@ -1292,7 +1292,7 @@ contains !! Wrapper functions local to this module {{{1 Cname = string2Cstring (name) Ccount = size(data) / natoms if ( Ccount /= 1 .and. Ccount /= 3 ) & - call lammps_error_all (ptr, FLERR, 'lammps_gather_atoms requires& + call lammps_error_all (ptr, FLERR, 'lammps_scatter_atoms requires& & count to be either 1 or 3') Fdata = data Cdata = C_loc (Fdata(1)) @@ -1355,7 +1355,7 @@ contains !! Wrapper functions local to this module {{{1 Cname = string2Cstring (name) Ccount = size(data) / ndata if ( Ccount /= 1 .and. Ccount /= 3 ) & - call lammps_error_all (ptr, FLERR, 'lammps_gather_atoms requires& + call lammps_error_all (ptr, FLERR, 'lammps_scatter_atoms requires& & count to be either 1 or 3') Fdata = data Cdata = C_loc (Fdata(1)) diff --git a/examples/COUPLE/fortran2/Makefile b/examples/COUPLE/fortran2/Makefile index fe66914892..f8a2126233 100644 --- a/examples/COUPLE/fortran2/Makefile +++ b/examples/COUPLE/fortran2/Makefile @@ -14,7 +14,7 @@ CXXLIB = -lstdc++ # replace with your C++ runtime libs # Flags for Fortran compiler, C++ compiler, and C preprocessor, respectively FFLAGS = -O2 -fPIC CXXFLAGS = -O2 -fPIC -CPPFLAGS = -DOMPI_SKIP_MPICXX=1 -DMPICH_SKIP_MPICXX -DLAMMPS_LIB_MPI +CPPFLAGS = -DOMPI_SKIP_MPICXX=1 -DMPICH_SKIP_MPICXX all : liblammps_fortran.a liblammps_fortran.so diff --git a/examples/COUPLE/fortran2/README b/examples/COUPLE/fortran2/README index 1710ca8ba6..0a7b862d86 100644 --- a/examples/COUPLE/fortran2/README +++ b/examples/COUPLE/fortran2/README @@ -11,9 +11,8 @@ This interface was created by Karl Hammond who you can contact with questions: Karl D. Hammond -University of Tennessee, Knoxville -karlh at ugcs.caltech.edu -karlh at utk.edu +University of Missouri +hammondkd at missouri.edu ------------------------------------- diff --git a/examples/COUPLE/fortran_dftb/LAMMPS-wrapper.cpp b/examples/COUPLE/fortran_dftb/LAMMPS-wrapper.cpp index 8b6a257755..eb6f421606 100644 --- a/examples/COUPLE/fortran_dftb/LAMMPS-wrapper.cpp +++ b/examples/COUPLE/fortran_dftb/LAMMPS-wrapper.cpp @@ -12,8 +12,8 @@ ------------------------------------------------------------------------- */ /* ------------------------------------------------------------------------ - Contributing author: Karl D. Hammond - University of Tennessee, Knoxville (USA), 2012 + Contributing author: Karl D. Hammond + University of Missouri (USA), 2012 ------------------------------------------------------------------------- */ /* This is set of "wrapper" functions to assist LAMMPS.F90, which itself @@ -23,6 +23,7 @@ #include #include "LAMMPS-wrapper.h" +#define LAMMPS_LIB_MPI 1 #include #include #include @@ -56,181 +57,40 @@ void lammps_error_all (void *ptr, const char *file, int line, const char *str) int lammps_extract_compute_vectorsize (void *ptr, char *id, int style) { - class LAMMPS *lmp = (class LAMMPS *) ptr; - int icompute = lmp->modify->find_compute(id); - if ( icompute < 0 ) return 0; - class Compute *compute = lmp->modify->compute[icompute]; - - if ( style == 0 ) - { - if ( !compute->vector_flag ) - return 0; - else - return compute->size_vector; - } - else if ( style == 1 ) - { - return lammps_get_natoms (ptr); - } - else if ( style == 2 ) - { - if ( !compute->local_flag ) - return 0; - else - return compute->size_local_rows; - } - else - return 0; + int *size; + size = (int *) lammps_extract_compute(ptr, id, style, LMP_SIZE_VECTOR); + if (size) return *size; + return 0; } void lammps_extract_compute_arraysize (void *ptr, char *id, int style, int *nrows, int *ncols) { - class LAMMPS *lmp = (class LAMMPS *) ptr; - int icompute = lmp->modify->find_compute(id); - if ( icompute < 0 ) - { - *nrows = 0; - *ncols = 0; - } - class Compute *compute = lmp->modify->compute[icompute]; - - if ( style == 0 ) - { - if ( !compute->array_flag ) - { - *nrows = 0; - *ncols = 0; - } - else - { - *nrows = compute->size_array_rows; - *ncols = compute->size_array_cols; - } - } - else if ( style == 1 ) - { - if ( !compute->peratom_flag ) - { - *nrows = 0; - *ncols = 0; - } - else - { - *nrows = lammps_get_natoms (ptr); - *ncols = compute->size_peratom_cols; - } - } - else if ( style == 2 ) - { - if ( !compute->local_flag ) - { - *nrows = 0; - *ncols = 0; - } - else - { - *nrows = compute->size_local_rows; - *ncols = compute->size_local_cols; - } - } - else - { - *nrows = 0; - *ncols = 0; - } - + int *tmp; + tmp = (int *) lammps_extract_compute(ptr, id, style, LMP_SIZE_ROWS); + if (tmp) *nrows = *tmp; + tmp = (int *) lammps_extract_compute(ptr, id, style, LMP_SIZE_COLS); + if (tmp) *ncols = *tmp; return; } int lammps_extract_fix_vectorsize (void *ptr, char *id, int style) { - class LAMMPS *lmp = (class LAMMPS *) ptr; - int ifix = lmp->modify->find_fix(id); - if ( ifix < 0 ) return 0; - class Fix *fix = lmp->modify->fix[ifix]; - - if ( style == 0 ) - { - if ( !fix->vector_flag ) - return 0; - else - return fix->size_vector; - } - else if ( style == 1 ) - { - return lammps_get_natoms (ptr); - } - else if ( style == 2 ) - { - if ( !fix->local_flag ) - return 0; - else - return fix->size_local_rows; - } - else - return 0; + int *size; + size = (int *) lammps_extract_fix(ptr, id, style, LMP_SIZE_VECTOR, 0, 0); + if (size) return *size; + return 0; } void lammps_extract_fix_arraysize (void *ptr, char *id, int style, int *nrows, int *ncols) { - class LAMMPS *lmp = (class LAMMPS *) ptr; - int ifix = lmp->modify->find_fix(id); - if ( ifix < 0 ) - { - *nrows = 0; - *ncols = 0; - } - class Fix *fix = lmp->modify->fix[ifix]; - - if ( style == 0 ) - { - if ( !fix->array_flag ) - { - *nrows = 0; - *ncols = 0; - } - else - { - *nrows = fix->size_array_rows; - *ncols = fix->size_array_cols; - } - } - else if ( style == 1 ) - { - if ( !fix->peratom_flag ) - { - *nrows = 0; - *ncols = 0; - } - else - { - *nrows = lammps_get_natoms (ptr); - *ncols = fix->size_peratom_cols; - } - } - else if ( style == 2 ) - { - if ( !fix->local_flag ) - { - *nrows = 0; - *ncols = 0; - } - else - { - *nrows = fix->size_local_rows; - *ncols = fix->size_local_cols; - } - } - else - { - *nrows = 0; - *ncols = 0; - } - + int *tmp; + tmp = (int *) lammps_extract_fix(ptr, id, style, LMP_SIZE_ROWS, 0, 0); + if (tmp) *nrows = *tmp; + tmp = (int *) lammps_extract_fix(ptr, id, style, LMP_SIZE_COLS, 0, 0); + if (tmp) *ncols = *tmp; return; - } /* vim: set ts=3 sts=3 expandtab: */ diff --git a/examples/COUPLE/fortran_dftb/LAMMPS-wrapper.h b/examples/COUPLE/fortran_dftb/LAMMPS-wrapper.h index 7d94362436..553466d138 100644 --- a/examples/COUPLE/fortran_dftb/LAMMPS-wrapper.h +++ b/examples/COUPLE/fortran_dftb/LAMMPS-wrapper.h @@ -12,8 +12,8 @@ ------------------------------------------------------------------------- */ /* ------------------------------------------------------------------------ - Contributing author: Karl D. Hammond - University of Tennessee, Knoxville (USA), 2012 + Contributing author: Karl D. Hammond + University of Missouri (USA), 2012 ------------------------------------------------------------------------- */ /* This is set of "wrapper" functions to assist LAMMPS.F90, which itself diff --git a/examples/COUPLE/fortran_dftb/LAMMPS-wrapper2.cpp b/examples/COUPLE/fortran_dftb/LAMMPS-wrapper2.cpp index 5f94363346..44d7b5bca4 100644 --- a/examples/COUPLE/fortran_dftb/LAMMPS-wrapper2.cpp +++ b/examples/COUPLE/fortran_dftb/LAMMPS-wrapper2.cpp @@ -12,8 +12,7 @@ ------------------------------------------------------------------------- */ /* ------------------------------------------------------------------------ - Contributing author: Karl D. Hammond - University of Tennessee, Knoxville (USA), 2012 + Contributing author: Nir Goldman, LLNL , 2016 ------------------------------------------------------------------------- */ /* This is set of "wrapper" functions to assist LAMMPS.F90, which itself diff --git a/examples/COUPLE/fortran_dftb/LAMMPS-wrapper2.h b/examples/COUPLE/fortran_dftb/LAMMPS-wrapper2.h index ed79015e78..d3705179ed 100644 --- a/examples/COUPLE/fortran_dftb/LAMMPS-wrapper2.h +++ b/examples/COUPLE/fortran_dftb/LAMMPS-wrapper2.h @@ -12,7 +12,7 @@ ------------------------------------------------------------------------- */ /* ------------------------------------------------------------------------ - Contributing author: Nir Goldman, ngoldman@llnl.gov, Oct. 19th, 2016 + Contributing author: Nir Goldman, LLNL , 2016 ------------------------------------------------------------------------- */ /* This is set of "wrapper" functions to assist LAMMPS.F90, which itself diff --git a/examples/COUPLE/fortran_dftb/LAMMPS.F90 b/examples/COUPLE/fortran_dftb/LAMMPS.F90 index 9b18bbfa5f..188fff9d60 100644 --- a/examples/COUPLE/fortran_dftb/LAMMPS.F90 +++ b/examples/COUPLE/fortran_dftb/LAMMPS.F90 @@ -12,8 +12,8 @@ !-------------------------------------------------------------------------- !! ------------------------------------------------------------------------ -! Contributing author: Karl D. Hammond -! University of Tennessee, Knoxville (USA), 2012 +! Contributing author: Karl D. Hammond +! University of Missouri (USA), 2012 !-------------------------------------------------------------------------- !! LAMMPS, a Fortran 2003 module containing an interface between Fortran diff --git a/examples/ELASTIC_T/BORN_MATRIX/Argon/Analytical/compute_born.py b/examples/ELASTIC_T/BORN_MATRIX/Argon/Analytical/compute_born.py index 0f3265faeb..88cc8dcca2 100644 --- a/examples/ELASTIC_T/BORN_MATRIX/Argon/Analytical/compute_born.py +++ b/examples/ELASTIC_T/BORN_MATRIX/Argon/Analytical/compute_born.py @@ -1,7 +1,6 @@ #!/usr/bin/env python3 # -*- coding: utf-8 -*- -import sys import numpy as np def reduce_Born(Cf): diff --git a/examples/ELASTIC_T/BORN_MATRIX/Argon/Numdiff/compute_born.py b/examples/ELASTIC_T/BORN_MATRIX/Argon/Numdiff/compute_born.py index 0f3265faeb..88cc8dcca2 100644 --- a/examples/ELASTIC_T/BORN_MATRIX/Argon/Numdiff/compute_born.py +++ b/examples/ELASTIC_T/BORN_MATRIX/Argon/Numdiff/compute_born.py @@ -1,7 +1,6 @@ #!/usr/bin/env python3 # -*- coding: utf-8 -*- -import sys import numpy as np def reduce_Born(Cf): diff --git a/examples/ELASTIC_T/BORN_MATRIX/Silicon/elastic_utils.py b/examples/ELASTIC_T/BORN_MATRIX/Silicon/elastic_utils.py index 0daf62b1b0..2307f3d413 100644 --- a/examples/ELASTIC_T/BORN_MATRIX/Silicon/elastic_utils.py +++ b/examples/ELASTIC_T/BORN_MATRIX/Silicon/elastic_utils.py @@ -1,5 +1,5 @@ import numpy as np -from lammps import lammps, LAMMPS_INT, LMP_STYLE_GLOBAL, LMP_VAR_EQUAL, LMP_VAR_ATOM +from lammps import lammps, LMP_VAR_EQUAL # method for rotating elastic constants diff --git a/examples/PACKAGES/pace/plugin/CMakeLists.txt b/examples/PACKAGES/pace/plugin/CMakeLists.txt new file mode 100644 index 0000000000..f4068a03c9 --- /dev/null +++ b/examples/PACKAGES/pace/plugin/CMakeLists.txt @@ -0,0 +1,36 @@ +########################################## +# CMake build system for plugin examples. +# The is meant to be used as a template for plugins that are +# distributed independent from the LAMMPS package. +########################################## + +cmake_minimum_required(VERSION 3.10) + +project(paceplugin VERSION 1.0 LANGUAGES CXX) + +set(CMAKE_MODULE_PATH ${CMAKE_CURRENT_SOURCE_DIR}) +include(CheckIncludeFileCXX) +include(LAMMPSInterfacePlugin) +include(ML-PACE) + +########################## +# building the plugins + +add_library(paceplugin MODULE paceplugin.cpp ${LAMMPS_SOURCE_DIR}/ML-PACE/pair_pace.cpp) +target_link_libraries(paceplugin PRIVATE pace) +target_link_libraries(paceplugin PRIVATE lammps) +target_include_directories(paceplugin PRIVATE ${LAMMPS_SOURCE_DIR}/ML-PACE) +set_target_properties(paceplugin PROPERTIES PREFIX "" SUFFIX ".so") + +# MacOS seems to need this +if(CMAKE_SYSTEM_NAME STREQUAL Darwin) + set_target_properties(paceplugin PROPERTIES LINK_FLAGS "-Wl,-undefined,dynamic_lookup") +elseif(CMAKE_SYSTEM_NAME STREQUAL "Windows") +# tell CMake to export all symbols to a .dll on Windows with special case for MinGW cross-compilers + set_target_properties(paceplugin PROPERTIES WINDOWS_EXPORT_ALL_SYMBOLS TRUE) + if(CMAKE_CROSSCOMPILING) + set_target_properties(paceplugin PROPERTIES LINK_FLAGS "-Wl,--export-all-symbols") + endif() +else() + set_target_properties(paceplugin PROPERTIES LINK_FLAGS "-rdynamic") +endif() diff --git a/examples/PACKAGES/pace/plugin/LAMMPSInterfacePlugin.cmake b/examples/PACKAGES/pace/plugin/LAMMPSInterfacePlugin.cmake new file mode 120000 index 0000000000..2ac6d20a54 --- /dev/null +++ b/examples/PACKAGES/pace/plugin/LAMMPSInterfacePlugin.cmake @@ -0,0 +1 @@ +../../../../cmake/Modules/LAMMPSInterfacePlugin.cmake \ No newline at end of file diff --git a/examples/PACKAGES/pace/plugin/ML-PACE.cmake b/examples/PACKAGES/pace/plugin/ML-PACE.cmake new file mode 120000 index 0000000000..3990c0d334 --- /dev/null +++ b/examples/PACKAGES/pace/plugin/ML-PACE.cmake @@ -0,0 +1 @@ +../../../../cmake/Modules/Packages/ML-PACE.cmake \ No newline at end of file diff --git a/examples/PACKAGES/pace/plugin/README.txt b/examples/PACKAGES/pace/plugin/README.txt new file mode 100644 index 0000000000..85490c3f6c --- /dev/null +++ b/examples/PACKAGES/pace/plugin/README.txt @@ -0,0 +1,2 @@ +This folder contains a loader and support files to build the ML-PACE package as plugin. +For more information please see: https://docs.lammps.org/Developer_plugins.html diff --git a/examples/PACKAGES/pace/plugin/lammps-text-logo-wide.bmp b/examples/PACKAGES/pace/plugin/lammps-text-logo-wide.bmp new file mode 100644 index 0000000000..b9ec4c35f2 Binary files /dev/null and b/examples/PACKAGES/pace/plugin/lammps-text-logo-wide.bmp differ diff --git a/examples/PACKAGES/pace/plugin/lammps.ico b/examples/PACKAGES/pace/plugin/lammps.ico new file mode 100644 index 0000000000..cce156bf79 Binary files /dev/null and b/examples/PACKAGES/pace/plugin/lammps.ico differ diff --git a/examples/PACKAGES/pace/plugin/paceplugin.cpp b/examples/PACKAGES/pace/plugin/paceplugin.cpp new file mode 100644 index 0000000000..adf1c168f9 --- /dev/null +++ b/examples/PACKAGES/pace/plugin/paceplugin.cpp @@ -0,0 +1,28 @@ + +#include "lammpsplugin.h" +#include "version.h" + +#include "pair_pace.h" + +using namespace LAMMPS_NS; + +static Pair *pair_pace_creator(LAMMPS *lmp) +{ + return new PairPACE(lmp); +} + +extern "C" void lammpsplugin_init(void *lmp, void *handle, void *regfunc) +{ + lammpsplugin_t plugin; + lammpsplugin_regfunc register_plugin = (lammpsplugin_regfunc) regfunc; + + // register pace pair style + plugin.version = LAMMPS_VERSION; + plugin.style = "pair"; + plugin.name = "pace"; + plugin.info = "PACE plugin pair style v1.0"; + plugin.author = "Axel Kohlmeyer (akohlmey@gmail.com)"; + plugin.creator.v1 = (lammpsplugin_factory1 *) &pair_pace_creator; + plugin.handle = handle; + (*register_plugin)(&plugin, lmp); +} diff --git a/examples/PACKAGES/pace/plugin/paceplugin.nsis b/examples/PACKAGES/pace/plugin/paceplugin.nsis new file mode 100644 index 0000000000..de8d1d8478 --- /dev/null +++ b/examples/PACKAGES/pace/plugin/paceplugin.nsis @@ -0,0 +1,157 @@ +#!Nsis Installer Command Script +# +# The following external defines are recognized: +# ${VERSION} = YYYYMMDD + +!include "MUI2.nsh" +!include "FileFunc.nsh" + +!define MUI_ICON "lammps.ico" +!define MUI_UNICON "lammps.ico" +!define MUI_HEADERIMAGE +!define MUI_HEADERIMAGE_BITMAP "lammps-text-logo-wide.bmp" +!define MUI_HEADERIMAGE_RIGHT + +Unicode true +XPStyle on + +!include "LogicLib.nsh" +!addplugindir "envvar/Plugins/x86-unicode" +!include "x64.nsh" + +RequestExecutionLevel user + +!macro VerifyUserIsAdmin +UserInfo::GetAccountType +pop $0 +${If} $0 != "admin" + messageBox mb_iconstop "Administrator rights required!" + setErrorLevel 740 ;ERROR_ELEVATION_REQUIRED + quit +${EndIf} +!macroend + +!define PACEPLUGIN "LAMMPS ML-PACE Plugin ${VERSION}" +OutFile "LAMMPS-ML-PACE-plugin-${VERSION}.exe" + +Name "${PACEPLUGIN}" +InstallDir "$LOCALAPPDATA\${PACEPLUGIN}" + +ShowInstDetails show +ShowUninstDetails show +SetCompressor lzma + +!define MUI_ABORTWARNING + +!insertmacro MUI_PAGE_DIRECTORY +!insertmacro MUI_PAGE_INSTFILES + +!insertmacro MUI_UNPAGE_CONFIRM +!insertmacro MUI_UNPAGE_INSTFILES + +!insertmacro MUI_LANGUAGE "English" + +function .onInit + # Determine if LAMMPS was already installed and check whether it was in 32-bit + # or 64-bit. Then look up path to uninstaller and offer to uninstall or quit + SetRegView 32 + ReadRegDWORD $0 HKCU "Software\LAMMPS-ML-PACE" "Bits" + SetRegView LastUsed + ${If} $0 == "32" + SetRegView 32 + ${ElseIf} $0 == "64" + SetRegView 64 + ${Else} + SetRegView 64 + ${EndIf} + ClearErrors + ReadRegStr $R0 HKCU "Software\Microsoft\Windows\CurrentVersion\Uninstall\LAMMPS-ML-PACE" "UninstallString" + SetRegView LastUsed + ${If} ${Errors} + DetailPrint "LAMMPS ML-PACE plugin not (yet) installed" + ${Else} + MessageBox MB_YESNO "LAMMPS ML-PACE plugin ($0 bit) is already installed. Uninstall existing version?" /SD IDYES IDNO Quit + Pop $R1 + StrCmp $R1 2 Quit +1 + Exec $R0 + Quit: + Quit + ${EndIf} + setShellVarContext all +functionEnd + +Section "${PACEPLUGIN}" SecPaceplugin + SectionIn RO + # Write LAMMPS installation bitness marker. Always use 32-bit registry view + SetRegView 32 + IntFmt $0 "0x%08X" 64 + WriteRegDWORD HKCU "Software\LAMMPS-ML-PACE" "Bits" $0 + + # Switch to "native" registry view + SetRegView 64 + SetShellVarContext current + + SetOutPath "$INSTDIR" + File lammps.ico + File paceplugin.so + + # Register Application and its uninstaller + WriteRegStr HKCU "Software\Microsoft\Windows\CurrentVersion\Uninstall\LAMMPS-ML-PACE" \ + "DisplayName" "${PACEPLUGIN}" + WriteRegStr HKCU "Software\Microsoft\Windows\CurrentVersion\Uninstall\LAMMPS-ML-PACE" \ + "Publisher" "The LAMMPS and PACE Developers" + WriteRegStr HKCU "Software\Microsoft\Windows\CurrentVersion\Uninstall\LAMMPS-ML-PACE" \ + "URLInfoAbout" "lammps.org" + WriteRegStr HKCU "Software\Microsoft\Windows\CurrentVersion\Uninstall\LAMMPS-ML-PACE" \ + "DisplayIcon" "$INSTDIR\lammps.ico" + WriteRegStr HKCU "Software\Microsoft\Windows\CurrentVersion\Uninstall\LAMMPS-ML-PACE" \ + "DisplayVersion" "${VERSION}" + WriteRegStr HKCU "Software\Microsoft\Windows\CurrentVersion\Uninstall\LAMMPS-ML-PACE" \ + "InstallLocation" "$INSTDIR" + WriteRegStr HKCU "Software\Microsoft\Windows\CurrentVersion\Uninstall\LAMMPS-ML-PACE" \ + "UninstallString" "$\"$INSTDIR\uninstall.exe$\"" + WriteRegStr HKCU "Software\Microsoft\Windows\CurrentVersion\Uninstall\LAMMPS-ML-PACE" \ + "QuietUninstallString" "$\"$INSTDIR\uninstall.exe$\" /S" + + ${GetSize} "$INSTDIR" "/S=0K" $0 $1 $2 + IntFmt $0 "0x%08X" $0 + WriteRegDWORD HKCU "Software\Microsoft\Windows\CurrentVersion\Uninstall\LAMMPS-ML-PACE" \ + "EstimatedSize" "$0" + + # update path variables + EnVar::SetHKCU + # add to LAMMPS plugin search path + EnVar::AddValue "LAMMPS_PLUGIN_PATH" "$INSTDIR" + + WriteUninstaller "$INSTDIR\Uninstall.exe" +SectionEnd + +function un.onInit + SetShellVarContext current +functionEnd + +Section "Uninstall" + # remove LAMMPS bitness/installation indicator always in 32-bit registry view + SetRegView 32 + DeleteRegKey HKCU "Software\LAMMPS-ML-PACE" + + # unregister extension, and uninstall info + SetRegView 64 + SetShellVarContext current + # unregister installation + DeleteRegKey HKCU "Software\Microsoft\Windows\CurrentVersion\Uninstall\LAMMPS-ML-PACE" + + # update path variables + EnVar::SetHKCU + # remove entry from LAMMPS plugin search path + EnVar::DeleteValue "LAMMPS_PLUGIN_PATH" "$INSTDIR" + + Delete /REBOOTOK "$INSTDIR\paceplugin.so" + Delete /REBOOTOK "$INSTDIR\Uninstall.exe" + Delete /REBOOTOK "$INSTDIR\lammps.ico" + RMDir /REBOOTOK "$INSTDIR" +SectionEnd + +# Local Variables: +# mode: sh +# End: diff --git a/examples/kim/plugin/CMakeLists.txt b/examples/kim/plugin/CMakeLists.txt index f4cb5f598d..78d117469d 100644 --- a/examples/kim/plugin/CMakeLists.txt +++ b/examples/kim/plugin/CMakeLists.txt @@ -6,46 +6,11 @@ cmake_minimum_required(VERSION 3.10) -# enforce out-of-source build -if(${CMAKE_SOURCE_DIR} STREQUAL ${CMAKE_BINARY_DIR}) - message(FATAL_ERROR "In-source builds are not allowed. You must create and use a build directory. " - "Please remove CMakeCache.txt and CMakeFiles first.") -endif() - project(kimplugin VERSION 1.0 LANGUAGES CXX) -set(LAMMPS_SOURCE_DIR "" CACHE PATH "Location of LAMMPS sources folder") -if(NOT LAMMPS_SOURCE_DIR) - message(FATAL_ERROR "Must set LAMMPS_SOURCE_DIR") -endif() - -# by default, install into $HOME/.local (not /usr/local), -# so that no root access (and sudo) is needed -if(CMAKE_INSTALL_PREFIX_INITIALIZED_TO_DEFAULT) - set(CMAKE_INSTALL_PREFIX "$ENV{HOME}/.local" CACHE PATH "Default install path" FORCE) -endif() - -# ugly hacks for MSVC which by default always reports an old C++ standard in the __cplusplus macro -# and prints lots of pointless warnings about "unsafe" functions -if(MSVC) - add_compile_options(/Zc:__cplusplus) - add_compile_options(/wd4244) - add_compile_options(/wd4267) - add_compile_definitions(_CRT_SECURE_NO_WARNINGS) -endif() - -# C++11 is required -set(CMAKE_CXX_STANDARD 11) -set(CMAKE_CXX_STANDARD_REQUIRED ON) - -# Need -restrict with Intel compilers -if(CMAKE_CXX_COMPILER_ID STREQUAL "Intel") - set(CMAKE_CXX_FLAGS "${CMAKE_CXX_FLAGS} -restrict") -endif() - set(CMAKE_MODULE_PATH ${CMAKE_CURRENT_SOURCE_DIR}) include(CheckIncludeFileCXX) -include(LAMMPSInterfaceCXX) +include(LAMMPSInterfacePlugin) ########################## # building the plugins @@ -90,9 +55,9 @@ if(CMAKE_SYSTEM_NAME STREQUAL Darwin) set_target_properties(kimplugin PROPERTIES LINK_FLAGS "-Wl,-undefined,dynamic_lookup") elseif(CMAKE_SYSTEM_NAME STREQUAL "Windows") # tell CMake to export all symbols to a .dll on Windows with special case for MinGW cross-compilers - set_target_properties(kimplugin.so PROPERTIES WINDOWS_EXPORT_ALL_SYMBOLS TRUE) + set_target_properties(kimplugin PROPERTIES WINDOWS_EXPORT_ALL_SYMBOLS TRUE) if(CMAKE_CROSSCOMPILING) - set_target_properties(kimplugin PROPERTIES LINK_FLAGS "-Wl,--export-all-symbols") + set_target_properties(kimplugin PROPERTIES LINK_FLAGS "-Wl,--export-all-symbols") endif() else() set_target_properties(kimplugin PROPERTIES LINK_FLAGS "-rdynamic") diff --git a/examples/kim/plugin/LAMMPSInterfaceCXX.cmake b/examples/kim/plugin/LAMMPSInterfaceCXX.cmake deleted file mode 100644 index dfbd77e28a..0000000000 --- a/examples/kim/plugin/LAMMPSInterfaceCXX.cmake +++ /dev/null @@ -1,88 +0,0 @@ -# Cmake script code to define the LAMMPS C++ interface -# settings required for building LAMMPS plugins - -################################################################################ -# helper function -function(validate_option name values) - string(TOLOWER ${${name}} needle_lower) - string(TOUPPER ${${name}} needle_upper) - list(FIND ${values} ${needle_lower} IDX_LOWER) - list(FIND ${values} ${needle_upper} IDX_UPPER) - if(${IDX_LOWER} LESS 0 AND ${IDX_UPPER} LESS 0) - list_to_bulletpoints(POSSIBLE_VALUE_LIST ${${values}}) - message(FATAL_ERROR "\n########################################################################\n" - "Invalid value '${${name}}' for option ${name}\n" - "\n" - "Possible values are:\n" - "${POSSIBLE_VALUE_LIST}" - "########################################################################") - endif() -endfunction(validate_option) - -################################################################################# -# LAMMPS C++ interface. We only need the header related parts. -add_library(lammps INTERFACE) -target_include_directories(lammps INTERFACE ${LAMMPS_SOURCE_DIR}) -if((CMAKE_SYSTEM_NAME STREQUAL "Windows") AND CMAKE_CROSSCOMPILING) - target_link_libraries(lammps INTERFACE ${CMAKE_BINARY_DIR}/../liblammps.dll.a) -endif() -################################################################################ -# MPI configuration -if(NOT CMAKE_CROSSCOMPILING) - set(MPI_CXX_SKIP_MPICXX TRUE) - find_package(MPI QUIET) - option(BUILD_MPI "Build MPI version" ${MPI_FOUND}) -else() - option(BUILD_MPI "Build MPI version" OFF) -endif() - -if(BUILD_MPI) - find_package(MPI REQUIRED) - option(LAMMPS_LONGLONG_TO_LONG "Workaround if your system or MPI version does not recognize 'long long' data types" OFF) - if(LAMMPS_LONGLONG_TO_LONG) - target_compile_definitions(lammps INTERFACE -DLAMMPS_LONGLONG_TO_LONG) - endif() - target_link_libraries(lammps INTERFACE MPI::MPI_CXX) -else() - target_include_directories(lammps INTERFACE "${LAMMPS_SOURCE_DIR}/STUBS") -endif() - -set(LAMMPS_SIZES "smallbig" CACHE STRING "LAMMPS integer sizes (smallsmall: all 32-bit, smallbig: 64-bit #atoms #timesteps, bigbig: also 64-bit imageint, 64-bit atom ids)") -set(LAMMPS_SIZES_VALUES smallbig bigbig smallsmall) -set_property(CACHE LAMMPS_SIZES PROPERTY STRINGS ${LAMMPS_SIZES_VALUES}) -validate_option(LAMMPS_SIZES LAMMPS_SIZES_VALUES) -string(TOUPPER ${LAMMPS_SIZES} LAMMPS_SIZES) -target_compile_definitions(lammps INTERFACE -DLAMMPS_${LAMMPS_SIZES}) - -################################################################################ -# detect if we may enable OpenMP support by default -set(BUILD_OMP_DEFAULT OFF) -find_package(OpenMP QUIET) -if(OpenMP_FOUND) - check_include_file_cxx(omp.h HAVE_OMP_H_INCLUDE) - if(HAVE_OMP_H_INCLUDE) - set(BUILD_OMP_DEFAULT ON) - endif() -endif() - -option(BUILD_OMP "Build with OpenMP support" ${BUILD_OMP_DEFAULT}) - -if(BUILD_OMP) - find_package(OpenMP REQUIRED) - check_include_file_cxx(omp.h HAVE_OMP_H_INCLUDE) - if(NOT HAVE_OMP_H_INCLUDE) - message(FATAL_ERROR "Cannot find the 'omp.h' header file required for full OpenMP support") - endif() - - if (((CMAKE_CXX_COMPILER_ID STREQUAL "GNU") AND (CMAKE_CXX_COMPILER_VERSION VERSION_GREATER_EQUAL 9.0)) OR - (CMAKE_CXX_COMPILER_ID STREQUAL "PGI") OR - ((CMAKE_CXX_COMPILER_ID STREQUAL "Clang") AND (CMAKE_CXX_COMPILER_VERSION VERSION_GREATER_EQUAL 10.0)) OR - ((CMAKE_CXX_COMPILER_ID STREQUAL "Intel") AND (CMAKE_CXX_COMPILER_VERSION VERSION_GREATER_EQUAL 19.0))) - # GCC 9.x and later plus Clang 10.x and later implement strict OpenMP 4.0 semantics for consts. - # Intel 18.0 was tested to support both, so we switch to OpenMP 4+ from 19.x onward to be safe. - target_compile_definitions(lammps INTERFACE -DLAMMPS_OMP_COMPAT=4) - else() - target_compile_definitions(lammps INTERFACE -DLAMMPS_OMP_COMPAT=3) - endif() - target_link_libraries(lammps INTERFACE OpenMP::OpenMP_CXX) -endif() diff --git a/examples/kim/plugin/LAMMPSInterfacePlugin.cmake b/examples/kim/plugin/LAMMPSInterfacePlugin.cmake new file mode 120000 index 0000000000..2b78a6ebcc --- /dev/null +++ b/examples/kim/plugin/LAMMPSInterfacePlugin.cmake @@ -0,0 +1 @@ +../../../cmake/Modules/LAMMPSInterfacePlugin.cmake \ No newline at end of file diff --git a/examples/kim/plugin/README.txt b/examples/kim/plugin/README.txt new file mode 100644 index 0000000000..eecc6cc7b4 --- /dev/null +++ b/examples/kim/plugin/README.txt @@ -0,0 +1,2 @@ +This folder contains a loader and support files to build the KIM package as plugin. +For more information please see: https://docs.lammps.org/Developer_plugins.html diff --git a/examples/mdi/aimd_driver.py b/examples/mdi/aimd_driver.py index cb01a85e42..062d680a1f 100644 --- a/examples/mdi/aimd_driver.py +++ b/examples/mdi/aimd_driver.py @@ -26,7 +26,7 @@ # -nsteps 10 # number of timesteps, default = 10 -import sys,math,random +import sys import mdi import numpy as np from mpi4py import MPI @@ -42,10 +42,9 @@ def error(txt=None): def perform_aimd(world,mm_comm,qm_comm): me = world.Get_rank() - nprocs = world.Get_size() # receive number of atoms from the MM engine - + mdi.MDI_Send_command(" 0: print("Build was successful") else: - error("Build of lib/%s/lib%s.so was NOT successful" % (lib,lib)) + error("Build of lib/%s/lib%s.a was NOT successful" % (lib,lib)) if has_extramake and not os.path.exists("Makefile.lammps"): print("lib/%s/Makefile.lammps was NOT created" % lib) diff --git a/lib/mdi/Makefile.g++ b/lib/mdi/Makefile.g++ index 15cc9c3a3b..41300f9992 100644 --- a/lib/mdi/Makefile.g++ +++ b/lib/mdi/Makefile.g++ @@ -8,7 +8,7 @@ EXTRAMAKE = Makefile.lammps.empty lib: $(OBJ) mkdir -p build - cd build; cmake -Dlibtype=SHARED -Dlanguage=CXX -D CMAKE_C_COMPILER=gcc -D CMAKE_CXX_COMPILER=g++ ../MDI_Library; make + cd build; cmake -Dlibtype=STATIC -Dlanguage=C -D CMAKE_POSITION_INDEPENDENT_CODE=yes -D CMAKE_C_COMPILER=gcc -D CMAKE_CXX_COMPILER=g++ ../MDI_Library; make @cp $(EXTRAMAKE) Makefile.lammps # ------ CLEAN ------ diff --git a/lib/mdi/Makefile.gcc b/lib/mdi/Makefile.gcc index a67209202b..023c0e20df 100644 --- a/lib/mdi/Makefile.gcc +++ b/lib/mdi/Makefile.gcc @@ -8,7 +8,7 @@ EXTRAMAKE = Makefile.lammps.empty lib: $(OBJ) mkdir -p build - cd build; cmake -Dlibtype=SHARED -D CMAKE_C_COMPILER=gcc ../MDI_Library; make + cd build; cmake -Dlibtype=STATIC -Dlanguage=C -D CMAKE_POSITION_INDEPENDENT_CODE=yes -D CMAKE_C_COMPILER=gcc ../MDI_Library; make @cp $(EXTRAMAKE) Makefile.lammps # ------ CLEAN ------ diff --git a/lib/mdi/Makefile.icc b/lib/mdi/Makefile.icc index eac6ba087f..8c5049561b 100644 --- a/lib/mdi/Makefile.icc +++ b/lib/mdi/Makefile.icc @@ -8,7 +8,7 @@ EXTRAMAKE = Makefile.lammps.empty lib: $(OBJ) mkdir -p build - cd build; cmake -Dlibtype=SHARED -D CMAKE_C_COMPILER=icc -D CMAKE_CXX_COMPILER=icc ../MDI_Library; make + cd build; cmake -Dlibtype=STATIC -Dlanguage=C -D CMAKE_POSITION_INDEPENDENT_CODE=yes -D CMAKE_C_COMPILER=icc -D CMAKE_CXX_COMPILER=icpc ../MDI_Library; make @cp $(EXTRAMAKE) Makefile.lammps # ------ CLEAN ------ diff --git a/lib/mdi/Makefile.mpi b/lib/mdi/Makefile.mpi index 0a40520c95..ce40bc9d4d 100644 --- a/lib/mdi/Makefile.mpi +++ b/lib/mdi/Makefile.mpi @@ -8,7 +8,7 @@ EXTRAMAKE = Makefile.lammps.empty lib: $(OBJ) mkdir -p build - cd build; cmake -Dlibtype=SHARED -D CMAKE_C_COMPILER=mpicc -D CMAKE_CXX_COMPILER=mpicxx ../MDI_Library; make + cd build; cmake -Dlibtype=STATIC -Dlanguage=C -D CMAKE_POSITION_INDEPENDENT_CODE=yes -D CMAKE_C_COMPILER=mpicc -D CMAKE_CXX_COMPILER=mpicxx ../MDI_Library; make @cp $(EXTRAMAKE) Makefile.lammps # ------ CLEAN ------ diff --git a/python/install.py b/python/install.py index 03b3366ba6..591e8525dc 100644 --- a/python/install.py +++ b/python/install.py @@ -11,7 +11,7 @@ independently and used to build the wheel without installing it. """ from __future__ import print_function -import sys,os,shutil,time,glob,subprocess +import sys,os,shutil,glob,subprocess from argparse import ArgumentParser parser = ArgumentParser(prog='install.py', @@ -23,6 +23,8 @@ parser.add_argument("-l", "--lib", required=True, help="path to the compiled LAMMPS shared library") parser.add_argument("-n", "--noinstall", action="store_true", default=False, help="only build a binary wheel. Don't attempt to install it") +parser.add_argument("-w", "--wheeldir", required=False, + help="path to a directory where the created wheel will be stored") args = parser.parse_args() @@ -30,7 +32,7 @@ args = parser.parse_args() if args.package: if not os.path.exists(args.package): - print( "ERROR: LAMMPS package %s does not exist" % args.package) + print("ERROR: LAMMPS package %s does not exist" % args.package) parser.print_help() sys.exit(1) else: @@ -38,12 +40,20 @@ if args.package: if args.lib: if not os.path.exists(args.lib): - print( "ERROR: LAMMPS shared library %s does not exist" % args.lib) + print("ERROR: LAMMPS shared library %s does not exist" % args.lib) parser.print_help() sys.exit(1) else: args.lib = os.path.abspath(args.lib) +if args.wheeldir: + if not os.path.exists(args.wheeldir): + print("ERROR: directory %s to store the wheel does not exist" % args.wheeldir) + parser.print_help() + sys.exit(1) + else: + args.wheeldir = os.path.abspath(args.wheeldir) + # we need to switch to the folder of the python package olddir = os.path.abspath('.') os.chdir(os.path.dirname(args.package)) @@ -80,10 +90,18 @@ os.remove(os.path.join('lammps',os.path.basename(args.lib))) # stop here if we were asked not to install the wheel we created if args.noinstall: - exit(0) + if args.wheeldir: + for wheel in glob.glob('lammps-*.whl'): + shutil.copy(wheel, args.wheeldir) + os.remove(wheel) + exit(0) # install the wheel with pip. first try to install in the default environment. # that will be a virtual environment, if active, or the system folder. +# if in a virtual environment, we must not use the python executable +# that is running this script (configured by cmake), but use "python" +# from the regular system path. The user may have changed to the virtual +# environment *after* running cmake. # recent versions of pip will automatically drop to use the user folder # in case the system folder is not writable. @@ -93,21 +111,35 @@ if args.noinstall: # must be uninstalled manually. We must not ignore this and drop # back to install into a (forced) user folder. -print("Installing wheel") +if "VIRTUAL_ENV" in os.environ: + print("Installing wheel into virtual environment") + py_exe = 'python' +else: + print("Installing wheel into system site-packages folder") + py_exe = sys.executable + for wheel in glob.glob('lammps-*.whl'): try: - txt = subprocess.check_output([sys.executable, '-m', 'pip', 'install', '--force-reinstall', wheel], stderr=subprocess.STDOUT, shell=False) + txt = subprocess.check_output([py_exe, '-m', 'pip', 'install', '--force-reinstall', wheel], stderr=subprocess.STDOUT, shell=False) print(txt.decode('UTF-8')) + if args.wheeldir: + shutil.copy(wheel, args.wheeldir) + else: + shutil.copy(wheel, olddir) + os.remove(wheel) continue except subprocess.CalledProcessError as err: errmsg = err.output.decode('UTF-8') if errmsg.find("distutils installed"): sys.exit(errmsg + "You need to uninstall the LAMMPS python module manually first.\n") try: - print('Installing wheel into standard site-packages folder failed. Trying user folder now') + print('Installing wheel into system site-packages folder failed. Trying user folder now') txt = subprocess.check_output([sys.executable, '-m', 'pip', 'install', '--user', '--force-reinstall', wheel], stderr=subprocess.STDOUT, shell=False) print(txt.decode('UTF-8')) + if args.wheeldir: + shutil.copy(wheel, args.wheeldir) + else: + shutil.copy(wheel, olddir) + os.remove(wheel) except: sys.exit('Failed to install wheel ' + wheel) - shutil.copy(wheel, olddir) - os.remove(wheel) diff --git a/python/lammps/core.py b/python/lammps/core.py index 23002b210c..930a40a4b0 100644 --- a/python/lammps/core.py +++ b/python/lammps/core.py @@ -188,20 +188,17 @@ class lammps(object): [c_void_p,POINTER(c_double),POINTER(c_double),c_double,c_double,c_double] self.lib.lammps_reset_box.restype = None - self.lib.lammps_gather_atoms.argtypes = \ - [c_void_p,c_char_p,c_int,c_int,c_void_p] + self.lib.lammps_gather_atoms.argtypes = [c_void_p,c_char_p,c_int,c_int,c_void_p] self.lib.lammps_gather_atoms.restype = None - self.lib.lammps_gather_atoms_concat.argtypes = \ - [c_void_p,c_char_p,c_int,c_int,c_void_p] + self.lib.lammps_gather_atoms_concat.argtypes = [c_void_p,c_char_p,c_int,c_int,c_void_p] self.lib.lammps_gather_atoms_concat.restype = None self.lib.lammps_gather_atoms_subset.argtypes = \ [c_void_p,c_char_p,c_int,c_int,c_int,POINTER(c_int),c_void_p] self.lib.lammps_gather_atoms_subset.restype = None - self.lib.lammps_scatter_atoms.argtypes = \ - [c_void_p,c_char_p,c_int,c_int,c_void_p] + self.lib.lammps_scatter_atoms.argtypes = [c_void_p,c_char_p,c_int,c_int,c_void_p] self.lib.lammps_scatter_atoms.restype = None self.lib.lammps_scatter_atoms_subset.argtypes = \ @@ -211,20 +208,17 @@ class lammps(object): self.lib.lammps_gather_bonds.argtypes = [c_void_p,c_void_p] self.lib.lammps_gather_bonds.restype = None - self.lib.lammps_gather.argtypes = \ - [c_void_p,c_char_p,c_int,c_int,c_void_p] + self.lib.lammps_gather.argtypes = [c_void_p,c_char_p,c_int,c_int,c_void_p] self.lib.lammps_gather.restype = None - self.lib.lammps_gather_concat.argtypes = \ - [c_void_p,c_char_p,c_int,c_int,c_void_p] + self.lib.lammps_gather_concat.argtypes = [c_void_p,c_char_p,c_int,c_int,c_void_p] self.lib.lammps_gather_concat.restype = None self.lib.lammps_gather_subset.argtypes = \ [c_void_p,c_char_p,c_int,c_int,c_int,POINTER(c_int),c_void_p] self.lib.lammps_gather_subset.restype = None - self.lib.lammps_scatter.argtypes = \ - [c_void_p,c_char_p,c_int,c_int,c_void_p] + self.lib.lammps_scatter.argtypes = [c_void_p,c_char_p,c_int,c_int,c_void_p] self.lib.lammps_scatter.restype = None self.lib.lammps_scatter_subset.argtypes = \ @@ -244,7 +238,8 @@ class lammps(object): self.lib.lammps_neighlist_num_elements.argtypes = [c_void_p, c_int] self.lib.lammps_neighlist_num_elements.restype = c_int - self.lib.lammps_neighlist_element_neighbors.argtypes = [c_void_p, c_int, c_int, POINTER(c_int), POINTER(c_int), POINTER(POINTER(c_int))] + self.lib.lammps_neighlist_element_neighbors.argtypes = \ + [c_void_p, c_int, c_int, POINTER(c_int), POINTER(c_int), POINTER(POINTER(c_int))] self.lib.lammps_neighlist_element_neighbors.restype = None self.lib.lammps_is_running.argtypes = [c_void_p] @@ -368,11 +363,9 @@ class lammps(object): if type(cmdargs[i]) is str: cmdargs[i] = cmdargs[i].encode() cargs = (c_char_p*narg)(*cmdargs) - self.lib.lammps_open.argtypes = [c_int, c_char_p*narg, \ - MPI_Comm, c_void_p] + self.lib.lammps_open.argtypes = [c_int, c_char_p*narg, MPI_Comm, c_void_p] else: - self.lib.lammps_open.argtypes = [c_int, c_char_p, \ - MPI_Comm, c_void_p] + self.lib.lammps_open.argtypes = [c_int, c_char_p, MPI_Comm, c_void_p] self.opened = 1 comm_ptr = self.MPI._addressof(comm) @@ -390,8 +383,7 @@ class lammps(object): if type(cmdargs[i]) is str: cmdargs[i] = cmdargs[i].encode() cargs = (c_char_p*narg)(*cmdargs) - self.lib.lammps_open_no_mpi.argtypes = [c_int, c_char_p*narg, \ - c_void_p] + self.lib.lammps_open_no_mpi.argtypes = [c_int, c_char_p*narg, c_void_p] self.lmp = c_void_p(self.lib.lammps_open_no_mpi(narg,cargs,None)) else: self.lib.lammps_open_no_mpi.argtypes = [c_int, c_char_p, c_void_p] @@ -963,17 +955,14 @@ class lammps(object): return ptr elif ctype == LMP_SIZE_COLS: - if cstyle == LMP_STYLE_GLOBAL \ - or cstyle == LMP_STYLE_ATOM \ - or cstyle == LMP_STYLE_LOCAL: + if cstyle == LMP_STYLE_GLOBAL or cstyle == LMP_STYLE_ATOM or cstyle == LMP_STYLE_LOCAL: self.lib.lammps_extract_compute.restype = POINTER(c_int) with ExceptionCheck(self): ptr = self.lib.lammps_extract_compute(self.lmp,cid,cstyle,ctype) return ptr[0] elif ctype == LMP_SIZE_VECTOR or ctype == LMP_SIZE_ROWS: - if cstyle == LMP_STYLE_GLOBAL \ - or cstyle == LMP_STYLE_LOCAL: + if cstyle == LMP_STYLE_GLOBAL or cstyle == LMP_STYLE_LOCAL: self.lib.lammps_extract_compute.restype = POINTER(c_int) with ExceptionCheck(self): ptr = self.lib.lammps_extract_compute(self.lmp,cid,cstyle,ctype) diff --git a/python/lammps/numpy_wrapper.py b/python/lammps/numpy_wrapper.py index 3619728081..ce0cb35e47 100644 --- a/python/lammps/numpy_wrapper.py +++ b/python/lammps/numpy_wrapper.py @@ -165,7 +165,7 @@ class numpy_wrapper: """ value = self.lmp.extract_compute(cid, cstyle, ctype) - if cstyle in (LMP_STYLE_GLOBAL, LMP_STYLE_LOCAL): + if cstyle == LMP_STYLE_GLOBAL: if ctype == LMP_TYPE_VECTOR: nrows = self.lmp.extract_compute(cid, cstyle, LMP_SIZE_VECTOR) return self.darray(value, nrows) @@ -173,6 +173,13 @@ class numpy_wrapper: nrows = self.lmp.extract_compute(cid, cstyle, LMP_SIZE_ROWS) ncols = self.lmp.extract_compute(cid, cstyle, LMP_SIZE_COLS) return self.darray(value, nrows, ncols) + elif cstyle == LMP_STYLE_LOCAL: + nrows = self.lmp.extract_compute(cid, cstyle, LMP_SIZE_ROWS) + ncols = self.lmp.extract_compute(cid, cstyle, LMP_SIZE_COLS) + if ncols == 0: + return self.darray(value, nrows) + else: + return self.darray(value, nrows, ncols) elif cstyle == LMP_STYLE_ATOM: if ctype == LMP_TYPE_VECTOR: nlocal = self.lmp.extract_global("nlocal") diff --git a/python/makewheel.py b/python/makewheel.py index 64ecbe2464..f13ad110ce 100644 --- a/python/makewheel.py +++ b/python/makewheel.py @@ -1,14 +1,26 @@ #!/usr/bin/env python -import sys,os,shutil +import sys,os,site -# find python script to activate the virtual environment and source it +base = os.path.abspath('buildwheel') if sys.platform == 'win32': - virtenv=os.path.join('buildwheel','Scripts','activate_this.py') + bin_dir=os.path.join(base,'Scripts') else: - virtenv=os.path.join('buildwheel','bin','activate_this.py') + bin_dir=os.path.join(base,'bin') -exec(open(virtenv).read(), {'__file__': virtenv}) +# prepend bin to PATH, set venv path +os.environ["PATH"] = os.pathsep.join([bin_dir] + os.environ.get("PATH", "").split(os.pathsep)) +os.environ["VIRTUAL_ENV"] = base + +# add the virtual environments libraries to the host python import mechanism +prev_length = len(sys.path) +for lib in "__LIB_FOLDERS__".split(os.pathsep): + path = os.path.realpath(os.path.join(bin_dir, lib)) + site.addsitedir(path) +sys.path[:] = sys.path[prev_length:] + sys.path[0:prev_length] + +sys.real_prefix = sys.prefix +sys.prefix = base # update pip and install all requirements to build the wheel os.system('python -m pip install --upgrade pip') diff --git a/python/setup.py b/python/setup.py index 0097e3d596..7794119930 100644 --- a/python/setup.py +++ b/python/setup.py @@ -3,7 +3,7 @@ from setuptools import setup from setuptools.dist import Distribution from sys import version_info -import os,time,shutil +import os,time LAMMPS_PYTHON_DIR = os.path.dirname(os.path.realpath(__file__)) LAMMPS_DIR = os.path.dirname(LAMMPS_PYTHON_DIR) LAMMPS_SOURCE_DIR = os.path.join(LAMMPS_DIR, 'src') @@ -24,7 +24,7 @@ def get_lammps_version(): class BinaryDistribution(Distribution): """Wrapper to enforce creating a binary package""" - def has_ext_modules(foo): + def has_ext_modules(self): return True if version_info.major >= 3: diff --git a/src/.gitignore b/src/.gitignore index de157734fc..4c1cf49b3b 100644 --- a/src/.gitignore +++ b/src/.gitignore @@ -173,12 +173,20 @@ /pair_tdpd.cpp /pair_tdpd.h +/compute_grid.cpp +/compute_grid.h +/compute_grid_local.cpp +/compute_grid_local.h /compute_sna_atom.cpp /compute_sna_atom.h /compute_snad_atom.cpp /compute_snad_atom.h /compute_snav_atom.cpp /compute_snav_atom.h +/compute_sna_grid.cpp +/compute_sna_grid.h +/compute_sna_grid_local.cpp +/compute_sna_grid_local.h /compute_snap.cpp /compute_snap.h /openmp_snap.h @@ -997,8 +1005,8 @@ /neb.h /netcdf_units.cpp /netcdf_units.h -/pair_3b_table.cpp -/pair_3b_table.h +/pair_threebody_table.cpp +/pair_threebody_table.h /pair_adp.cpp /pair_adp.h /pair_agni.cpp @@ -1293,8 +1301,8 @@ /pair_sph_taitwater_morris.h /pair_sw.cpp /pair_sw.h -/pair_sw_3b_table.cpp -/pair_sw_3b_table.h +/pair_sw_angle_table.cpp +/pair_sw_angle_table.h /pair_sw_mod.cpp /pair_sw_mod.h /pair_tersoff.cpp diff --git a/src/ATC/fix_atc.cpp b/src/ATC/fix_atc.cpp index a7c0cf5ade..356aa2596c 100644 --- a/src/ATC/fix_atc.cpp +++ b/src/ATC/fix_atc.cpp @@ -284,7 +284,7 @@ FixATC::FixATC(LAMMPS *lmp, int narg, char **arg) : Fix(lmp, narg, arg), int me = ATC::LammpsInterface::instance()->comm_rank(); string groupName(arg[1]); - int igroup = group->find(groupName.c_str()); + int igroup = group->find(groupName); int atomCount = group->count(igroup); try { diff --git a/src/DIELECTRIC/atom_vec_dielectric.cpp b/src/DIELECTRIC/atom_vec_dielectric.cpp index b15b233991..67bb9f7dc3 100644 --- a/src/DIELECTRIC/atom_vec_dielectric.cpp +++ b/src/DIELECTRIC/atom_vec_dielectric.cpp @@ -78,9 +78,9 @@ AtomVecDielectric::AtomVecDielectric(LAMMPS *_lmp) : AtomVec(_lmp) "mu", "area", "ed", "em", "epsilon", "curvature", "q_unscaled"}; fields_create = {"q", "molecule", "num_bond", "num_angle", "num_dihedral", "num_improper", "nspecial", "mu", "area", "ed", "em", "epsilon", "curvature", "q_unscaled"}; - fields_data_atom = { "id", "molecule", "type", "q", "x", "mu3", "area", "ed", "em", "epsilon", + fields_data_atom = {"id", "molecule", "type", "q", "x", "mu3", "area", "ed", "em", "epsilon", "curvature"}; - fields_data_vel = {"id v"}; + fields_data_vel = {"id", "v"}; // clang-format on setup_fields(); diff --git a/src/DPD-MESO/fix_mvv_dpd.cpp b/src/DPD-MESO/fix_mvv_dpd.cpp index 73290535b1..d96031c9f9 100644 --- a/src/DPD-MESO/fix_mvv_dpd.cpp +++ b/src/DPD-MESO/fix_mvv_dpd.cpp @@ -17,8 +17,8 @@ modified velocity-Verlet (MVV) algorithm. Setting verlet = 0.5 recovers the standard velocity-Verlet algorithm. - Contributing author: Zhen Li (Brown University) - Email: zhen_li@brown.edu + Contributing author: Zhen Li (Clemson University) + Email: zli7@clemson.edu ------------------------------------------------------------------------- */ #include "fix_mvv_dpd.h" diff --git a/src/DPD-MESO/fix_mvv_edpd.cpp b/src/DPD-MESO/fix_mvv_edpd.cpp index 856647d792..9b31bc4c5d 100644 --- a/src/DPD-MESO/fix_mvv_edpd.cpp +++ b/src/DPD-MESO/fix_mvv_edpd.cpp @@ -17,8 +17,8 @@ v and edpd_T) using the modified velocity-Verlet (MVV) algorithm. Setting verlet = 0.5 recovers the standard velocity-Verlet algorithm. - Contributing author: Zhen Li (Brown University) - Email: zhen_li@brown.edu + Contributing author: Zhen Li (Clemson University) + Email: zli7@clemson.edu Please cite the related publication: Z. Li, Y.-H. Tang, H. Lei, B. Caswell and G.E. Karniadakis. "Energy- diff --git a/src/DPD-MESO/fix_mvv_tdpd.cpp b/src/DPD-MESO/fix_mvv_tdpd.cpp index 00c6e29968..1a0dc5d520 100644 --- a/src/DPD-MESO/fix_mvv_tdpd.cpp +++ b/src/DPD-MESO/fix_mvv_tdpd.cpp @@ -17,8 +17,8 @@ v and cc) using the modified velocity-Verlet (MVV) algorithm. Setting verlet = 0.5 recovers the standard velocity-Verlet algorithm. - Contributing author: Zhen Li (Brown University) - Email: zhen_li@brown.edu + Contributing author: Zhen Li (Clemson University) + Email: zli7@clemson.edu Please cite the related publication: Z. Li, A. Yazdani, A. Tartakovsky and G.E. Karniadakis. "Transport diff --git a/src/DPD-MESO/pair_edpd.cpp b/src/DPD-MESO/pair_edpd.cpp index ddbbd05085..b05f588b7c 100644 --- a/src/DPD-MESO/pair_edpd.cpp +++ b/src/DPD-MESO/pair_edpd.cpp @@ -13,8 +13,8 @@ ------------------------------------------------------------------------- */ /* ---------------------------------------------------------------------- - Contributing author: Zhen Li (Brown University) - Email: zhen_li@brown.edu + Contributing author: Zhen Li (Clemson University) + Email: zli7@clemson.edu ------------------------------------------------------------------------- */ #include "pair_edpd.h" diff --git a/src/DPD-MESO/pair_mdpd.cpp b/src/DPD-MESO/pair_mdpd.cpp index 053d322f00..ec0a57be15 100644 --- a/src/DPD-MESO/pair_mdpd.cpp +++ b/src/DPD-MESO/pair_mdpd.cpp @@ -13,8 +13,8 @@ ------------------------------------------------------------------------- */ /* ---------------------------------------------------------------------- - Contributing author: Zhen Li (Brown University) - Email: zhen_li@brown.edu + Contributing author: Zhen Li (Clemson University) + Email: zli7@clemson.edu ------------------------------------------------------------------------- */ #include "pair_mdpd.h" diff --git a/src/DPD-MESO/pair_mdpd_rhosum.cpp b/src/DPD-MESO/pair_mdpd_rhosum.cpp index b44ca94c4e..773248a212 100644 --- a/src/DPD-MESO/pair_mdpd_rhosum.cpp +++ b/src/DPD-MESO/pair_mdpd_rhosum.cpp @@ -17,7 +17,7 @@ before the force calculation. The code uses 3D Lucy kernel, it can be modified for other kernels. - Contributing author: Zhen Li (Brown University) + Contributing author: Zhen Li (Clemson University) ------------------------------------------------------------------------- */ #include "pair_mdpd_rhosum.h" diff --git a/src/DPD-MESO/pair_tdpd.cpp b/src/DPD-MESO/pair_tdpd.cpp index 76f4b59108..39d9a151d9 100644 --- a/src/DPD-MESO/pair_tdpd.cpp +++ b/src/DPD-MESO/pair_tdpd.cpp @@ -13,8 +13,8 @@ ------------------------------------------------------------------------- */ /* ---------------------------------------------------------------------- - Contributing author: Zhen Li (Brown University) - Email: zhen_li@brown.edu + Contributing author: Zhen Li (Clemson University) + Email: zli7@clemson.edu ------------------------------------------------------------------------- */ #include "pair_tdpd.h" diff --git a/src/DPD-SMOOTH/pair_sdpd_taitwater_isothermal.cpp b/src/DPD-SMOOTH/pair_sdpd_taitwater_isothermal.cpp index f5ccb15eac..7a36cdfc5c 100644 --- a/src/DPD-SMOOTH/pair_sdpd_taitwater_isothermal.cpp +++ b/src/DPD-SMOOTH/pair_sdpd_taitwater_isothermal.cpp @@ -44,7 +44,7 @@ static const double sqrt_2_inv = std::sqrt(0.5); /* ---------------------------------------------------------------------- */ PairSDPDTaitwaterIsothermal::PairSDPDTaitwaterIsothermal (LAMMPS *lmp) -: Pair (lmp) { +: Pair (lmp), random(nullptr) { restartinfo = 0; single_enable =0; } @@ -61,6 +61,7 @@ PairSDPDTaitwaterIsothermal::~PairSDPDTaitwaterIsothermal () { memory->destroy (soundspeed); memory->destroy (B); } + delete random; } /* ---------------------------------------------------------------------- */ diff --git a/src/EXTRA-COMPUTE/compute_ave_sphere_atom.cpp b/src/EXTRA-COMPUTE/compute_ave_sphere_atom.cpp index 07803aca20..010051029a 100644 --- a/src/EXTRA-COMPUTE/compute_ave_sphere_atom.cpp +++ b/src/EXTRA-COMPUTE/compute_ave_sphere_atom.cpp @@ -11,10 +11,15 @@ See the README file in the top-level LAMMPS directory. ------------------------------------------------------------------------- */ +/* ---------------------------------------------------------------------- + Contributing author: Stan Moore (SNL) +------------------------------------------------------------------------- */ + #include "compute_ave_sphere_atom.h" #include "atom.h" #include "comm.h" +#include "domain.h" #include "error.h" #include "force.h" #include "math_const.h" @@ -98,7 +103,10 @@ void ComputeAveSphereAtom::init() } cutsq = cutoff * cutoff; - sphere_vol = 4.0 / 3.0 * MY_PI * cutsq * cutoff; + if (domain->dimension == 3) + volume = 4.0 / 3.0 * MY_PI * cutsq * cutoff; + else + volume = MY_PI * cutsq; // need an occasional full neighbor list @@ -121,7 +129,7 @@ void ComputeAveSphereAtom::compute_peratom() double xtmp, ytmp, ztmp, delx, dely, delz, rsq; int *ilist, *jlist, *numneigh, **firstneigh; int count; - double vsum[3], vavg[3], vnet[3]; + double p[3], vcom[3], vnet[3]; invoked_peratom = update->ntimestep; @@ -152,12 +160,26 @@ void ComputeAveSphereAtom::compute_peratom() double **x = atom->x; double **v = atom->v; + double *mass = atom->mass; + double *rmass = atom->rmass; + int *type = atom->type; int *mask = atom->mask; + double massone_i, massone_j, totalmass; + + double adof = domain->dimension; + double mvv2e = force->mvv2e; + double mv2d = force->mv2d; + double boltz = force->boltz; for (ii = 0; ii < inum; ii++) { i = ilist[ii]; if (mask[i] & groupbit) { + if (rmass) + massone_i = rmass[i]; + else + massone_i = mass[type[i]]; + xtmp = x[i][0]; ytmp = x[i][1]; ztmp = x[i][2]; @@ -167,13 +189,18 @@ void ComputeAveSphereAtom::compute_peratom() // i atom contribution count = 1; - vsum[0] = v[i][0]; - vsum[1] = v[i][1]; - vsum[2] = v[i][2]; + totalmass = massone_i; + p[0] = v[i][0] * massone_i; + p[1] = v[i][1] * massone_i; + p[2] = v[i][2] * massone_i; for (jj = 0; jj < jnum; jj++) { j = jlist[jj]; j &= NEIGHMASK; + if (rmass) + massone_j = rmass[j]; + else + massone_j = mass[type[j]]; delx = xtmp - x[j][0]; dely = ytmp - x[j][1]; @@ -181,42 +208,45 @@ void ComputeAveSphereAtom::compute_peratom() rsq = delx * delx + dely * dely + delz * delz; if (rsq < cutsq) { count++; - vsum[0] += v[j][0]; - vsum[1] += v[j][1]; - vsum[2] += v[j][2]; + totalmass += massone_j; + p[0] += v[j][0] * massone_j; + p[1] += v[j][1] * massone_j; + p[2] += v[j][2] * massone_j; } } - vavg[0] = vsum[0] / count; - vavg[1] = vsum[1] / count; - vavg[2] = vsum[2] / count; + vcom[0] = p[0] / totalmass; + vcom[1] = p[1] / totalmass; + vcom[2] = p[2] / totalmass; // i atom contribution - count = 1; - vnet[0] = v[i][0] - vavg[0]; - vnet[1] = v[i][1] - vavg[1]; - vnet[2] = v[i][2] - vavg[2]; - double ke_sum = vnet[0] * vnet[0] + vnet[1] * vnet[1] + vnet[2] * vnet[2]; + vnet[0] = v[i][0] - vcom[0]; + vnet[1] = v[i][1] - vcom[1]; + vnet[2] = v[i][2] - vcom[2]; + double ke_sum = massone_i * (vnet[0] * vnet[0] + vnet[1] * vnet[1] + vnet[2] * vnet[2]); for (jj = 0; jj < jnum; jj++) { j = jlist[jj]; j &= NEIGHMASK; + if (rmass) + massone_j = rmass[j]; + else + massone_j = mass[type[j]]; delx = xtmp - x[j][0]; dely = ytmp - x[j][1]; delz = ztmp - x[j][2]; rsq = delx * delx + dely * dely + delz * delz; if (rsq < cutsq) { - count++; - vnet[0] = v[j][0] - vavg[0]; - vnet[1] = v[j][1] - vavg[1]; - vnet[2] = v[j][2] - vavg[2]; - ke_sum += vnet[0] * vnet[0] + vnet[1] * vnet[1] + vnet[2] * vnet[2]; + vnet[0] = v[j][0] - vcom[0]; + vnet[1] = v[j][1] - vcom[1]; + vnet[2] = v[j][2] - vcom[2]; + ke_sum += massone_j * (vnet[0] * vnet[0] + vnet[1] * vnet[1] + vnet[2] * vnet[2]); } } - double density = count / sphere_vol; - double temp = ke_sum / 3.0 / count; + double density = mv2d * totalmass / volume; + double temp = mvv2e * ke_sum / (adof * count * boltz); result[i][0] = density; result[i][1] = temp; } diff --git a/src/EXTRA-COMPUTE/compute_ave_sphere_atom.h b/src/EXTRA-COMPUTE/compute_ave_sphere_atom.h index ffed09bae5..76350997f9 100644 --- a/src/EXTRA-COMPUTE/compute_ave_sphere_atom.h +++ b/src/EXTRA-COMPUTE/compute_ave_sphere_atom.h @@ -37,7 +37,7 @@ class ComputeAveSphereAtom : public Compute { protected: int nmax; - double cutoff, cutsq, sphere_vol; + double cutoff, cutsq, volume; class NeighList *list; double **result; diff --git a/src/EXTRA-COMPUTE/compute_stress_cartesian.cpp b/src/EXTRA-COMPUTE/compute_stress_cartesian.cpp index 25c3117b0f..2dde2397db 100644 --- a/src/EXTRA-COMPUTE/compute_stress_cartesian.cpp +++ b/src/EXTRA-COMPUTE/compute_stress_cartesian.cpp @@ -266,7 +266,7 @@ void ComputeStressCartesian::compute_array() Pair *pair = force->pair; double **cutsq = force->pair->cutsq; - double xi1, xi2, xj1, xj2; + double xi1, xi2; for (ii = 0; ii < inum; ii++) { i = ilist[ii]; @@ -305,9 +305,6 @@ void ComputeStressCartesian::compute_array() } } } - xj1 = x[j][dir1]; - xj2 = x[j][dir2]; - delx = x[j][0] - xtmp; dely = x[j][1] - ytmp; delz = x[j][2] - ztmp; @@ -318,7 +315,7 @@ void ComputeStressCartesian::compute_array() // Check if inside cut-off if (rsq >= cutsq[itype][jtype]) continue; pair->single(i, j, itype, jtype, rsq, factor_coul, factor_lj, fpair); - compute_pressure(fpair, xi1, xi2, xj1, xj2, delx, dely, delz); + compute_pressure(fpair, xi1, xi2, delx, dely, delz); } } @@ -356,8 +353,8 @@ void ComputeStressCartesian::compute_array() } } -void ComputeStressCartesian::compute_pressure(double fpair, double xi, double yi, double xj, - double yj, double delx, double dely, double delz) +void ComputeStressCartesian::compute_pressure(double fpair, double xi, double yi, double delx, + double dely, double delz) { int bin1, bin2, next_bin1, next_bin2; double la = 0.0, lb = 0.0, l_sum = 0.0; diff --git a/src/EXTRA-COMPUTE/compute_stress_cartesian.h b/src/EXTRA-COMPUTE/compute_stress_cartesian.h index 63e33a5407..d4505a1e8e 100644 --- a/src/EXTRA-COMPUTE/compute_stress_cartesian.h +++ b/src/EXTRA-COMPUTE/compute_stress_cartesian.h @@ -41,7 +41,7 @@ class ComputeStressCartesian : public Compute { double *dens, *pkxx, *pkyy, *pkzz, *pcxx, *pcyy, *pczz; double *tdens, *tpkxx, *tpkyy, *tpkzz, *tpcxx, *tpcyy, *tpczz; class NeighList *list; - void compute_pressure(double, double, double, double, double, double, double, double); + void compute_pressure(double, double, double, double, double, double); }; } // namespace LAMMPS_NS diff --git a/src/EXTRA-DUMP/dump_xtc.cpp b/src/EXTRA-DUMP/dump_xtc.cpp index 3c5be6b9be..8e0bb4a0d7 100644 --- a/src/EXTRA-DUMP/dump_xtc.cpp +++ b/src/EXTRA-DUMP/dump_xtc.cpp @@ -433,10 +433,10 @@ int xdropen(XDR *xdrs, const char *filename, const char *type) return 0; } if (*type == 'w' || *type == 'W') { - type = (char *) "w+"; + type = (char *) "wb+"; lmode = XDR_ENCODE; } else { - type = (char *) "r"; + type = (char *) "rb"; lmode = XDR_DECODE; } xdrfiles[xdrid] = fopen(filename, type); diff --git a/src/EXTRA-DUMP/dump_yaml.cpp b/src/EXTRA-DUMP/dump_yaml.cpp index d4f3208ffb..be1c9768bf 100644 --- a/src/EXTRA-DUMP/dump_yaml.cpp +++ b/src/EXTRA-DUMP/dump_yaml.cpp @@ -124,6 +124,12 @@ void DumpYAML::write_data(int n, double *mybuf) } fputs("]\n", fp); } +} + +/* ---------------------------------------------------------------------- */ + +void DumpYAML::write_footer() +{ fputs("...\n", fp); } diff --git a/src/EXTRA-DUMP/dump_yaml.h b/src/EXTRA-DUMP/dump_yaml.h index bf621dcb51..60ab894de4 100644 --- a/src/EXTRA-DUMP/dump_yaml.h +++ b/src/EXTRA-DUMP/dump_yaml.h @@ -35,6 +35,7 @@ class DumpYAML : public DumpCustom { void write() override; void write_header(bigint) override; void write_data(int, double *) override; + void write_footer() override; int modify_param(int, char **) override; }; diff --git a/src/INTERLAYER/pair_ilp_tmd.cpp b/src/INTERLAYER/pair_ilp_tmd.cpp index 15e1beede7..52119cbf12 100644 --- a/src/INTERLAYER/pair_ilp_tmd.cpp +++ b/src/INTERLAYER/pair_ilp_tmd.cpp @@ -483,7 +483,7 @@ void PairILPTMD::calc_normal() } } //############################ For the edge atoms of TMD ################################ - else if (cont > 1 && cont < Nnei) { + else if (cont < Nnei) { if (strcmp(elements[itype], "Mo") == 0 || strcmp(elements[itype], "W") == 0 || strcmp(elements[itype], "S") == 0 || strcmp(elements[itype], "Se") == 0) { // derivatives of Ni[l] respect to the cont neighbors diff --git a/src/KOKKOS/compute_ave_sphere_atom_kokkos.cpp b/src/KOKKOS/compute_ave_sphere_atom_kokkos.cpp index d2cb6682a7..7873bd0d92 100644 --- a/src/KOKKOS/compute_ave_sphere_atom_kokkos.cpp +++ b/src/KOKKOS/compute_ave_sphere_atom_kokkos.cpp @@ -11,11 +11,16 @@ See the README file in the top-level LAMMPS directory. ------------------------------------------------------------------------- */ +/* ---------------------------------------------------------------------- + Contributing author: Stan Moore (SNL) +------------------------------------------------------------------------- */ + #include "compute_ave_sphere_atom_kokkos.h" #include "atom_kokkos.h" #include "atom_masks.h" #include "comm.h" +#include "domain.h" #include "error.h" #include "force.h" #include "memory_kokkos.h" @@ -105,11 +110,19 @@ void ComputeAveSphereAtomKokkos::compute_peratom() // compute properties for each atom in group // use full neighbor list to count atoms less than cutoff - atomKK->sync(execution_space,X_MASK|V_MASK|TYPE_MASK|MASK_MASK); + atomKK->sync(execution_space,X_MASK|V_MASK|RMASS_MASK|TYPE_MASK|MASK_MASK); x = atomKK->k_x.view(); v = atomKK->k_v.view(); + rmass = atomKK->k_rmass.view(); + mass = atomKK->k_mass.view(); + type = atomKK->k_type.view(); mask = atomKK->k_mask.view(); + adof = domain->dimension; + mvv2e = force->mvv2e; + mv2d = force->mv2d; + boltz = force->boltz; + Kokkos::deep_copy(d_result,0.0); copymode = 1; @@ -125,8 +138,13 @@ template KOKKOS_INLINE_FUNCTION void ComputeAveSphereAtomKokkos::operator()(TagComputeAveSphereAtom, const int &ii) const { + double massone_i,massone_j; + const int i = d_ilist[ii]; if (mask[i] & groupbit) { + if (rmass.data()) massone_i = rmass[i]; + else massone_i = mass[type[i]]; + const X_FLOAT xtmp = x(i,0); const X_FLOAT ytmp = x(i,1); const X_FLOAT ztmp = x(i,2); @@ -135,14 +153,17 @@ void ComputeAveSphereAtomKokkos::operator()(TagComputeAveSphereAtom, // i atom contribution int count = 1; - double vsum[3]; - vsum[0] = v(i,0); - vsum[1] = v(i,1); - vsum[2] = v(i,2); + double totalmass = massone_i; + double p[3]; + p[0] = v(i,0)*massone_i; + p[1] = v(i,1)*massone_i; + p[2] = v(i,2)*massone_i; for (int jj = 0; jj < jnum; jj++) { int j = d_neighbors(i,jj); j &= NEIGHMASK; + if (rmass.data()) massone_j = rmass[j]; + else massone_j = mass[type[j]]; const F_FLOAT delx = x(j,0) - xtmp; const F_FLOAT dely = x(j,1) - ytmp; @@ -150,44 +171,45 @@ void ComputeAveSphereAtomKokkos::operator()(TagComputeAveSphereAtom, const F_FLOAT rsq = delx*delx + dely*dely + delz*delz; if (rsq < cutsq) { count++; - vsum[0] += v(j,0); - vsum[1] += v(j,1); - vsum[2] += v(j,2); + totalmass += massone_j; + p[0] += v(j,0)*massone_j; + p[1] += v(j,1)*massone_j; + p[2] += v(j,2)*massone_j; } } - double vavg[3]; - vavg[0] = vsum[0]/count; - vavg[1] = vsum[1]/count; - vavg[2] = vsum[2]/count; + double vcom[3]; + vcom[0] = p[0]/totalmass; + vcom[1] = p[1]/totalmass; + vcom[2] = p[2]/totalmass; // i atom contribution - count = 1; double vnet[3]; - vnet[0] = v(i,0) - vavg[0]; - vnet[1] = v(i,1) - vavg[1]; - vnet[2] = v(i,2) - vavg[2]; - double ke_sum = vnet[0]*vnet[0] + vnet[1]*vnet[1] + vnet[2]*vnet[2]; + vnet[0] = v(i,0) - vcom[0]; + vnet[1] = v(i,1) - vcom[1]; + vnet[2] = v(i,2) - vcom[2]; + double ke_sum = massone_i * (vnet[0]*vnet[0] + vnet[1]*vnet[1] + vnet[2]*vnet[2]); for (int jj = 0; jj < jnum; jj++) { int j = d_neighbors(i,jj); j &= NEIGHMASK; + if (rmass.data()) massone_j = rmass[j]; + else massone_j = mass[type[j]]; const F_FLOAT delx = x(j,0) - xtmp; const F_FLOAT dely = x(j,1) - ytmp; const F_FLOAT delz = x(j,2) - ztmp; const F_FLOAT rsq = delx*delx + dely*dely + delz*delz; if (rsq < cutsq) { - count++; - vnet[0] = v(j,0) - vavg[0]; - vnet[1] = v(j,1) - vavg[1]; - vnet[2] = v(j,2) - vavg[2]; - ke_sum += vnet[0]*vnet[0] + vnet[1]*vnet[1] + vnet[2]*vnet[2]; + vnet[0] = v(j,0) - vcom[0]; + vnet[1] = v(j,1) - vcom[1]; + vnet[2] = v(j,2) - vcom[2]; + ke_sum += massone_j * (vnet[0]*vnet[0] + vnet[1]*vnet[1] + vnet[2]*vnet[2]); } } - double density = count/sphere_vol; - double temp = ke_sum/3.0/count; + double density = mv2d*totalmass/volume; + double temp = mvv2e*ke_sum/(adof*count*boltz); d_result(i,0) = density; d_result(i,1) = temp; } diff --git a/src/KOKKOS/compute_ave_sphere_atom_kokkos.h b/src/KOKKOS/compute_ave_sphere_atom_kokkos.h index 75b5ca3aba..1ddf943e7d 100644 --- a/src/KOKKOS/compute_ave_sphere_atom_kokkos.h +++ b/src/KOKKOS/compute_ave_sphere_atom_kokkos.h @@ -46,13 +46,18 @@ template class ComputeAveSphereAtomKokkos : public ComputeAve void operator()(TagComputeAveSphereAtom, const int &) const; private: - typename AT::t_x_array_randomread x; - typename AT::t_v_array_randomread v; + double adof,mvv2e,mv2d,boltz; + + typename AT::t_x_array x; + typename AT::t_v_array v; + typename ArrayTypes::t_float_1d rmass; + typename ArrayTypes::t_float_1d mass; + typename ArrayTypes::t_int_1d type; typename ArrayTypes::t_int_1d mask; typename AT::t_neighbors_2d d_neighbors; - typename AT::t_int_1d_randomread d_ilist; - typename AT::t_int_1d_randomread d_numneigh; + typename AT::t_int_1d d_ilist; + typename AT::t_int_1d d_numneigh; DAT::tdual_float_2d k_result; typename AT::t_float_2d d_result; diff --git a/src/KOKKOS/pair_pace_kokkos.cpp b/src/KOKKOS/pair_pace_kokkos.cpp index 9d78b63167..4ffc971cea 100644 --- a/src/KOKKOS/pair_pace_kokkos.cpp +++ b/src/KOKKOS/pair_pace_kokkos.cpp @@ -534,8 +534,7 @@ void PairPACEKokkos::compute(int eflag_in, int vflag_in) } copymode = 1; - int newton_pair = force->newton_pair; - if (newton_pair == false) + if (!force->newton_pair) error->all(FLERR,"PairPACEKokkos requires 'newton on'"); if (recursive) diff --git a/src/KOKKOS/pair_reaxff_kokkos.cpp b/src/KOKKOS/pair_reaxff_kokkos.cpp index 7f2477adba..bb6ee0c1f1 100644 --- a/src/KOKKOS/pair_reaxff_kokkos.cpp +++ b/src/KOKKOS/pair_reaxff_kokkos.cpp @@ -2628,7 +2628,7 @@ int PairReaxFFKokkos::preprocess_angular(int i, int itype, int j_sta template template KOKKOS_INLINE_FUNCTION -int PairReaxFFKokkos::preprocess_torsion(int i, int /*itype*/, int itag, +int PairReaxFFKokkos::preprocess_torsion(int i, int /*itype*/, tagint itag, F_FLOAT xtmp, F_FLOAT ytmp, F_FLOAT ztmp, int j_start, int j_end, int location_torsion) const { // in reaxff_torsion_angles: j = i, k = j, i = k; diff --git a/src/KOKKOS/pair_reaxff_kokkos.h b/src/KOKKOS/pair_reaxff_kokkos.h index 39b323a0fe..836a2de731 100644 --- a/src/KOKKOS/pair_reaxff_kokkos.h +++ b/src/KOKKOS/pair_reaxff_kokkos.h @@ -257,7 +257,7 @@ class PairReaxFFKokkos : public PairReaxFF { // Abstraction for counting and populating torsion intermediated template KOKKOS_INLINE_FUNCTION - int preprocess_torsion(int, int, int, F_FLOAT, F_FLOAT, F_FLOAT, int, int, int) const; + int preprocess_torsion(int, int, tagint, F_FLOAT, F_FLOAT, F_FLOAT, int, int, int) const; template KOKKOS_INLINE_FUNCTION diff --git a/src/MANYBODY/pair_threebody_table.cpp b/src/MANYBODY/pair_threebody_table.cpp index 14f3c261d9..2f4bc83f5a 100644 --- a/src/MANYBODY/pair_threebody_table.cpp +++ b/src/MANYBODY/pair_threebody_table.cpp @@ -450,7 +450,7 @@ void PairThreebodyTable::read_table(Table *tb, char *file, char *keyword, bool s param_extract(tb, line); // if it is a symmetric threebody interaction, less table entries are required - if (symmetric == true) { + if (symmetric) { memory->create(tb->r12file, tb->ninput * tb->ninput * (tb->ninput + 1), "mltable:r12file"); memory->create(tb->r13file, tb->ninput * tb->ninput * (tb->ninput + 1), "mltable:r13file"); memory->create(tb->thetafile, tb->ninput * tb->ninput * (tb->ninput + 1), "mltable:thetafile"); @@ -481,7 +481,7 @@ void PairThreebodyTable::read_table(Table *tb, char *file, char *keyword, bool s int cerror = 0; reader.skip_line(); // if it is a symmetric threebody interaction, less table entries are required - if (symmetric == true) { + if (symmetric) { for (int i = 0; i < tb->ninput * tb->ninput * (tb->ninput + 1); i++) { line = reader.next_line(11); try { @@ -583,7 +583,7 @@ void PairThreebodyTable::bcast_table(Table *tb, bool symmetric) MPI_Comm_rank(world, &me); if (me > 0) { // if it is a symmetric threebody interaction, less table entries are required - if (symmetric == true) { + if (symmetric) { memory->create(tb->r12file, tb->ninput * tb->ninput * (tb->ninput + 1), "mltable:r12file"); memory->create(tb->r13file, tb->ninput * tb->ninput * (tb->ninput + 1), "mltable:r13file"); memory->create(tb->thetafile, tb->ninput * tb->ninput * (tb->ninput + 1), @@ -612,7 +612,7 @@ void PairThreebodyTable::bcast_table(Table *tb, bool symmetric) } // if it is a symmetric threebody interaction, less table entries are required - if (symmetric == true) { + if (symmetric) { MPI_Bcast(tb->r12file, tb->ninput * tb->ninput * (tb->ninput + 1), MPI_DOUBLE, 0, world); MPI_Bcast(tb->r13file, tb->ninput * tb->ninput * (tb->ninput + 1), MPI_DOUBLE, 0, world); MPI_Bcast(tb->thetafile, tb->ninput * tb->ninput * (tb->ninput + 1), MPI_DOUBLE, 0, world); @@ -697,7 +697,7 @@ void PairThreebodyTable::uf_lookup(Param *pm, double r12, double r13, double the //lookup scheme // if it is a symmetric threebody interaction, less table entries are required - if (pm->symmetric == true) { + if (pm->symmetric) { nr12 = (r12 - pm->mltable->rmin + 0.5 * dr - 0.00000001) / dr; if (r12 == (pm->mltable->rmin - 0.5 * dr)) { nr12 = 0; } nr13 = (r13 - pm->mltable->rmin + 0.5 * dr - 0.00000001) / dr; @@ -778,7 +778,7 @@ void PairThreebodyTable::threebody(Param *paramijk, double rsq1, double rsq2, do } // if the indices have been swapped, swap them back - if (swapped == true) { + if (swapped) { temp = r12; r12 = r13; r13 = temp; diff --git a/src/MC/fix_widom.cpp b/src/MC/fix_widom.cpp index 0adabe5eae..0a20e7adf3 100644 --- a/src/MC/fix_widom.cpp +++ b/src/MC/fix_widom.cpp @@ -273,7 +273,7 @@ void FixWidom::init() triclinic = domain->triclinic; - ave_widom_chemical_potential = 0; + ave_widom_chemical_potential = 0.0; if (region) volume = region_volume; else volume = domain->xprd * domain->yprd * domain->zprd; diff --git a/src/MDI/fix_mdi_qm.cpp b/src/MDI/fix_mdi_qm.cpp index 1c725a0871..db1c59ada9 100644 --- a/src/MDI/fix_mdi_qm.cpp +++ b/src/MDI/fix_mdi_qm.cpp @@ -33,7 +33,7 @@ FixMDIQM::FixMDIQM(LAMMPS *lmp, int narg, char **arg) : Fix(lmp, narg, arg) { // check requirements for LAMMPS to work with MDI as an engine - if (atom->tag_enable == 0) + if (atom->tag_enable == 0) error->all(FLERR, "Cannot use MDI engine without atom IDs"); if (atom->natoms && atom->tag_consecutive() == 0) error->all(FLERR, "MDI engine requires consecutive atom IDs"); @@ -201,7 +201,7 @@ void FixMDIQM::init() // this fix makes one-time connection to engine if (connectflag) { - + // if MDI's mdicomm not set, need to Accept_comm() with stand-alone engine // othewise are already connected to plugin engine @@ -210,24 +210,24 @@ void FixMDIQM::init() if (mdicomm == MDI_COMM_NULL) { plugin = 0; MDI_Accept_communicator(&mdicomm); - if (mdicomm == MDI_COMM_NULL) + if (mdicomm == MDI_COMM_NULL) error->all(FLERR, "MDI unable to connect to stand-alone engine"); } else { plugin = 1; int method; MDI_Get_method(&method, mdicomm); - if (method != MDI_PLUGIN) + if (method != MDI_PLUGIN) error->all(FLERR, "MDI internal error for plugin engine"); } - + // connection should have been already made by "mdi connect" command // only works for stand-alone engines } else { plugin = 0; - if (lmp->mdicomm == nullptr) + if (lmp->mdicomm == nullptr) error->all(FLERR,"Fix mdi/qm is not connected to engine via mdi connect"); int nbytes = sizeof(MDI_Comm); @@ -424,7 +424,7 @@ void FixMDIQM::reallocate() bigint nsize = atom->natoms * 3; if (nsize > MAXSMALLINT) error->all(FLERR, "Natoms too large to use with fix mdi/qm"); - + maxbuf = static_cast (atom->natoms); memory->destroy(ibuf1); memory->destroy(buf3); @@ -447,10 +447,10 @@ void FixMDIQM::send_types() // use local atomID to index into ordered ibuf1 - tagint *tag = atom->tag; + tagint *tag = atom->tag; int *type = atom->type; int nlocal = atom->nlocal; - + int index; for (int i = 0; i < nlocal; i++) { index = static_cast(tag[i]) - 1; @@ -476,10 +476,10 @@ void FixMDIQM::send_elements() // use local atomID to index into ordered ibuf1 - tagint *tag = atom->tag; + tagint *tag = atom->tag; int *type = atom->type; int nlocal = atom->nlocal; - + int index; for (int i = 0; i < nlocal; i++) { index = static_cast(tag[i]) - 1; diff --git a/src/MDI/mdi_command.cpp b/src/MDI/mdi_command.cpp index 96ba839dfc..42d8b2691a 100644 --- a/src/MDI/mdi_command.cpp +++ b/src/MDI/mdi_command.cpp @@ -26,7 +26,7 @@ using namespace LAMMPS_NS; mdi command: engine or plugin or connect or exit engine is used when LAMMPS is an MDI engine, to start listening for requests plugin is used when LAMMPS is an MDI driver to load a plugin library - connect and exit are used when LAMMPS is an MDI driver to + connect and exit are used when LAMMPS is an MDI driver to (a) connect = setup comm with a stand-alone MDI engine (b) exit = terminate comm with a stand-alone MDI engine ---------------------------------------------------------------------- */ @@ -51,7 +51,7 @@ void MDICommand::command(int narg, char **arg) if (mdicomm == MDI_COMM_NULL) { MDI_Accept_communicator(&mdicomm); - if (mdicomm == MDI_COMM_NULL) + if (mdicomm == MDI_COMM_NULL) error->all(FLERR, "MDI unable to connect to stand-alone engine"); } else error->all(FLERR, "Cannot use mdi connect with plugin engine"); diff --git a/src/MDI/mdi_engine.cpp b/src/MDI/mdi_engine.cpp index 79bfbe097e..91774659ff 100644 --- a/src/MDI/mdi_engine.cpp +++ b/src/MDI/mdi_engine.cpp @@ -101,7 +101,7 @@ MDIEngine::MDIEngine(LAMMPS *_lmp, int narg, char ** arg) : Pointers(_lmp) for (int i = 1; i < ntypes; i++) for (int j = i+1; j <= ntypes; j++) { if (elements[i] == 0 || elements[j] == 0) continue; - if (elements[i] == elements[j]) + if (elements[i] == elements[j]) error->all(FLERR,"MDI engine element cannot map to multiple types"); } } @@ -1227,7 +1227,7 @@ void MDIEngine::receive_elements() break; } } - if (itype > ntypes) + if (itype > ntypes) error->all(FLERR,"MDI element not found in element list"); } } diff --git a/src/ML-SNAP/compute_grid.cpp b/src/ML-SNAP/compute_grid.cpp new file mode 100644 index 0000000000..7a5ffe1e24 --- /dev/null +++ b/src/ML-SNAP/compute_grid.cpp @@ -0,0 +1,242 @@ +/* ---------------------------------------------------------------------- + LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator + https://www.lammps.org/ Sandia National Laboratories + Steve Plimpton, sjplimp@sandia.gov + + Copyright (2003) Sandia Corporation. Under the terms of Contract + DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains + certain rights in this software. This software is distributed under + the GNU General Public License. + + See the README file in the top-level LAMMPS directory. +------------------------------------------------------------------------- */ + +#include "compute_grid.h" + +#include "atom.h" +#include "comm.h" +#include "domain.h" +#include "error.h" +#include "force.h" +#include "memory.h" +#include "modify.h" +#include "update.h" + +#include + +using namespace LAMMPS_NS; + +/* ---------------------------------------------------------------------- */ + +ComputeGrid::ComputeGrid(LAMMPS *lmp, int narg, char **arg) : + Compute(lmp, narg, arg), grid(nullptr), gridall(nullptr), gridlocal(nullptr) +{ + if (narg < 6) error->all(FLERR, "Illegal compute grid command"); + + array_flag = 1; + size_array_cols = 0; + size_array_rows = 0; + extarray = 0; + + int iarg0 = 3; + int iarg = iarg0; + if (strcmp(arg[iarg], "grid") == 0) { + if (iarg + 4 > narg) error->all(FLERR, "Illegal compute grid command"); + nx = utils::inumeric(FLERR, arg[iarg + 1], false, lmp); + ny = utils::inumeric(FLERR, arg[iarg + 2], false, lmp); + nz = utils::inumeric(FLERR, arg[iarg + 3], false, lmp); + if (nx <= 0 || ny <= 0 || nz <= 0) error->all(FLERR, "All grid dimensions must be positive"); + iarg += 4; + } else + error->all(FLERR, "Illegal compute grid command"); + + nargbase = iarg - iarg0; + + size_array_rows = nx * ny * nz; + size_array_cols_base = 3; + gridlocal_allocated = 0; +} + +/* ---------------------------------------------------------------------- */ + +ComputeGrid::~ComputeGrid() +{ + deallocate(); +} + +/* ---------------------------------------------------------------------- */ + +void ComputeGrid::setup() +{ + deallocate(); + set_grid_global(); + set_grid_local(); + allocate(); +} + +/* ---------------------------------------------------------------------- + convert global array index to box coords +------------------------------------------------------------------------- */ + +void ComputeGrid::grid2x(int igrid, double *x) +{ + int iz = igrid / (nx * ny); + igrid -= iz * (nx * ny); + int iy = igrid / nx; + igrid -= iy * nx; + int ix = igrid; + + x[0] = ix * delx; + x[1] = iy * dely; + x[2] = iz * delz; + + if (triclinic) domain->lamda2x(x, x); +} + +/* ---------------------------------------------------------------------- + copy coords to global array +------------------------------------------------------------------------- */ + +void ComputeGrid::assign_coords_all() +{ + double x[3]; + for (int igrid = 0; igrid < size_array_rows; igrid++) { + grid2x(igrid, x); + gridall[igrid][0] = x[0]; + gridall[igrid][1] = x[1]; + gridall[igrid][2] = x[2]; + } +} + +/* ---------------------------------------------------------------------- + create arrays +------------------------------------------------------------------------- */ + +void ComputeGrid::allocate() +{ + // allocate arrays + + memory->create(grid, size_array_rows, size_array_cols, "grid:grid"); + memory->create(gridall, size_array_rows, size_array_cols, "grid:gridall"); + if (nxlo <= nxhi && nylo <= nyhi && nzlo <= nzhi) { + gridlocal_allocated = 1; + memory->create4d_offset(gridlocal, size_array_cols, nzlo, nzhi, nylo, nyhi, nxlo, nxhi, + "grid:gridlocal"); + } + array = gridall; +} + +/* ---------------------------------------------------------------------- + free arrays +------------------------------------------------------------------------- */ + +void ComputeGrid::deallocate() +{ + memory->destroy(grid); + memory->destroy(gridall); + if (gridlocal_allocated) { + gridlocal_allocated = 0; + memory->destroy4d_offset(gridlocal, nzlo, nylo, nxlo); + } + array = nullptr; +} + +/* ---------------------------------------------------------------------- + set global grid +------------------------------------------------------------------------- */ + +void ComputeGrid::set_grid_global() +{ + // calculate grid layout + + triclinic = domain->triclinic; + + if (triclinic == 0) { + prd = domain->prd; + boxlo = domain->boxlo; + sublo = domain->sublo; + subhi = domain->subhi; + } else { + prd = domain->prd_lamda; + boxlo = domain->boxlo_lamda; + sublo = domain->sublo_lamda; + subhi = domain->subhi_lamda; + } + + double xprd = prd[0]; + double yprd = prd[1]; + double zprd = prd[2]; + + delxinv = nx / xprd; + delyinv = ny / yprd; + delzinv = nz / zprd; + + delx = 1.0 / delxinv; + dely = 1.0 / delyinv; + delz = 1.0 / delzinv; +} + +/* ---------------------------------------------------------------------- + set local subset of grid that I own + n xyz lo/hi = 3d brick that I own (inclusive) +------------------------------------------------------------------------- */ + +void ComputeGrid::set_grid_local() +{ + // nx,ny,nz = extent of global grid + // indices into the global grid range from 0 to N-1 in each dim + // if grid point is inside my sub-domain I own it, + // this includes sub-domain lo boundary but excludes hi boundary + // ixyz lo/hi = inclusive lo/hi bounds of global grid sub-brick I own + // if proc owns no grid cells in a dim, then ilo > ihi + // if 2 procs share a boundary a grid point is exactly on, + // the 2 equality if tests insure a consistent decision + // as to which proc owns it + + double xfraclo, xfrachi, yfraclo, yfrachi, zfraclo, zfrachi; + + if (comm->layout != Comm::LAYOUT_TILED) { + xfraclo = comm->xsplit[comm->myloc[0]]; + xfrachi = comm->xsplit[comm->myloc[0] + 1]; + yfraclo = comm->ysplit[comm->myloc[1]]; + yfrachi = comm->ysplit[comm->myloc[1] + 1]; + zfraclo = comm->zsplit[comm->myloc[2]]; + zfrachi = comm->zsplit[comm->myloc[2] + 1]; + } else { + xfraclo = comm->mysplit[0][0]; + xfrachi = comm->mysplit[0][1]; + yfraclo = comm->mysplit[1][0]; + yfrachi = comm->mysplit[1][1]; + zfraclo = comm->mysplit[2][0]; + zfrachi = comm->mysplit[2][1]; + } + + nxlo = static_cast(xfraclo * nx); + if (1.0 * nxlo != xfraclo * nx) nxlo++; + nxhi = static_cast(xfrachi * nx); + if (1.0 * nxhi == xfrachi * nx) nxhi--; + + nylo = static_cast(yfraclo * ny); + if (1.0 * nylo != yfraclo * ny) nylo++; + nyhi = static_cast(yfrachi * ny); + if (1.0 * nyhi == yfrachi * ny) nyhi--; + + nzlo = static_cast(zfraclo * nz); + if (1.0 * nzlo != zfraclo * nz) nzlo++; + nzhi = static_cast(zfrachi * nz); + if (1.0 * nzhi == zfrachi * nz) nzhi--; + + ngridlocal = (nxhi - nxlo + 1) * (nyhi - nylo + 1) * (nzhi - nzlo + 1); +} + +/* ---------------------------------------------------------------------- + memory usage of local data +------------------------------------------------------------------------- */ + +double ComputeGrid::memory_usage() +{ + double nbytes = size_array_rows * size_array_cols * sizeof(double); // grid + nbytes += size_array_rows * size_array_cols * sizeof(double); // gridall + nbytes += size_array_cols * ngridlocal * sizeof(double); // gridlocal + return nbytes; +} diff --git a/src/ML-SNAP/compute_grid.h b/src/ML-SNAP/compute_grid.h new file mode 100644 index 0000000000..84b3bc4e98 --- /dev/null +++ b/src/ML-SNAP/compute_grid.h @@ -0,0 +1,58 @@ +/* -*- c++ -*- ---------------------------------------------------------- + LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator + https://www.lammps.org/ Sandia National Laboratories + Steve Plimpton, sjplimp@sandia.gov + + Copyright (2003) Sandia Corporation. Under the terms of Contract + DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains + certain rights in this software. This software is distributed under + the GNU General Public License. + + See the README file in the top-level LAMMPS directory. +------------------------------------------------------------------------- */ + +#ifndef LMP_COMPUTE_GRID_H +#define LMP_COMPUTE_GRID_H + +#include "compute.h" + +namespace LAMMPS_NS { + +class ComputeGrid : public Compute { + public: + ComputeGrid(class LAMMPS *, int, char **); + ~ComputeGrid() override; + void setup() override; + void compute_array() override = 0; + + double memory_usage() override; + + protected: + int nx, ny, nz; // global grid dimensions + int nxlo, nxhi, nylo, nyhi, nzlo, nzhi; // local grid bounds, inclusive + int ngridlocal; // number of local grid points + int nvalues; // number of values per grid point + double **grid; // global grid + double **gridall; // global grid summed over procs + double ****gridlocal; // local grid + int triclinic; // triclinic flag + double *boxlo, *prd; // box info (units real/ortho or reduced/tri) + double *sublo, *subhi; // subdomain info (units real/ortho or reduced/tri) + double delxinv, delyinv, delzinv; // inverse grid spacing + double delx, dely, delz; // grid spacing + int nargbase; // number of base class args + double cutmax; // largest cutoff distance + int size_array_cols_base; // number of columns used for coords, etc. + int gridlocal_allocated; // shows if gridlocal allocated + + void allocate(); // create arrays + void deallocate(); // free arrays + void grid2x(int, double *); // convert grid point to coord + void assign_coords_all(); // assign coords for global grid + void set_grid_global(); // set global grid + void set_grid_local(); // set bounds for local grid +}; + +} // namespace LAMMPS_NS + +#endif diff --git a/src/ML-SNAP/compute_grid_local.cpp b/src/ML-SNAP/compute_grid_local.cpp new file mode 100644 index 0000000000..48714962ac --- /dev/null +++ b/src/ML-SNAP/compute_grid_local.cpp @@ -0,0 +1,270 @@ +/* ---------------------------------------------------------------------- + LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator + https://www.lammps.org/ Sandia National Laboratories + Steve Plimpton, sjplimp@sandia.gov + + Copyright (2003) Sandia Corporation. Under the terms of Contract + DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains + certain rights in this software. This software is distributed under + the GNU General Public License. + + See the README file in the top-level LAMMPS directory. +------------------------------------------------------------------------- */ + +#include "compute_grid_local.h" + +#include "atom.h" +#include "comm.h" +#include "domain.h" +#include "error.h" +#include "force.h" +#include "memory.h" +#include "modify.h" +#include "update.h" + +#include + +// For the subdomain test below; grid-points and subdomain boundaries +// sometimes differ by minimal amounts (in the order of 2e-17). +static constexpr double EPSILON = 1.0e-10; + +using namespace LAMMPS_NS; + +/* ---------------------------------------------------------------------- */ + +ComputeGridLocal::ComputeGridLocal(LAMMPS *lmp, int narg, char **arg) : + Compute(lmp, narg, arg), alocal(nullptr) +{ + if (narg < 6) error->all(FLERR, "Illegal compute grid/local command"); + + local_flag = 1; + size_local_cols = 0; + size_local_rows = 0; + extarray = 0; + + int iarg0 = 3; + int iarg = iarg0; + if (strcmp(arg[iarg], "grid") == 0) { + if (iarg + 4 > narg) error->all(FLERR, "Illegal compute grid/local command"); + nx = utils::inumeric(FLERR, arg[iarg + 1], false, lmp); + ny = utils::inumeric(FLERR, arg[iarg + 2], false, lmp); + nz = utils::inumeric(FLERR, arg[iarg + 3], false, lmp); + if (nx <= 0 || ny <= 0 || nz <= 0) + error->all(FLERR, "All grid/local dimensions must be positive"); + iarg += 4; + } else + error->all(FLERR, "Illegal compute grid/local command"); + + nargbase = iarg - iarg0; + + size_local_cols_base = 6; + gridlocal_allocated = 0; +} + +/* ---------------------------------------------------------------------- */ + +ComputeGridLocal::~ComputeGridLocal() +{ + deallocate(); +} + +/* ---------------------------------------------------------------------- */ + +void ComputeGridLocal::setup() +{ + deallocate(); + set_grid_global(); + set_grid_local(); + allocate(); + assign_coords(); +} + +/* ---------------------------------------------------------------------- + convert global array indexes to box coords +------------------------------------------------------------------------- */ + +void ComputeGridLocal::grid2x(int ix, int iy, int iz, double *x) +{ + x[0] = ix * delx; + x[1] = iy * dely; + x[2] = iz * delz; + + if (triclinic) domain->lamda2x(x, x); +} + +/* ---------------------------------------------------------------------- + convert global array indexes to lamda coords; for orthorombic + cells defaults to grid2x. +------------------------------------------------------------------------- */ + +void ComputeGridLocal::grid2lamda(int ix, int iy, int iz, double *x) +{ + x[0] = ix * delx; + x[1] = iy * dely; + x[2] = iz * delz; +} + +/* ---------------------------------------------------------------------- + create arrays +------------------------------------------------------------------------- */ + +void ComputeGridLocal::allocate() +{ + if (nxlo <= nxhi && nylo <= nyhi && nzlo <= nzhi) { + gridlocal_allocated = 1; + memory->create(alocal, size_local_rows, size_local_cols, "compute/grid/local:alocal"); + array_local = alocal; + } +} + +/* ---------------------------------------------------------------------- + free arrays +------------------------------------------------------------------------- */ + +void ComputeGridLocal::deallocate() +{ + if (gridlocal_allocated) { + gridlocal_allocated = 0; + memory->destroy(alocal); + } + array_local = nullptr; +} + +/* ---------------------------------------------------------------------- + set global grid +------------------------------------------------------------------------- */ + +void ComputeGridLocal::set_grid_global() +{ + // calculate grid layout + + triclinic = domain->triclinic; + + if (triclinic == 0) { + prd = domain->prd; + boxlo = domain->boxlo; + sublo = domain->sublo; + subhi = domain->subhi; + } else { + prd = domain->prd_lamda; + boxlo = domain->boxlo_lamda; + sublo = domain->sublo_lamda; + subhi = domain->subhi_lamda; + } + + double xprd = prd[0]; + double yprd = prd[1]; + double zprd = prd[2]; + + delxinv = nx / xprd; + delyinv = ny / yprd; + delzinv = nz / zprd; + + delx = 1.0 / delxinv; + dely = 1.0 / delyinv; + delz = 1.0 / delzinv; +} + +/* ---------------------------------------------------------------------- + set local subset of grid that I own + n xyz lo/hi = 3d brick that I own (inclusive) +------------------------------------------------------------------------- */ + +void ComputeGridLocal::set_grid_local() +{ + // nx,ny,nz = extent of global grid + // indices into the global grid range from 0 to N-1 in each dim + // if grid point is inside my sub-domain I own it, + // this includes sub-domain lo boundary but excludes hi boundary + // ixyz lo/hi = inclusive lo/hi bounds of global grid sub-brick I own + // if proc owns no grid cells in a dim, then ilo > ihi + // if 2 procs share a boundary a grid point is exactly on, + // the 2 equality if tests insure a consistent decision + // as to which proc owns it + + double xfraclo, xfrachi, yfraclo, yfrachi, zfraclo, zfrachi; + + if (comm->layout != Comm::LAYOUT_TILED) { + xfraclo = comm->xsplit[comm->myloc[0]]; + xfrachi = comm->xsplit[comm->myloc[0] + 1]; + yfraclo = comm->ysplit[comm->myloc[1]]; + yfrachi = comm->ysplit[comm->myloc[1] + 1]; + zfraclo = comm->zsplit[comm->myloc[2]]; + zfrachi = comm->zsplit[comm->myloc[2] + 1]; + } else { + xfraclo = comm->mysplit[0][0]; + xfrachi = comm->mysplit[0][1]; + yfraclo = comm->mysplit[1][0]; + yfrachi = comm->mysplit[1][1]; + zfraclo = comm->mysplit[2][0]; + zfrachi = comm->mysplit[2][1]; + } + + nxlo = static_cast(xfraclo * nx); + if (1.0 * nxlo != xfraclo * nx) nxlo++; + nxhi = static_cast(xfrachi * nx); + if (1.0 * nxhi == xfrachi * nx) nxhi--; + + nylo = static_cast(yfraclo * ny); + if (1.0 * nylo != yfraclo * ny) nylo++; + nyhi = static_cast(yfrachi * ny); + if (1.0 * nyhi == yfrachi * ny) nyhi--; + + nzlo = static_cast(zfraclo * nz); + if (1.0 * nzlo != zfraclo * nz) nzlo++; + nzhi = static_cast(zfrachi * nz); + if (1.0 * nzhi == zfrachi * nz) nzhi--; + + size_local_rows = (nxhi - nxlo + 1) * (nyhi - nylo + 1) * (nzhi - nzlo + 1); +} + +/* ---------------------------------------------------------------------- + copy coords to local array +------------------------------------------------------------------------- */ + +void ComputeGridLocal::assign_coords() +{ + int igrid = 0; + for (int iz = nzlo; iz <= nzhi; iz++) + for (int iy = nylo; iy <= nyhi; iy++) + for (int ix = nxlo; ix <= nxhi; ix++) { + alocal[igrid][0] = ix; + alocal[igrid][1] = iy; + alocal[igrid][2] = iz; + double xgrid[3]; + + // for triclinic: create gridpoint in lamda coordinates and transform after check. + // for orthorombic: create gridpoint in box coordinates. + + if (triclinic) + grid2lamda(ix, iy, iz, xgrid); + else + grid2x(ix, iy, iz, xgrid); + + // ensure gridpoint is not strictly outside subdomain + + if ((sublo[0] - xgrid[0]) > EPSILON || (xgrid[0] - subhi[0]) > EPSILON || + (sublo[1] - xgrid[1]) > EPSILON || (xgrid[1] - subhi[1]) > EPSILON || + (sublo[2] - xgrid[2]) > EPSILON || (xgrid[2] - subhi[2]) > EPSILON) + error->one(FLERR, "Invalid gridpoint position in compute grid/local"); + + // convert lamda to x, y, z, after sudomain check + + if (triclinic) domain->lamda2x(xgrid, xgrid); + + alocal[igrid][3] = xgrid[0]; + alocal[igrid][4] = xgrid[1]; + alocal[igrid][5] = xgrid[2]; + igrid++; + } +} + +/* ---------------------------------------------------------------------- + memory usage of local data +------------------------------------------------------------------------- */ + +double ComputeGridLocal::memory_usage() +{ + int nbytes = size_local_rows * size_local_cols * sizeof(double); // gridlocal + return nbytes; +} diff --git a/src/ML-SNAP/compute_grid_local.h b/src/ML-SNAP/compute_grid_local.h new file mode 100644 index 0000000000..5bc6f2082a --- /dev/null +++ b/src/ML-SNAP/compute_grid_local.h @@ -0,0 +1,56 @@ +/* -*- c++ -*- ---------------------------------------------------------- + LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator + https://www.lammps.org/ Sandia National Laboratories + Steve Plimpton, sjplimp@sandia.gov + + Copyright (2003) Sandia Corporation. Under the terms of Contract + DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains + certain rights in this software. This software is distributed under + the GNU General Public License. + + See the README file in the top-level LAMMPS directory. +------------------------------------------------------------------------- */ + +#ifndef LMP_COMPUTE_GRID_LOCAL_H +#define LMP_COMPUTE_GRID_LOCAL_H + +#include "compute.h" + +namespace LAMMPS_NS { + +class ComputeGridLocal : public Compute { + public: + ComputeGridLocal(class LAMMPS *, int, char **); + ~ComputeGridLocal() override; + void setup() override; + void compute_local() override = 0; + + double memory_usage() override; + + protected: + int nx, ny, nz; // global grid dimensions + int nxlo, nxhi, nylo, nyhi, nzlo, nzhi; // local grid bounds, inclusive + int nvalues; // number of values per grid point + double **alocal; // pointer to Compute::array_local + int triclinic; // triclinic flag + double *boxlo, *prd; // box info (units real/ortho or reduced/tri) + double *sublo, *subhi; // subdomain info (units real/ortho or reduced/tri) + double delxinv, delyinv, delzinv; // inverse grid spacing + double delx, dely, delz; // grid spacing + int nargbase; // number of base class args + double cutmax; // largest cutoff distance + int size_local_cols_base; // number of columns used for coords, etc. + int gridlocal_allocated; // shows if gridlocal allocated + + void allocate(); // create arrays + void deallocate(); // free arrays + void grid2x(int, int, int, double *); // convert global indices to coordinates + void grid2lamda(int, int, int, double *); // convert global indices to lamda coordinates + void set_grid_global(); // set global grid + void set_grid_local(); // set bounds for local grid + void assign_coords(); // assign coords for grid +}; + +} // namespace LAMMPS_NS + +#endif diff --git a/src/ML-SNAP/compute_sna_atom.cpp b/src/ML-SNAP/compute_sna_atom.cpp index 5f24ebeaa6..1d2190d50d 100644 --- a/src/ML-SNAP/compute_sna_atom.cpp +++ b/src/ML-SNAP/compute_sna_atom.cpp @@ -35,20 +35,21 @@ ComputeSNAAtom::ComputeSNAAtom(LAMMPS *lmp, int narg, char **arg) : radelem(nullptr), wjelem(nullptr), sinnerelem(nullptr), dinnerelem(nullptr) { - double rmin0, rfac0; + // begin code common to all SNAP computes + + double rfac0, rmin0; int twojmax, switchflag, bzeroflag, bnormflag, wselfallflag; int ntypes = atom->ntypes; - int nargmin = 6+2*ntypes; + int nargmin = 6 + 2 * ntypes; - if (narg < nargmin) error->all(FLERR,"Illegal compute sna/atom command"); + if (narg < nargmin) error->all(FLERR, "Illegal compute {} command", style); // default values rmin0 = 0.0; switchflag = 1; bzeroflag = 1; - bnormflag = 0; quadraticflag = 0; chemflag = 0; bnormflag = 0; @@ -56,32 +57,34 @@ ComputeSNAAtom::ComputeSNAAtom(LAMMPS *lmp, int narg, char **arg) : switchinnerflag = 0; nelements = 1; - // offset by 1 to match up with types + // process required arguments - memory->create(radelem,ntypes+1,"sna/atom:radelem"); - memory->create(wjelem,ntypes+1,"sna/atom:wjelem"); + memory->create(radelem, ntypes + 1, "sna/atom:radelem"); // offset by 1 to match up with types + memory->create(wjelem, ntypes + 1, "sna/atom:wjelem"); - rcutfac = atof(arg[3]); - rfac0 = atof(arg[4]); - twojmax = atoi(arg[5]); + rcutfac = utils::numeric(FLERR, arg[3], false, lmp); + rfac0 = utils::numeric(FLERR, arg[4], false, lmp); + twojmax = utils::inumeric(FLERR, arg[5], false, lmp); for (int i = 0; i < ntypes; i++) - radelem[i+1] = atof(arg[6+i]); + radelem[i + 1] = + utils::numeric(FLERR, arg[6 + i], false, lmp); for (int i = 0; i < ntypes; i++) - wjelem[i+1] = atof(arg[6+ntypes+i]); + wjelem[i + 1] = + utils::numeric(FLERR, arg[6 + ntypes + i], false, lmp); // construct cutsq double cut; cutmax = 0.0; - memory->create(cutsq,ntypes+1,ntypes+1,"sna/atom:cutsq"); + memory->create(cutsq, ntypes + 1, ntypes + 1, "sna/atom:cutsq"); for (int i = 1; i <= ntypes; i++) { - cut = 2.0*radelem[i]*rcutfac; + cut = 2.0 * radelem[i] * rcutfac; if (cut > cutmax) cutmax = cut; - cutsq[i][i] = cut*cut; - for (int j = i+1; j <= ntypes; j++) { - cut = (radelem[i]+radelem[j])*rcutfac; - cutsq[i][j] = cutsq[j][i] = cut*cut; + cutsq[i][i] = cut * cut; + for (int j = i + 1; j <= ntypes; j++) { + cut = (radelem[i] + radelem[j]) * rcutfac; + cutsq[i][j] = cutsq[j][i] = cut * cut; } } @@ -95,89 +98,87 @@ ComputeSNAAtom::ComputeSNAAtom(LAMMPS *lmp, int narg, char **arg) : int iarg = nargmin; while (iarg < narg) { - if (strcmp(arg[iarg],"rmin0") == 0) { - if (iarg+2 > narg) - error->all(FLERR,"Illegal compute sna/atom command"); - rmin0 = atof(arg[iarg+1]); + if (strcmp(arg[iarg], "rmin0") == 0) { + if (iarg + 2 > narg) error->all(FLERR, "Illegal compute {} command", style); + rmin0 = utils::numeric(FLERR, arg[iarg + 1], false, lmp); iarg += 2; - } else if (strcmp(arg[iarg],"switchflag") == 0) { - if (iarg+2 > narg) - error->all(FLERR,"Illegal compute sna/atom command"); - switchflag = atoi(arg[iarg+1]); + } else if (strcmp(arg[iarg], "switchflag") == 0) { + if (iarg + 2 > narg) error->all(FLERR, "Illegal compute {} command", style); + switchflag = utils::inumeric(FLERR, arg[iarg + 1], false, lmp); iarg += 2; - } else if (strcmp(arg[iarg],"bzeroflag") == 0) { - if (iarg+2 > narg) - error->all(FLERR,"Illegal compute sna/atom command"); - bzeroflag = atoi(arg[iarg+1]); + } else if (strcmp(arg[iarg], "bzeroflag") == 0) { + if (iarg + 2 > narg) error->all(FLERR, "Illegal compute {} command", style); + bzeroflag = utils::inumeric(FLERR, arg[iarg + 1], false, lmp); iarg += 2; - } else if (strcmp(arg[iarg],"quadraticflag") == 0) { - if (iarg+2 > narg) - error->all(FLERR,"Illegal compute sna/atom command"); - quadraticflag = atoi(arg[iarg+1]); + } else if (strcmp(arg[iarg], "quadraticflag") == 0) { + if (iarg + 2 > narg) error->all(FLERR, "Illegal compute {} command", style); + quadraticflag = utils::inumeric(FLERR, arg[iarg + 1], false, lmp); iarg += 2; - } else if (strcmp(arg[iarg],"chem") == 0) { - if (iarg+2+ntypes > narg) - error->all(FLERR,"Illegal compute sna/atom command"); + } else if (strcmp(arg[iarg], "chem") == 0) { + if (iarg + 2 > narg) error->all(FLERR, "Illegal compute {} command", style); chemflag = 1; - memory->create(map,ntypes+1,"compute_sna_atom:map"); - nelements = utils::inumeric(FLERR,arg[iarg+1],false,lmp); + memory->create(map, ntypes + 1, "compute_sna_grid:map"); + nelements = utils::inumeric(FLERR, arg[iarg + 1], false, lmp); for (int i = 0; i < ntypes; i++) { - int jelem = utils::inumeric(FLERR,arg[iarg+2+i],false,lmp); - if (jelem < 0 || jelem >= nelements) - error->all(FLERR,"Illegal compute sna/atom command"); - map[i+1] = jelem; + int jelem = utils::inumeric(FLERR, arg[iarg + 2 + i], false, lmp); + if (jelem < 0 || jelem >= nelements) error->all(FLERR, "Illegal compute {} command", style); + map[i + 1] = jelem; } - iarg += 2+ntypes; - } else if (strcmp(arg[iarg],"bnormflag") == 0) { - if (iarg+2 > narg) - error->all(FLERR,"Illegal compute sna/atom command"); - bnormflag = atoi(arg[iarg+1]); + iarg += 2 + ntypes; + } else if (strcmp(arg[iarg], "bnormflag") == 0) { + if (iarg + 2 > narg) error->all(FLERR, "Illegal compute {} command", style); + bnormflag = utils::inumeric(FLERR, arg[iarg + 1], false, lmp); iarg += 2; - } else if (strcmp(arg[iarg],"wselfallflag") == 0) { - if (iarg+2 > narg) - error->all(FLERR,"Illegal compute sna/atom command"); - wselfallflag = atoi(arg[iarg+1]); + } else if (strcmp(arg[iarg], "wselfallflag") == 0) { + if (iarg + 2 > narg) error->all(FLERR, "Illegal compute {} command", style); + wselfallflag = utils::inumeric(FLERR, arg[iarg + 1], false, lmp); iarg += 2; - } else if (strcmp(arg[iarg],"switchinnerflag") == 0) { - if (iarg+2 > narg) - error->all(FLERR,"Illegal compute sna/atom command"); - switchinnerflag = atoi(arg[iarg+1]); + } else if (strcmp(arg[iarg], "switchinnerflag") == 0) { + if (iarg + 2 > narg) error->all(FLERR, "Illegal compute {} command", style); + switchinnerflag = utils::inumeric(FLERR, arg[iarg + 1], false, lmp); iarg += 2; - } else if (strcmp(arg[iarg],"sinner") == 0) { + } else if (strcmp(arg[iarg], "sinner") == 0) { iarg++; - if (iarg+ntypes > narg) - error->all(FLERR,"Illegal compute sna/atom command"); - memory->create(sinnerelem,ntypes+1,"sna/atom:sinnerelem"); + if (iarg + ntypes > narg) error->all(FLERR, "Illegal compute {} command", style); + memory->create(sinnerelem, ntypes + 1, "snap:sinnerelem"); for (int i = 0; i < ntypes; i++) - sinnerelem[i+1] = utils::numeric(FLERR,arg[iarg+i],false,lmp); + sinnerelem[i + 1] = utils::numeric(FLERR, arg[iarg + i], false, lmp); sinnerflag = 1; iarg += ntypes; - } else if (strcmp(arg[iarg],"dinner") == 0) { + } else if (strcmp(arg[iarg], "dinner") == 0) { iarg++; - if (iarg+ntypes > narg) - error->all(FLERR,"Illegal compute sna/atom command"); - memory->create(dinnerelem,ntypes+1,"sna/atom:dinnerelem"); + if (iarg + ntypes > narg) error->all(FLERR, "Illegal compute {} command", style); + memory->create(dinnerelem, ntypes + 1, "snap:dinnerelem"); for (int i = 0; i < ntypes; i++) - dinnerelem[i+1] = utils::numeric(FLERR,arg[iarg+i],false,lmp); + dinnerelem[i + 1] = utils::numeric(FLERR, arg[iarg + i], false, lmp); dinnerflag = 1; iarg += ntypes; - } else error->all(FLERR,"Illegal compute sna/atom command"); + } else + error->all(FLERR, "Illegal compute {} command", style); } if (switchinnerflag && !(sinnerflag && dinnerflag)) - error->all(FLERR,"Illegal compute sna/atom command: switchinnerflag = 1, missing sinner/dinner keyword"); + error->all( + FLERR, + "Illegal compute {} command: switchinnerflag = 1, missing sinner/dinner keyword", + style); if (!switchinnerflag && (sinnerflag || dinnerflag)) - error->all(FLERR,"Illegal compute sna/atom command: switchinnerflag = 0, unexpected sinner/dinner keyword"); + error->all( + FLERR, + "Illegal compute {} command: switchinnerflag = 0, unexpected sinner/dinner keyword", + style); - snaptr = new SNA(lmp, rfac0, twojmax, - rmin0, switchflag, bzeroflag, - chemflag, bnormflag, wselfallflag, - nelements, switchinnerflag); + snaptr = new SNA(lmp, rfac0, twojmax, rmin0, switchflag, bzeroflag, chemflag, bnormflag, + wselfallflag, nelements, switchinnerflag); ncoeff = snaptr->ncoeff; - size_peratom_cols = ncoeff; - if (quadraticflag) size_peratom_cols += (ncoeff*(ncoeff+1))/2; + nvalues = ncoeff; + if (quadraticflag) nvalues += (ncoeff * (ncoeff + 1)) / 2; + + // end code common to all SNAP computes + + size_peratom_cols = nvalues; peratom_flag = 1; nmax = 0; diff --git a/src/ML-SNAP/compute_sna_atom.h b/src/ML-SNAP/compute_sna_atom.h index 7d1ebfa2f8..059062670f 100644 --- a/src/ML-SNAP/compute_sna_atom.h +++ b/src/ML-SNAP/compute_sna_atom.h @@ -50,6 +50,7 @@ class ComputeSNAAtom : public Compute { class SNA *snaptr; double cutmax; int quadraticflag; + int nvalues; }; } // namespace LAMMPS_NS diff --git a/src/ML-SNAP/compute_sna_grid.cpp b/src/ML-SNAP/compute_sna_grid.cpp new file mode 100644 index 0000000000..497390f99a --- /dev/null +++ b/src/ML-SNAP/compute_sna_grid.cpp @@ -0,0 +1,320 @@ +/* ---------------------------------------------------------------------- + LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator + https://www.lammps.org/ Sandia National Laboratories + Steve Plimpton, sjplimp@sandia.gov + + Copyright (2003) Sandia Corporation. Under the terms of Contract + DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains + certain rights in this software. This software is distributed under + the GNU General Public License. + + See the README file in the top-level LAMMPS directory. +------------------------------------------------------------------------- */ + +#include "compute_sna_grid.h" + +#include "atom.h" +#include "comm.h" +#include "error.h" +#include "force.h" +#include "memory.h" +#include "modify.h" +#include "sna.h" +#include "update.h" + +#include +#include + +using namespace LAMMPS_NS; + +ComputeSNAGrid::ComputeSNAGrid(LAMMPS *lmp, int narg, char **arg) : + ComputeGrid(lmp, narg, arg), cutsq(nullptr), radelem(nullptr), wjelem(nullptr) +{ + // skip over arguments used by base class + // so that argument positions are identical to + // regular per-atom compute + + arg += nargbase; + narg -= nargbase; + + // begin code common to all SNAP computes + + double rfac0, rmin0; + int twojmax, switchflag, bzeroflag, bnormflag, wselfallflag; + + int ntypes = atom->ntypes; + int nargmin = 6 + 2 * ntypes; + + if (narg < nargmin) error->all(FLERR, "Illegal compute {} command", style); + + // default values + + rmin0 = 0.0; + switchflag = 1; + bzeroflag = 1; + quadraticflag = 0; + chemflag = 0; + bnormflag = 0; + wselfallflag = 0; + switchinnerflag = 0; + nelements = 1; + + // process required arguments + + memory->create(radelem, ntypes + 1, "sna/atom:radelem"); // offset by 1 to match up with types + memory->create(wjelem, ntypes + 1, "sna/atom:wjelem"); + + rcutfac = utils::numeric(FLERR, arg[3], false, lmp); + rfac0 = utils::numeric(FLERR, arg[4], false, lmp); + twojmax = utils::inumeric(FLERR, arg[5], false, lmp); + + for (int i = 0; i < ntypes; i++) radelem[i + 1] = utils::numeric(FLERR, arg[6 + i], false, lmp); + for (int i = 0; i < ntypes; i++) + wjelem[i + 1] = utils::numeric(FLERR, arg[6 + ntypes + i], false, lmp); + + // construct cutsq + + double cut; + cutmax = 0.0; + memory->create(cutsq, ntypes + 1, ntypes + 1, "sna/atom:cutsq"); + for (int i = 1; i <= ntypes; i++) { + cut = 2.0 * radelem[i] * rcutfac; + if (cut > cutmax) cutmax = cut; + cutsq[i][i] = cut * cut; + for (int j = i + 1; j <= ntypes; j++) { + cut = (radelem[i] + radelem[j]) * rcutfac; + cutsq[i][j] = cutsq[j][i] = cut * cut; + } + } + + // set local input checks + + int sinnerflag = 0; + int dinnerflag = 0; + + // process optional args + + int iarg = nargmin; + + while (iarg < narg) { + if (strcmp(arg[iarg], "rmin0") == 0) { + if (iarg + 2 > narg) error->all(FLERR, "Illegal compute {} command", style); + rmin0 = utils::numeric(FLERR, arg[iarg + 1], false, lmp); + iarg += 2; + } else if (strcmp(arg[iarg], "switchflag") == 0) { + if (iarg + 2 > narg) error->all(FLERR, "Illegal compute {} command", style); + switchflag = utils::inumeric(FLERR, arg[iarg + 1], false, lmp); + iarg += 2; + } else if (strcmp(arg[iarg], "bzeroflag") == 0) { + if (iarg + 2 > narg) error->all(FLERR, "Illegal compute {} command", style); + bzeroflag = utils::inumeric(FLERR, arg[iarg + 1], false, lmp); + iarg += 2; + } else if (strcmp(arg[iarg], "quadraticflag") == 0) { + if (iarg + 2 > narg) error->all(FLERR, "Illegal compute {} command", style); + quadraticflag = utils::inumeric(FLERR, arg[iarg + 1], false, lmp); + iarg += 2; + } else if (strcmp(arg[iarg], "chem") == 0) { + if (iarg + 2 > narg) error->all(FLERR, "Illegal compute {} command", style); + chemflag = 1; + memory->create(map, ntypes + 1, "compute_sna_grid:map"); + nelements = utils::inumeric(FLERR, arg[iarg + 1], false, lmp); + for (int i = 0; i < ntypes; i++) { + int jelem = utils::inumeric(FLERR, arg[iarg + 2 + i], false, lmp); + if (jelem < 0 || jelem >= nelements) error->all(FLERR, "Illegal compute {} command", style); + map[i + 1] = jelem; + } + iarg += 2 + ntypes; + } else if (strcmp(arg[iarg], "bnormflag") == 0) { + if (iarg + 2 > narg) error->all(FLERR, "Illegal compute {} command", style); + bnormflag = utils::inumeric(FLERR, arg[iarg + 1], false, lmp); + iarg += 2; + } else if (strcmp(arg[iarg], "wselfallflag") == 0) { + if (iarg + 2 > narg) error->all(FLERR, "Illegal compute {} command", style); + wselfallflag = utils::inumeric(FLERR, arg[iarg + 1], false, lmp); + iarg += 2; + } else if (strcmp(arg[iarg], "switchinnerflag") == 0) { + if (iarg + 2 > narg) error->all(FLERR, "Illegal compute {} command", style); + switchinnerflag = utils::inumeric(FLERR, arg[iarg + 1], false, lmp); + iarg += 2; + } else if (strcmp(arg[iarg], "sinner") == 0) { + iarg++; + if (iarg + ntypes > narg) error->all(FLERR, "Illegal compute {} command", style); + memory->create(sinnerelem, ntypes + 1, "snap:sinnerelem"); + for (int i = 0; i < ntypes; i++) + sinnerelem[i + 1] = utils::numeric(FLERR, arg[iarg + i], false, lmp); + sinnerflag = 1; + iarg += ntypes; + } else if (strcmp(arg[iarg], "dinner") == 0) { + iarg++; + if (iarg + ntypes > narg) error->all(FLERR, "Illegal compute {} command", style); + memory->create(dinnerelem, ntypes + 1, "snap:dinnerelem"); + for (int i = 0; i < ntypes; i++) + dinnerelem[i + 1] = utils::numeric(FLERR, arg[iarg + i], false, lmp); + dinnerflag = 1; + iarg += ntypes; + } else + error->all(FLERR, "Illegal compute {} command", style); + } + + if (switchinnerflag && !(sinnerflag && dinnerflag)) + error->all(FLERR, + "Illegal compute {} command: switchinnerflag = 1, missing sinner/dinner keyword", + style); + + if (!switchinnerflag && (sinnerflag || dinnerflag)) + error->all(FLERR, + "Illegal compute {} command: switchinnerflag = 0, unexpected sinner/dinner keyword", + style); + + snaptr = new SNA(lmp, rfac0, twojmax, rmin0, switchflag, bzeroflag, chemflag, bnormflag, + wselfallflag, nelements, switchinnerflag); + + ncoeff = snaptr->ncoeff; + nvalues = ncoeff; + if (quadraticflag) nvalues += (ncoeff * (ncoeff + 1)) / 2; + + // end code common to all SNAP computes + + size_array_cols = size_array_cols_base + nvalues; + array_flag = 1; +} + +/* ---------------------------------------------------------------------- */ + +ComputeSNAGrid::~ComputeSNAGrid() +{ + memory->destroy(radelem); + memory->destroy(wjelem); + memory->destroy(cutsq); + delete snaptr; + + if (chemflag) memory->destroy(map); +} + +/* ---------------------------------------------------------------------- */ + +void ComputeSNAGrid::init() +{ + if ((modify->get_compute_by_style("^sna/grid$").size() > 1) && (comm->me == 0)) + error->warning(FLERR, "More than one instance of compute sna/grid"); + snaptr->init(); +} + +/* ---------------------------------------------------------------------- */ + +void ComputeSNAGrid::compute_array() +{ + invoked_array = update->ntimestep; + + // compute sna for each gridpoint + + double **const x = atom->x; + const int *const mask = atom->mask; + int *const type = atom->type; + const int ntotal = atom->nlocal + atom->nghost; + + // insure rij, inside, and typej are of size jnum + + snaptr->grow_rij(ntotal); + + for (int iz = nzlo; iz <= nzhi; iz++) + for (int iy = nylo; iy <= nyhi; iy++) + for (int ix = nxlo; ix <= nxhi; ix++) { + double xgrid[3]; + const int igrid = iz * (nx * ny) + iy * nx + ix; + grid2x(igrid, xgrid); + const double xtmp = xgrid[0]; + const double ytmp = xgrid[1]; + const double ztmp = xgrid[2]; + + // currently, all grid points are type 1 + // not clear what a better choice would be + + const int itype = 1; + int ielem = 0; + if (chemflag) ielem = map[itype]; + + // rij[][3] = displacements between atom I and those neighbors + // inside = indices of neighbors of I within cutoff + // typej = types of neighbors of I within cutoff + + int ninside = 0; + for (int j = 0; j < ntotal; j++) { + + // check that j is in compute group + + if (!(mask[j] & groupbit)) continue; + + const double delx = xtmp - x[j][0]; + const double dely = ytmp - x[j][1]; + const double delz = ztmp - x[j][2]; + const double rsq = delx * delx + dely * dely + delz * delz; + int jtype = type[j]; + int jelem = 0; + if (chemflag) jelem = map[jtype]; + + if (rsq < cutsq[jtype][jtype] && rsq > 1e-20) { + snaptr->rij[ninside][0] = delx; + snaptr->rij[ninside][1] = dely; + snaptr->rij[ninside][2] = delz; + snaptr->inside[ninside] = j; + snaptr->wj[ninside] = wjelem[jtype]; + snaptr->rcutij[ninside] = 2.0 * radelem[jtype] * rcutfac; + if (switchinnerflag) { + snaptr->sinnerij[ninside] = sinnerelem[jelem]; + snaptr->dinnerij[ninside] = dinnerelem[jelem]; + } + if (chemflag) snaptr->element[ninside] = jelem; + ninside++; + } + } + + snaptr->compute_ui(ninside, ielem); + snaptr->compute_zi(); + snaptr->compute_bi(ielem); + + // linear contributions + + for (int icoeff = 0; icoeff < ncoeff; icoeff++) + gridlocal[size_array_cols_base + icoeff][iz][iy][ix] = snaptr->blist[icoeff]; + + // quadratic contributions + + if (quadraticflag) { + int ncount = ncoeff; + for (int icoeff = 0; icoeff < ncoeff; icoeff++) { + double bveci = snaptr->blist[icoeff]; + gridlocal[size_array_cols_base + ncount++][iz][iy][ix] = 0.5 * bveci * bveci; + for (int jcoeff = icoeff + 1; jcoeff < ncoeff; jcoeff++) + gridlocal[size_array_cols_base + ncount++][iz][iy][ix] = + bveci * snaptr->blist[jcoeff]; + } + } + } + + memset(&grid[0][0], 0, size_array_rows * size_array_cols * sizeof(double)); + + for (int iz = nzlo; iz <= nzhi; iz++) + for (int iy = nylo; iy <= nyhi; iy++) + for (int ix = nxlo; ix <= nxhi; ix++) { + const int igrid = iz * (nx * ny) + iy * nx + ix; + for (int j = 0; j < nvalues; j++) + grid[igrid][size_array_cols_base + j] = gridlocal[size_array_cols_base + j][iz][iy][ix]; + } + MPI_Allreduce(&grid[0][0], &gridall[0][0], size_array_rows * size_array_cols, MPI_DOUBLE, MPI_SUM, + world); + assign_coords_all(); +} + +/* ---------------------------------------------------------------------- + memory usage +------------------------------------------------------------------------- */ + +double ComputeSNAGrid::memory_usage() +{ + double nbytes = snaptr->memory_usage(); // SNA object + int n = atom->ntypes + 1; + nbytes += (double) n * sizeof(int); // map + + return nbytes; +} diff --git a/src/ML-SNAP/compute_sna_grid.h b/src/ML-SNAP/compute_sna_grid.h new file mode 100644 index 0000000000..839003dc2e --- /dev/null +++ b/src/ML-SNAP/compute_sna_grid.h @@ -0,0 +1,54 @@ +/* -*- c++ -*- ---------------------------------------------------------- + LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator + https://www.lammps.org/ Sandia National Laboratories + Steve Plimpton, sjplimp@sandia.gov + + Copyright (2003) Sandia Corporation. Under the terms of Contract + DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains + certain rights in this software. This software is distributed under + the GNU General Public License. + + See the README file in the top-level LAMMPS directory. +------------------------------------------------------------------------- */ + +#ifdef COMPUTE_CLASS +// clang-format off +ComputeStyle(sna/grid,ComputeSNAGrid); +// clang-format on +#else + +#ifndef LMP_COMPUTE_SNA_GRID_H +#define LMP_COMPUTE_SNA_GRID_H + +#include "compute_grid.h" + +namespace LAMMPS_NS { + +class ComputeSNAGrid : public ComputeGrid { + public: + ComputeSNAGrid(class LAMMPS *, int, char **); + ~ComputeSNAGrid() override; + void init() override; + void compute_array() override; + double memory_usage() override; + + private: + int ncoeff; + double **cutsq; + double rcutfac; + double *radelem; + double *wjelem; + int *map; // map types to [0,nelements) + int nelements, chemflag; + int switchinnerflag; + double *sinnerelem; + double *dinnerelem; + class SNA *snaptr; + double cutmax; + int quadraticflag; +}; + +} // namespace LAMMPS_NS + +#endif +#endif diff --git a/src/ML-SNAP/compute_sna_grid_local.cpp b/src/ML-SNAP/compute_sna_grid_local.cpp new file mode 100644 index 0000000000..80a1baddab --- /dev/null +++ b/src/ML-SNAP/compute_sna_grid_local.cpp @@ -0,0 +1,306 @@ +/* ---------------------------------------------------------------------- + LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator + https://www.lammps.org/ Sandia National Laboratories + Steve Plimpton, sjplimp@sandia.gov + + Copyright (2003) Sandia Corporation. Under the terms of Contract + DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains + certain rights in this software. This software is distributed under + the GNU General Public License. + + See the README file in the top-level LAMMPS directory. +------------------------------------------------------------------------- */ + +#include "compute_sna_grid_local.h" + +#include "atom.h" +#include "comm.h" +#include "error.h" +#include "force.h" +#include "memory.h" +#include "modify.h" +#include "sna.h" +#include "update.h" + +#include +#include + +using namespace LAMMPS_NS; + +ComputeSNAGridLocal::ComputeSNAGridLocal(LAMMPS *lmp, int narg, char **arg) : + ComputeGridLocal(lmp, narg, arg), cutsq(nullptr), radelem(nullptr), wjelem(nullptr) +{ + // skip over arguments used by base class + // so that argument positions are identical to + // regular per-atom compute + + arg += nargbase; + narg -= nargbase; + + // begin code common to all SNAP computes + + double rfac0, rmin0; + int twojmax, switchflag, bzeroflag, bnormflag, wselfallflag; + + int ntypes = atom->ntypes; + int nargmin = 6 + 2 * ntypes; + + if (narg < nargmin) error->all(FLERR, "Illegal compute {} command", style); + + // default values + + rmin0 = 0.0; + switchflag = 1; + bzeroflag = 1; + quadraticflag = 0; + chemflag = 0; + bnormflag = 0; + wselfallflag = 0; + switchinnerflag = 0; + nelements = 1; + + // process required arguments + + memory->create(radelem, ntypes + 1, "sna/atom:radelem"); // offset by 1 to match up with types + memory->create(wjelem, ntypes + 1, "sna/atom:wjelem"); + + rcutfac = utils::numeric(FLERR, arg[3], false, lmp); + rfac0 = utils::numeric(FLERR, arg[4], false, lmp); + twojmax = utils::inumeric(FLERR, arg[5], false, lmp); + + for (int i = 0; i < ntypes; i++) radelem[i + 1] = utils::numeric(FLERR, arg[6 + i], false, lmp); + for (int i = 0; i < ntypes; i++) + wjelem[i + 1] = utils::numeric(FLERR, arg[6 + ntypes + i], false, lmp); + + // construct cutsq + + double cut; + cutmax = 0.0; + memory->create(cutsq, ntypes + 1, ntypes + 1, "sna/atom:cutsq"); + for (int i = 1; i <= ntypes; i++) { + cut = 2.0 * radelem[i] * rcutfac; + if (cut > cutmax) cutmax = cut; + cutsq[i][i] = cut * cut; + for (int j = i + 1; j <= ntypes; j++) { + cut = (radelem[i] + radelem[j]) * rcutfac; + cutsq[i][j] = cutsq[j][i] = cut * cut; + } + } + + // set local input checks + + int sinnerflag = 0; + int dinnerflag = 0; + + // process optional args + + int iarg = nargmin; + + while (iarg < narg) { + if (strcmp(arg[iarg], "rmin0") == 0) { + if (iarg + 2 > narg) error->all(FLERR, "Illegal compute {} command", style); + rmin0 = utils::numeric(FLERR, arg[iarg + 1], false, lmp); + iarg += 2; + } else if (strcmp(arg[iarg], "switchflag") == 0) { + if (iarg + 2 > narg) error->all(FLERR, "Illegal compute {} command", style); + switchflag = utils::inumeric(FLERR, arg[iarg + 1], false, lmp); + iarg += 2; + } else if (strcmp(arg[iarg], "bzeroflag") == 0) { + if (iarg + 2 > narg) error->all(FLERR, "Illegal compute {} command", style); + bzeroflag = utils::inumeric(FLERR, arg[iarg + 1], false, lmp); + iarg += 2; + } else if (strcmp(arg[iarg], "quadraticflag") == 0) { + if (iarg + 2 > narg) error->all(FLERR, "Illegal compute {} command", style); + quadraticflag = utils::inumeric(FLERR, arg[iarg + 1], false, lmp); + iarg += 2; + } else if (strcmp(arg[iarg], "chem") == 0) { + if (iarg + 2 > narg) error->all(FLERR, "Illegal compute {} command", style); + chemflag = 1; + memory->create(map, ntypes + 1, "compute_sna_grid:map"); + nelements = utils::inumeric(FLERR, arg[iarg + 1], false, lmp); + for (int i = 0; i < ntypes; i++) { + int jelem = utils::inumeric(FLERR, arg[iarg + 2 + i], false, lmp); + if (jelem < 0 || jelem >= nelements) error->all(FLERR, "Illegal compute {} command", style); + map[i + 1] = jelem; + } + iarg += 2 + ntypes; + } else if (strcmp(arg[iarg], "bnormflag") == 0) { + if (iarg + 2 > narg) error->all(FLERR, "Illegal compute {} command", style); + bnormflag = utils::inumeric(FLERR, arg[iarg + 1], false, lmp); + iarg += 2; + } else if (strcmp(arg[iarg], "wselfallflag") == 0) { + if (iarg + 2 > narg) error->all(FLERR, "Illegal compute {} command", style); + wselfallflag = utils::inumeric(FLERR, arg[iarg + 1], false, lmp); + iarg += 2; + } else if (strcmp(arg[iarg], "switchinnerflag") == 0) { + if (iarg + 2 > narg) error->all(FLERR, "Illegal compute {} command", style); + switchinnerflag = utils::inumeric(FLERR, arg[iarg + 1], false, lmp); + iarg += 2; + } else if (strcmp(arg[iarg], "sinner") == 0) { + iarg++; + if (iarg + ntypes > narg) error->all(FLERR, "Illegal compute {} command", style); + memory->create(sinnerelem, ntypes + 1, "snap:sinnerelem"); + for (int i = 0; i < ntypes; i++) + sinnerelem[i + 1] = utils::numeric(FLERR, arg[iarg + i], false, lmp); + sinnerflag = 1; + iarg += ntypes; + } else if (strcmp(arg[iarg], "dinner") == 0) { + iarg++; + if (iarg + ntypes > narg) error->all(FLERR, "Illegal compute {} command", style); + memory->create(dinnerelem, ntypes + 1, "snap:dinnerelem"); + for (int i = 0; i < ntypes; i++) + dinnerelem[i + 1] = utils::numeric(FLERR, arg[iarg + i], false, lmp); + dinnerflag = 1; + iarg += ntypes; + } else + error->all(FLERR, "Illegal compute {} command", style); + } + + if (switchinnerflag && !(sinnerflag && dinnerflag)) + error->all(FLERR, + "Illegal compute {} command: switchinnerflag = 1, missing sinner/dinner keyword", + style); + + if (!switchinnerflag && (sinnerflag || dinnerflag)) + error->all(FLERR, + "Illegal compute {} command: switchinnerflag = 0, unexpected sinner/dinner keyword", + style); + + snaptr = new SNA(lmp, rfac0, twojmax, rmin0, switchflag, bzeroflag, chemflag, bnormflag, + wselfallflag, nelements, switchinnerflag); + + ncoeff = snaptr->ncoeff; + nvalues = ncoeff; + if (quadraticflag) nvalues += (ncoeff * (ncoeff + 1)) / 2; + + // end code common to all SNAP computes + + size_local_cols = size_local_cols_base + nvalues; +} + +/* ---------------------------------------------------------------------- */ + +ComputeSNAGridLocal::~ComputeSNAGridLocal() +{ + memory->destroy(radelem); + memory->destroy(wjelem); + memory->destroy(cutsq); + delete snaptr; + + if (chemflag) memory->destroy(map); +} + +/* ---------------------------------------------------------------------- */ + +void ComputeSNAGridLocal::init() +{ + if ((modify->get_compute_by_style("^sna/grid/local$").size() > 1) && (comm->me == 0)) + error->warning(FLERR, "More than one instance of compute sna/grid/local"); + snaptr->init(); +} + +/* ---------------------------------------------------------------------- */ + +void ComputeSNAGridLocal::compute_local() +{ + invoked_array = update->ntimestep; + + // compute sna for each gridpoint + + double **const x = atom->x; + const int *const mask = atom->mask; + int *const type = atom->type; + const int ntotal = atom->nlocal + atom->nghost; + + // insure rij, inside, and typej are of size jnum + + snaptr->grow_rij(ntotal); + + int igrid = 0; + for (int iz = nzlo; iz <= nzhi; iz++) + for (int iy = nylo; iy <= nyhi; iy++) + for (int ix = nxlo; ix <= nxhi; ix++) { + double xgrid[3]; + grid2x(ix, iy, iz, xgrid); + const double xtmp = xgrid[0]; + const double ytmp = xgrid[1]; + const double ztmp = xgrid[2]; + + // currently, all grid points are type 1 + // not clear what a better choice would be + + const int itype = 1; + int ielem = 0; + if (chemflag) ielem = map[itype]; + + // rij[][3] = displacements between atom I and those neighbors + // inside = indices of neighbors of I within cutoff + // typej = types of neighbors of I within cutoff + + int ninside = 0; + for (int j = 0; j < ntotal; j++) { + + // check that j is in compute group + + if (!(mask[j] & groupbit)) continue; + + const double delx = xtmp - x[j][0]; + const double dely = ytmp - x[j][1]; + const double delz = ztmp - x[j][2]; + const double rsq = delx * delx + dely * dely + delz * delz; + int jtype = type[j]; + int jelem = 0; + if (chemflag) jelem = map[jtype]; + if (rsq < cutsq[jtype][jtype] && rsq > 1e-20) { + snaptr->rij[ninside][0] = delx; + snaptr->rij[ninside][1] = dely; + snaptr->rij[ninside][2] = delz; + snaptr->inside[ninside] = j; + snaptr->wj[ninside] = wjelem[jtype]; + snaptr->rcutij[ninside] = 2.0 * radelem[jtype] * rcutfac; + if (switchinnerflag) { + snaptr->sinnerij[ninside] = sinnerelem[jelem]; + snaptr->dinnerij[ninside] = dinnerelem[jelem]; + } + if (chemflag) + snaptr->element[ninside] = jelem; // element index for multi-element snap + ninside++; + } + } + + snaptr->compute_ui(ninside, ielem); + snaptr->compute_zi(); + snaptr->compute_bi(ielem); + + // linear contributions + + for (int icoeff = 0; icoeff < ncoeff; icoeff++) + alocal[igrid][size_local_cols_base + icoeff] = snaptr->blist[icoeff]; + + // quadratic contributions + + if (quadraticflag) { + int ncount = ncoeff; + for (int icoeff = 0; icoeff < ncoeff; icoeff++) { + double bveci = snaptr->blist[icoeff]; + alocal[igrid][size_local_cols_base + ncount++] = 0.5 * bveci * bveci; + for (int jcoeff = icoeff + 1; jcoeff < ncoeff; jcoeff++) + alocal[igrid][size_local_cols_base + ncount++] = bveci * snaptr->blist[jcoeff]; + } + } + igrid++; + } +} + +/* ---------------------------------------------------------------------- + memory usage +------------------------------------------------------------------------- */ + +double ComputeSNAGridLocal::memory_usage() +{ + double nbytes = snaptr->memory_usage(); // SNA object + int n = atom->ntypes + 1; + nbytes += (double) n * sizeof(int); // map + + return nbytes; +} diff --git a/src/ML-SNAP/compute_sna_grid_local.h b/src/ML-SNAP/compute_sna_grid_local.h new file mode 100644 index 0000000000..60f93b92b0 --- /dev/null +++ b/src/ML-SNAP/compute_sna_grid_local.h @@ -0,0 +1,54 @@ +/* -*- c++ -*- ---------------------------------------------------------- + LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator + https://www.lammps.org/ Sandia National Laboratories + Steve Plimpton, sjplimp@sandia.gov + + Copyright (2003) Sandia Corporation. Under the terms of Contract + DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains + certain rights in this software. This software is distributed under + the GNU General Public License. + + See the README file in the top-level LAMMPS directory. +------------------------------------------------------------------------- */ + +#ifdef COMPUTE_CLASS +// clang-format off +ComputeStyle(sna/grid/local,ComputeSNAGridLocal); +// clang-format on +#else + +#ifndef LMP_COMPUTE_SNA_GRID_LOCAL_H +#define LMP_COMPUTE_SNA_GRID_LOCAL_H + +#include "compute_grid_local.h" + +namespace LAMMPS_NS { + +class ComputeSNAGridLocal : public ComputeGridLocal { + public: + ComputeSNAGridLocal(class LAMMPS *, int, char **); + ~ComputeSNAGridLocal() override; + void init() override; + void compute_local() override; + double memory_usage() override; + + private: + int ncoeff; + double **cutsq; + double rcutfac; + double *radelem; + double *wjelem; + int *map; // map types to [0,nelements) + int nelements, chemflag; + int switchinnerflag; + double *sinnerelem; + double *dinnerelem; + class SNA *snaptr; + double cutmax; + int quadraticflag; +}; + +} // namespace LAMMPS_NS + +#endif +#endif diff --git a/src/ML-SNAP/compute_snad_atom.cpp b/src/ML-SNAP/compute_snad_atom.cpp index 838d8a85c9..12839629a9 100644 --- a/src/ML-SNAP/compute_snad_atom.cpp +++ b/src/ML-SNAP/compute_snad_atom.cpp @@ -34,20 +34,22 @@ ComputeSNADAtom::ComputeSNADAtom(LAMMPS *lmp, int narg, char **arg) : Compute(lmp, narg, arg), cutsq(nullptr), list(nullptr), snad(nullptr), radelem(nullptr), wjelem(nullptr), sinnerelem(nullptr), dinnerelem(nullptr) { + + // begin code common to all SNAP computes + double rfac0, rmin0; int twojmax, switchflag, bzeroflag, bnormflag, wselfallflag; int ntypes = atom->ntypes; - int nargmin = 6+2*ntypes; + int nargmin = 6 + 2 * ntypes; - if (narg < nargmin) error->all(FLERR,"Illegal compute snad/atom command"); + if (narg < nargmin) error->all(FLERR, "Illegal compute {} command", style); // default values rmin0 = 0.0; switchflag = 1; bzeroflag = 1; - bnormflag = 0; quadraticflag = 0; chemflag = 0; bnormflag = 0; @@ -57,28 +59,32 @@ ComputeSNADAtom::ComputeSNADAtom(LAMMPS *lmp, int narg, char **arg) : // process required arguments - memory->create(radelem,ntypes+1,"snad/atom:radelem"); // offset by 1 to match up with types - memory->create(wjelem,ntypes+1,"snad/atom:wjelem"); - rcutfac = atof(arg[3]); - rfac0 = atof(arg[4]); - twojmax = atoi(arg[5]); + memory->create(radelem, ntypes + 1, "sna/atom:radelem"); // offset by 1 to match up with types + memory->create(wjelem, ntypes + 1, "sna/atom:wjelem"); + + rcutfac = utils::numeric(FLERR, arg[3], false, lmp); + rfac0 = utils::numeric(FLERR, arg[4], false, lmp); + twojmax = utils::inumeric(FLERR, arg[5], false, lmp); + for (int i = 0; i < ntypes; i++) - radelem[i+1] = atof(arg[6+i]); + radelem[i + 1] = + utils::numeric(FLERR, arg[6 + i], false, lmp); for (int i = 0; i < ntypes; i++) - wjelem[i+1] = atof(arg[6+ntypes+i]); + wjelem[i + 1] = + utils::numeric(FLERR, arg[6 + ntypes + i], false, lmp); // construct cutsq double cut; cutmax = 0.0; - memory->create(cutsq,ntypes+1,ntypes+1,"snad/atom:cutsq"); + memory->create(cutsq, ntypes + 1, ntypes + 1, "sna/atom:cutsq"); for (int i = 1; i <= ntypes; i++) { - cut = 2.0*radelem[i]*rcutfac; + cut = 2.0 * radelem[i] * rcutfac; if (cut > cutmax) cutmax = cut; - cutsq[i][i] = cut*cut; - for (int j = i+1; j <= ntypes; j++) { - cut = (radelem[i]+radelem[j])*rcutfac; - cutsq[i][j] = cutsq[j][i] = cut*cut; + cutsq[i][i] = cut * cut; + for (int j = i + 1; j <= ntypes; j++) { + cut = (radelem[i] + radelem[j]) * rcutfac; + cutsq[i][j] = cutsq[j][i] = cut * cut; } } @@ -92,93 +98,89 @@ ComputeSNADAtom::ComputeSNADAtom(LAMMPS *lmp, int narg, char **arg) : int iarg = nargmin; while (iarg < narg) { - if (strcmp(arg[iarg],"rmin0") == 0) { - if (iarg+2 > narg) - error->all(FLERR,"Illegal compute snad/atom command"); - rmin0 = atof(arg[iarg+1]); + if (strcmp(arg[iarg], "rmin0") == 0) { + if (iarg + 2 > narg) error->all(FLERR, "Illegal compute {} command", style); + rmin0 = utils::numeric(FLERR, arg[iarg + 1], false, lmp); iarg += 2; - } else if (strcmp(arg[iarg],"bzeroflag") == 0) { - if (iarg+2 > narg) - error->all(FLERR,"Illegal compute snad/atom command"); - bzeroflag = atoi(arg[iarg+1]); + } else if (strcmp(arg[iarg], "switchflag") == 0) { + if (iarg + 2 > narg) error->all(FLERR, "Illegal compute {} command", style); + switchflag = utils::inumeric(FLERR, arg[iarg + 1], false, lmp); iarg += 2; - } else if (strcmp(arg[iarg],"switchflag") == 0) { - if (iarg+2 > narg) - error->all(FLERR,"Illegal compute snad/atom command"); - switchflag = atoi(arg[iarg+1]); + } else if (strcmp(arg[iarg], "bzeroflag") == 0) { + if (iarg + 2 > narg) error->all(FLERR, "Illegal compute {} command", style); + bzeroflag = utils::inumeric(FLERR, arg[iarg + 1], false, lmp); iarg += 2; - } else if (strcmp(arg[iarg],"quadraticflag") == 0) { - if (iarg+2 > narg) - error->all(FLERR,"Illegal compute snad/atom command"); - quadraticflag = atoi(arg[iarg+1]); + } else if (strcmp(arg[iarg], "quadraticflag") == 0) { + if (iarg + 2 > narg) error->all(FLERR, "Illegal compute {} command", style); + quadraticflag = utils::inumeric(FLERR, arg[iarg + 1], false, lmp); iarg += 2; - } else if (strcmp(arg[iarg],"chem") == 0) { - if (iarg+2+ntypes > narg) - error->all(FLERR,"Illegal compute snad/atom command"); + } else if (strcmp(arg[iarg], "chem") == 0) { + if (iarg + 2 > narg) error->all(FLERR, "Illegal compute {} command", style); chemflag = 1; - memory->create(map,ntypes+1,"compute_snad_atom:map"); - nelements = utils::inumeric(FLERR,arg[iarg+1],false,lmp); + memory->create(map, ntypes + 1, "compute_sna_grid:map"); + nelements = utils::inumeric(FLERR, arg[iarg + 1], false, lmp); for (int i = 0; i < ntypes; i++) { - int jelem = utils::inumeric(FLERR,arg[iarg+2+i],false,lmp); - if (jelem < 0 || jelem >= nelements) - error->all(FLERR,"Illegal compute snad/atom command"); - map[i+1] = jelem; + int jelem = utils::inumeric(FLERR, arg[iarg + 2 + i], false, lmp); + if (jelem < 0 || jelem >= nelements) error->all(FLERR, "Illegal compute {} command", style); + map[i + 1] = jelem; } - iarg += 2+ntypes; - } else if (strcmp(arg[iarg],"bnormflag") == 0) { - if (iarg+2 > narg) - error->all(FLERR,"Illegal compute snad/atom command"); - bnormflag = atoi(arg[iarg+1]); + iarg += 2 + ntypes; + } else if (strcmp(arg[iarg], "bnormflag") == 0) { + if (iarg + 2 > narg) error->all(FLERR, "Illegal compute {} command", style); + bnormflag = utils::inumeric(FLERR, arg[iarg + 1], false, lmp); iarg += 2; - } else if (strcmp(arg[iarg],"wselfallflag") == 0) { - if (iarg+2 > narg) - error->all(FLERR,"Illegal compute snad/atom command"); - wselfallflag = atoi(arg[iarg+1]); + } else if (strcmp(arg[iarg], "wselfallflag") == 0) { + if (iarg + 2 > narg) error->all(FLERR, "Illegal compute {} command", style); + wselfallflag = utils::inumeric(FLERR, arg[iarg + 1], false, lmp); iarg += 2; - } else if (strcmp(arg[iarg],"switchinnerflag") == 0) { - if (iarg+2 > narg) - error->all(FLERR,"Illegal compute snad/atom command"); - switchinnerflag = atoi(arg[iarg+1]); + } else if (strcmp(arg[iarg], "switchinnerflag") == 0) { + if (iarg + 2 > narg) error->all(FLERR, "Illegal compute {} command", style); + switchinnerflag = utils::inumeric(FLERR, arg[iarg + 1], false, lmp); iarg += 2; - } else if (strcmp(arg[iarg],"sinner") == 0) { + } else if (strcmp(arg[iarg], "sinner") == 0) { iarg++; - if (iarg+ntypes > narg) - error->all(FLERR,"Illegal compute snad/atom command"); - memory->create(sinnerelem,ntypes+1,"snad/atom:sinnerelem"); + if (iarg + ntypes > narg) error->all(FLERR, "Illegal compute {} command", style); + memory->create(sinnerelem, ntypes + 1, "snap:sinnerelem"); for (int i = 0; i < ntypes; i++) - sinnerelem[i+1] = utils::numeric(FLERR,arg[iarg+i],false,lmp); + sinnerelem[i + 1] = utils::numeric(FLERR, arg[iarg + i], false, lmp); sinnerflag = 1; iarg += ntypes; - } else if (strcmp(arg[iarg],"dinner") == 0) { + } else if (strcmp(arg[iarg], "dinner") == 0) { iarg++; - if (iarg+ntypes > narg) - error->all(FLERR,"Illegal compute snad/atom command"); - memory->create(dinnerelem,ntypes+1,"snad/atom:dinnerelem"); + if (iarg + ntypes > narg) error->all(FLERR, "Illegal compute {} command", style); + memory->create(dinnerelem, ntypes + 1, "snap:dinnerelem"); for (int i = 0; i < ntypes; i++) - dinnerelem[i+1] = utils::numeric(FLERR,arg[iarg+i],false,lmp); + dinnerelem[i + 1] = utils::numeric(FLERR, arg[iarg + i], false, lmp); dinnerflag = 1; iarg += ntypes; - } else error->all(FLERR,"Illegal compute snad/atom command"); + } else + error->all(FLERR, "Illegal compute {} command", style); } if (switchinnerflag && !(sinnerflag && dinnerflag)) - error->all(FLERR,"Illegal compute snad/atom command: switchinnerflag = 1, missing sinner/dinner keyword"); + error->all( + FLERR, + "Illegal compute {} command: switchinnerflag = 1, missing sinner/dinner keyword", + style); if (!switchinnerflag && (sinnerflag || dinnerflag)) - error->all(FLERR,"Illegal compute snad/atom command: switchinnerflag = 0, unexpected sinner/dinner keyword"); + error->all( + FLERR, + "Illegal compute {} command: switchinnerflag = 0, unexpected sinner/dinner keyword", + style); - - snaptr = new SNA(lmp, rfac0, twojmax, - rmin0, switchflag, bzeroflag, - chemflag, bnormflag, wselfallflag, - nelements, switchinnerflag); + snaptr = new SNA(lmp, rfac0, twojmax, rmin0, switchflag, bzeroflag, chemflag, bnormflag, + wselfallflag, nelements, switchinnerflag); ncoeff = snaptr->ncoeff; - nperdim = ncoeff; - if (quadraticflag) nperdim += (ncoeff*(ncoeff+1))/2; - yoffset = nperdim; - zoffset = 2*nperdim; - size_peratom_cols = 3*nperdim*atom->ntypes; + nvalues = ncoeff; + if (quadraticflag) nvalues += (ncoeff * (ncoeff + 1)) / 2; + + // end code common to all SNAP computes + + yoffset = nvalues; + zoffset = 2*nvalues; + size_peratom_cols = 3*nvalues*atom->ntypes; comm_reverse = size_peratom_cols; peratom_flag = 1; @@ -289,7 +291,7 @@ void ComputeSNADAtom::compute_peratom() // const int typeoffset = threencoeff*(atom->type[i]-1); // const int quadraticoffset = threencoeff*atom->ntypes + // threencoeffq*(atom->type[i]-1); - const int typeoffset = 3*nperdim*(atom->type[i]-1); + const int typeoffset = 3*nvalues*(atom->type[i]-1); // insure rij, inside, and typej are of size jnum diff --git a/src/ML-SNAP/compute_snad_atom.h b/src/ML-SNAP/compute_snad_atom.h index 0cdb075daf..ac27e10769 100644 --- a/src/ML-SNAP/compute_snad_atom.h +++ b/src/ML-SNAP/compute_snad_atom.h @@ -37,7 +37,7 @@ class ComputeSNADAtom : public Compute { private: int nmax; - int ncoeff, nperdim, yoffset, zoffset; + int ncoeff, nvalues, yoffset, zoffset; double **cutsq; class NeighList *list; double **snad; diff --git a/src/ML-SNAP/compute_snap.cpp b/src/ML-SNAP/compute_snap.cpp index 83f56e338c..d1f75f8382 100644 --- a/src/ML-SNAP/compute_snap.cpp +++ b/src/ML-SNAP/compute_snap.cpp @@ -41,13 +41,15 @@ ComputeSnap::ComputeSnap(LAMMPS *lmp, int narg, char **arg) : array_flag = 1; extarray = 0; + // begin code common to all SNAP computes + double rfac0, rmin0; int twojmax, switchflag, bzeroflag, bnormflag, wselfallflag; int ntypes = atom->ntypes; - int nargmin = 6+2*ntypes; + int nargmin = 6 + 2 * ntypes; - if (narg < nargmin) error->all(FLERR,"Illegal compute snap command"); + if (narg < nargmin) error->all(FLERR, "Illegal compute {} command", style); // default values @@ -55,7 +57,6 @@ ComputeSnap::ComputeSnap(LAMMPS *lmp, int narg, char **arg) : switchflag = 1; bzeroflag = 1; quadraticflag = 0; - bikflag = 0; chemflag = 0; bnormflag = 0; wselfallflag = 0; @@ -64,28 +65,32 @@ ComputeSnap::ComputeSnap(LAMMPS *lmp, int narg, char **arg) : // process required arguments - memory->create(radelem,ntypes+1,"snap:radelem"); // offset by 1 to match up with types - memory->create(wjelem,ntypes+1,"snap:wjelem"); - rcutfac = atof(arg[3]); - rfac0 = atof(arg[4]); - twojmax = atoi(arg[5]); + memory->create(radelem, ntypes + 1, "sna/atom:radelem"); // offset by 1 to match up with types + memory->create(wjelem, ntypes + 1, "sna/atom:wjelem"); + + rcutfac = utils::numeric(FLERR, arg[3], false, lmp); + rfac0 = utils::numeric(FLERR, arg[4], false, lmp); + twojmax = utils::inumeric(FLERR, arg[5], false, lmp); + for (int i = 0; i < ntypes; i++) - radelem[i+1] = atof(arg[6+i]); + radelem[i + 1] = + utils::numeric(FLERR, arg[6 + i], false, lmp); for (int i = 0; i < ntypes; i++) - wjelem[i+1] = atof(arg[6+ntypes+i]); + wjelem[i + 1] = + utils::numeric(FLERR, arg[6 + ntypes + i], false, lmp); // construct cutsq double cut; cutmax = 0.0; - memory->create(cutsq,ntypes+1,ntypes+1,"snap:cutsq"); + memory->create(cutsq, ntypes + 1, ntypes + 1, "sna/atom:cutsq"); for (int i = 1; i <= ntypes; i++) { - cut = 2.0*radelem[i]*rcutfac; + cut = 2.0 * radelem[i] * rcutfac; if (cut > cutmax) cutmax = cut; - cutsq[i][i] = cut*cut; - for (int j = i+1; j <= ntypes; j++) { - cut = (radelem[i]+radelem[j])*rcutfac; - cutsq[i][j] = cutsq[j][i] = cut*cut; + cutsq[i][i] = cut * cut; + for (int j = i + 1; j <= ntypes; j++) { + cut = (radelem[i] + radelem[j]) * rcutfac; + cutsq[i][j] = cutsq[j][i] = cut * cut; } } @@ -99,107 +104,99 @@ ComputeSnap::ComputeSnap(LAMMPS *lmp, int narg, char **arg) : int iarg = nargmin; while (iarg < narg) { - if (strcmp(arg[iarg],"rmin0") == 0) { - if (iarg+2 > narg) - error->all(FLERR,"Illegal compute snap command"); - rmin0 = atof(arg[iarg+1]); + if (strcmp(arg[iarg], "rmin0") == 0) { + if (iarg + 2 > narg) error->all(FLERR, "Illegal compute {} command", style); + rmin0 = utils::numeric(FLERR, arg[iarg + 1], false, lmp); iarg += 2; - } else if (strcmp(arg[iarg],"bzeroflag") == 0) { - if (iarg+2 > narg) - error->all(FLERR,"Illegal compute snap command"); - bzeroflag = atoi(arg[iarg+1]); + } else if (strcmp(arg[iarg], "switchflag") == 0) { + if (iarg + 2 > narg) error->all(FLERR, "Illegal compute {} command", style); + switchflag = utils::inumeric(FLERR, arg[iarg + 1], false, lmp); iarg += 2; - } else if (strcmp(arg[iarg],"switchflag") == 0) { - if (iarg+2 > narg) - error->all(FLERR,"Illegal compute snap command"); - switchflag = atoi(arg[iarg+1]); + } else if (strcmp(arg[iarg], "bzeroflag") == 0) { + if (iarg + 2 > narg) error->all(FLERR, "Illegal compute {} command", style); + bzeroflag = utils::inumeric(FLERR, arg[iarg + 1], false, lmp); iarg += 2; - } else if (strcmp(arg[iarg],"quadraticflag") == 0) { - if (iarg+2 > narg) - error->all(FLERR,"Illegal compute snap command"); - quadraticflag = atoi(arg[iarg+1]); + } else if (strcmp(arg[iarg], "quadraticflag") == 0) { + if (iarg + 2 > narg) error->all(FLERR, "Illegal compute {} command", style); + quadraticflag = utils::inumeric(FLERR, arg[iarg + 1], false, lmp); iarg += 2; - } else if (strcmp(arg[iarg],"chem") == 0) { - if (iarg+2+ntypes > narg) - error->all(FLERR,"Illegal compute snap command"); + } else if (strcmp(arg[iarg], "chem") == 0) { + if (iarg + 2 > narg) error->all(FLERR, "Illegal compute {} command", style); chemflag = 1; - memory->create(map,ntypes+1,"compute_snap:map"); - nelements = utils::inumeric(FLERR,arg[iarg+1],false,lmp); + memory->create(map, ntypes + 1, "compute_sna_grid:map"); + nelements = utils::inumeric(FLERR, arg[iarg + 1], false, lmp); for (int i = 0; i < ntypes; i++) { - int jelem = utils::inumeric(FLERR,arg[iarg+2+i],false,lmp); - if (jelem < 0 || jelem >= nelements) - error->all(FLERR,"Illegal compute snap command"); - map[i+1] = jelem; + int jelem = utils::inumeric(FLERR, arg[iarg + 2 + i], false, lmp); + if (jelem < 0 || jelem >= nelements) error->all(FLERR, "Illegal compute {} command", style); + map[i + 1] = jelem; } - iarg += 2+ntypes; - } else if (strcmp(arg[iarg],"bnormflag") == 0) { - if (iarg+2 > narg) - error->all(FLERR,"Illegal compute snap command"); - bnormflag = atoi(arg[iarg+1]); + iarg += 2 + ntypes; + } else if (strcmp(arg[iarg], "bnormflag") == 0) { + if (iarg + 2 > narg) error->all(FLERR, "Illegal compute {} command", style); + bnormflag = utils::inumeric(FLERR, arg[iarg + 1], false, lmp); iarg += 2; - } else if (strcmp(arg[iarg],"wselfallflag") == 0) { - if (iarg+2 > narg) - error->all(FLERR,"Illegal compute snap command"); - wselfallflag = atoi(arg[iarg+1]); + } else if (strcmp(arg[iarg], "wselfallflag") == 0) { + if (iarg + 2 > narg) error->all(FLERR, "Illegal compute {} command", style); + wselfallflag = utils::inumeric(FLERR, arg[iarg + 1], false, lmp); iarg += 2; - } else if (strcmp(arg[iarg],"bikflag") == 0) { - if (iarg+2 > narg) - error->all(FLERR,"Illegal compute snap command"); - bikflag = atoi(arg[iarg+1]); + } else if (strcmp(arg[iarg], "switchinnerflag") == 0) { + if (iarg + 2 > narg) error->all(FLERR, "Illegal compute {} command", style); + switchinnerflag = utils::inumeric(FLERR, arg[iarg + 1], false, lmp); iarg += 2; - } else if (strcmp(arg[iarg],"switchinnerflag") == 0) { - if (iarg+2 > narg) - error->all(FLERR,"Illegal compute snap command"); - switchinnerflag = atoi(arg[iarg+1]); - iarg += 2; - } else if (strcmp(arg[iarg],"sinner") == 0) { + } else if (strcmp(arg[iarg], "sinner") == 0) { iarg++; - if (iarg+ntypes > narg) - error->all(FLERR,"Illegal compute snap command"); - memory->create(sinnerelem,ntypes+1,"snap:sinnerelem"); + if (iarg + ntypes > narg) error->all(FLERR, "Illegal compute {} command", style); + memory->create(sinnerelem, ntypes + 1, "snap:sinnerelem"); for (int i = 0; i < ntypes; i++) - sinnerelem[i+1] = utils::numeric(FLERR,arg[iarg+i],false,lmp); + sinnerelem[i + 1] = utils::numeric(FLERR, arg[iarg + i], false, lmp); sinnerflag = 1; iarg += ntypes; - } else if (strcmp(arg[iarg],"dinner") == 0) { + } else if (strcmp(arg[iarg], "dinner") == 0) { iarg++; - if (iarg+ntypes > narg) - error->all(FLERR,"Illegal compute snap command"); - memory->create(dinnerelem,ntypes+1,"snap:dinnerelem"); + if (iarg + ntypes > narg) error->all(FLERR, "Illegal compute {} command", style); + memory->create(dinnerelem, ntypes + 1, "snap:dinnerelem"); for (int i = 0; i < ntypes; i++) - dinnerelem[i+1] = utils::numeric(FLERR,arg[iarg+i],false,lmp); + dinnerelem[i + 1] = utils::numeric(FLERR, arg[iarg + i], false, lmp); dinnerflag = 1; iarg += ntypes; - } else error->all(FLERR,"Illegal compute snap command"); + } else + error->all(FLERR, "Illegal compute {} command", style); } if (switchinnerflag && !(sinnerflag && dinnerflag)) - error->all(FLERR,"Illegal compute snap command: switchinnerflag = 1, missing sinner/dinner keyword"); + error->all( + FLERR, + "Illegal compute {} command: switchinnerflag = 1, missing sinner/dinner keyword", + style); if (!switchinnerflag && (sinnerflag || dinnerflag)) - error->all(FLERR,"Illegal compute snap command: switchinnerflag = 0, unexpected sinner/dinner keyword"); + error->all( + FLERR, + "Illegal compute {} command: switchinnerflag = 0, unexpected sinner/dinner keyword", + style); - snaptr = new SNA(lmp, rfac0, twojmax, - rmin0, switchflag, bzeroflag, - chemflag, bnormflag, wselfallflag, - nelements, switchinnerflag); + snaptr = new SNA(lmp, rfac0, twojmax, rmin0, switchflag, bzeroflag, chemflag, bnormflag, + wselfallflag, nelements, switchinnerflag); ncoeff = snaptr->ncoeff; - nperdim = ncoeff; - if (quadraticflag) nperdim += (ncoeff*(ncoeff+1))/2; + nvalues = ncoeff; + if (quadraticflag) nvalues += (ncoeff * (ncoeff + 1)) / 2; + + // end code common to all SNAP computes + ndims_force = 3; ndims_virial = 6; - yoffset = nperdim; - zoffset = 2*nperdim; + yoffset = nvalues; + zoffset = 2*nvalues; natoms = atom->natoms; bik_rows = 1; if (bikflag) bik_rows = natoms; size_array_rows = bik_rows+ndims_force*natoms+ndims_virial; - size_array_cols = nperdim*atom->ntypes+1; + size_array_cols = nvalues*atom->ntypes+1; lastcol = size_array_cols-1; ndims_peratom = ndims_force; - size_peratom = ndims_peratom*nperdim*atom->ntypes; + size_peratom = ndims_peratom*nvalues*atom->ntypes; nmax = 0; } @@ -341,8 +338,8 @@ void ComputeSnap::compute_array() const double radi = radelem[itype]; const int* const jlist = firstneigh[i]; const int jnum = numneigh[i]; - const int typeoffset_local = ndims_peratom*nperdim*(itype-1); - const int typeoffset_global = nperdim*(itype-1); + const int typeoffset_local = ndims_peratom*nvalues*(itype-1); + const int typeoffset_global = nvalues*(itype-1); // insure rij, inside, and typej are of size jnum @@ -481,9 +478,9 @@ void ComputeSnap::compute_array() // accumulate bispectrum force contributions to global array for (int itype = 0; itype < atom->ntypes; itype++) { - const int typeoffset_local = ndims_peratom*nperdim*itype; - const int typeoffset_global = nperdim*itype; - for (int icoeff = 0; icoeff < nperdim; icoeff++) { + const int typeoffset_local = ndims_peratom*nvalues*itype; + const int typeoffset_global = nvalues*itype; + for (int icoeff = 0; icoeff < nvalues; icoeff++) { for (int i = 0; i < ntotal; i++) { double *snadi = snap_peratom[i]+typeoffset_local; int iglobal = atom->tag[i]; @@ -549,10 +546,10 @@ void ComputeSnap::dbdotr_compute() int nall = atom->nlocal + atom->nghost; for (int i = 0; i < nall; i++) for (int itype = 0; itype < atom->ntypes; itype++) { - const int typeoffset_local = ndims_peratom*nperdim*itype; - const int typeoffset_global = nperdim*itype; + const int typeoffset_local = ndims_peratom*nvalues*itype; + const int typeoffset_global = nvalues*itype; double *snadi = snap_peratom[i]+typeoffset_local; - for (int icoeff = 0; icoeff < nperdim; icoeff++) { + for (int icoeff = 0; icoeff < nvalues; icoeff++) { double dbdx = snadi[icoeff]; double dbdy = snadi[icoeff+yoffset]; double dbdz = snadi[icoeff+zoffset]; diff --git a/src/ML-SNAP/compute_snap.h b/src/ML-SNAP/compute_snap.h index bc0670e2c7..9f4162d6c4 100644 --- a/src/ML-SNAP/compute_snap.h +++ b/src/ML-SNAP/compute_snap.h @@ -35,7 +35,7 @@ class ComputeSnap : public Compute { private: int natoms, nmax, size_peratom, lastcol; - int ncoeff, nperdim, yoffset, zoffset; + int ncoeff, nvalues, yoffset, zoffset; int ndims_peratom, ndims_force, ndims_virial; double **cutsq; class NeighList *list; diff --git a/src/ML-SNAP/compute_snav_atom.cpp b/src/ML-SNAP/compute_snav_atom.cpp index 0baa4bc110..949dda8e5c 100644 --- a/src/ML-SNAP/compute_snav_atom.cpp +++ b/src/ML-SNAP/compute_snav_atom.cpp @@ -21,6 +21,7 @@ #include "neighbor.h" #include "neigh_list.h" #include "force.h" +#include "pair.h" #include "comm.h" #include "memory.h" #include "error.h" @@ -33,20 +34,22 @@ ComputeSNAVAtom::ComputeSNAVAtom(LAMMPS *lmp, int narg, char **arg) : Compute(lmp, narg, arg), cutsq(nullptr), list(nullptr), snav(nullptr), radelem(nullptr), wjelem(nullptr), sinnerelem(nullptr), dinnerelem(nullptr) { + + // begin code common to all SNAP computes + double rfac0, rmin0; int twojmax, switchflag, bzeroflag, bnormflag, wselfallflag; int ntypes = atom->ntypes; - int nargmin = 6+2*ntypes; + int nargmin = 6 + 2 * ntypes; - if (narg < nargmin) error->all(FLERR,"Illegal compute snav/atom command"); + if (narg < nargmin) error->all(FLERR, "Illegal compute {} command", style); // default values rmin0 = 0.0; switchflag = 1; bzeroflag = 1; - bnormflag = 0; quadraticflag = 0; chemflag = 0; bnormflag = 0; @@ -56,24 +59,32 @@ ComputeSNAVAtom::ComputeSNAVAtom(LAMMPS *lmp, int narg, char **arg) : // process required arguments - memory->create(radelem,ntypes+1,"snav/atom:radelem"); // offset by 1 to match up with types - memory->create(wjelem,ntypes+1,"snav/atom:wjelem"); - rcutfac = atof(arg[3]); - rfac0 = atof(arg[4]); - twojmax = atoi(arg[5]); + memory->create(radelem, ntypes + 1, "sna/atom:radelem"); // offset by 1 to match up with types + memory->create(wjelem, ntypes + 1, "sna/atom:wjelem"); + + rcutfac = utils::numeric(FLERR, arg[3], false, lmp); + rfac0 = utils::numeric(FLERR, arg[4], false, lmp); + twojmax = utils::inumeric(FLERR, arg[5], false, lmp); + for (int i = 0; i < ntypes; i++) - radelem[i+1] = atof(arg[6+i]); + radelem[i + 1] = + utils::numeric(FLERR, arg[6 + i], false, lmp); for (int i = 0; i < ntypes; i++) - wjelem[i+1] = atof(arg[6+ntypes+i]); + wjelem[i + 1] = + utils::numeric(FLERR, arg[6 + ntypes + i], false, lmp); + // construct cutsq + double cut; - memory->create(cutsq,ntypes+1,ntypes+1,"snav/atom:cutsq"); + cutmax = 0.0; + memory->create(cutsq, ntypes + 1, ntypes + 1, "sna/atom:cutsq"); for (int i = 1; i <= ntypes; i++) { - cut = 2.0*radelem[i]*rcutfac; - cutsq[i][i] = cut*cut; - for (int j = i+1; j <= ntypes; j++) { - cut = (radelem[i]+radelem[j])*rcutfac; - cutsq[i][j] = cutsq[j][i] = cut*cut; + cut = 2.0 * radelem[i] * rcutfac; + if (cut > cutmax) cutmax = cut; + cutsq[i][i] = cut * cut; + for (int j = i + 1; j <= ntypes; j++) { + cut = (radelem[i] + radelem[j]) * rcutfac; + cutsq[i][j] = cutsq[j][i] = cut * cut; } } @@ -87,90 +98,87 @@ ComputeSNAVAtom::ComputeSNAVAtom(LAMMPS *lmp, int narg, char **arg) : int iarg = nargmin; while (iarg < narg) { - if (strcmp(arg[iarg],"rmin0") == 0) { - if (iarg+2 > narg) - error->all(FLERR,"Illegal compute snav/atom command"); - rmin0 = atof(arg[iarg+1]); + if (strcmp(arg[iarg], "rmin0") == 0) { + if (iarg + 2 > narg) error->all(FLERR, "Illegal compute {} command", style); + rmin0 = utils::numeric(FLERR, arg[iarg + 1], false, lmp); iarg += 2; - } else if (strcmp(arg[iarg],"switchflag") == 0) { - if (iarg+2 > narg) - error->all(FLERR,"Illegal compute snav/atom command"); - switchflag = atoi(arg[iarg+1]); + } else if (strcmp(arg[iarg], "switchflag") == 0) { + if (iarg + 2 > narg) error->all(FLERR, "Illegal compute {} command", style); + switchflag = utils::inumeric(FLERR, arg[iarg + 1], false, lmp); iarg += 2; - } else if (strcmp(arg[iarg],"bzeroflag") == 0) { - if (iarg+2 > narg) - error->all(FLERR,"Illegal compute snav/atom command"); - bzeroflag = atoi(arg[iarg+1]); + } else if (strcmp(arg[iarg], "bzeroflag") == 0) { + if (iarg + 2 > narg) error->all(FLERR, "Illegal compute {} command", style); + bzeroflag = utils::inumeric(FLERR, arg[iarg + 1], false, lmp); iarg += 2; - } else if (strcmp(arg[iarg],"quadraticflag") == 0) { - if (iarg+2 > narg) - error->all(FLERR,"Illegal compute snav/atom command"); - quadraticflag = atoi(arg[iarg+1]); + } else if (strcmp(arg[iarg], "quadraticflag") == 0) { + if (iarg + 2 > narg) error->all(FLERR, "Illegal compute {} command", style); + quadraticflag = utils::inumeric(FLERR, arg[iarg + 1], false, lmp); iarg += 2; - } else if (strcmp(arg[iarg],"chem") == 0) { - if (iarg+2+ntypes > narg) - error->all(FLERR,"Illegal compute snav/atom command"); + } else if (strcmp(arg[iarg], "chem") == 0) { + if (iarg + 2 > narg) error->all(FLERR, "Illegal compute {} command", style); chemflag = 1; - memory->create(map,ntypes+1,"compute_sna_atom:map"); - nelements = utils::inumeric(FLERR,arg[iarg+1],false,lmp); + memory->create(map, ntypes + 1, "compute_sna_grid:map"); + nelements = utils::inumeric(FLERR, arg[iarg + 1], false, lmp); for (int i = 0; i < ntypes; i++) { - int jelem = utils::inumeric(FLERR,arg[iarg+2+i],false,lmp); - if (jelem < 0 || jelem >= nelements) - error->all(FLERR,"Illegal compute snav/atom command"); - map[i+1] = jelem; + int jelem = utils::inumeric(FLERR, arg[iarg + 2 + i], false, lmp); + if (jelem < 0 || jelem >= nelements) error->all(FLERR, "Illegal compute {} command", style); + map[i + 1] = jelem; } - iarg += 2+ntypes; - } else if (strcmp(arg[iarg],"bnormflag") == 0) { - if (iarg+2 > narg) - error->all(FLERR,"Illegal compute snav/atom command"); - bnormflag = atoi(arg[iarg+1]); + iarg += 2 + ntypes; + } else if (strcmp(arg[iarg], "bnormflag") == 0) { + if (iarg + 2 > narg) error->all(FLERR, "Illegal compute {} command", style); + bnormflag = utils::inumeric(FLERR, arg[iarg + 1], false, lmp); iarg += 2; - } else if (strcmp(arg[iarg],"wselfallflag") == 0) { - if (iarg+2 > narg) - error->all(FLERR,"Illegal compute snav/atom command"); - wselfallflag = atoi(arg[iarg+1]); + } else if (strcmp(arg[iarg], "wselfallflag") == 0) { + if (iarg + 2 > narg) error->all(FLERR, "Illegal compute {} command", style); + wselfallflag = utils::inumeric(FLERR, arg[iarg + 1], false, lmp); iarg += 2; - } else if (strcmp(arg[iarg],"switchinnerflag") == 0) { - if (iarg+2 > narg) - error->all(FLERR,"Illegal compute snav/atom command"); - switchinnerflag = atoi(arg[iarg+1]); + } else if (strcmp(arg[iarg], "switchinnerflag") == 0) { + if (iarg + 2 > narg) error->all(FLERR, "Illegal compute {} command", style); + switchinnerflag = utils::inumeric(FLERR, arg[iarg + 1], false, lmp); iarg += 2; - } else if (strcmp(arg[iarg],"sinner") == 0) { + } else if (strcmp(arg[iarg], "sinner") == 0) { iarg++; - if (iarg+ntypes > narg) - error->all(FLERR,"Illegal compute snav/atom command"); - memory->create(sinnerelem,ntypes+1,"snav/atom:sinnerelem"); + if (iarg + ntypes > narg) error->all(FLERR, "Illegal compute {} command", style); + memory->create(sinnerelem, ntypes + 1, "snap:sinnerelem"); for (int i = 0; i < ntypes; i++) - sinnerelem[i+1] = utils::numeric(FLERR,arg[iarg+i],false,lmp); + sinnerelem[i + 1] = utils::numeric(FLERR, arg[iarg + i], false, lmp); sinnerflag = 1; iarg += ntypes; - } else if (strcmp(arg[iarg],"dinner") == 0) { + } else if (strcmp(arg[iarg], "dinner") == 0) { iarg++; - if (iarg+ntypes > narg) - error->all(FLERR,"Illegal compute snav/atom command"); - memory->create(dinnerelem,ntypes+1,"snav/atom:dinnerelem"); + if (iarg + ntypes > narg) error->all(FLERR, "Illegal compute {} command", style); + memory->create(dinnerelem, ntypes + 1, "snap:dinnerelem"); for (int i = 0; i < ntypes; i++) - dinnerelem[i+1] = utils::numeric(FLERR,arg[iarg+i],false,lmp); + dinnerelem[i + 1] = utils::numeric(FLERR, arg[iarg + i], false, lmp); dinnerflag = 1; iarg += ntypes; - } else error->all(FLERR,"Illegal compute snav/atom command"); + } else + error->all(FLERR, "Illegal compute {} command", style); } if (switchinnerflag && !(sinnerflag && dinnerflag)) - error->all(FLERR,"Illegal compute snav/atom command: switchinnerflag = 1, missing sinner/dinner keyword"); + error->all( + FLERR, + "Illegal compute {} command: switchinnerflag = 1, missing sinner/dinner keyword", + style); if (!switchinnerflag && (sinnerflag || dinnerflag)) - error->all(FLERR,"Illegal compute snav/atom command: switchinnerflag = 0, unexpected sinner/dinner keyword"); + error->all( + FLERR, + "Illegal compute {} command: switchinnerflag = 0, unexpected sinner/dinner keyword", + style); - snaptr = new SNA(lmp, rfac0, twojmax, - rmin0, switchflag, bzeroflag, - chemflag, bnormflag, wselfallflag, - nelements, switchinnerflag); + snaptr = new SNA(lmp, rfac0, twojmax, rmin0, switchflag, bzeroflag, chemflag, bnormflag, + wselfallflag, nelements, switchinnerflag); ncoeff = snaptr->ncoeff; - nperdim = ncoeff; - if (quadraticflag) nperdim += (ncoeff*(ncoeff+1))/2; - size_peratom_cols = 6*nperdim*atom->ntypes; + nvalues = ncoeff; + if (quadraticflag) nvalues += (ncoeff * (ncoeff + 1)) / 2; + + // end code common to all SNAP computes + + size_peratom_cols = 6*nvalues*atom->ntypes; comm_reverse = size_peratom_cols; peratom_flag = 1; @@ -203,10 +211,9 @@ void ComputeSNAVAtom::init() { if (force->pair == nullptr) error->all(FLERR,"Compute snav/atom requires a pair style be defined"); - // TODO: Not sure what to do with this error check since cutoff radius is not - // a single number - //if (sqrt(cutsq) > force->pair->cutforce) - // error->all(FLERR,"Compute snav/atom cutoff is longer than pairwise cutoff"); + + if (cutmax > force->pair->cutforce) + error->all(FLERR,"Compute snav/atom cutoff is longer than pairwise cutoff"); // need an occasional full neighbor list @@ -280,7 +287,7 @@ void ComputeSNAVAtom::compute_peratom() const int* const jlist = firstneigh[i]; const int jnum = numneigh[i]; - const int typeoffset = 6*nperdim*(atom->type[i]-1); + const int typeoffset = 6*nvalues*(atom->type[i]-1); // insure rij, inside, and typej are of size jnum @@ -339,17 +346,17 @@ void ComputeSNAVAtom::compute_peratom() for (int icoeff = 0; icoeff < ncoeff; icoeff++) { snavi[icoeff] += snaptr->dblist[icoeff][0]*xtmp; - snavi[icoeff+nperdim] += snaptr->dblist[icoeff][1]*ytmp; - snavi[icoeff+2*nperdim] += snaptr->dblist[icoeff][2]*ztmp; - snavi[icoeff+3*nperdim] += snaptr->dblist[icoeff][1]*ztmp; - snavi[icoeff+4*nperdim] += snaptr->dblist[icoeff][0]*ztmp; - snavi[icoeff+5*nperdim] += snaptr->dblist[icoeff][0]*ytmp; + snavi[icoeff+nvalues] += snaptr->dblist[icoeff][1]*ytmp; + snavi[icoeff+2*nvalues] += snaptr->dblist[icoeff][2]*ztmp; + snavi[icoeff+3*nvalues] += snaptr->dblist[icoeff][1]*ztmp; + snavi[icoeff+4*nvalues] += snaptr->dblist[icoeff][0]*ztmp; + snavi[icoeff+5*nvalues] += snaptr->dblist[icoeff][0]*ytmp; snavj[icoeff] -= snaptr->dblist[icoeff][0]*x[j][0]; - snavj[icoeff+nperdim] -= snaptr->dblist[icoeff][1]*x[j][1]; - snavj[icoeff+2*nperdim] -= snaptr->dblist[icoeff][2]*x[j][2]; - snavj[icoeff+3*nperdim] -= snaptr->dblist[icoeff][1]*x[j][2]; - snavj[icoeff+4*nperdim] -= snaptr->dblist[icoeff][0]*x[j][2]; - snavj[icoeff+5*nperdim] -= snaptr->dblist[icoeff][0]*x[j][1]; + snavj[icoeff+nvalues] -= snaptr->dblist[icoeff][1]*x[j][1]; + snavj[icoeff+2*nvalues] -= snaptr->dblist[icoeff][2]*x[j][2]; + snavj[icoeff+3*nvalues] -= snaptr->dblist[icoeff][1]*x[j][2]; + snavj[icoeff+4*nvalues] -= snaptr->dblist[icoeff][0]*x[j][2]; + snavj[icoeff+5*nvalues] -= snaptr->dblist[icoeff][0]*x[j][1]; } if (quadraticflag) { @@ -369,17 +376,17 @@ void ComputeSNAVAtom::compute_peratom() double dbytmp = bi*biy; double dbztmp = bi*biz; snavi[ncount] += dbxtmp*xtmp; - snavi[ncount+nperdim] += dbytmp*ytmp; - snavi[ncount+2*nperdim] += dbztmp*ztmp; - snavi[ncount+3*nperdim] += dbytmp*ztmp; - snavi[ncount+4*nperdim] += dbxtmp*ztmp; - snavi[ncount+5*nperdim] += dbxtmp*ytmp; + snavi[ncount+nvalues] += dbytmp*ytmp; + snavi[ncount+2*nvalues] += dbztmp*ztmp; + snavi[ncount+3*nvalues] += dbytmp*ztmp; + snavi[ncount+4*nvalues] += dbxtmp*ztmp; + snavi[ncount+5*nvalues] += dbxtmp*ytmp; snavj[ncount] -= dbxtmp*x[j][0]; - snavj[ncount+nperdim] -= dbytmp*x[j][1]; - snavj[ncount+2*nperdim] -= dbztmp*x[j][2]; - snavj[ncount+3*nperdim] -= dbytmp*x[j][2]; - snavj[ncount+4*nperdim] -= dbxtmp*x[j][2]; - snavj[ncount+5*nperdim] -= dbxtmp*x[j][1]; + snavj[ncount+nvalues] -= dbytmp*x[j][1]; + snavj[ncount+2*nvalues] -= dbztmp*x[j][2]; + snavj[ncount+3*nvalues] -= dbytmp*x[j][2]; + snavj[ncount+4*nvalues] -= dbxtmp*x[j][2]; + snavj[ncount+5*nvalues] -= dbxtmp*x[j][1]; ncount++; // upper-triangular elements of quadratic matrix @@ -392,17 +399,17 @@ void ComputeSNAVAtom::compute_peratom() double dbztmp = bi*snaptr->dblist[jcoeff][2] + biz*snaptr->blist[jcoeff]; snavi[ncount] += dbxtmp*xtmp; - snavi[ncount+nperdim] += dbytmp*ytmp; - snavi[ncount+2*nperdim] += dbztmp*ztmp; - snavi[ncount+3*nperdim] += dbytmp*ztmp; - snavi[ncount+4*nperdim] += dbxtmp*ztmp; - snavi[ncount+5*nperdim] += dbxtmp*ytmp; + snavi[ncount+nvalues] += dbytmp*ytmp; + snavi[ncount+2*nvalues] += dbztmp*ztmp; + snavi[ncount+3*nvalues] += dbytmp*ztmp; + snavi[ncount+4*nvalues] += dbxtmp*ztmp; + snavi[ncount+5*nvalues] += dbxtmp*ytmp; snavj[ncount] -= dbxtmp*x[j][0]; - snavj[ncount+nperdim] -= dbytmp*x[j][1]; - snavj[ncount+2*nperdim] -= dbztmp*x[j][2]; - snavj[ncount+3*nperdim] -= dbytmp*x[j][2]; - snavj[ncount+4*nperdim] -= dbxtmp*x[j][2]; - snavj[ncount+5*nperdim] -= dbxtmp*x[j][1]; + snavj[ncount+nvalues] -= dbytmp*x[j][1]; + snavj[ncount+2*nvalues] -= dbztmp*x[j][2]; + snavj[ncount+3*nvalues] -= dbytmp*x[j][2]; + snavj[ncount+4*nvalues] -= dbxtmp*x[j][2]; + snavj[ncount+5*nvalues] -= dbxtmp*x[j][1]; ncount++; } } diff --git a/src/ML-SNAP/compute_snav_atom.h b/src/ML-SNAP/compute_snav_atom.h index 1eb84d8df7..2eac0b7b28 100644 --- a/src/ML-SNAP/compute_snav_atom.h +++ b/src/ML-SNAP/compute_snav_atom.h @@ -37,7 +37,7 @@ class ComputeSNAVAtom : public Compute { private: int nmax; - int ncoeff, nperdim; + int ncoeff, nvalues; double **cutsq; class NeighList *list; double **snav; @@ -50,6 +50,7 @@ class ComputeSNAVAtom : public Compute { double *sinnerelem; double *dinnerelem; class SNA *snaptr; + double cutmax; int quadraticflag; }; diff --git a/src/Makefile b/src/Makefile index ab8e5c4fea..649a1ad002 100644 --- a/src/Makefile +++ b/src/Makefile @@ -461,7 +461,7 @@ mpi-stubs: sinclude ../lib/python/Makefile.lammps install-python: @rm -rf ../python/build - @$(PYTHON) ../python/install.py -p ../python/lammps -l ../src/liblammps.so + @$(PYTHON) ../python/install.py -p ../python/lammps -l ../src/liblammps.so -w $(PWD) # Create a tarball of src dir and packages diff --git a/src/REAXFF/fix_acks2_reaxff.cpp b/src/REAXFF/fix_acks2_reaxff.cpp index bd93dec0b7..be71bc761f 100644 --- a/src/REAXFF/fix_acks2_reaxff.cpp +++ b/src/REAXFF/fix_acks2_reaxff.cpp @@ -94,8 +94,8 @@ FixACKS2ReaxFF::~FixACKS2ReaxFF() memory->destroy(s_hist_X); memory->destroy(s_hist_last); - deallocate_storage(); - deallocate_matrix(); + FixACKS2ReaxFF::deallocate_storage(); + FixACKS2ReaxFF::deallocate_matrix(); } /* ---------------------------------------------------------------------- */ diff --git a/src/REAXFF/fix_qeq_reaxff.cpp b/src/REAXFF/fix_qeq_reaxff.cpp index aeeee7b71a..48e93f682a 100644 --- a/src/REAXFF/fix_qeq_reaxff.cpp +++ b/src/REAXFF/fix_qeq_reaxff.cpp @@ -163,7 +163,7 @@ FixQEqReaxFF::~FixQEqReaxFF() memory->destroy(t_hist); FixQEqReaxFF::deallocate_storage(); - deallocate_matrix(); + FixQEqReaxFF::deallocate_matrix(); memory->destroy(shld); @@ -640,7 +640,7 @@ void FixQEqReaxFF::compute_H() int jnum; int i, j, ii, jj, flag; double dx, dy, dz, r_sqr; - const double SMALL = 0.0001; + constexpr double EPSILON = 0.0001; int *type = atom->type; tagint *tag = atom->tag; @@ -671,10 +671,10 @@ void FixQEqReaxFF::compute_H() if (j < atom->nlocal) flag = 1; else if (tag[i] < tag[j]) flag = 1; else if (tag[i] == tag[j]) { - if (dz > SMALL) flag = 1; - else if (fabs(dz) < SMALL) { - if (dy > SMALL) flag = 1; - else if (fabs(dy) < SMALL && dx > SMALL) + if (dz > EPSILON) flag = 1; + else if (fabs(dz) < EPSILON) { + if (dy > EPSILON) flag = 1; + else if (fabs(dy) < EPSILON && dx > EPSILON) flag = 1; } } diff --git a/src/create_bonds.cpp b/src/create_bonds.cpp index 7b3929db62..d316644e15 100644 --- a/src/create_bonds.cpp +++ b/src/create_bonds.cpp @@ -73,7 +73,6 @@ void CreateBonds::command(int narg, char **arg) iarg = 6; } else if (strcmp(arg[0], "single/bond") == 0) { style = SBOND; - if (narg < 4) error->all(FLERR, "Illegal create_bonds command"); btype = utils::inumeric(FLERR, arg[1], false, lmp); batom1 = utils::tnumeric(FLERR, arg[2], false, lmp); batom2 = utils::tnumeric(FLERR, arg[3], false, lmp); diff --git a/src/dump.cpp b/src/dump.cpp index 3569d32165..a354be4f04 100644 --- a/src/dump.cpp +++ b/src/dump.cpp @@ -517,6 +517,8 @@ void Dump::write() if (refreshflag) modify->compute[irefresh]->refresh(); + if (filewriter && fp != nullptr) write_footer(); + // if file per timestep, close file if I am filewriter if (multifile) { diff --git a/src/dump.h b/src/dump.h index a95029ada8..90866a3567 100644 --- a/src/dump.h +++ b/src/dump.h @@ -153,6 +153,8 @@ class Dump : protected Pointers { virtual void pack(tagint *) = 0; virtual int convert_string(int, double *) { return 0; } virtual void write_data(int, double *) = 0; + virtual void write_footer() {} + void pbc_allocate(); double compute_time(); diff --git a/src/fmt/core.h b/src/fmt/core.h index 8444cd9546..2fafa777ba 100644 --- a/src/fmt/core.h +++ b/src/fmt/core.h @@ -315,7 +315,9 @@ // Enable minimal optimizations for more compact code in debug mode. FMT_GCC_PRAGMA("GCC push_options") -#ifndef __OPTIMIZE__ +// LAMMPS CUSTOMIZATION: suppress warning about pragma with KOKKOS +#if !defined(__OPTIMIZE__) && !defined(LMP_KOKKOS) +// END LAMMPS CUSTOMIZATION FMT_GCC_PRAGMA("GCC optimize(\"Og\")") #endif diff --git a/src/library.cpp b/src/library.cpp index 079534d663..cdc0a812c5 100644 --- a/src/library.cpp +++ b/src/library.cpp @@ -1659,14 +1659,18 @@ lists the available options. - LMP_TYPE_ARRAY - ``double **`` - Local data array + * - LMP_STYLE_LOCAL + - LMP_SIZE_VECTOR + - ``int *`` + - Alias for using LMP_SIZE_ROWS * - LMP_STYLE_LOCAL - LMP_SIZE_ROWS - ``int *`` - - Number of local data rows + - Number of local array rows or length of vector * - LMP_STYLE_LOCAL - LMP_SIZE_COLS - ``int *`` - - Number of local data columns + - Number of local array columns, 0 if vector .. warning:: @@ -1749,6 +1753,7 @@ void *lammps_extract_compute(void *handle, const char *id, int style, int type) if (type == LMP_TYPE_SCALAR) return (void *) &compute->size_local_rows; /* for backward compatibility */ if (type == LMP_TYPE_VECTOR) return (void *) compute->vector_local; if (type == LMP_TYPE_ARRAY) return (void *) compute->array_local; + if (type == LMP_SIZE_VECTOR) return (void *) &compute->size_local_rows; /* alias for LMP_SIZE_ROWS */ if (type == LMP_SIZE_ROWS) return (void *) &compute->size_local_rows; if (type == LMP_SIZE_COLS) return (void *) &compute->size_local_cols; } diff --git a/src/read_restart.cpp b/src/read_restart.cpp index b42dcce1ca..a7fb5a828d 100644 --- a/src/read_restart.cpp +++ b/src/read_restart.cpp @@ -61,9 +61,10 @@ void ReadRestart::command(int narg, char **arg) // check for remap option - int remapflag = 0; + int remapflag = 1; if (narg == 2) { - if (strcmp(arg[1],"remap") == 0) remapflag = 1; + if (strcmp(arg[1],"noremap") == 0) remapflag = 0; + else if (strcmp(arg[1],"remap") == 0) remapflag = 1; // for backward compatibility else error->all(FLERR,"Illegal read_restart command"); } @@ -643,8 +644,7 @@ void ReadRestart::header() int dimension = read_int(); domain->dimension = dimension; if (domain->dimension == 2 && domain->zperiodic == 0) - error->all(FLERR, - "Cannot run 2d simulation with non-periodic Z dimension"); + error->all(FLERR, "Cannot run 2d simulation with non-periodic Z dimension"); // read nprocs from restart file, warn if different diff --git a/src/reader.cpp b/src/reader.cpp index c7b99260e7..025445ca24 100644 --- a/src/reader.cpp +++ b/src/reader.cpp @@ -28,6 +28,12 @@ Reader::Reader(LAMMPS *lmp) : Pointers(lmp) compressed = false; } +// avoid resource leak +Reader::~Reader() +{ + if (fp != nullptr) close_file(); +} + /* ---------------------------------------------------------------------- try to open given file generic version for ASCII files with optional compression or for native binary dumps diff --git a/src/reader.h b/src/reader.h index 93f9197b7a..753b6956f7 100644 --- a/src/reader.h +++ b/src/reader.h @@ -23,6 +23,7 @@ namespace LAMMPS_NS { class Reader : protected Pointers { public: Reader(class LAMMPS *); + ~Reader() override; virtual void settings(int, char **); diff --git a/src/version.h b/src/version.h index 60ed7ad75f..4b95d1c910 100644 --- a/src/version.h +++ b/src/version.h @@ -1 +1 @@ -#define LAMMPS_VERSION "2 Jun 2022" +#define LAMMPS_VERSION "23 Jun 2022" diff --git a/tools/msi2lmp/src/msi2lmp.c b/tools/msi2lmp/src/msi2lmp.c index 87cffb4e6e..3136678023 100644 --- a/tools/msi2lmp/src/msi2lmp.c +++ b/tools/msi2lmp/src/msi2lmp.c @@ -66,9 +66,9 @@ * The program is started by supplying information at the command prompt * according to the usage described below. * -* USAGE: msi2lmp3 ROOTNAME {-print #} {-class #} {-frc FRC_FILE} {-ignore} {-nocenter} {-oldstyle} +* USAGE: msi2lmp ROOTNAME {-print #} {-class #} {-frc FRC_FILE} {-ignore} {-nocenter} {-oldstyle} * -* -- msi2lmp3 is the name of the executable +* -- msi2lmp is the name of the executable * -- ROOTNAME is the base name of the .car and .mdf files * -- all opther flags are optional and can be abbreviated (e.g. -p instead of -print) * diff --git a/tools/singularity/README.md b/tools/singularity/README.md index db7aa9e3b0..4700dac6ec 100644 --- a/tools/singularity/README.md +++ b/tools/singularity/README.md @@ -1,17 +1,19 @@ -# Singularity container definitions for compiling/testing LAMMPS +# Apptainer (aka Singularity) container definitions for compiling/testing LAMMPS The *.def files in this folder can be used to build container images -for [Singularity](https://sylabs.io), suitable for compiling and testing +for [Apptainer](https://apptainer.org) (previously called +[Singularity](https://sylabs.io)), suitable for compiling and testing LAMMPS on a variety of OS variants with support for most standard packages and - for some of them - also building/spellchecking the manual -in all supported formats. This allows to test and debug LAMMPS code on +in all supported formats. This allows to test and debug LAMMPS code on different OS variants without doing a full installation on your development workstation, e.g. when bugs are reported that can only be reproduced on a specific OS or with specific (mostly older) versions of tools, compilers, or libraries. Here is a workflow for testing a compilation of LAMMPS with a locally -built CentOS 7.x singularity container. +built CentOS 7.x Singularity container. For Apptainer replace the +`singularity` command with `apptainer`. ``` cd some/work/directory diff --git a/tools/singularity/ubuntu18.04_amd_rocm.def b/tools/singularity/ubuntu18.04_amd_rocm.def index ceedfd8144..febdf1172e 100644 --- a/tools/singularity/ubuntu18.04_amd_rocm.def +++ b/tools/singularity/ubuntu18.04_amd_rocm.def @@ -3,7 +3,7 @@ From: ubuntu:18.04 %environment export PATH=/usr/lib/ccache:${PATH}:/opt/rocm/bin:/opt/rocm/profiler/bin:/opt/rocm/opencl/bin/x86_64 - export LD_LIBRARY_PATH=$LD_LIBRARY_PATH:/opt/rocm/lib:/opt/rocm-5.1.2/llvm/lib + export LD_LIBRARY_PATH=$LD_LIBRARY_PATH:/opt/rocm/lib:/opt/rocm-5.1.3/llvm/lib %post export DEBIAN_FRONTEND=noninteractive apt-get update @@ -22,10 +22,10 @@ From: ubuntu:18.04 apt install -y cmake ########################################################################### - # ROCm 5.1.2 + # ROCm 5.1.3 ########################################################################### - wget https://repo.radeon.com/amdgpu-install/22.10.2/ubuntu/bionic/amdgpu-install_22.10.2.50102-1_all.deb - apt-get install -y ./amdgpu-install_22.10.2.50102-1_all.deb + wget https://repo.radeon.com/amdgpu-install/22.10.3/ubuntu/bionic/amdgpu-install_22.10.3.50103-1_all.deb + apt-get install -y ./amdgpu-install_22.10.3.50103-1_all.deb apt-get update apt-get install --no-install-recommends -y \ diff --git a/tools/singularity/ubuntu18.04_gpu.def b/tools/singularity/ubuntu18.04_gpu.def index 6aa37ccf84..90034285ff 100644 --- a/tools/singularity/ubuntu18.04_gpu.def +++ b/tools/singularity/ubuntu18.04_gpu.def @@ -2,11 +2,11 @@ BootStrap: docker From: ubuntu:18.04 %environment - export PATH=/usr/lib/ccache:/usr/local/cuda-11.5/bin:${PATH}:/opt/rocm/bin:/opt/rocm/profiler/bin:/opt/rocm/opencl/bin/x86_64 - export CUDADIR=/usr/local/cuda-11.5 - export CUDA_PATH=/usr/local/cuda-11.5 - export LD_LIBRARY_PATH=$LD_LIBRARY_PATH:/usr/local/cuda-11.5/lib64:/opt/rocm/lib:/opt/rocm-4.5.0/llvm/lib - export LIBRARY_PATH=/usr/local/cuda-11.5/lib64/stubs + export PATH=/usr/lib/ccache:/usr/local/cuda-11.7/bin:${PATH}:/opt/rocm/bin:/opt/rocm/profiler/bin:/opt/rocm/opencl/bin/x86_64 + export CUDADIR=/usr/local/cuda-11.7 + export CUDA_PATH=/usr/local/cuda-11.7 + export LD_LIBRARY_PATH=$LD_LIBRARY_PATH:/usr/local/cuda-11.7/lib64:/opt/rocm/lib:/opt/rocm-5.1.3/llvm/lib + export LIBRARY_PATH=/usr/local/cuda-11.7/lib64/stubs %post export DEBIAN_FRONTEND=noninteractive apt-get update @@ -27,8 +27,8 @@ From: ubuntu:18.04 ########################################################################### # ROCm 4.5 ########################################################################### - wget https://repo.radeon.com/amdgpu-install/21.40/ubuntu/focal/amdgpu-install-21.40.40500-1_all.deb - apt-get install -y ./amdgpu-install-21.40.40500-1_all.deb + wget https://repo.radeon.com/amdgpu-install/22.10.3/ubuntu/focal/amdgpu-install_22.10.3.50103-1_all.deb + apt-get install -y ./amdgpu-install_22.10.3.50103-1_all.deb apt-get update apt-get install --no-install-recommends -y \ @@ -122,11 +122,11 @@ From: ubuntu:18.04 wget https://developer.download.nvidia.com/compute/cuda/repos/ubuntu1804/x86_64/cuda-ubuntu1804.pin mv cuda-ubuntu1804.pin /etc/apt/preferences.d/cuda-repository-pin-600 - apt-key adv --fetch-keys https://developer.download.nvidia.com/compute/cuda/repos/ubuntu1804/x86_64/7fa2af80.pub + apt-key adv --fetch-keys https://developer.download.nvidia.com/compute/cuda/repos/ubuntu1804/x86_64/3bf863cc.pub add-apt-repository "deb http://developer.download.nvidia.com/compute/cuda/repos/ubuntu1804/x86_64/ /" apt-get update - export CUDA_PKG_VERSION=11.5 + export CUDA_PKG_VERSION=11.7 apt-get install -y --no-install-recommends \ cuda-libraries-${CUDA_PKG_VERSION} \ diff --git a/tools/singularity/ubuntu20.04_amd_rocm.def b/tools/singularity/ubuntu20.04_amd_rocm.def index 5e351e49a8..6034014370 100644 --- a/tools/singularity/ubuntu20.04_amd_rocm.def +++ b/tools/singularity/ubuntu20.04_amd_rocm.def @@ -3,7 +3,7 @@ From: ubuntu:20.04 %environment export PATH=/usr/lib/ccache:${PATH}:/opt/rocm/bin:/opt/rocm/profiler/bin:/opt/rocm/opencl/bin/x86_64 - export LD_LIBRARY_PATH=$LD_LIBRARY_PATH:/opt/rocm/lib:/opt/rocm-5.1.2/llvm/lib + export LD_LIBRARY_PATH=$LD_LIBRARY_PATH:/opt/rocm/lib:/opt/rocm-5.1.3/llvm/lib %post export DEBIAN_FRONTEND=noninteractive apt-get update @@ -13,10 +13,10 @@ From: ubuntu:20.04 apt-get install --no-install-recommends -y software-properties-common ########################################################################### - # ROCm 5.1.2 + # ROCm 5.1.3 ########################################################################### - wget https://repo.radeon.com/amdgpu-install/22.10.2/ubuntu/focal/amdgpu-install_22.10.2.50102-1_all.deb - apt-get install -y ./amdgpu-install_22.10.2.50102-1_all.deb + wget https://repo.radeon.com/amdgpu-install/22.10.3/ubuntu/focal/amdgpu-install_22.10.3.50103-1_all.deb + apt-get install -y ./amdgpu-install_22.10.3.50103-1_all.deb apt-get update apt-get install --no-install-recommends -y \ diff --git a/tools/singularity/ubuntu20.04_gpu.def b/tools/singularity/ubuntu20.04_gpu.def index 23bddeb14f..f7bab1ee9d 100644 --- a/tools/singularity/ubuntu20.04_gpu.def +++ b/tools/singularity/ubuntu20.04_gpu.def @@ -2,11 +2,11 @@ BootStrap: docker From: ubuntu:20.04 %environment - export PATH=/usr/lib/ccache:/usr/local/cuda-11.5/bin:${PATH}:/opt/rocm/bin:/opt/rocm/profiler/bin:/opt/rocm/opencl/bin/x86_64 - export CUDADIR=/usr/local/cuda-11.5 - export CUDA_PATH=/usr/local/cuda-11.5 - export LD_LIBRARY_PATH=$LD_LIBRARY_PATH:/usr/local/cuda-11.5/lib64:/opt/rocm/lib:/opt/rocm-4.5.0/llvm/lib - export LIBRARY_PATH=/usr/local/cuda-11.5/lib64/stubs + export PATH=/usr/lib/ccache:/usr/local/cuda-11.7/bin:${PATH}:/opt/rocm/bin:/opt/rocm/profiler/bin:/opt/rocm/opencl/bin/x86_64 + export CUDADIR=/usr/local/cuda-11.7 + export CUDA_PATH=/usr/local/cuda-11.7 + export LD_LIBRARY_PATH=$LD_LIBRARY_PATH:/usr/local/cuda-11.7/lib64:/opt/rocm/lib:/opt/rocm-5.1.3/llvm/lib + export LIBRARY_PATH=/usr/local/cuda-11.7/lib64/stubs %post export DEBIAN_FRONTEND=noninteractive apt-get update @@ -15,10 +15,10 @@ From: ubuntu:20.04 apt-get install -y --no-install-recommends curl wget libnuma-dev gnupg ca-certificates ########################################################################### - # ROCm 4.5 + # ROCm 5.1.3 ########################################################################### - wget https://repo.radeon.com/amdgpu-install/21.40/ubuntu/focal/amdgpu-install-21.40.40500-1_all.deb - apt-get install -y ./amdgpu-install-21.40.40500-1_all.deb + wget https://repo.radeon.com/amdgpu-install/22.10.3/ubuntu/focal/amdgpu-install_22.10.3.50103-1_all.deb + apt-get install -y ./amdgpu-install_22.10.3.50103-1_all.deb apt-get update apt-get install --no-install-recommends -y \ @@ -109,7 +109,7 @@ From: ubuntu:20.04 wget https://developer.download.nvidia.com/compute/cuda/repos/ubuntu2004/x86_64/cuda-ubuntu2004.pin mv cuda-ubuntu2004.pin /etc/apt/preferences.d/cuda-repository-pin-600 - apt-key adv --fetch-keys https://developer.download.nvidia.com/compute/cuda/repos/ubuntu2004/x86_64/7fa2af80.pub + apt-key adv --fetch-keys https://developer.download.nvidia.com/compute/cuda/repos/ubuntu2004/x86_64/3bf863cc.pub add-apt-repository "deb http://developer.download.nvidia.com/compute/cuda/repos/ubuntu2004/x86_64/ /" apt-get update diff --git a/tools/valgrind/OpenMP.supp b/tools/valgrind/OpenMP.supp index e1870e668c..15531fa8ec 100644 --- a/tools/valgrind/OpenMP.supp +++ b/tools/valgrind/OpenMP.supp @@ -134,3 +134,31 @@ fun:GOMP_parallel obj:* } +{ + OpnMP_open_part1 + Memcheck:Leak + match-leak-kinds: reachable + fun:malloc + ... + fun:openaux + ... + fun:dl_open_worker_begin + ... + fun:dl_open_worker + ... + fun:_dl_open +} +{ + OpnMP_open_part2 + Memcheck:Leak + match-leak-kinds: reachable + fun:calloc + ... + fun:openaux + ... + fun:dl_open_worker_begin + ... + fun:dl_open_worker + ... + fun:_dl_open +} diff --git a/unittest/cplusplus/test_lammps_class.cpp b/unittest/cplusplus/test_lammps_class.cpp index d25bf43656..64ce1eefb1 100644 --- a/unittest/cplusplus/test_lammps_class.cpp +++ b/unittest/cplusplus/test_lammps_class.cpp @@ -83,8 +83,8 @@ TEST_F(LAMMPS_plain, InitMembers) EXPECT_STREQ(lmp->exename, "LAMMPS_test"); EXPECT_EQ(lmp->num_package, 0); - EXPECT_EQ(lmp->clientserver, 0); + EXPECT_EQ(lmp->mdicomm, nullptr); EXPECT_EQ(lmp->kokkos, nullptr); EXPECT_EQ(lmp->atomKK, nullptr); EXPECT_EQ(lmp->memoryKK, nullptr); @@ -218,8 +218,8 @@ TEST_F(LAMMPS_omp, InitMembers) EXPECT_STREQ(lmp->exename, "LAMMPS_test"); EXPECT_EQ(lmp->num_package, 1); - EXPECT_EQ(lmp->clientserver, 0); + EXPECT_EQ(lmp->mdicomm, nullptr); EXPECT_EQ(lmp->kokkos, nullptr); EXPECT_EQ(lmp->atomKK, nullptr); EXPECT_EQ(lmp->memoryKK, nullptr); @@ -302,12 +302,12 @@ TEST_F(LAMMPS_kokkos, InitMembers) EXPECT_STREQ(lmp->exename, "LAMMPS_test"); EXPECT_EQ(lmp->num_package, 0); - EXPECT_EQ(lmp->clientserver, 0); if (Info::has_accelerator_feature("KOKKOS", "api", "openmp")) EXPECT_EQ(lmp->comm->nthreads, 2); else EXPECT_EQ(lmp->comm->nthreads, 1); + EXPECT_EQ(lmp->mdicomm, nullptr); EXPECT_NE(lmp->kokkos, nullptr); EXPECT_NE(lmp->atomKK, nullptr); EXPECT_NE(lmp->memoryKK, nullptr); diff --git a/unittest/force-styles/tests/atomic-pair-sw_angle_table.yaml b/unittest/force-styles/tests/atomic-pair-sw_angle_table.yaml index b8eae7848b..7229a3cd26 100644 --- a/unittest/force-styles/tests/atomic-pair-sw_angle_table.yaml +++ b/unittest/force-styles/tests/atomic-pair-sw_angle_table.yaml @@ -1,6 +1,5 @@ --- lammps_version: 4 May 2022 -tags: generated date_generated: Wed Jun 1 15:17:22 2022 epsilon: 1e-12 skip_tests: single diff --git a/unittest/force-styles/tests/atomic-pair-threebody_table.yaml b/unittest/force-styles/tests/atomic-pair-threebody_table.yaml index 605f18a9a6..bb3c8c1c1b 100644 --- a/unittest/force-styles/tests/atomic-pair-threebody_table.yaml +++ b/unittest/force-styles/tests/atomic-pair-threebody_table.yaml @@ -1,6 +1,5 @@ --- lammps_version: 4 May 2022 -tags: generated date_generated: Wed Jun 1 15:28:13 2022 epsilon: 1e-05 skip_tests: single diff --git a/unittest/force-styles/tests/manybody-pair-pace_product.yaml b/unittest/force-styles/tests/manybody-pair-pace_product.yaml index 9aac9ddcae..b178586ab2 100644 --- a/unittest/force-styles/tests/manybody-pair-pace_product.yaml +++ b/unittest/force-styles/tests/manybody-pair-pace_product.yaml @@ -1,7 +1,7 @@ --- lammps_version: 17 Feb 2022 date_generated: Fri Mar 18 22:17:48 2022 -epsilon: 5e-13 +epsilon: 5e-12 skip_tests: prerequisites: ! | pair pace diff --git a/unittest/force-styles/tests/manybody-pair-pace_recursive.yaml b/unittest/force-styles/tests/manybody-pair-pace_recursive.yaml index eef7509606..6105debb67 100644 --- a/unittest/force-styles/tests/manybody-pair-pace_recursive.yaml +++ b/unittest/force-styles/tests/manybody-pair-pace_recursive.yaml @@ -1,7 +1,7 @@ --- lammps_version: 10 Mar 2021 date_generated: Wed Apr 7 19:30:07 2021 -epsilon: 5e-13 +epsilon: 5e-12 prerequisites: ! | pair pace pre_commands: ! | diff --git a/unittest/python/python-commands.py b/unittest/python/python-commands.py index 58a0772510..cd77fa21a1 100644 --- a/unittest/python/python-commands.py +++ b/unittest/python/python-commands.py @@ -1,6 +1,8 @@ import sys,os,unittest,ctypes -from lammps import lammps, LMP_VAR_ATOM, LMP_STYLE_GLOBAL, LMP_TYPE_VECTOR, LAMMPS_DOUBLE_2D, LAMMPS_AUTODETECT +from lammps import lammps, LMP_VAR_ATOM, LMP_STYLE_GLOBAL, LMP_STYLE_LOCAL +from lammps import LMP_TYPE_VECTOR, LMP_SIZE_VECTOR, LMP_SIZE_ROWS, LMP_SIZE_COLS +from lammps import LAMMPS_DOUBLE_2D, LAMMPS_AUTODETECT has_manybody=False try: @@ -311,6 +313,14 @@ create_atoms 1 single & self.assertEqual(minval,1.0) self.assertEqual(maxval,2.1) + ndist1 = self.lmp.extract_compute("dist",LMP_STYLE_LOCAL,LMP_SIZE_VECTOR) + ndist2 = self.lmp.extract_compute("dist",LMP_STYLE_LOCAL,LMP_SIZE_ROWS) + ndist3 = self.lmp.extract_compute("dist",LMP_STYLE_LOCAL,LMP_SIZE_COLS) + + self.assertEqual(ndist1,21) + self.assertEqual(ndist2,21) + self.assertEqual(ndist3,0) + self.assertNotEqual(self.lmp.find_pair_neighlist("lj/cut"),-1) self.assertNotEqual(self.lmp.find_compute_neighlist("dist"),-1) self.assertEqual(self.lmp.find_compute_neighlist("xxx"),-1) diff --git a/unittest/python/python-numpy.py b/unittest/python/python-numpy.py index 010255a81e..7cbae8e48d 100644 --- a/unittest/python/python-numpy.py +++ b/unittest/python/python-numpy.py @@ -93,17 +93,39 @@ class PythonNumpy(unittest.TestCase): # TODO pass - def testExtractComputeLocalScalar(self): - # TODO - pass - def testExtractComputeLocalVector(self): - # TODO - pass + self.lmp.command("region box block 0 2 0 2 0 2") + self.lmp.command("create_box 1 box") + self.lmp.command("create_atoms 1 single 1.0 1.0 1.0") + self.lmp.command("create_atoms 1 single 1.0 1.0 1.5") + self.lmp.command("mass 1 1.0") + self.lmp.command("pair_style lj/cut 1.9") + self.lmp.command("pair_coeff 1 1 1.0 1.0") + self.lmp.command("compute r0 all pair/local dist") + self.lmp.command("run 0 post no") + values = self.lmp.numpy.extract_compute("r0", LMP_STYLE_LOCAL, LMP_TYPE_VECTOR) + self.assertEqual(values.ndim, 1) + self.assertEqual(values.size, 2) + self.assertEqual(values[0], 0.5) + self.assertEqual(values[1], 1.5) def testExtractComputeLocalArray(self): - # TODO - pass + self.lmp.command("region box block 0 2 0 2 0 2") + self.lmp.command("create_box 1 box") + self.lmp.command("create_atoms 1 single 1.0 1.0 1.0") + self.lmp.command("create_atoms 1 single 1.0 1.0 1.5") + self.lmp.command("mass 1 1.0") + self.lmp.command("pair_style lj/cut 1.9") + self.lmp.command("pair_coeff 1 1 1.0 1.0") + self.lmp.command("compute r0 all pair/local dist dx dy dz") + self.lmp.command("run 0 post no") + values = self.lmp.numpy.extract_compute("r0", LMP_STYLE_LOCAL, LMP_TYPE_ARRAY) + self.assertEqual(values.ndim, 2) + self.assertEqual(values.size, 8) + self.assertEqual(values[0,0], 0.5) + self.assertEqual(values[0,3], -0.5) + self.assertEqual(values[1,0], 1.5) + self.assertEqual(values[1,3], 1.5) def testExtractAtomDeprecated(self): self.lmp.command("units lj")