Merge remote-tracking branch 'origin/master' into latte

cmake/CMakeLists.txt

@@ -23,13 +23,11 @@ if(NOT CMAKE_BUILD_TYPE AND NOT CMAKE_CXX_FLAGS)
   set(CMAKE_BUILD_TYPE Release CACHE STRING "Choose the type of build, options are: None Debug Release RelWithDebInfo MinSizeRel." FORCE)
 endif(NOT CMAKE_BUILD_TYPE AND NOT CMAKE_CXX_FLAGS)
 
-foreach(STYLE_FILE style_angle.h style_atom.h style_body.h style_bond.h style_command.h style_compute.h style_dihedral.h style_dump.h
+# remove any style headers in the src dir
-             style_fix.h style_improper.h style_integrate.h style_kspace.h style_minimize.h style_nbin.h style_npair.h style_nstencil.h
+file(GLOB SRC_STYLE_FILES ${LAMMPS_SOURCE_DIR}/style_*.h)
-             style_ntopo.h style_pair.h style_reader.h style_region.h)
+if(SRC_STYLE_FILES)
-  if(EXISTS ${LAMMPS_SOURCE_DIR}/${STYLE_FILE})
+  file(REMOVE ${SRC_STYLE_FILES})
-    message(FATAL_ERROR "There is a ${STYLE_FILE} in ${LAMMPS_SOURCE_DIR}, please clean up the source directory first")
+endif()
-  endif()
-endforeach()
 
 enable_language(CXX)
 
@@ -43,10 +41,13 @@ include(CheckCCompilerFlag)
 # User input options                                                   #
 ########################################################################
 option(BUILD_SHARED_LIBS "Build shared libs" OFF)
-option(INSTALL_LIB "Install lammps library and header" ON)
+if(BUILD_SHARED_LIBS) # for all pkg libs, mpi_stubs and linalg
+  set(CMAKE_POSITION_INDEPENDENT_CODE ON)
+endif()
 include(GNUInstallDirs)
 
 set(LAMMPS_LINK_LIBS)
+set(LAMMPS_API_DEFINES)
 option(ENABLE_MPI "Build MPI version" OFF)
 if(ENABLE_MPI)
   find_package(MPI REQUIRED)
@@ -58,13 +59,15 @@ if(ENABLE_MPI)
   endif()
 else()
   file(GLOB MPI_SOURCES ${LAMMPS_SOURCE_DIR}/STUBS/mpi.c)
-  list(APPEND LIB_SOURCES ${MPI_SOURCES})
+  add_library(mpi_stubs STATIC ${MPI_SOURCES})
   include_directories(${LAMMPS_SOURCE_DIR}/STUBS)
+  list(APPEND LAMMPS_LINK_LIBS mpi_stubs)
 endif()
 
 set(LAMMPS_SIZE_LIMIT "LAMMPS_SMALLBIG" CACHE STRING "Lammps size limit")
 set_property(CACHE LAMMPS_SIZE_LIMIT PROPERTY STRINGS LAMMPS_SMALLBIG LAMMPS_BIGBIG LAMMPS_SMALLSMALL)
 add_definitions(-D${LAMMPS_SIZE_LIMIT})
+set(LAMMPS_API_DEFINES "${LAMMPS_API_DEFINES} -D${LAMMPS_SIZE_LIMIT}")
 
 set(LAMMPS_MEMALIGN "64" CACHE STRING "enables the use of the posix_memalign() call instead of malloc() when large chunks or memory are allocated by LAMMPS")
 add_definitions(-DLAMMPS_MEMALIGN=${LAMMPS_MEMALIGN})
@@ -72,6 +75,13 @@ add_definitions(-DLAMMPS_MEMALIGN=${LAMMPS_MEMALIGN})
 option(LAMMPS_EXCEPTIONS "enable the use of C++ exceptions for error messages (useful for library interface)" OFF)
 if(LAMMPS_EXCEPTIONS)
   add_definitions(-DLAMMPS_EXCEPTIONS)
+  set(LAMMPS_API_DEFINES "${LAMMPS_API_DEFINES -DLAMMPS_EXCEPTIONS")
+endif()
+
+set(LAMMPS_MACHINE "" CACHE STRING "Suffix to append to lmp binary and liblammps (WON'T enable any features automatically")
+mark_as_advanced(LAMMPS_MACHINE)
+if(LAMMPS_MACHINE)
+  set(LAMMPS_MACHINE "_${LAMMPS_MACHINE}")
 endif()
 
 option(CMAKE_VERBOSE_MAKEFILE "Verbose makefile" OFF)
@@ -86,10 +96,10 @@ set(DEFAULT_PACKAGES ASPHERE BODY CLASS2 COLLOID COMPRESS CORESHELL DIPOLE GRANU
   KSPACE MANYBODY MC MEAM MISC MOLECULE PERI QEQ
   REAX REPLICA RIGID SHOCK SNAP SRD)
 set(OTHER_PACKAGES KIM PYTHON MSCG MPIIO VORONOI POEMS LATTE
-   USER-ATC USER-AWPMD USER-CGDNA
+  USER-ATC USER-AWPMD USER-CGDNA USER-MESO
   USER-CGSDK USER-COLVARS USER-DIFFRACTION USER-DPD USER-DRUDE USER-EFF
   USER-FEP USER-H5MD USER-LB USER-MANIFOLD USER-MEAMC USER-MGPT USER-MISC
-  USER-MOLFILE USER-NETCDF USER-PHONON USER-QTB USER-REAXC USER-SMD 
+  USER-MOLFILE USER-NETCDF USER-PHONON USER-QTB USER-REAXC USER-SMD
   USER-SMTBQ USER-SPH USER-TALLY USER-VTK USER-QUIP USER-QMMM)
 set(ACCEL_PACKAGES USER-OMP KOKKOS OPT USER-INTEL GPU)
 foreach(PKG ${DEFAULT_PACKAGES})
@@ -122,6 +132,8 @@ endif()
 ######################################################
 if(ENABLE_REAX OR ENABLE_MEAM OR ENABLE_USER-QUIP OR ENABLE_USER-QMMM OR ENABLE_LATTE)
   enable_language(Fortran)
+  include(CheckFortranCompilerFlag)
+  check_Fortran_compiler_flag("-fno-second-underscore" FC_HAS_NO_SECOND_UNDERSCORE)
 endif()
 
 if(ENABLE_KOKKOS OR ENABLE_MSCG)
@@ -156,18 +168,22 @@ endif()
 if(ENABLE_MISC)
   option(LAMMPS_XDR "include XDR compatibility files for doing particle dumps in XTC format" OFF)
   if(LAMMPS_XDR)
-    add_definitions(-DLAMMPS_XDR)
+    add_definitions(-DLAMMPS_XDR) # for liblammps
   endif()
 endif()
 
 if(ENABLE_MSCG OR ENABLE_USER-ATC OR ENABLE_USER-AWPMD OR ENABLE_USER-QUIP)
   find_package(LAPACK)
-  if(LAPACK_FOUND)
+  if(NOT LAPACK_FOUND)
-    list(APPEND LAMMPS_LINK_LIBS ${LAPACK_LIBRARIES})
-  else()
     enable_language(Fortran)
     file(GLOB LAPACK_SOURCES ${LAMMPS_LIB_SOURCE_DIR}/linalg/*.f)
-    list(APPEND LIB_SOURCES ${LAPACK_SOURCES})
+    add_library(linalg STATIC ${LAPACK_SOURCES})
+    include(CheckFortranCompilerFlag)
+    check_Fortran_compiler_flag("-fno-second-underscore" FC_HAS_NO_SECOND_UNDERSCORE)
+    if(FC_HAS_NO_SECOND_UNDERSCORE)
+      target_compile_options(linalg PRIVATE -fno-second-underscore)
+    endif()
+    set(LAPACK_LIBRARIES linalg)
   endif()
 endif()
 
@@ -227,7 +243,10 @@ if(ENABLE_LATTE)
 endif()
 
 if(ENABLE_USER-MOLFILE)
-  list(APPEND LAMMPS_LINK_LIBS ${CMAKE_DL_LIBS})
+  add_library(molfile INTERFACE)
+  target_include_directories(molfile INTERFACE ${LAMMPS_LIB_SOURCE_DIR}/molfile)
+  target_link_libraries(molfile INTERFACE ${CMAKE_DL_LIBS})
+  list(APPEND LAMMPS_LINK_LIBS molfile)
 endif()
 
 if(ENABLE_USER-NETCDF)
@@ -244,7 +263,7 @@ endif()
 
 if(ENABLE_USER-QUIP)
   find_package(QUIP REQUIRED)
-  list(APPEND LAMMPS_LINK_LIBS ${QUIP_LIBRARIES} ${CMAKE_Fortran_IMPLICIT_LINK_LIBRARIES})
+  list(APPEND LAMMPS_LINK_LIBS ${QUIP_LIBRARIES} ${LAPACK_LIBRARIES} ${CMAKE_Fortran_IMPLICIT_LINK_LIBRARIES})
 endif()
 
 if(ENABLE_USER-QMMM)
@@ -253,17 +272,6 @@ if(ENABLE_USER-QMMM)
   list(APPEND LAMMPS_LINK_LIBS ${QE_LIBRARIES} ${CMAKE_Fortran_IMPLICIT_LINK_LIBRARIES})
 endif()
 
-if(ENABLE_USER-AWPMD)
-  include_directories(${LAMMPS_LIB_SOURCE_DIR}/awpmd/systems/interact
-    ${LAMMPS_LIB_SOURCE_DIR}/awpmd/ivutils/include)
-endif()
-
-if(ENABLE_USER-H5MD)
-  find_package(HDF5 REQUIRED)
-  list(APPEND LAMMPS_LINK_LIBS ${HDF5_LIBRARIES})
-  include_directories(${HDF5_INCLUDE_DIRS} ${LAMMPS_LIB_SOURCE_DIR}/h5md/include)
-endif()
-
 if(ENABLE_USER-VTK)
   find_package(VTK REQUIRED NO_MODULE)
   include(${VTK_USE_FILE})
@@ -294,13 +302,11 @@ if(ENABLE_MSCG)
       WORKING_DIRECTORY ${LAMMPS_LIB_BINARY_DIR}/mscg)
   endif()
   file(GLOB MSCG_SOURCES ${LAMMPS_LIB_MSCG_BIN_DIR}/*.cpp)
-  list(APPEND LIB_SOURCES ${MSCG_SOURCES})
+  add_library(mscg STATIC ${MSCG_SOURCES})
-  foreach(MSCG_SOURCE ${MSCG_SOURCES})
+  list(APPEND LAMMPS_LINK_LIBS mscg)
-    set_property(SOURCE ${MSCG_SOURCE} APPEND PROPERTY COMPILE_DEFINITIONS
+  target_compile_options(mscg PRIVATE -DDIMENSION=3 -D_exclude_gromacs=1)
-      DIMENSION=3 _exclude_gromacs=1)
+  target_include_directories(mscg PUBLIC ${LAMMPS_LIB_MSCG_BIN_DIR})
-  endforeach()
+  target_link_libraries(mscg ${GSL_LIBRARIES} ${LAPACK_LIBRARIES})
-  include_directories(${LAMMPS_LIB_MSCG_BIN_DIR} ${GSL_INCLUDE_DIRS})
-  list(APPEND LAMMPS_LINK_LIBS ${GSL_LIBRARIES})
 endif()
 
 ########################################################################
@@ -326,7 +332,7 @@ endforeach(FUNC)
 list(APPEND LAMMPS_LINK_LIBS ${MATH_LIBRARIES})
 
 ######################################
-# Generate Basic Style files 
+# Generate Basic Style files
 ######################################
 include(StyleHeaderUtils)
 RegisterStyles(${LAMMPS_SOURCE_DIR})
@@ -335,13 +341,27 @@ RegisterStyles(${LAMMPS_SOURCE_DIR})
 # add sources of enabled packages
 ############################################
 foreach(PKG ${DEFAULT_PACKAGES} ${OTHER_PACKAGES})
-  if(ENABLE_${PKG})
+  set(${PKG}_SOURCES_DIR ${LAMMPS_SOURCE_DIR}/${PKG})
-    set(${PKG}_SOURCES_DIR ${LAMMPS_SOURCE_DIR}/${PKG})
 
+  # ignore PKG files which were manually installed in src folder
+  # headers are ignored during RegisterStyles
+  file(GLOB ${PKG}_SOURCES ${${PKG}_SOURCES_DIR}/*.cpp)
+  file(GLOB ${PKG}_HEADERS ${${PKG}_SOURCES_DIR}/*.h)
+
+  foreach(PKG_FILE in ${${PKG}_SOURCES})
+      get_filename_component(FNAME ${PKG_FILE} NAME)
+      list(REMOVE_ITEM LIB_SOURCES ${LAMMPS_SOURCE_DIR}/${FNAME})
+  endforeach()
+
+  foreach(PKG_FILE in ${${PKG}_HEADERS})
+      get_filename_component(FNAME ${PKG_FILE} NAME)
+      DetectAndRemovePackageHeader(${LAMMPS_SOURCE_DIR}/${FNAME})
+  endforeach()
+
+  if(ENABLE_${PKG})
     # detects styles in package and adds them to global list
     RegisterStyles(${${PKG}_SOURCES_DIR})
 
-    file(GLOB ${PKG}_SOURCES ${${PKG}_SOURCES_DIR}/*.cpp)
     list(APPEND LIB_SOURCES ${${PKG}_SOURCES})
     include_directories(${${PKG}_SOURCES_DIR})
   endif()
@@ -351,17 +371,53 @@ endforeach()
 # add lib sources of (simple) enabled packages
 ############################################
 foreach(SIMPLE_LIB REAX MEAM POEMS USER-ATC USER-AWPMD USER-COLVARS USER-H5MD
-  USER-MOLFILE USER-QMMM)
+  USER-QMMM)
   if(ENABLE_${SIMPLE_LIB})
-    string(REGEX REPLACE "^USER-" "" SIMPLE_LIB "${SIMPLE_LIB}")
+    string(REGEX REPLACE "^USER-" "" PKG_LIB "${SIMPLE_LIB}")
-    string(TOLOWER "${SIMPLE_LIB}" INC_DIR)
+    string(TOLOWER "${PKG_LIB}" PKG_LIB)
-    file(GLOB_RECURSE ${SIMPLE_LIB}_SOURCES ${LAMMPS_LIB_SOURCE_DIR}/${INC_DIR}/*.F
+    file(GLOB_RECURSE ${PKG_LIB}_SOURCES ${LAMMPS_LIB_SOURCE_DIR}/${PKG_LIB}/*.F
-      ${LAMMPS_LIB_SOURCE_DIR}/${INC_DIR}/*.c ${LAMMPS_LIB_SOURCE_DIR}/${INC_DIR}/*.cpp)
+      ${LAMMPS_LIB_SOURCE_DIR}/${PKG_LIB}/*.c ${LAMMPS_LIB_SOURCE_DIR}/${PKG_LIB}/*.cpp)
-    list(APPEND LIB_SOURCES ${${SIMPLE_LIB}_SOURCES})
+    add_library(${PKG_LIB} STATIC ${${PKG_LIB}_SOURCES})
-    include_directories(${LAMMPS_LIB_SOURCE_DIR}/${INC_DIR})
+    list(APPEND LAMMPS_LINK_LIBS ${PKG_LIB})
+    if(PKG_LIB STREQUAL awpmd)
+      target_include_directories(awpmd PUBLIC ${LAMMPS_LIB_SOURCE_DIR}/awpmd/systems/interact ${LAMMPS_LIB_SOURCE_DIR}/awpmd/ivutils/include)
+    elseif(PKG_LIB STREQUAL h5md)
+      target_include_directories(h5md PUBLIC ${LAMMPS_LIB_SOURCE_DIR}/h5md/include)
+    else()
+      target_include_directories(${PKG_LIB} PUBLIC ${LAMMPS_LIB_SOURCE_DIR}/${PKG_LIB})
+    endif()
   endif()
 endforeach()
 
+if(ENABLE_USER-AWPMD)
+  target_link_libraries(awpmd ${LAPACK_LIBRARIES})
+endif()
+
+if(ENABLE_USER-ATC)
+  target_link_libraries(atc ${LAPACK_LIBRARIES})
+endif()
+
+if(ENABLE_USER-H5MD)
+  find_package(HDF5 REQUIRED)
+  target_link_libraries(h5md ${HDF5_LIBRARIES})
+  target_include_directories(h5md PRIVATE ${HDF5_INCLUDE_DIRS})
+endif()
+
+if(ENABLE_MEAM AND FC_HAS_NO_SECOND_UNDERSCORE)
+  foreach(FSRC ${meam_SOURCES})
+    string(REGEX REPLACE "^.*\\." "" FEXT "${FSRC}")
+    list(FIND CMAKE_Fortran_SOURCE_FILE_EXTENSIONS "${FEXT}" FINDEX)
+    if(FINDEX GREATER -1)
+      set_property(SOURCE ${FSRC} APPEND PROPERTY COMPILE_FLAGS "-fno-second-underscore")
+    endif()
+  endforeach()
+endif()
+
+if(ENABLE_REAX AND FC_HAS_NO_SECOND_UNDERSCORE)
+  target_compile_options(reax PRIVATE -fno-second-underscore)
+endif()
+
+
 ######################################################################
 # packages which selectively include variants based on enabled styles
 # e.g. accelerator packages
@@ -458,62 +514,130 @@ if(ENABLE_USER-INTEL)
 endif()
 
 if(ENABLE_GPU)
-    find_package(CUDA REQUIRED)
-    find_program(BIN2C bin2c)
-    if(NOT BIN2C)
-      message(FATAL_ERROR "Couldn't find bin2c, use -DBIN2C helping cmake to find it.")
-    endif()
-    include_directories(${CUDA_INCLUDE_DIRS})
-    list(APPEND LAMMPS_LINK_LIBS ${CUDA_LIBRARIES} ${CUDA_CUDA_LIBRARY})
-    set(GPU_PREC "SINGLE_DOUBLE" CACHE STRING "Lammps gpu precision size")
-    set_property(CACHE GPU_PREC PROPERTY STRINGS SINGLE_DOUBLE SINGLE_SINGLE DOUBLE_DOUBLE)
-    add_definitions(-D_${GPU_PREC})
-    add_definitions(-DNV_KERNEL -DUCL_CUDADR)
-    option(CUDPP_OPT "Enable CUDPP_OPT" ON)
-
     set(GPU_SOURCES_DIR ${LAMMPS_SOURCE_DIR}/GPU)
-    set(GPU_SOURCES ${GPU_SOURCES_DIR}/gpu_extra.h)
+    set(GPU_SOURCES ${GPU_SOURCES_DIR}/gpu_extra.h
+                    ${GPU_SOURCES_DIR}/fix_gpu.h
+                    ${GPU_SOURCES_DIR}/fix_gpu.cpp)
+
+    set(GPU_API "OpenCL" CACHE STRING "API used by GPU package")
+    set_property(CACHE GPU_API PROPERTY STRINGS OpenCL CUDA)
+
+    set(GPU_PREC "SINGLE_DOUBLE" CACHE STRING "LAMMPS GPU precision size")
+    set_property(CACHE GPU_PREC PROPERTY STRINGS SINGLE_DOUBLE SINGLE_SINGLE DOUBLE_DOUBLE)
+
+    file(GLOB GPU_LIB_SOURCES ${LAMMPS_LIB_SOURCE_DIR}/gpu/*.cpp)
+    file(MAKE_DIRECTORY ${LAMMPS_LIB_BINARY_DIR}/gpu)
+
+    if(GPU_API STREQUAL "CUDA")
+      find_package(CUDA REQUIRED)
+      find_program(BIN2C bin2c)
+      if(NOT BIN2C)
+        message(FATAL_ERROR "Couldn't find bin2c, use -DBIN2C helping cmake to find it.")
+      endif()
+      option(CUDPP_OPT "Enable CUDPP_OPT" ON)
+
+      set(GPU_ARCH "sm_30" CACHE STRING "LAMMPS GPU CUDA SM architecture")
+      set_property(CACHE GPU_ARCH PROPERTY STRINGS sm_10 sm_20 sm_30 sm_60)
+
+      file(GLOB GPU_LIB_CU ${LAMMPS_LIB_SOURCE_DIR}/gpu/*.cu ${CMAKE_CURRENT_SOURCE_DIR}/gpu/*.cu)
+      list(REMOVE_ITEM GPU_LIB_CU ${LAMMPS_LIB_SOURCE_DIR}/gpu/lal_pppm.cu)
+
+      cuda_include_directories(${LAMMPS_LIB_SOURCE_DIR}/gpu ${LAMMPS_LIB_BINARY_DIR}/gpu)
+
+      if(CUDPP_OPT)
+        cuda_include_directories(${LAMMPS_LIB_SOURCE_DIR}/gpu/cudpp_mini)
+        file(GLOB GPU_LIB_CUDPP_SOURCES ${LAMMPS_LIB_SOURCE_DIR}/gpu/cudpp_mini/*.cpp)
+        file(GLOB GPU_LIB_CUDPP_CU ${LAMMPS_LIB_SOURCE_DIR}/gpu/cudpp_mini/*.cu)
+      endif()
+
+      cuda_compile_cubin(GPU_GEN_OBJS ${GPU_LIB_CU} OPTIONS
+                   -DUNIX -O3 -Xptxas -v --use_fast_math -DNV_KERNEL -DUCL_CUDADR -arch=${GPU_ARCH} -D_${GPU_PREC})
+
+      cuda_compile(GPU_OBJS ${GPU_LIB_CUDPP_CU} OPTIONS $<$<BOOL:${BUILD_SHARED_LIBS}>:-Xcompiler=-fPIC>
+                   -DUNIX -O3 -Xptxas -v --use_fast_math -DUCL_CUDADR -arch=${GPU_ARCH} -D_${GPU_PREC})
+
+      foreach(CU_OBJ ${GPU_GEN_OBJS})
+        get_filename_component(CU_NAME ${CU_OBJ} NAME_WE)
+        string(REGEX REPLACE "^.*_lal_" "" CU_NAME "${CU_NAME}")
+        add_custom_command(OUTPUT ${LAMMPS_LIB_BINARY_DIR}/gpu/${CU_NAME}_cubin.h
+          COMMAND ${BIN2C} -c -n ${CU_NAME} ${CU_OBJ} > ${LAMMPS_LIB_BINARY_DIR}/gpu/${CU_NAME}_cubin.h
+          DEPENDS ${CU_OBJ}
+          COMMENT "Generating ${CU_NAME}_cubin.h")
+        list(APPEND GPU_LIB_SOURCES ${LAMMPS_LIB_BINARY_DIR}/gpu/${CU_NAME}_cubin.h)
+      endforeach()
+      set_directory_properties(PROPERTIES ADDITIONAL_MAKE_CLEAN_FILES "${LAMMPS_LIB_BINARY_DIR}/gpu/*_cubin.h")
+
+
+      add_library(gpu STATIC ${GPU_LIB_SOURCES} ${GPU_LIB_CUDPP_SOURCES} ${GPU_OBJS})
+      target_link_libraries(gpu ${CUDA_LIBRARIES} ${CUDA_CUDA_LIBRARY})
+      target_include_directories(gpu PRIVATE ${LAMMPS_LIB_BINARY_DIR}/gpu ${CUDA_INCLUDE_DIRS})
+      target_compile_definitions(gpu PRIVATE -D_${GPU_PREC} -DMPI_GERYON -DUCL_NO_EXIT)
+      if(CUDPP_OPT)
+        target_include_directories(gpu PRIVATE ${LAMMPS_LIB_SOURCE_DIR}/gpu/cudpp_mini)
+        target_compile_definitions(gpu PRIVATE -DUSE_CUDPP)
+      endif()
+
+      list(APPEND LAMMPS_LINK_LIBS gpu)
+
+      add_executable(nvc_get_devices ${LAMMPS_LIB_SOURCE_DIR}/gpu/geryon/ucl_get_devices.cpp)
+      target_compile_definitions(nvc_get_devices PRIVATE -DUCL_CUDADR)
+      target_link_libraries(nvc_get_devices PRIVATE ${CUDA_LIBRARIES} ${CUDA_CUDA_LIBRARY})
+      target_include_directories(nvc_get_devices PRIVATE ${CUDA_INCLUDE_DIRS})
+
+
+    elseif(GPU_API STREQUAL "OpenCL")
+      find_package(OpenCL REQUIRED)
+      set(OCL_TUNE "GENERIC" CACHE STRING "OpenCL Device Tuning")
+      set_property(CACHE OCL_TUNE PROPERTY STRINGS INTEL FERMI KEPLER CYPRESS GENERIC)
+
+      include(OpenCLUtils)
+      set(OCL_COMMON_HEADERS ${LAMMPS_LIB_SOURCE_DIR}/gpu/lal_preprocessor.h ${LAMMPS_LIB_SOURCE_DIR}/gpu/lal_aux_fun1.h)
+
+      file(GLOB GPU_LIB_CU ${LAMMPS_LIB_SOURCE_DIR}/gpu/*.cu)
+      list(REMOVE_ITEM GPU_LIB_CU ${LAMMPS_LIB_SOURCE_DIR}/gpu/lal_gayberne.cu ${LAMMPS_LIB_SOURCE_DIR}/gpu/lal_gayberne_lj.cu)
+
+      foreach(GPU_KERNEL ${GPU_LIB_CU})
+          get_filename_component(basename ${GPU_KERNEL} NAME_WE)
+          string(SUBSTRING ${basename} 4 -1 KERNEL_NAME)
+          GenerateOpenCLHeader(${KERNEL_NAME} ${CMAKE_CURRENT_BINARY_DIR}/gpu/${KERNEL_NAME}_cl.h ${OCL_COMMON_HEADERS} ${GPU_KERNEL})
+          list(APPEND GPU_LIB_SOURCES ${CMAKE_CURRENT_BINARY_DIR}/gpu/${KERNEL_NAME}_cl.h)
+      endforeach()
+
+      GenerateOpenCLHeader(gayberne ${CMAKE_CURRENT_BINARY_DIR}/gpu/gayberne_cl.h ${OCL_COMMON_HEADERS} ${LAMMPS_LIB_SOURCE_DIR}/gpu/lal_ellipsoid_extra.h ${LAMMPS_LIB_SOURCE_DIR}/gpu/lal_gayberne.cu)
+      GenerateOpenCLHeader(gayberne_lj ${CMAKE_CURRENT_BINARY_DIR}/gpu/gayberne_lj_cl.h ${OCL_COMMON_HEADERS} ${LAMMPS_LIB_SOURCE_DIR}/gpu/lal_ellipsoid_extra.h ${LAMMPS_LIB_SOURCE_DIR}/gpu/lal_gayberne_lj.cu)
+      list(APPEND GPU_LIB_SOURCES ${CMAKE_CURRENT_BINARY_DIR}/gpu/gayberne_cl.h ${CMAKE_CURRENT_BINARY_DIR}/gpu/gayberne_lj_cl.h)
+
+      add_library(gpu STATIC ${GPU_LIB_SOURCES})
+      target_link_libraries(gpu ${OpenCL_LIBRARIES})
+      target_include_directories(gpu PRIVATE ${CMAKE_CURRENT_BINARY_DIR}/gpu ${OpenCL_INCLUDE_DIRS})
+      target_compile_definitions(gpu PRIVATE -D_${GPU_PREC} -DMPI_GERYON -DUCL_NO_EXIT)
+      target_compile_definitions(gpu PRIVATE -DUSE_OPENCL)
+
+      list(APPEND LAMMPS_LINK_LIBS gpu)
+
+      add_executable(ocl_get_devices ${LAMMPS_LIB_SOURCE_DIR}/gpu/geryon/ucl_get_devices.cpp)
+      target_compile_definitions(ocl_get_devices PRIVATE -DUCL_OPENCL)
+      target_link_libraries(ocl_get_devices PRIVATE ${OpenCL_LIBRARIES})
+      target_include_directories(ocl_get_devices PRIVATE ${OpenCL_INCLUDE_DIRS})
+    endif()
+
+    # GPU package
+    FindStyleHeaders(${GPU_SOURCES_DIR} FIX_CLASS fix_ FIX)
 
     set_property(GLOBAL PROPERTY "GPU_SOURCES" "${GPU_SOURCES}")
 
     # detects styles which have GPU version
-    RegisterStylesExt(${GPU_SOURCES_DIR} opt GPU_SOURCES)
+    RegisterStylesExt(${GPU_SOURCES_DIR} gpu GPU_SOURCES)
 
     get_property(GPU_SOURCES GLOBAL PROPERTY GPU_SOURCES)
-    
+
-    file(GLOB GPU_LIB_SOURCES ${LAMMPS_LIB_SOURCE_DIR}/gpu/*.cpp)
+    list(APPEND LIB_SOURCES ${GPU_SOURCES})
-    file(GLOB GPU_LIB_CU ${LAMMPS_LIB_SOURCE_DIR}/gpu/*.cu ${CMAKE_CURRENT_SOURCE_DIR}/gpu/*.cu)
+    include_directories(${GPU_SOURCES_DIR})
-    file(GLOB_RECURSE GPU_NOT_LIB_CU ${LAMMPS_LIB_SOURCE_DIR}/gpu/lal_pppm.cu)
-    list(REMOVE_ITEM GPU_LIB_CU ${GPU_NOT_LIB_CU})
-    include_directories(${GPU_SOURCES_DIR} ${LAMMPS_LIB_SOURCE_DIR}/gpu ${LAMMPS_LIB_BINARY_DIR}/gpu)
-    if(CUDPP_OPT)
-      include_directories(${LAMMPS_LIB_SOURCE_DIR}/gpu/cudpp_mini)
-      add_definitions(-DCUDPP_OPT)
-      file(GLOB GPU_LIB_CUDPP_SOURCES ${LAMMPS_LIB_SOURCE_DIR}/gpu/cudpp_mini/*.cpp)
-      file(GLOB GPU_LIB_CUDPP_CU ${LAMMPS_LIB_SOURCE_DIR}/gpu/cudpp_mini/*.cu)
-    endif()
-    cuda_compile(GPU_OBJS ${GPU_LIB_CU} ${GPU_LIB_CUDPP_CU} OPTIONS $<$<BOOL:${BUILD_SHARED_LIBS}>:-Xcompiler=-fPIC>)    
-    file(MAKE_DIRECTORY ${LAMMPS_LIB_BINARY_DIR}/gpu)
-    foreach(CU_OBJ ${GPU_OBJS})
-      get_filename_component(CU_NAME ${CU_OBJ} NAME_WE)
-      string(REGEX REPLACE "^.*_lal_" "" CU_NAME "${CU_NAME}")
-      add_custom_command(OUTPUT ${LAMMPS_LIB_BINARY_DIR}/gpu/${CU_NAME}_cubin.h
-        COMMAND ${BIN2C} -c -n ${CU_NAME} ${CU_OBJ} > ${LAMMPS_LIB_BINARY_DIR}/gpu/${CU_NAME}_cubin.h
-        DEPENDS ${CU_OBJ} 
-        COMMENT "Generating ${CU_NAME}_cubin.h")
-      list(APPEND LIB_SOURCES ${LAMMPS_LIB_BINARY_DIR}/gpu/${CU_NAME}_cubin.h)
-      if(${CU_NAME} STREQUAL "pppm_d") #pppm_d doesn't get linked into the lib
-        set(CU_FORBIDDEN_OBJ "${CU_OBJ}") 
-      endif()
-    endforeach()
-    list(REMOVE_ITEM GPU_OBJS "${CU_FORBIDDEN_OBJ}")
-    list(APPEND LIB_SOURCES ${GPU_SOURCES} ${GPU_LIB_SOURCES} ${GPU_LIB_CUDPP_SOURCES} ${GPU_OBJS})
-    set_directory_properties(PROPERTIES ADDITIONAL_MAKE_CLEAN_FILES "${LAMMPS_LIB_BINARY_DIR}/gpu/*_cubin.h")
 endif()
 
 ######################################################
 # Generate style headers based on global list of
-# styles registered during package selection           
+# styles registered during package selection
 ######################################################
 set(LAMMPS_STYLE_HEADERS_DIR ${CMAKE_CURRENT_BINARY_DIR}/styles)
 
@@ -527,19 +651,21 @@ include_directories(${LAMMPS_STYLE_HEADERS_DIR})
 ############################################
 add_library(lammps ${LIB_SOURCES})
 target_link_libraries(lammps ${LAMMPS_LINK_LIBS})
-set_target_properties(lammps PROPERTIES SOVERSION ${SOVERSION})
+set_target_properties(lammps PROPERTIES OUTPUT_NAME lammps${LAMMPS_MACHINE})
-if(INSTALL_LIB)
+if(BUILD_SHARED_LIBS)
+  set_target_properties(lammps PROPERTIES SOVERSION ${SOVERSION})
   install(TARGETS lammps LIBRARY DESTINATION ${CMAKE_INSTALL_LIBDIR} ARCHIVE DESTINATION ${CMAKE_INSTALL_LIBDIR})
-  install(FILES ${LAMMPS_SOURCE_DIR}/lammps.h DESTINATION ${CMAKE_INSTALL_INCLUDEDIR})
+  install(FILES ${LAMMPS_SOURCE_DIR}/library.h DESTINATION ${CMAKE_INSTALL_INCLUDEDIR}/lammps)
-elseif(BUILD_SHARED_LIBS)
+  configure_file(pkgconfig/liblammps.pc.in ${CMAKE_CURRENT_BINARY_DIR}/liblammps${LAMMPS_MACHINE}.pc @ONLY)
-  message(FATAL_ERROR "Shared library has to be installed, use -DINSTALL_LIB=ON to install lammps with a library")
+  install(FILES ${CMAKE_CURRENT_BINARY_DIR}/liblammps${LAMMPS_MACHINE}.pc DESTINATION ${CMAKE_INSTALL_LIBDIR}/pkgconfig)
 endif()
 
 add_executable(lmp ${LMP_SOURCES})
 target_link_libraries(lmp lammps)
+set_target_properties(lmp PROPERTIES OUTPUT_NAME lmp${LAMMPS_MACHINE})
 install(TARGETS lmp DESTINATION ${CMAKE_INSTALL_BINDIR})
 if(ENABLE_TESTING)
-  add_test(ShowHelp ${CMAKE_CURRENT_BINARY_DIR}/lmp -help)
+  add_test(ShowHelp lmp${LAMMPS_MACHINE} -help)
 endif()
 
 ##################################
@@ -550,3 +676,27 @@ foreach(PKG ${DEFAULT_PACKAGES} ${OTHER_PACKAGES} ${ACCEL_PACKAGES})
     message(STATUS "Building package: ${PKG}")
   endif()
 endforeach()
+
+string(TOUPPER "${CMAKE_BUILD_TYPE}" BTYPE)
+message(STATUS "<<< Build configuration >>>
+   Build type       ${CMAKE_BUILD_TYPE}
+   Install path     ${CMAKE_INSTALL_PREFIX}
+   Compilers and Flags:
+   C++ Compiler     ${CMAKE_CXX_COMPILER}
+       Type         ${CMAKE_CXX_COMPILER_ID}
+   C++ Flags        ${CMAKE_CXX_FLAGS} ${CMAKE_CXX_FLAGS_${BTYPE}}")
+get_property(LANGUAGES GLOBAL PROPERTY ENABLED_LANGUAGES)
+if(LANGUAGES MATCHES ".*Fortran.*")
+  message(STATUS "Fortran Compiler ${CMAKE_Fortran_COMPILER} 
+           Type     ${CMAKE_Fortran_COMPILER_ID}
+   Fortran Flags    ${CMAKE_Fortran_FLAGS} ${CMAKE_Fortran_FLAGS_${BTYPE}}")
+endif()
+message(STATUS "Linker flags:
+   Executable      ${CMAKE_EXE_LINKER_FLAGS}")
+if(BUILD_SHARED_LIBS)
+  message(STATUS "Shared libries  ${CMAKE_SHARED_LINKER_FLAGS}")
+else()
+  message(STATUS "Static libries  ${CMAKE_STATIC_LINKER_FLAGS}")
+endif()
+message(STATUS "Link libraries: ${LAMMPS_LINK_LIBS}")
+

cmake/Modules/OpenCLUtils.cmake

@@ -0,0 +1,18 @@
+function(GenerateOpenCLHeader varname outfile files)
+    message("Creating ${outfile}...")
+    file(WRITE ${outfile} "const char * ${varname} = \n")
+    math(EXPR ARG_END   "${ARGC}-1")
+
+    foreach(IDX RANGE 2 ${ARG_END})
+        list(GET ARGV ${IDX} filename)
+        file(READ ${filename} content)
+        string(REGEX REPLACE "\\s*//[^\n]*\n" "" content "${content}")
+        string(REGEX REPLACE "\\\\" "\\\\\\\\" content "${content}")
+        string(REGEX REPLACE "\"" "\\\\\"" content "${content}")
+        string(REGEX REPLACE "([^\n]+)\n" "\"\\1\\\\n\"\n" content "${content}")
+        string(REGEX REPLACE "\n+" "\n" content "${content}")
+        file(APPEND ${outfile} "${content}")
+    endforeach()
+
+    file(APPEND ${outfile} ";\n")
+endfunction(GenerateOpenCLHeader)

cmake/Modules/StyleHeaderUtils.cmake

@@ -85,6 +85,35 @@ function(RegisterStyles search_path)
     FindStyleHeaders(${search_path} REGION_CLASS    region_    REGION    ) # region    ) # domain
 endfunction(RegisterStyles)
 
+function(RemovePackageHeader headers pkg_header)
+    get_property(hlist GLOBAL PROPERTY ${headers})
+    list(REMOVE_ITEM hlist ${pkg_header})
+    set_property(GLOBAL PROPERTY ${headers} "${hlist}")
+endfunction(RemovePackageHeader)
+
+function(DetectAndRemovePackageHeader fname)
+    RemovePackageHeader(ANGLE     ${fname})
+    RemovePackageHeader(ATOM_VEC  ${fname})
+    RemovePackageHeader(BODY      ${fname})
+    RemovePackageHeader(BOND      ${fname})
+    RemovePackageHeader(COMMAND   ${fname})
+    RemovePackageHeader(COMPUTE   ${fname})
+    RemovePackageHeader(DIHEDRAL  ${fname})
+    RemovePackageHeader(DUMP      ${fname})
+    RemovePackageHeader(FIX       ${fname})
+    RemovePackageHeader(IMPROPER  ${fname})
+    RemovePackageHeader(INTEGRATE ${fname})
+    RemovePackageHeader(KSPACE    ${fname})
+    RemovePackageHeader(MINIMIZE  ${fname})
+    RemovePackageHeader(NBIN      ${fname})
+    RemovePackageHeader(NPAIR     ${fname})
+    RemovePackageHeader(NSTENCIL  ${fname})
+    RemovePackageHeader(NTOPO     ${fname})
+    RemovePackageHeader(PAIR      ${fname})
+    RemovePackageHeader(READER    ${fname})
+    RemovePackageHeader(REGION    ${fname})
+endfunction(DetectAndRemovePackageHeader)
+
 function(RegisterStylesExt search_path extension sources)
     FindStyleHeadersExt(${search_path} ANGLE_CLASS     ${extension}  ANGLE     ${sources})
     FindStyleHeadersExt(${search_path} ATOM_CLASS      ${extension}  ATOM_VEC  ${sources})

cmake/pkgconfig/liblammps.pc.in

@@ -0,0 +1,18 @@
+# pkg-config file for lammps
+# https://people.freedesktop.org/~dbn/pkg-config-guide.html
+# Usage: cc `pkg-config --cflags --libs liblammps` -o myapp myapp.c
+# after you added @CMAKE_INSTALL_FULL_LIBDIR@/pkg-config to PKG_CONFIG_PATH,
+# e.g. export PKG_CONFIG_PATH=@CMAKE_INSTALL_FULL_LIBDIR@/pkgconfig
+
+prefix=@CMAKE_INSTALL_FULL_PREFIX@
+libdir=@CMAKE_INSTALL_FULL_LIBDIR@
+includedir=@CMAKE_INSTALL_FULL_INCLUDEDIR@
+
+Name: liblammps@LAMMPS_MACHINE@
+Description: Large-scale Atomic/Molecular Massively Parallel Simulator Library
+URL: http://lammps.sandia.gov
+Version:
+Requires:
+Libs: -L${libdir} -llammps@LAMMPS_MACHINE@
+Libs.private: -lm
+Cflags: -I${includedir} @LAMMPS_API_DEFINES@

doc/src/Section_commands.txt

@@ -532,7 +532,7 @@ package"_Section_start.html#start_3.
 "dump vtk"_dump_vtk.html,
 "group2ndx"_group2ndx.html,
 "ndx2group"_group2ndx.html,
-"temper/grem"_temper_grem.html
+"temper/grem"_temper_grem.html,
 "temper/npt"_temper_npt.html :tb(c=3,ea=c)
 
 :line

doc/src/Section_packages.txt

@@ -2191,7 +2191,7 @@ src/USER-MESO/README
 "pair_style edpd"_pair_meso.html
 "pair_style mdpd"_pair_meso.html
 "pair_style tdpd"_pair_meso.html
-"fix mvv/dpd"_fix_mvv.html
+"fix mvv/dpd"_fix_mvv_dpd.html
 examples/USER/meso
 http://lammps.sandia.gov/movies.html#mesodpd :ul
 

doc/src/fix_addforce.txt

@@ -139,6 +139,11 @@ forces added by this fix in a consistent manner.  I.e. there is a
 decrease in potential energy when atoms move in the direction of the
 added force.
 
+The "fix_modify"_fix_modify.html {virial} option is supported by this
+fix to add the contribution due to the added forces on atoms to the
+system's virial as part of "thermodynamic output"_thermo_style.html.
+The default is {virial no}
+
 The "fix_modify"_fix_modify.html {respa} option is supported by this
 fix. This allows to set at which level of the "r-RESPA"_run_style.html
 integrator the fix is adding its forces. Default is the outermost

doc/src/fix_cmap.txt

@@ -97,6 +97,11 @@ The "fix_modify"_fix_modify.html {energy} option is supported by this
 fix to add the potential "energy" of the CMAP interactions system's
 potential energy as part of "thermodynamic output"_thermo_style.html.
 
+The "fix_modify"_fix_modify.html {virial} option is supported by this
+fix to add the contribution due to the interaction between atoms to
+the system's virial as part of "thermodynamic output"_thermo_style.html.
+The default is {virial yes}
+
 This fix computes a global scalar which can be accessed by various
 "output commands"_Section_howto.html#howto_15.  The scalar is the
 potential energy discussed above.  The scalar value calculated by this

doc/src/fix_efield.txt

@@ -124,6 +124,11 @@ can include the forces added by this fix in a consistent manner.
 I.e. there is a decrease in potential energy when atoms move in the
 direction of the added force due to the electric field.
 
+The "fix_modify"_fix_modify.html {virial} option is supported by this
+fix to add the contribution due to the added forces on atoms to the
+system's virial as part of "thermodynamic output"_thermo_style.html.
+The default is {virial no}
+
 The "fix_modify"_fix_modify.html {respa} option is supported by this
 fix. This allows to set at which level of the "r-RESPA"_run_style.html
 integrator the fix adding its forces. Default is the outermost level.

doc/src/fix_external.txt

@@ -131,6 +131,11 @@ forces added by this fix in a consistent manner.  I.e. there is a
 decrease in potential energy when atoms move in the direction of the
 added force.
 
+The "fix_modify"_fix_modify.html {virial} option is supported by this
+fix to add the contribution due to the interactions computed by the
+external program to the system's virial as part of "thermodynamic
+output"_thermo_style.html. The default is {virial yes}
+
 This fix computes a global scalar which can be accessed by various
 "output commands"_Section_howto.html#howto_15.  The scalar is the
 potential energy discussed above.  The scalar stored by this fix

doc/src/fix_modify.txt

@@ -14,10 +14,11 @@ fix_modify fix-ID keyword value ... :pre
 
 fix-ID = ID of the fix to modify :ulb,l
 one or more keyword/value pairs may be appended :l
-keyword = {temp} or {press} or {energy} or {respa} or {dynamic/dof} :l
+keyword = {temp} or {press} or {energy} or {virial} or {respa} or {dynamic/dof} :l
   {temp} value = compute ID that calculates a temperature
   {press} value = compute ID that calculates a pressure
   {energy} value = {yes} or {no}
+  {virial} value = {yes} or {no}
   {respa} value = {1} to {max respa level} or {0} (for outermost level)
   {dynamic/dof} value = {yes} or {no}
     yes/no = do or do not recompute the number of degrees of freedom (DOF) contributing to the temperature :pre
@@ -52,11 +53,10 @@ define their own compute by default, as described in their
 documentation.  Thus this option allows the user to override the
 default method for computing P.
 
-For fixes that calculate a contribution to the potential energy of the
+The {energy} keyword can be used with fixes that support it.
-system, the {energy} keyword will include that contribution in
+{energy yes} adds a contribution to the potential energy of the
-thermodynamic output of potential energy.  This is because the {energy
+system. The fix's global and per-atom
-yes} setting must be specified to include the fix's global or per-atom
+energy is included in the calculation performed by the "compute
-energy in the calculation performed by the "compute
 pe"_compute_pe.html or "compute pe/atom"_compute_pe_atom.html
 commands.  See the "thermo_style"_thermo_style.html command for info
 on how potential energy is output.  For fixes that tally a global
@@ -69,6 +69,25 @@ are using it when performing an "energy minimization"_minimize.html
 and if you want the energy and forces it produces to be part of the
 optimization criteria.
 
+The {virial} keyword can be used with fixes that support it.
+{virial yes} adds a contribution to the virial of the
+system. The fix's global and per-atom
+virial is included in the calculation performed by the "compute
+pressure"_compute_pressure.html or
+"compute stress/atom"_compute_stress_atom.html
+commands.  See the "thermo_style"_thermo_style.html command for info
+on how pressure is output. 
+
+NOTE: You must specify the {virial yes} setting for a fix if you
+are doing "box relaxation"_fix_box_relax.html and
+if you want virial contribution of the fix to be part of the
+relaxation criteria, although this seems unlikely.
+
+NOTE: This option is only supported by fixes that explicitly say
+so. For some of these (e.g. the
+"fix shake"_fix_shake.html command) the default setting is
+{virial yes}, for others it is {virial no}.
+
 For fixes that set or modify forces, it may be possible to select at
 which "r-RESPA"_run_style.html level the fix operates via the {respa}
 keyword. The RESPA level at which the fix is active can be selected.
@@ -111,4 +130,4 @@ pressure"_compute_pressure.html, "thermo_style"_thermo_style.html
 [Default:]
 
 The option defaults are temp = ID defined by fix, press = ID defined
-by fix, energy = no, respa = 0.
+by fix, energy = no, virial = different for each fix style, respa = 0.

doc/src/fix_rigid.txt

@@ -703,6 +703,11 @@ NVT, NPT, NPH rigid styles to add the energy change induced by the
 thermostatting to the system's potential energy as part of
 "thermodynamic output"_thermo_style.html.
 
+The "fix_modify"_fix_modify.html {virial} option is supported by this
+fix to add the contribution due to keeping the objects rigid to the
+system's virial as part of "thermodynamic output"_thermo_style.html.
+The default is {virial yes}
+
 The "fix_modify"_fix_modify.html {temp} and {press} options are
 supported by the 4 NPT and NPH rigid styles to change the computes
 used to calculate the instantaneous pressure tensor.  Note that the 2

doc/src/fix_shake.txt

@@ -186,6 +186,11 @@ to 1 and recompiling LAMMPS.
 
 [Restart, fix_modify, output, run start/stop, minimize info:]
 
+The "fix_modify"_fix_modify.html {virial} option is supported by this
+fix to add the contribution due to keeping the constraints to the
+system's virial as part of "thermodynamic output"_thermo_style.html.
+The default is {virial yes}
+
 No information about these fixes is written to "binary restart
 files"_restart.html.  None of the "fix_modify"_fix_modify.html options
 are relevant to these fixes.  No global or per-atom quantities are

doc/src/fix_smd.txt

@@ -101,6 +101,11 @@ See the "read_restart"_read_restart.html command for info on how to
 re-specify a fix in an input script that reads a restart file, so that
 the operation of the fix continues in an uninterrupted fashion.
 
+The "fix_modify"_fix_modify.html {virial} option is supported by this
+fix to add the contribution due to the added forces on atoms to the
+system's virial as part of "thermodynamic output"_thermo_style.html.
+The default is {virial no}
+
 The "fix_modify"_fix_modify.html {respa} option is supported by
 this fix. This allows to set at which level of the "r-RESPA"_run_style.html
 integrator the fix is adding its forces. Default is the outermost level.

doc/src/fix_wall.txt

@@ -252,6 +252,11 @@ fix to add the energy of interaction between atoms and each wall to
 the system's potential energy as part of "thermodynamic
 output"_thermo_style.html.
 
+The "fix_modify"_fix_modify.html {virial} option is supported by this
+fix to add the contribution due to the interaction between
+atoms and each wall to the system's virial as part of "thermodynamic
+output"_thermo_style.html. The default is {virial no}
+
 The "fix_modify"_fix_modify.html {respa} option is supported by this
 fix. This allows to set at which level of the "r-RESPA"_run_style.html
 integrator the fix is adding its forces. Default is the outermost level.

doc/src/fix_wall_region.txt

@@ -15,7 +15,7 @@ fix ID group-ID wall/region region-ID style epsilon sigma cutoff :pre
 ID, group-ID are documented in "fix"_fix.html command
 wall/region = style name of this fix command
 region-ID = region whose boundary will act as wall
-style = {lj93} or {lj126} or {colloid} or {harmonic}
+style = {lj93} or {lj126} or {lj1043} or {colloid} or {harmonic}
 epsilon = strength factor for wall-particle interaction (energy or energy/distance^2 units)
 sigma = size factor for wall-particle interaction (distance units)
 cutoff = distance from wall at which wall-particle interaction is cut off (distance units) :ul
@@ -112,6 +112,10 @@ For style {lj126}, the energy E is given by the 12/6 potential:
 
 :c,image(Eqs/pair_lj.jpg)
 
+For style {wall/lj1043}, the energy E is given by the 10/4/3 potential:
+
+:c,image(Eqs/fix_wall_lj1043.jpg)
+
 For style {colloid}, the energy E is given by an integrated form of
 the "pair_style colloid"_pair_colloid.html potential:
 
@@ -128,49 +132,8 @@ surface no longer interact.  The energy of the wall potential is
 shifted so that the wall-particle interaction energy is 0.0 at the
 cutoff distance.
 
-For the {lj93} and {lj126} styles, {epsilon} and {sigma} are the usual
+For a full description of these wall styles, see fix_style
-Lennard-Jones parameters, which determine the strength and size of the
+"wall"_fix_wall.html
-particle as it interacts with the wall.  Epsilon has energy units.
-Note that this {epsilon} and {sigma} may be different than any
-{epsilon} or {sigma} values defined for a pair style that computes
-particle-particle interactions.
-
-The {lj93} interaction is derived by integrating over a 3d
-half-lattice of Lennard-Jones 12/6 particles.  The {lj126} interaction
-is effectively a harder, more repulsive wall interaction.
-
-For the {colloid} style, {epsilon} is effectively a Hamaker constant
-with energy units for the colloid-wall interaction, {R} is the radius
-of the colloid particle, {D} is the distance from the surface of the
-colloid particle to the wall (r-R), and {sigma} is the size of a
-constituent LJ particle inside the colloid particle.  Note that the
-cutoff distance Rc in this case is the distance from the colloid
-particle center to the wall.
-
-The {colloid} interaction is derived by integrating over constituent
-LJ particles of size {sigma} within the colloid particle and a 3d
-half-lattice of Lennard-Jones 12/6 particles of size {sigma} in the
-wall.
-
-For the {wall/harmonic} style, {epsilon} is effectively the spring
-constant K, and has units (energy/distance^2).  The input parameter
-{sigma} is ignored.  The minimum energy position of the harmonic
-spring is at the {cutoff}.  This is a repulsive-only spring since the
-interaction is truncated at the {cutoff}
-
-NOTE: For all of the styles, you must insure that r is always > 0 for
-all particles in the group, or LAMMPS will generate an error.  This
-means you cannot start your simulation with particles on the region
-surface (r = 0) or with particles on the wrong side of the region
-surface (r < 0).  For the {wall/lj93} and {wall/lj126} styles, the
-energy of the wall/particle interaction (and hence the force on the
-particle) blows up as r -> 0.  The {wall/colloid} style is even more
-restrictive, since the energy blows up as D = r-R -> 0.  This means
-the finite-size particles of radius R must be a distance larger than R
-from the region surface.  The {harmonic} style is a softer potential
-and does not blow up as r -> 0, but you must use a large enough
-{epsilon} that particles always reamin on the correct side of the
-region surface (r > 0).
 
 [Restart, fix_modify, output, run start/stop, minimize info:]
 
@@ -182,6 +145,11 @@ fix to add the energy of interaction between atoms and the wall to the
 system's potential energy as part of "thermodynamic
 output"_thermo_style.html.
 
+The "fix_modify"_fix_modify.html {virial} option is supported by this
+fix to add the contribution due to the interaction between
+atoms and each wall to the system's virial as part of "thermodynamic
+output"_thermo_style.html. The default is {virial no}
+
 The "fix_modify"_fix_modify.html {respa} option is supported by this
 fix. This allows to set at which level of the "r-RESPA"_run_style.html
 integrator the fix is adding its forces. Default is the outermost level.

lib/colvars/colvarbias_restraint.cpp

@@ -308,7 +308,7 @@ int colvarbias_restraint_centers_moving::update_centers(cvm::real lambda)
     colvarvalue const c_new = colvarvalue::interpolate(initial_centers[i],
                                                        target_centers[i],
                                                        lambda);
-    centers_incr[i] = (c_new).dist2_grad(colvar_centers[i]);
+    centers_incr[i] = 0.5 * c_new.dist2_grad(colvar_centers[i]);
     colvar_centers[i] = c_new;
     variables(i)->wrap(colvar_centers[i]);
   }

lib/gpu/geryon/ocl_device.h

@@ -140,8 +140,13 @@ class UCL_Device {
   inline void push_command_queue() {
     cl_int errorv;
     _cq.push_back(cl_command_queue());
-    _cq.back()=clCreateCommandQueue(_context,_cl_device,
+
-                                    CL_QUEUE_PROFILING_ENABLE,&errorv);
+#ifdef CL_VERSION_2_0
+    cl_queue_properties props[] = {CL_QUEUE_PROPERTIES, CL_QUEUE_PROFILING_ENABLE, 0};
+    _cq.back()=clCreateCommandQueueWithProperties(_context, _cl_device, props, &errorv);
+#else
+    _cq.back()=clCreateCommandQueue(_context, _cl_device, CL_QUEUE_PROFILING_ENABLE, &errorv);
+#endif
     if (errorv!=CL_SUCCESS) {
       std::cerr << "Could not create command queue on device: " << name()
                 << std::endl;

lib/gpu/lal_sw.cu

@@ -217,7 +217,7 @@ __kernel void k_sw(const __global numtyp4 *restrict x_,
 
   if (ii<inum) {
     int nbor, nbor_end, i, numj;
-    const int* nbor_mem = dev_packed;
+    const __global int* nbor_mem = dev_packed;
     nbor_info(dev_nbor,dev_packed,nbor_pitch,t_per_atom,ii,offset,i,numj,
               n_stride,nbor_end,nbor);
 
@@ -428,7 +428,7 @@ __kernel void k_sw_three_center(const __global numtyp4 *restrict x_,
 
   if (ii<inum) {
     int i, numj, nbor_j, nbor_end;
-    const int* nbor_mem = dev_packed;
+    const __global int* nbor_mem = dev_packed;
     int offset_j=offset/t_per_atom;
     nbor_info(dev_nbor,dev_packed,nbor_pitch,t_per_atom,ii,offset_j,i,numj,
               n_stride,nbor_end,nbor_j);
@@ -570,7 +570,7 @@ __kernel void k_sw_three_end(const __global numtyp4 *restrict x_,
 
   if (ii<inum) {
     int i, numj, nbor_j, nbor_end, k_end;
-    const int* nbor_mem = dev_packed;
+    const __global int* nbor_mem = dev_packed;
     int offset_j=offset/t_per_atom;
     nbor_info(dev_nbor,dev_packed,nbor_pitch,t_per_atom,ii,offset_j,i,numj,
               n_stride,nbor_end,nbor_j);
@@ -724,7 +724,7 @@ __kernel void k_sw_three_end_vatom(const __global numtyp4 *restrict x_,
 
   if (ii<inum) {
     int i, numj, nbor_j, nbor_end, k_end;
-    const int* nbor_mem = dev_packed;
+    const __global int* nbor_mem = dev_packed;
     int offset_j=offset/t_per_atom;
     nbor_info(dev_nbor,dev_packed,nbor_pitch,t_per_atom,ii,offset_j,i,numj,
               n_stride,nbor_end,nbor_j);

lib/gpu/lal_tersoff.cu

@@ -272,7 +272,7 @@ __kernel void k_tersoff_zeta(const __global numtyp4 *restrict x_,
 
   if (ii<inum) {
     int nbor_j, nbor_end, i, numj;
-    const int* nbor_mem=dev_packed;
+    const __global int* nbor_mem=dev_packed;
     int offset_j=offset/t_per_atom;
     nbor_info(dev_nbor,dev_packed,nbor_pitch,t_per_atom,ii,offset_j,i,numj,
               n_stride,nbor_end,nbor_j);
@@ -433,7 +433,7 @@ __kernel void k_tersoff_repulsive(const __global numtyp4 *restrict x_,
 
   if (ii<inum) {
     int nbor, nbor_end, i, numj;
-    const int* nbor_mem=dev_packed;
+    const __global int* nbor_mem=dev_packed;
     nbor_info(dev_nbor,dev_packed,nbor_pitch,t_per_atom,ii,offset,i,numj,
               n_stride,nbor_end,nbor);
 
@@ -545,7 +545,7 @@ __kernel void k_tersoff_three_center(const __global numtyp4 *restrict x_,
 
   if (ii<inum) {
     int i, numj, nbor_j, nbor_end;
-    const int* nbor_mem=dev_packed;
+    const __global int* nbor_mem=dev_packed;
     int offset_j=offset/t_per_atom;
     nbor_info(dev_nbor,dev_packed,nbor_pitch,t_per_atom,ii,offset_j,i,numj,
               n_stride,nbor_end,nbor_j);
@@ -732,7 +732,7 @@ __kernel void k_tersoff_three_end(const __global numtyp4 *restrict x_,
 
   if (ii<inum) {
     int i, numj, nbor_j, nbor_end, k_end;
-    const int* nbor_mem=dev_packed;
+    const __global int* nbor_mem=dev_packed;
     int offset_j=offset/t_per_atom;
     nbor_info(dev_nbor,dev_packed,nbor_pitch,t_per_atom,ii,offset_j,i,numj,
               n_stride,nbor_end,nbor_j);
@@ -977,7 +977,7 @@ __kernel void k_tersoff_three_end_vatom(const __global numtyp4 *restrict x_,
 
   if (ii<inum) {
     int i, numj, nbor_j, nbor_end, k_end;
-    const int* nbor_mem = dev_packed;
+    const __global int* nbor_mem = dev_packed;
     int offset_j=offset/t_per_atom;
     nbor_info(dev_nbor,dev_packed,nbor_pitch,t_per_atom,ii,offset_j,i,numj,
               n_stride,nbor_end,nbor_j);

lib/gpu/lal_vashishta.cu

@@ -436,7 +436,7 @@ __kernel void k_vashishta_three_center(const __global numtyp4 *restrict x_,
 
   if (ii<inum) {
     int i, numj, nbor_j, nbor_end;
-    const int* nbor_mem = dev_packed;
+    const __global int* nbor_mem = dev_packed;
     int offset_j=offset/t_per_atom;
     nbor_info(dev_nbor,dev_packed,nbor_pitch,t_per_atom,ii,offset_j,i,numj,
               n_stride,nbor_end,nbor_j);
@@ -580,7 +580,7 @@ __kernel void k_vashishta_three_end(const __global numtyp4 *restrict x_,
 
   if (ii<inum) {
     int i, numj, nbor_j, nbor_end, k_end;
-    const int* nbor_mem = dev_packed;
+    const __global int* nbor_mem = dev_packed;
     int offset_j=offset/t_per_atom;
     nbor_info(dev_nbor,dev_packed,nbor_pitch,t_per_atom,ii,offset_j,i,numj,
               n_stride,nbor_end,nbor_j);
@@ -735,7 +735,7 @@ __kernel void k_vashishta_three_end_vatom(const __global numtyp4 *restrict x_,
 
   if (ii<inum) {
     int i, numj, nbor_j, nbor_end, k_end;
-    const int* nbor_mem = dev_packed;
+    const __global int* nbor_mem = dev_packed;
     int offset_j=offset/t_per_atom;
     nbor_info(dev_nbor,dev_packed,nbor_pitch,t_per_atom,ii,offset_j,i,numj,
               n_stride,nbor_end,nbor_j);

src/.gitignore

@@ -46,6 +46,7 @@
 
 /colvarproxy_lammps.cpp
 /colvarproxy_lammps.h
+/colvarproxy_lammps_version.h
 /fix_colvars.cpp
 /fix_colvars.h
 /dump_molfile.cpp

src/COLLOID/fix_wall_colloid.cpp

@@ -81,6 +81,7 @@ void FixWallColloid::wall_particle(int m, int which, double coord)
   double r3,rinv3,r2inv3,r4inv3;
   double rad,rad2,rad4,rad8,diam,new_coeff2;
   double eoffset;
+  double vn;
 
   double **x = atom->x;
   double **f = atom->f;
@@ -151,6 +152,12 @@ void FixWallColloid::wall_particle(int m, int which, double coord)
       ewall[0] -= eoffset;
 
       ewall[m+1] += fwall;
+
+      if (evflag) {
+        if (side < 0) vn = -fwall*delta;
+        else vn = fwall*delta;
+        v_tally(dim, i, vn);
+      }
     }
 
   if (onflag) error->one(FLERR,"Particle on or inside fix wall surface");

src/KOKKOS/fix_wall_lj93_kokkos.cpp

@@ -31,6 +31,7 @@ FixWallLJ93Kokkos<DeviceType>::FixWallLJ93Kokkos(LAMMPS *lmp, int narg, char **a
   execution_space = ExecutionSpaceFromDevice<DeviceType>::space;
   datamask_read = EMPTY_MASK;
   datamask_modify = EMPTY_MASK;
+  virial_flag = 0
 }
 
 /* ----------------------------------------------------------------------

src/MISC/fix_efield.cpp

@@ -55,6 +55,7 @@ FixEfield::FixEfield(LAMMPS *lmp, int narg, char **arg) :
   extscalar = 1;
   respa_level_support = 1;
   ilevel_respa = 0;
+  virial_flag = 1;
 
   qe2f = force->qe2f;
   xstr = ystr = zstr = NULL;
@@ -257,6 +258,11 @@ void FixEfield::post_force(int vflag)
   imageint *image = atom->image;
   int nlocal = atom->nlocal;
 
+  // energy and virial setup
+
+  if (vflag) v_setup(vflag);
+  else evflag = 0;
+
   // reallocate efield array if necessary
 
   if (varflag == ATOM && atom->nmax > maxatom) {
@@ -281,6 +287,7 @@ void FixEfield::post_force(int vflag)
 
   double **x = atom->x;
   double fx,fy,fz;
+  double v[6];
 
   // constant efield
 
@@ -306,6 +313,15 @@ void FixEfield::post_force(int vflag)
           fsum[1] += fx;
           fsum[2] += fy;
           fsum[3] += fz;
+          if (evflag) {
+            v[0] = fx*unwrap[0];
+            v[1] = fy*unwrap[1];
+            v[2] = fz*unwrap[2];
+            v[3] = fx*unwrap[1];
+            v[4] = fx*unwrap[2];
+            v[5] = fy*unwrap[2];
+            v_tally(i, v);
+          }
         }
     }
 

src/MOLECULE/fix_cmap.cpp

@@ -76,6 +76,7 @@ FixCMAP::FixCMAP(LAMMPS *lmp, int narg, char **arg) :
   restart_peratom = 1;
   peatom_flag = 1;
   virial_flag = 1;
+  thermo_virial = 1;
   peratom_freq = 1;
   scalar_flag = 1;
   global_freq = 1;

src/Purge.list

@@ -423,9 +423,6 @@ atom_vec_granular.h
 # deleted on Fri Nov 19 12:53:07 2010 -0500
 fix_pour_omp.cpp
 fix_pour_omp.h
-# deleted on Thu Aug 19 23:20:14 2010 +0000
-fix_qeq.cpp
-fix_qeq.h
 # deleted on Thu Jun 17 01:34:38 2010 +0000
 compute_vsum.cpp
 compute_vsum.h

src/RIGID/fix_rigid.cpp

@@ -74,6 +74,7 @@ FixRigid::FixRigid(LAMMPS *lmp, int narg, char **arg) :
   time_integrate = 1;
   rigid_flag = 1;
   virial_flag = 1;
+  thermo_virial = 1;
   create_attribute = 1;
   dof_flag = 1;
   enforce2d_flag = 1;

src/RIGID/fix_rigid_small.cpp

@@ -78,6 +78,7 @@ FixRigidSmall::FixRigidSmall(LAMMPS *lmp, int narg, char **arg) :
   time_integrate = 1;
   rigid_flag = 1;
   virial_flag = 1;
+  thermo_virial = 1;
   create_attribute = 1;
   dof_flag = 1;
   enforce2d_flag = 1;

src/RIGID/fix_shake.cpp

@@ -60,6 +60,7 @@ FixShake::FixShake(LAMMPS *lmp, int narg, char **arg) :
   MPI_Comm_size(world,&nprocs);
 
   virial_flag = 1;
+  thermo_virial = 1;
   create_attribute = 1;
   dof_flag = 1;
 

src/USER-INTEL/intel_preprocess.h

@@ -255,11 +255,10 @@ enum {TIME_PACK, TIME_HOST_NEIGHBOR, TIME_HOST_PAIR, TIME_OFFLOAD_NEIGHBOR,
     ito = inum;                                                 \
   }
 
-#define IP_PRE_omp_range(ifrom, ip, ito, tid, inum, nthreads)   \
+#define IP_PRE_omp_range(ifrom, ito, tid, inum, nthreads)       \
   {                                                             \
     ifrom = 0;                                                  \
     ito = inum;                                                 \
-    ip = 1;                                                     \
   }
 
 #define IP_PRE_omp_stride_id(ifrom, ip, ito, tid, inum, nthr)   \

src/USER-MISC/fix_smd.cpp

@@ -58,6 +58,7 @@ FixSMD::FixSMD(LAMMPS *lmp, int narg, char **arg) :
   extvector = 1;
   respa_level_support = 1;
   ilevel_respa = 0;
+  virial_flag = 1;
 
   int argoffs=3;
   if (strcmp(arg[argoffs],"cvel") == 0) {
@@ -181,6 +182,11 @@ void FixSMD::setup(int vflag)
 
 void FixSMD::post_force(int vflag)
 {
+  // energy and virial setup
+
+  if (vflag) v_setup(vflag);
+  else evflag = 0;
+
   if (styleflag & SMD_TETHER) smd_tether();
   else smd_couple();
 
@@ -238,12 +244,15 @@ void FixSMD::smd_tether()
   // apply restoring force to atoms in group
   // f = -k*(r-r0)*mass/masstotal
 
+  double **x = atom->x;
   double **f = atom->f;
+  imageint *image = atom->image;
   int *mask = atom->mask;
   int *type = atom->type;
   double *mass = atom->mass;
   double *rmass = atom->rmass;
   double massfrac;
+  double unwrap[3],v[6];
   int nlocal = atom->nlocal;
 
   ftotal[0] = ftotal[1] = ftotal[2] = 0.0;
@@ -259,6 +268,16 @@ void FixSMD::smd_tether()
         ftotal[0] -= fx*massfrac;
         ftotal[1] -= fy*massfrac;
         ftotal[2] -= fz*massfrac;
+        if (evflag) {
+          domain->unmap(x[i],image[i],unwrap);
+          v[0] = fx*massfrac*unwrap[0];
+          v[1] = fy*massfrac*unwrap[1];
+          v[2] = fz*massfrac*unwrap[2];
+          v[3] = fx*massfrac*unwrap[1];
+          v[4] = fx*massfrac*unwrap[2];
+          v[5] = fy*massfrac*unwrap[2];
+          v_tally(i, v);
+        }
       }
   } else {
     for (int i = 0; i < nlocal; i++)
@@ -270,6 +289,16 @@ void FixSMD::smd_tether()
         ftotal[0] -= fx*massfrac;
         ftotal[1] -= fy*massfrac;
         ftotal[2] -= fz*massfrac;
+        if (evflag) {
+          domain->unmap(x[i],image[i],unwrap);
+          v[0] = fx*massfrac*unwrap[0];
+          v[1] = fy*massfrac*unwrap[1];
+          v[2] = fz*massfrac*unwrap[2];
+          v[3] = fx*massfrac*unwrap[1];
+          v[4] = fx*massfrac*unwrap[2];
+          v[5] = fy*massfrac*unwrap[2];
+          v_tally(i, v);
+        }
       }
   }
 }

src/compute_chunk_atom.cpp

@@ -218,10 +218,10 @@ ComputeChunkAtom::ComputeChunkAtom(LAMMPS *lmp, int narg, char **arg) :
       if (limit && !compress) limitfirst = 1;
       iarg += 2;
       if (limit) {
-        if (iarg+1 > narg)
+        if (iarg > narg)
           error->all(FLERR,"Illegal compute chunk/atom command");
-        if (strcmp(arg[iarg+1],"max") == 0) limitstyle = LIMITMAX;
+        if (strcmp(arg[iarg],"max") == 0) limitstyle = LIMITMAX;
-        else if (strcmp(arg[iarg+1],"exact") == 0) limitstyle = LIMITEXACT;
+        else if (strcmp(arg[iarg],"exact") == 0) limitstyle = LIMITEXACT;
         else error->all(FLERR,"Illegal compute chunk/atom command");
         iarg++;
       }

src/fix.cpp

@@ -31,7 +31,7 @@ int Fix::instance_total = 0;
 
 /* ---------------------------------------------------------------------- */
 
-Fix::Fix(LAMMPS *lmp, int narg, char **arg) : 
+Fix::Fix(LAMMPS *lmp, int narg, char **arg) :
   Pointers(lmp),
   id(NULL), style(NULL), extlist(NULL), vector_atom(NULL), array_atom(NULL),
   vector_local(NULL), array_local(NULL), eatom(NULL), vatom(NULL)
@@ -61,6 +61,7 @@ Fix::Fix(LAMMPS *lmp, int narg, char **arg) :
   force_reneighbor = 0;
   box_change_size = box_change_shape = box_change_domain = 0;
   thermo_energy = 0;
+  thermo_virial = 0;
   rigid_flag = 0;
   peatom_flag = 0;
   virial_flag = 0;
@@ -140,8 +141,20 @@ void Fix::modify_params(int narg, char **arg)
     } else if (strcmp(arg[iarg],"energy") == 0) {
       if (iarg+2 > narg) error->all(FLERR,"Illegal fix_modify command");
       if (strcmp(arg[iarg+1],"no") == 0) thermo_energy = 0;
-      else if (strcmp(arg[iarg+1],"yes") == 0) thermo_energy = 1;
+      else if (strcmp(arg[iarg+1],"yes") == 0) {
-      else error->all(FLERR,"Illegal fix_modify command");
+        if (!(THERMO_ENERGY & setmask()))
+          error->all(FLERR,"Illegal fix_modify command");
+        thermo_energy = 1;
+      } else error->all(FLERR,"Illegal fix_modify command");
+      iarg += 2;
+    } else if (strcmp(arg[iarg],"virial") == 0) {
+      if (iarg+2 > narg) error->all(FLERR,"Illegal fix_modify command");
+      if (strcmp(arg[iarg+1],"no") == 0) thermo_virial = 0;
+      else if (strcmp(arg[iarg+1],"yes") == 0) {
+        if (virial_flag == 0)
+          error->all(FLERR,"Illegal fix_modify command");
+        thermo_virial = 1;
+      } else error->all(FLERR,"Illegal fix_modify command");
       iarg += 2;
     } else if (strcmp(arg[iarg],"respa") == 0) {
       if (iarg+2 > narg) error->all(FLERR,"Illegal fix_modify command");
@@ -214,15 +227,22 @@ void Fix::ev_setup(int eflag, int vflag)
 }
 
 /* ----------------------------------------------------------------------
-   setup for virial computation
+   if thermo_virial is on:
-   see integrate::ev_set() for values of vflag (0-6)
+     setup for virial computation
-   fixes call this if use v_tally()
+     see integrate::ev_set() for values of vflag (0-6)
+     fixes call this if use v_tally()
+   else: set evflag=0
 ------------------------------------------------------------------------- */
 
 void Fix::v_setup(int vflag)
 {
   int i,n;
 
+  if (!thermo_virial) {
+    evflag = 0;
+    return;
+  }
+
   evflag = 1;
 
   vflag_global = vflag % 4;
@@ -316,3 +336,58 @@ void Fix::v_tally(int n, int *list, double total, double *v)
     }
   }
 }
+
+/* ----------------------------------------------------------------------
+   tally virial into global and per-atom accumulators
+   i = local index of atom
+   v = total virial for the interaction
+   increment global virial by v
+   increment per-atom virial by v
+   this method can be used when fix computes forces in post_force()
+   and the force depends on a distance to some external object
+     e.g. fix wall/lj93: compute virial only on owned atoms
+------------------------------------------------------------------------- */
+
+void Fix::v_tally(int i, double *v)
+{
+
+  if (vflag_global) {
+    virial[0] += v[0];
+    virial[1] += v[1];
+    virial[2] += v[2];
+    virial[3] += v[3];
+    virial[4] += v[4];
+    virial[5] += v[5];
+  }
+
+  if (vflag_atom) {
+    vatom[i][0] += v[0];
+    vatom[i][1] += v[1];
+    vatom[i][2] += v[2];
+    vatom[i][3] += v[3];
+    vatom[i][4] += v[4];
+    vatom[i][5] += v[5];
+  }
+}
+
+/* ----------------------------------------------------------------------
+   tally virial component into global and per-atom accumulators
+   n = index of virial component (0-5)
+   i = local index of atom
+   vn = nth component of virial for the interaction
+   increment nth component of global virial by vn
+   increment nth component of per-atom virial by vn
+   this method can be used when fix computes forces in post_force()
+   and the force depends on a distance to some external object
+     e.g. fix wall/lj93: compute virial only on owned atoms
+------------------------------------------------------------------------- */
+
+void Fix::v_tally(int n, int i, double vn)
+{
+
+  if (vflag_global)
+    virial[n] += vn;
+
+  if (vflag_atom)
+    vatom[i][n] += vn;
+}

src/fix.h

@@ -36,6 +36,7 @@ class Fix : protected Pointers {
 
   bigint next_reneighbor;        // next timestep to force a reneighboring
   int thermo_energy;             // 1 if fix_modify enabled ThEng, 0 if not
+  int thermo_virial;             // 1 if fix_modify enabled ThVir, 0 if not
   int nevery;                    // how often to call an end_of_step fix
   int rigid_flag;                // 1 if Fix integrates rigid bodies, 0 if not
   int peatom_flag;               // 1 if Fix contributes per-atom eng, 0 if not
@@ -223,6 +224,8 @@ class Fix : protected Pointers {
   void ev_tally(int, int *, double, double, double *);
   void v_setup(int);
   void v_tally(int, int *, double, double *);
+  void v_tally(int, double *);
+  void v_tally(int, int, double);
 
   // union data struct for packing 32-bit and 64-bit ints into double bufs
   // see atom_vec.h for documentation

src/fix_addforce.cpp

@@ -51,6 +51,7 @@ FixAddForce::FixAddForce(LAMMPS *lmp, int narg, char **arg) :
   extvector = 1;
   respa_level_support = 1;
   ilevel_respa = 0;
+  virial_flag = 1;
 
   xstr = ystr = zstr = NULL;
 
@@ -237,10 +238,16 @@ void FixAddForce::post_force(int vflag)
   double **f = atom->f;
   int *mask = atom->mask;
   imageint *image = atom->image;
+  double v[6];
   int nlocal = atom->nlocal;
 
   if (update->ntimestep % nevery) return;
 
+  // energy and virial setup
+
+  if (vflag) v_setup(vflag);
+  else evflag = 0;
+
   if (lmp->kokkos)
     atom->sync_modify(Host, (unsigned int) (F_MASK | MASK_MASK),
                       (unsigned int) F_MASK);
@@ -283,6 +290,15 @@ void FixAddForce::post_force(int vflag)
         f[i][0] += xvalue;
         f[i][1] += yvalue;
         f[i][2] += zvalue;
+        if (evflag) {
+          v[0] = xvalue * unwrap[0];
+          v[1] = yvalue * unwrap[1];
+          v[2] = zvalue * unwrap[2];
+          v[3] = xvalue * unwrap[1];
+          v[4] = xvalue * unwrap[2];
+          v[5] = yvalue * unwrap[2];
+          v_tally(i,v);
+        }
       }
 
   // variable force, wrap with clear/add
@@ -290,6 +306,7 @@ void FixAddForce::post_force(int vflag)
   // wrap with clear/add
 
   } else {
+    double unwrap[3];
 
     modify->clearstep_compute();
 
@@ -307,20 +324,39 @@ void FixAddForce::post_force(int vflag)
 
     modify->addstep_compute(update->ntimestep + 1);
 
-    for (int i = 0; i < nlocal; i++)
+    for (int i = 0; i < nlocal; i++) {
       if (mask[i] & groupbit) {
         if (region && !region->match(x[i][0],x[i][1],x[i][2])) continue;
-        if (estyle == ATOM) foriginal[0] += sforce[i][3];
+        domain->unmap(x[i],image[i],unwrap);
+        if (xstyle == ATOM) xvalue = sforce[i][0];
+        if (ystyle == ATOM) yvalue = sforce[i][1];
+        if (zstyle == ATOM) zvalue = sforce[i][2];
+
+        if (estyle == ATOM) {
+          foriginal[0] += sforce[i][3];
+        } else {
+          if (xstyle) foriginal[0] -= xvalue*unwrap[0];
+          if (ystyle) foriginal[0] -= yvalue*unwrap[1];
+          if (zstyle) foriginal[0] -= zvalue*unwrap[2];
+        }
         foriginal[1] += f[i][0];
         foriginal[2] += f[i][1];
         foriginal[3] += f[i][2];
-        if (xstyle == ATOM) f[i][0] += sforce[i][0];
+
-        else if (xstyle) f[i][0] += xvalue;
+        if (xstyle) f[i][0] += xvalue;
-        if (ystyle == ATOM) f[i][1] += sforce[i][1];
+        if (ystyle) f[i][1] += yvalue;
-        else if (ystyle) f[i][1] += yvalue;
+        if (zstyle) f[i][2] += zvalue;
-        if (zstyle == ATOM) f[i][2] += sforce[i][2];
+        if (evflag) {
-        else if (zstyle) f[i][2] += zvalue;
+          v[0] = xstyle ? xvalue*unwrap[0] : 0.0;
+          v[1] = ystyle ? yvalue*unwrap[1] : 0.0;
+          v[2] = zstyle ? zvalue*unwrap[2] : 0.0;
+          v[3] = xstyle ? xvalue*unwrap[1] : 0.0;
+          v[4] = xstyle ? xvalue*unwrap[2] : 0.0;
+          v[5] = ystyle ? yvalue*unwrap[2] : 0.0;
+          v_tally(i, v);
+        }
       }
+    }
   }
 }
 

src/fix_external.cpp

@@ -37,6 +37,7 @@ FixExternal::FixExternal(LAMMPS *lmp, int narg, char **arg) :
   scalar_flag = 1;
   global_freq = 1;
   virial_flag = 1;
+  thermo_virial = 1;
   extscalar = 1;
 
   if (strcmp(arg[3],"pf/callback") == 0) {

src/fix_wall.cpp

@@ -45,6 +45,7 @@ FixWall::FixWall(LAMMPS *lmp, int narg, char **arg) :
   extvector = 1;
   respa_level_support = 1;
   ilevel_respa = 0;
+  virial_flag = 1;
 
   // parse args
 
@@ -298,7 +299,12 @@ void FixWall::pre_force(int vflag)
 
 void FixWall::post_force(int vflag)
 {
+
+  // energy and virial setup
+
   eflag = 0;
+  if (vflag) v_setup(vflag);
+  else evflag = 0;
   for (int m = 0; m <= nwall; m++) ewall[m] = 0.0;
 
   // coord = current position of wall

src/fix_wall_harmonic.cpp

@@ -34,6 +34,7 @@ FixWallHarmonic::FixWallHarmonic(LAMMPS *lmp, int narg, char **arg) :
 void FixWallHarmonic::wall_particle(int m, int which, double coord)
 {
   double delta,dr,fwall;
+  double vn;
 
   double **x = atom->x;
   double **f = atom->f;
@@ -60,6 +61,12 @@ void FixWallHarmonic::wall_particle(int m, int which, double coord)
       f[i][dim] -= fwall;
       ewall[0] += epsilon[m]*dr*dr;
       ewall[m+1] += fwall;
+
+      if (evflag) {
+        if (side < 0) vn = -fwall*delta;
+        else vn = fwall*delta;
+        v_tally(dim, i, vn);
+      }
     }
 
   if (onflag) error->one(FLERR,"Particle on or inside fix wall surface");

src/fix_wall_lj1043.cpp

@@ -52,6 +52,7 @@ void FixWallLJ1043::precompute(int m)
 void FixWallLJ1043::wall_particle(int m, int which, double coord)
 {
   double delta,rinv,r2inv,r4inv,r10inv,fwall;
+  double vn;
 
   double **x = atom->x;
   double **f = atom->f;
@@ -79,5 +80,11 @@ void FixWallLJ1043::wall_particle(int m, int which, double coord)
       ewall[0] += coeff1[m]*r10inv - coeff2[m]*r4inv -
 	coeff3[m]*pow(delta+coeff4[m],-3.0) - offset[m];
       ewall[m+1] += fwall;
+
+      if (evflag) {
+        if (side < 0) vn = -fwall*delta;
+        else vn = fwall*delta;
+        v_tally(dim, i, vn);
+      }
     }
 }

src/fix_wall_lj126.cpp

@@ -48,6 +48,7 @@ void FixWallLJ126::precompute(int m)
 void FixWallLJ126::wall_particle(int m, int which, double coord)
 {
   double delta,rinv,r2inv,r6inv,fwall;
+  double vn;
 
   double **x = atom->x;
   double **f = atom->f;
@@ -76,6 +77,12 @@ void FixWallLJ126::wall_particle(int m, int which, double coord)
       f[i][dim] -= fwall;
       ewall[0] += r6inv*(coeff3[m]*r6inv - coeff4[m]) - offset[m];
       ewall[m+1] += fwall;
+
+      if (evflag) {
+        if (side < 0) vn = -fwall*delta;
+        else vn = fwall*delta;
+        v_tally(dim, i, vn);
+      }
     }
 
   if (onflag) error->one(FLERR,"Particle on or inside fix wall surface");

src/fix_wall_lj93.cpp

@@ -49,6 +49,7 @@ void FixWallLJ93::precompute(int m)
 void FixWallLJ93::wall_particle(int m, int which, double coord)
 {
   double delta,rinv,r2inv,r4inv,r10inv,fwall;
+  double vn;
 
   double **x = atom->x;
   double **f = atom->f;
@@ -79,6 +80,12 @@ void FixWallLJ93::wall_particle(int m, int which, double coord)
       ewall[0] += coeff3[m]*r4inv*r4inv*rinv -
         coeff4[m]*r2inv*rinv - offset[m];
       ewall[m+1] += fwall;
+
+      if (evflag) {
+        if (side < 0) vn = -fwall*delta;
+        else vn = fwall*delta;
+        v_tally(dim, i, vn);
+      }
     }
 
   if (onflag) error->one(FLERR,"Particle on or inside fix wall surface");

src/fix_wall_region.cpp

@@ -25,11 +25,13 @@
 #include "output.h"
 #include "respa.h"
 #include "error.h"
+#include "math_const.h"
 
 using namespace LAMMPS_NS;
 using namespace FixConst;
+using namespace MathConst;
 
-enum{LJ93,LJ126,COLLOID,HARMONIC};
+enum{LJ93,LJ126,LJ1043,COLLOID,HARMONIC};
 
 /* ---------------------------------------------------------------------- */
 
@@ -47,6 +49,7 @@ FixWallRegion::FixWallRegion(LAMMPS *lmp, int narg, char **arg) :
   extvector = 1;
   respa_level_support = 1;
   ilevel_respa = 0;
+  virial_flag = 1;
 
   // parse args
 
@@ -59,6 +62,7 @@ FixWallRegion::FixWallRegion(LAMMPS *lmp, int narg, char **arg) :
 
   if (strcmp(arg[4],"lj93") == 0) style = LJ93;
   else if (strcmp(arg[4],"lj126") == 0) style = LJ126;
+  else if (strcmp(arg[4],"lj1043") == 0) style = LJ1043;
   else if (strcmp(arg[4],"colloid") == 0) style = COLLOID;
   else if (strcmp(arg[4],"harmonic") == 0) style = HARMONIC;
   else error->all(FLERR,"Illegal fix wall/region command");
@@ -143,6 +147,20 @@ void FixWallRegion::init()
     double r2inv = 1.0/(cutoff*cutoff);
     double r6inv = r2inv*r2inv*r2inv;
     offset = r6inv*(coeff3*r6inv - coeff4);
+  } else if (style == LJ1043) {
+    coeff1 = MY_2PI * 2.0/5.0 * epsilon * pow(sigma,10.0);
+    coeff2 = MY_2PI * epsilon * pow(sigma,4.0);
+    coeff3 = MY_2PI * pow(2.0,1/2.0) / 3 * epsilon * pow(sigma,3.0);
+    coeff4 = 0.61 / pow(2.0,1/2.0) * sigma;
+    coeff5 = coeff1 * 10.0;
+    coeff6 = coeff2 * 4.0;
+    coeff7 = coeff3 * 3.0;
+
+    double rinv = 1.0/cutoff;
+    double r2inv = rinv*rinv;
+    double r4inv = r2inv*r2inv;
+    offset = coeff1*r4inv*r4inv*r2inv - coeff2*r4inv -
+      coeff3*pow(cutoff+coeff4,-3.0);
   } else if (style == COLLOID) {
     coeff1 = -4.0/315.0 * epsilon * pow(sigma,6.0);
     coeff2 = -2.0/3.0 * epsilon;
@@ -186,6 +204,7 @@ void FixWallRegion::post_force(int vflag)
 {
   int i,m,n;
   double rinv,fx,fy,fz,tooclose;
+  double delx, dely, delz, v[6];
 
   double **x = atom->x;
   double **f = atom->f;
@@ -198,13 +217,18 @@ void FixWallRegion::post_force(int vflag)
 
   int onflag = 0;
 
+  // energy and virial setup
+
+  eflag = 0;
+  if (vflag) v_setup(vflag);
+  else evflag = 0;
+
   // region->match() insures particle is in region or on surface, else error
   // if returned contact dist r = 0, is on surface, also an error
   // in COLLOID case, r <= radius is an error
   // initilize ewall after region->prematch(),
   //   so a dynamic region can access last timestep values
 
-  eflag = 0;
   ewall[0] = ewall[1] = ewall[2] = ewall[3] = 0.0;
 
   for (i = 0; i < nlocal; i++)
@@ -226,19 +250,32 @@ void FixWallRegion::post_force(int vflag)
 
         if (style == LJ93) lj93(region->contact[m].r);
         else if (style == LJ126) lj126(region->contact[m].r);
+        else if (style == LJ1043) lj1043(region->contact[m].r);
         else if (style == COLLOID) colloid(region->contact[m].r,radius[i]);
         else harmonic(region->contact[m].r);
 
-        ewall[0] += eng;
+        delx = region->contact[m].delx;
-        fx = fwall * region->contact[m].delx * rinv;
+        dely = region->contact[m].dely;
-        fy = fwall * region->contact[m].dely * rinv;
+        delz = region->contact[m].delz;
-        fz = fwall * region->contact[m].delz * rinv;
+        fx = fwall * delx * rinv;
+        fy = fwall * dely * rinv;
+        fz = fwall * delz * rinv;
         f[i][0] += fx;
         f[i][1] += fy;
         f[i][2] += fz;
         ewall[1] -= fx;
         ewall[2] -= fy;
         ewall[3] -= fz;
+        ewall[0] += eng;
+        if (evflag) {
+          v[0] = fx*delx;
+          v[1] = fy*dely;
+          v[2] = fz*delz;
+          v[3] = fx*dely;
+          v[4] = fx*delz;
+          v[5] = fy*delz;
+          v_tally(i, v);
+        }
       }
     }
 
@@ -319,6 +356,23 @@ void FixWallRegion::lj126(double r)
   eng = r6inv*(coeff3*r6inv - coeff4) - offset;
 }
 
+/* ----------------------------------------------------------------------
+   LJ 10/4/3 interaction for particle with wall
+   compute eng and fwall = magnitude of wall force
+------------------------------------------------------------------------- */
+
+void FixWallRegion::lj1043(double r)
+{
+  double rinv = 1.0/r;
+  double r2inv = rinv*rinv;
+  double r4inv = r2inv*r2inv;
+  double r10inv = r4inv*r4inv*r2inv;
+  fwall = coeff5*r10inv*rinv - coeff6*r4inv*rinv -
+    coeff7*pow(r+coeff4,-4.0);
+  eng = coeff1*r10inv - coeff2*r4inv -
+    coeff3*pow(r+coeff4,-3.0) - offset;
+}
+
 /* ----------------------------------------------------------------------
    colloid interaction for finite-size particle of rad with wall
    compute eng and fwall = magnitude of wall force

src/fix_wall_region.h

@@ -47,10 +47,12 @@ class FixWallRegion : public Fix {
   char *idregion;
 
   double coeff1,coeff2,coeff3,coeff4,offset;
+  double coeff5,coeff6,coeff7;
   double eng,fwall;
 
   void lj93(double);
   void lj126(double);
+  void lj1043(double);
   void colloid(double, double);
   void harmonic(double);
 };