Merge release-5-0 into master
authorRoland Schulz <roland@utk.edu>
Wed, 26 Nov 2014 01:10:03 +0000 (20:10 -0500)
committerMark Abraham <mark.j.abraham@gmail.com>
Wed, 26 Nov 2014 10:05:54 +0000 (11:05 +0100)
Conflicts:
src/gromacs/mdlib/domdec.cpp (std::max)

Change-Id: Iaa77c561f52cf5815c35cb08008cf578d1be2314

1  2 
CMakeLists.txt
src/gromacs/domdec/domdec.cpp
src/gromacs/gmxana/gmx_trjcat.c
src/gromacs/gmxlib/txtdump.c
src/gromacs/gmxpreprocess/calc_verletbuf.c
src/gromacs/gmxpreprocess/readrot.c
src/gromacs/mdlib/nbnxn_internal.h
src/gromacs/mdlib/nbnxn_search.c
src/gromacs/mdlib/update.cpp
src/gromacs/pulling/pull_rotation.c

diff --combined CMakeLists.txt
index df7f82a21f70f8fff67b5f8aa2fa2d2e51f8da49,ce2320ca7a5998e93254350c92f5d9151166a260..59a8c3fa40e9164ff62b3b7842ddb66b0da5034a
@@@ -48,28 -48,75 +48,28 @@@ set(CMAKE_LIBRARY_OUTPUT_DIRECTORY ${CM
  set(CMAKE_ARCHIVE_OUTPUT_DIRECTORY ${CMAKE_BINARY_DIR}/lib)
  set(CMAKE_RUNTIME_OUTPUT_DIRECTORY ${CMAKE_BINARY_DIR}/bin)
  
 -# PROJECT_VERSION should have the following structure:
 -# VERSION-dev[-SUFFIX] where the VERSION should have the for: vMajor.vMinor.vPatch
 -#
 -# The "-dev" suffix is important to keep because it makes possible to distinguish
 -# between a build from official release and a build from git release branch on a
 -# machine with no git.
 -#
 -# NOTE: when releasing the "-dev" suffix needs to be stripped off!
 -# REGRESSIONTEST_VERSION and REGRESSIONTEST_BRANCH should always be
 -# defined.
 -set(PROJECT_VERSION "5.0.3-dev")
 -# If this is a released tarball, "-dev" will not be present in
 -# PROJECT_VERSION, and REGRESSIONTEST_VERSION specifies the version
 -# number of the regressiontest tarball against which the code tarball
 -# can be tested. This will be the version of the last patch release.
 -set(REGRESSIONTEST_VERSION "5.0.3-dev")
 -# The MD5 checksum of the regressiontest tarball. Only used if "-dev"
 -# is not present in the PROJECT_VERSION
 -set(REGRESSIONTEST_MD5SUM "004f2377adf252e2caa9c241390d4b9c")
 -# If this is not a released tarball, "-dev" will be present in
 -# PROJECT_VERSION, and REGRESSIONTEST_BRANCH specifies the name of the
 -# gerrit.gromacs.org branch whose HEAD can test this code, *if* this
 -# code contains all recent fixes from the corresponding code branch.
 -set(REGRESSIONTEST_BRANCH "refs/heads/release-5-0")
 -
 -set(CUSTOM_VERSION_STRING ""
 -    CACHE STRING "Custom version string (if empty, use hard-coded default)")
 -mark_as_advanced(CUSTOM_VERSION_STRING)
 -if (CUSTOM_VERSION_STRING)
 -    set(PROJECT_VERSION ${CUSTOM_VERSION_STRING})
 -endif()
 -set(LIBRARY_SOVERSION 0)
 -set(LIBRARY_VERSION ${LIBRARY_SOVERSION}.0.0)
 -# It is a bit irritating, but this has to be set separately for now!
 -SET(CPACK_PACKAGE_VERSION_MAJOR "5")
 -SET(CPACK_PACKAGE_VERSION_MINOR "0")
 -#SET(CPACK_PACKAGE_VERSION_PATCH "0")
 -
 -# The numerical gromacs version. It is 40600 for 4.6.0.
 -# The #define GMX_VERSION in gromacs/version.h is set to this value.
 -math(EXPR NUM_VERSION
 -    "${CPACK_PACKAGE_VERSION_MAJOR}*10000 + ${CPACK_PACKAGE_VERSION_MINOR}*100")
 -if(CPACK_PACKAGE_VERSION_PATCH)
 -    math(EXPR NUM_VERSION
 -         "${NUM_VERSION} + ${CPACK_PACKAGE_VERSION_PATCH}")
 -endif()
 -
 -# The API version tracks the numerical Gromacs version (for now).
 -# It is potentially different from the Gromacs version in the future, if
 -# the programs/libraries diverge from the presumably more stable API.
 -# The #define GMX_API_VERSION in version.h is set to this value to
 -# provide backward compatibility of software written against the Gromacs API.
 -set(API_VERSION ${NUM_VERSION})
 +# Set up common version variables, as well as general information about
 +# the build tree (whether the build is from a source package or from a git
 +# repository).  Also declares a few functions that will be used for generating
 +# version info files later.
 +include(gmxVersionInfo)
  
  if(CMAKE_INSTALL_PREFIX_INITIALIZED_TO_DEFAULT AND UNIX)
      set(CMAKE_INSTALL_PREFIX "/usr/local/gromacs" CACHE STRING "Installation prefix (installation will need write permissions here)" FORCE)
  endif()
  
  include(gmxBuildTypeReference)
 +include(gmxBuildTypeProfile)
  include(gmxBuildTypeTSAN)
  include(gmxBuildTypeASAN)
  include(gmxBuildTypeReleaseWithAssert)
  
  if(NOT CMAKE_BUILD_TYPE)
 -    set(CMAKE_BUILD_TYPE "Release" CACHE STRING "Choose the type of build, options are: Debug Release RelWithDebInfo MinSizeRel Reference RelWithAssert." FORCE)
 +    set(CMAKE_BUILD_TYPE "Release" CACHE STRING "Choose the type of build, options are: Debug Release RelWithDebInfo MinSizeRel Reference RelWithAssert Profile." FORCE)
      # There's no need to offer a user the choice of ThreadSanitizer
      # Set the possible values of build type for cmake-gui
      set_property(CACHE CMAKE_BUILD_TYPE PROPERTY STRINGS "Debug" "Release"
 -        "MinSizeRel" "RelWithDebInfo" "Reference" "RelWithAssert")
 +        "MinSizeRel" "RelWithDebInfo" "Reference" "RelWithAssert" "Profile")
  endif()
  if(CMAKE_CONFIGURATION_TYPES)
      # Add appropriate GROMACS-specific build types for the Visual
          "List of configuration types"
          FORCE)
  endif()
 -set(build_types_with_explicit_flags RELEASE DEBUG RELWITHDEBINFO RELWITHASSERT MINSIZEREL)
 +set(build_types_with_explicit_flags RELEASE DEBUG RELWITHDEBINFO RELWITHASSERT MINSIZEREL PROFILE)
  
 -enable_language(C)
 -enable_language(CXX)
  set_property(GLOBAL PROPERTY FIND_LIBRARY_USE_LIB64_PATHS ON)
  
 -set(CPACK_PACKAGE_NAME "gromacs")
 -set(CPACK_PACKAGE_VERSION ${PROJECT_VERSION})
 -set(CPACK_SOURCE_PACKAGE_FILE_NAME "${CPACK_PACKAGE_NAME}-${CPACK_PACKAGE_VERSION}")
 -set(CPACK_PACKAGE_VENDOR "gromacs.org")
 -set(CPACK_PACKAGE_DESCRIPTION_SUMMARY "Gromacs - a toolkit for high-performance molecular simulation")
 -set(CPACK_RESOURCE_FILE_WELCOME "${CMAKE_SOURCE_DIR}/admin/InstallWelcome.txt")
 -# Its GPL/LGPL, so they do not have to agree to a license for mere usage, but some installers require this...
 -set(CPACK_RESOURCE_FILE_LICENSE "${CMAKE_SOURCE_DIR}/COPYING")
 -set(CPACK_RESOURCE_FILE_README "${CMAKE_SOURCE_DIR}/admin/InstallInfo.txt")
 -set(CPACK_SOURCE_IGNORE_FILES "\\\\.isreposource$;\\\\.git/;\\\\.gitignore$;\\\\.gitattributes;")
 -set(CPACK_PROJECT_CONFIG_FILE "${CMAKE_SOURCE_DIR}/CPackInit.cmake")
 -# CPack source archives include only the directories we list here.
 -# This variable is a list of pairs of names of source and destination
 -# directories. Most of these are used for content GROMACS generates as
 -# part of the configuration or build.
 -set(CPACK_SOURCE_INSTALLED_DIRECTORIES "${CMAKE_SOURCE_DIR};/;${CMAKE_BINARY_DIR}/src/programs/completion;src/programs/completion;${CMAKE_BINARY_DIR}/docs/man/man1;docs/man/man1;${CMAKE_BINARY_DIR}/docs/man/man7;docs/man/man7;${CMAKE_BINARY_DIR}/docs/old-html/final;docs/old-html/final;${CMAKE_BINARY_DIR}/docs/install-guide/final;/")
 -set(CPACK_PACKAGE_CONTACT "gmx-users@gromacs.org")
 -set(CPACK_GMX_BUILD_HELP "${GMX_BUILD_HELP}") #Works even though GMX_BUILD_HELP is defined later because it is off by default.
 -
 -#must come after all cpack settings!
 -include(CPack)
 -
  # Set a default valgrind suppression file.
  # This unfortunately needs to duplicate information from CTest to work as
  # expected...
@@@ -94,12 -165,14 +94,12 @@@ set(MEMORYCHECK_SUPPRESSIONS_FIL
      "File that contains suppressions for the memory checker")
  include(CTest)
  
 -set(SOURCE_IS_GIT_REPOSITORY OFF)
 -set(SOURCE_IS_SOURCE_DISTRIBUTION OFF)
 -if(EXISTS "${CMAKE_SOURCE_DIR}/.git")
 -    set(SOURCE_IS_GIT_REPOSITORY ON)
 -endif()
 -if(NOT EXISTS "${CMAKE_SOURCE_DIR}/admin/.isreposource")
 -    set(SOURCE_IS_SOURCE_DISTRIBUTION ON)
 -endif()
 +include(gmxCPackUtilities)
 +gmx_cpack_init()
 +
 +# Variables that accumulate stuff influencing the installed headers
 +set(INSTALLED_HEADER_INCLUDE_DIRS "")
 +set(INSTALLED_HEADER_DEFINITIONS "")
  
  ########################################################################
  # Check and warn if cache generated on a different host is being reused
@@@ -111,7 -184,7 +111,7 @@@ if(CMAKE_HOST_UNIX
      if(GMX_BUILD_HOSTNAME AND NOT "${GMX_BUILD_HOSTNAME}" STREQUAL "${TMP_HOSTNAME}")
          message(WARNING "
              The CMake cache, probably generated on a different host (${GMX_BUILD_HOSTNAME}),
-             is being reused! This could lead to inconsitencies; therefore, it is
+             is being reused! This could lead to inconsistencies; therefore, it is
              recommended to regenerate the cache!")
      endif()
      set(GMX_BUILD_HOSTNAME "${TMP_HOSTNAME}" CACHE INTERNAL
@@@ -179,14 -252,18 +179,14 @@@ if(GMX_CPU_ACCELERATION
      # TODO remove all references to GMX_CPU_ACCELERATION in master branch
  endif()
  
 -if(NOT GMX_TARGET_MIC)
 -    include(gmxDetectSimd)
 -    gmx_detect_simd(GMX_SUGGESTED_SIMD)
 -else()
 -    set(GMX_SUGGESTED_SIMD "None")
 -endif()
 +include(gmxDetectSimd)
 +gmx_detect_simd(GMX_SUGGESTED_SIMD)
  
  gmx_option_multichoice(
      GMX_SIMD
      "SIMD instruction set for CPU kernels and compiler optimization"
      "${GMX_SUGGESTED_SIMD}"
 -    None SSE2 SSE4.1 AVX_128_FMA AVX_256 AVX2_256 IBM_QPX Sparc64_HPC_ACE Reference)
 +    None SSE2 SSE4.1 AVX_128_FMA AVX_256 AVX2_256 MIC ARM_NEON ARM_NEON_ASIMD IBM_QPX IBM_VMX IBM_VSX Sparc64_HPC_ACE Reference)
  
  gmx_option_multichoice(
      GMX_FFT_LIBRARY
@@@ -227,6 -304,9 +227,6 @@@ option(GMX_OPENMP "Enable OpenMP-based 
  
  option(GMX_USE_TNG "Use the TNG library for trajectory I/O" ON)
  
 -option(GMX_GIT_VERSION_INFO "Generate git version information" ${SOURCE_IS_GIT_REPOSITORY})
 -mark_as_advanced(GMX_GIT_VERSION_INFO)
 -
  if(UNIX)
      option(GMX_SYMLINK_OLD_BINARY_NAMES "Create symbolic links for pre-5.0 binary names" ON)
  endif()
@@@ -300,18 -380,20 +300,18 @@@ if(GMX_SIMD STREQUAL "AVX_256
  endif()
  
  
 -
 -set(PKG_CFLAGS "")
  if(GMX_DOUBLE)
      add_definitions(-DGMX_DOUBLE)
 -    set(PKG_CFLAGS "${PKG_CFLAGS} -DGMX_DOUBLE")
 +    list(APPEND INSTALLED_HEADER_DEFINITIONS "-DGMX_DOUBLE")
      if(GMX_RELAXED_DOUBLE_PRECISION)
          add_definitions(-DGMX_RELAXED_DOUBLE_PRECISION)
      endif()
  endif()
  if(GMX_SOFTWARE_INVSQRT)
 -  set(PKG_CFLAGS "${PKG_CFLAGS} -DGMX_SOFTWARE_INVSQRT")
 +    list(APPEND INSTALLED_HEADER_DEFINITIONS "-DGMX_SOFTWARE_INVSQRT")
  endif()
  
 -if(WIN32 AND NOT CYGWIN)
 +if(WIN32)
      list(APPEND GMX_EXTRA_LIBRARIES "wsock32")
      add_definitions(-DGMX_HAVE_WINSOCK)
  endif()
@@@ -350,7 -432,6 +350,7 @@@ check_function_exists(syscon
  check_function_exists(rsqrt             HAVE_RSQRT)
  check_function_exists(rsqrtf            HAVE_RSQRTF)
  check_function_exists(sqrtf             HAVE_SQRTF)
 +check_function_exists(nice              HAVE_NICE)
  
  include(CheckLibraryExists)
  check_library_exists(m sqrt "" HAVE_LIBM)
@@@ -362,11 -443,6 +362,11 @@@ test_sched_affinity(HAVE_SCHED_AFFINITY
  include(TestBigEndian)
  test_big_endian(GMX_INTEGER_BIG_ENDIAN)
  
 +set(GMX_USE_NICE 0)
 +if (HAVE_UNISTD_H AND HAVE_NICE)
 +    set(GMX_USE_NICE 1)
 +endif()
 +
  # Management of GROMACS options for specific toolchains should go
  # here. Because the initial settings for some of the main options have
  # already happened, but things like library detection and MPI compiler
@@@ -472,13 -548,32 +472,13 @@@ if(CYGWIN
      set(GMX_CYGWIN 1)
  endif()
  
 -if(WIN32 AND NOT CYGWIN)
 +if(WIN32)
      set(GMX_NATIVE_WINDOWS 1)
      # This makes windows.h not declare min/max as macros that would break
      # C++ code using std::min/std::max.
      add_definitions(-DNOMINMAX)
  endif()
  
 -# only bother with finding git and using version.h if the source is a git repo
 -if(GMX_GIT_VERSION_INFO)
 -    if (NOT SOURCE_IS_GIT_REPOSITORY)
 -        message(FATAL_ERROR
 -            "Cannot generate git version information from source tree not under git. "
 -            "Set GMX_GIT_VERSION_INFO=OFF to proceed.")
 -    endif()
 -    # We need at least git v1.5.3 be able to parse git's date output. If not
 -    # found or the version is too small, we can't generate version information.
 -    find_package(Git)
 -
 -    if(NOT GIT_FOUND OR GIT_VERSION_STRING VERSION_LESS "1.5.3")
 -        message(FATAL_ERROR
 -            "No compatible git version found (>= 1.5.3 required). "
 -            "Won't be able to generate development version information. "
 -            "Set GMX_GIT_VERSION_INFO=OFF to proceed.")
 -    endif()
 -endif()
 -
  # Detect boost unless GMX_EXTERNAL_BOOST is explicitly OFF
  # Used for default if GMX_EXTERNAL_BOOST is not defined (first CMake pass)
  if(NOT DEFINED GMX_EXTERNAL_BOOST OR GMX_EXTERNAL_BOOST)
@@@ -521,6 -616,12 +521,6 @@@ if(GMX_USE_TNG AND NOT GMX_EXTERNAL_TNG
      gmx_test_zlib(HAVE_ZLIB)
  endif()
  
 -########################################################################
 -# Generate development version info for cache
 -########################################################################
 -# set(GEN_VERSION_INFO_INTERNAL "ON")
 -# include(gmxGenerateVersionString)
 -
  ########################################################################
  # Our own GROMACS tests
  ########################################################################
@@@ -530,6 -631,10 +530,6 @@@ include_directories(BEFORE ${CMAKE_SOUR
  include_directories(BEFORE ${CMAKE_SOURCE_DIR}/src/external/thread_mpi/include)
  # Required for config.h, maybe should only be set in src/CMakeLists.txt
  include_directories(BEFORE ${CMAKE_BINARY_DIR}/src)
 -# Required for gmx_header_config_gen.h to be found before installation
 -include_directories(BEFORE ${CMAKE_BINARY_DIR}/src/gromacs/utility)
 -# Required for now to make old code compile
 -include_directories(BEFORE ${CMAKE_SOURCE_DIR}/src/gromacs/legacyheaders)
  
  include(gmxTestInlineASM)
  gmx_test_inline_asm_gcc_x86(GMX_X86_GCC_INLINE_ASM)
@@@ -634,14 -739,15 +634,14 @@@ include(gmxManageLinearAlgebraLibraries
  set(GMX_USE_PLUGINS OFF)
  
  if(GMX_LOAD_PLUGINS)
 -  if(CYGWIN OR NOT WIN32)
 +  if(NOT WIN32)
      # Native Windows does not have, nor need dlopen
 -    # Note that WIN32 is set with Cygwin, but Cygwin needs dlopen to use plug-ins
      include(gmxTestdlopen)
      gmx_test_dlopen(HAVE_DLOPEN)
    endif()
  
    # so, should we use plug-ins?
 -  if((WIN32 AND NOT CYGWIN) OR (HAVE_DLOPEN AND BUILD_SHARED_LIBS))
 +  if(WIN32 OR (HAVE_DLOPEN AND BUILD_SHARED_LIBS))
      if(NOT VMD_QUIETLY)
        MESSAGE(STATUS "Using dynamic plugins (e.g VMD-supported file formats)")
      endif()
@@@ -672,8 -778,7 +672,8 @@@ option(GMX_NACL "Configure for Native C
  if (GMX_NACL)
    list(APPEND GMX_EXTRA_LIBRARIES nosys)
    set(GMX_EXE_LINKER_FLAGS "${CMAKE_EXE_LINKER_FLAGS} -lnosys")
 -  set(GMX_NO_NICE 1)
 +  # TODO: Is this still necessary with the check for its presence?
 +  set(GMX_USE_NICE 0)
    set(GMX_NO_RENAME 1)
  endif()
  mark_as_advanced(GMX_NACL)
@@@ -684,18 -789,14 +684,18 @@@ if(GMX_FAHCORE
    include_directories(${COREWRAP_INCLUDE_DIR})
  endif()
  
 -option(GMX_BUILD_HELP "Build man pages, HTML help, and completions automatically (requires that compiled binaries can be executed on the build host)" OFF)
 +# Value of AUTO tries to generate things, but will only produce warnings if
 +# that fails.
 +set(build_help_default AUTO)
 +if (SOURCE_IS_SOURCE_DISTRIBUTION)
 +    set(build_help_default OFF)
 +endif()
 +gmx_option_trivalue(GMX_BUILD_HELP "Build man pages, HTML help, and completions automatically (requires that compiled binaries can be executed on the build host)" ${build_help_default})
  mark_as_advanced(GMX_BUILD_HELP)
 -if (GMX_BUILD_HELP AND SOURCE_IS_SOURCE_DISTRIBUTION AND
 -    "${CMAKE_CURRENT_SOURCE_DIR}" STREQUAL "${CMAKE_CURRENT_BINARY_DIR}")
 -
 +if (GMX_BUILD_HELP AND SOURCE_IS_SOURCE_DISTRIBUTION AND BUILD_IS_INSOURCE)
      message(FATAL_ERROR
 -        "Rebuilding HTML and man pages is not supported for in-source "
 -        "builds from a source distribution. "
 +        "Rebuilding HTML and man pages or shell completions is not supported "
 +        "for in-source builds from a source distribution. "
          "Set GMX_BUILD_HELP=OFF or do an out-of-source build to proceed.")
  endif()
  
@@@ -727,6 -828,7 +727,6 @@@ else(
      set(GMX_EXE_LINKER_FLAGS ${GMX_EXE_LINKER_FLAGS} ${OpenMP_LINKER_FLAGS})
      set(GMX_SHARED_LINKER_FLAGS ${GMX_SHARED_LINKER_FLAGS} ${OpenMP_SHARED_LINKER_FLAGS})
  endif()
 -set(PKG_CFLAGS "${PKG_CFLAGS} ${OpenMP_C_FLAGS}")
  
  ########################################################################
  # Specify install locations
@@@ -743,24 -845,17 +743,24 @@@ mark_as_advanced(GMX_LIB_INSTALL_DIR GM
  
  # These variables are used internally to provide a central location for
  # customizing the install locations.
 -set(LIB_INSTALL_DIR  ${GMX_LIB_INSTALL_DIR})
 -set(BIN_INSTALL_DIR  bin)
 -set(DATA_INSTALL_DIR share/${GMX_DATA_INSTALL_DIR})
 -set(MAN_INSTALL_DIR  share/man)
 -set(INCL_INSTALL_DIR include)
 +set(LIB_INSTALL_DIR       ${GMX_LIB_INSTALL_DIR})
 +set(BIN_INSTALL_DIR       bin)
 +set(DATA_INSTALL_DIR      share/${GMX_DATA_INSTALL_DIR})
 +set(MAN_INSTALL_DIR       share/man)
 +# If the nesting level wrt. the installation root is changed,
 +# gromacs-config.cmake.cmakein needs to be adapted.
 +set(CMAKE_INSTALL_DIR     share/cmake)
 +# TODO: Make GMXRC adapt if this is changed
 +set(PKGCONFIG_INSTALL_DIR ${LIB_INSTALL_DIR}/pkgconfig)
 +set(INCL_INSTALL_DIR      include)
  
  # These variables get written into config.h for use in finding the data
  # directories.
  set(GMXLIB_SEARCH_DIR share/${GMX_DATA_INSTALL_DIR}/top)
  set(GMXLIB_FALLBACK   ${CMAKE_INSTALL_PREFIX}/${DATA_INSTALL_DIR}/top)
  
 +list(APPEND INSTALLED_HEADER_INCLUDE_DIRS ${INCL_INSTALL_DIR})
 +
  # Binary and library suffix options
  include(gmxManageSuffixes)
  
@@@ -806,14 -901,15 +806,14 @@@ install(FILES COPYING DESTINATION ${DAT
  
  if(GMX_EXTERNAL_BOOST)
      include_directories(${Boost_INCLUDE_DIRS})
 -    set(PKG_CFLAGS "${PKG_CFLAGS} -I${Boost_INCLUDE_DIRS}")
 +    list(APPEND INSTALLED_HEADER_INCLUDE_DIRS ${Boost_INCLUDE_DIRS})
  else()
      include_directories(BEFORE ${CMAKE_SOURCE_DIR}/src/external/boost)
 +    list(APPEND INSTALLED_HEADER_INCLUDE_DIRS ${INCL_INSTALL_DIR}/gromacs/external/boost)
 +    list(APPEND INSTALLED_HEADED_DEFINITIONS "-DBOOST_NO_TYPEID")
      # typeid not supported for minimal internal version
      # (would add significant amount of code)
      add_definitions(-DBOOST_NO_TYPEID)
 -    # TODO: Propagate the above settings to the installed CMakeFiles.txt template
 -    # (from share/template/)
 -    set(PKG_CFLAGS "${PKG_CFLAGS} -DBOOST_NO_TYPEID -I${CMAKE_INSTALL_PREFIX}/${INCL_INSTALL_DIR}/gromacs/external/boost")
      if (NOT GMX_BUILD_MDRUN_ONLY)
          install(DIRECTORY ${CMAKE_SOURCE_DIR}/src/external/boost/boost
                  DESTINATION ${INCL_INSTALL_DIR}/gromacs/external/boost
@@@ -856,8 -952,6 +856,8 @@@ if (BUILD_TESTING
      add_subdirectory(tests)
  endif()
  
 +gmx_cpack_write_config()
 +
  # Issue a warning if NVIDIA GPUs were detected, but CUDA was not found.
  # Don't bother the user after the first configure pass.
  if ((CUDA_NOTFOUND_AUTO AND GMX_DETECT_GPU_AVAILABLE) AND NOT GMX_GPU_DETECTION_DONE)
index 6074bcc9ab5ed6b640e121d08bb7d3cca9b5b9e6,307276d49dc73c35e2c6c5320c3311d54330bb52..f70b5c902d483c1398212b4ff11c54ce97efa2bd
   * the research papers on the package. Check out http://www.gromacs.org.
   */
  
 -#ifdef HAVE_CONFIG_H
 -#include <config.h>
 -#endif
 +#include "gmxpre.h"
  
 -#include <stdio.h>
 -#include <time.h>
 +#include "domdec.h"
 +
 +#include "config.h"
 +
 +#include <assert.h>
 +#include <limits.h>
  #include <math.h>
 -#include <string.h>
 +#include <stdio.h>
  #include <stdlib.h>
 -#include <assert.h>
 +#include <string.h>
  
 -#include "typedefs.h"
 -#include "gromacs/utility/smalloc.h"
 -#include "gmx_fatal.h"
 -#include "gmx_fatal_collective.h"
 -#include "vec.h"
 -#include "domdec.h"
 -#include "domdec_network.h"
 -#include "nrnb.h"
 -#include "pbc.h"
 -#include "chargegroup.h"
 -#include "constr.h"
 -#include "mdatoms.h"
 -#include "names.h"
 -#include "force.h"
 -#include "pme.h"
 -#include "mdrun.h"
 -#include "nsgrid.h"
 -#include "shellfc.h"
 -#include "mtop_util.h"
 -#include "gmx_ga2la.h"
 -#include "macros.h"
 -#include "nbnxn_search.h"
 -#include "bondf.h"
 -#include "gmx_omp_nthreads.h"
 -#include "gpu_utils.h"
 -
 -#include "gromacs/fileio/futil.h"
 +#include <algorithm>
 +
 +#include "gromacs/domdec/domdec_network.h"
 +#include "gromacs/ewald/pme.h"
  #include "gromacs/fileio/gmxfio.h"
  #include "gromacs/fileio/pdbio.h"
 +#include "gromacs/imd/imd.h"
 +#include "gromacs/legacyheaders/bonded-threading.h"
 +#include "gromacs/legacyheaders/chargegroup.h"
 +#include "gromacs/legacyheaders/constr.h"
 +#include "gromacs/legacyheaders/force.h"
 +#include "gromacs/legacyheaders/genborn.h"
 +#include "gromacs/legacyheaders/gmx_ga2la.h"
 +#include "gromacs/legacyheaders/gmx_omp_nthreads.h"
 +#include "gromacs/legacyheaders/gpu_utils.h"
 +#include "gromacs/legacyheaders/mdatoms.h"
 +#include "gromacs/legacyheaders/mdrun.h"
 +#include "gromacs/legacyheaders/names.h"
 +#include "gromacs/legacyheaders/network.h"
 +#include "gromacs/legacyheaders/nrnb.h"
 +#include "gromacs/legacyheaders/nsgrid.h"
 +#include "gromacs/legacyheaders/shellfc.h"
 +#include "gromacs/legacyheaders/typedefs.h"
 +#include "gromacs/legacyheaders/vsite.h"
 +#include "gromacs/legacyheaders/types/commrec.h"
 +#include "gromacs/legacyheaders/types/constr.h"
 +#include "gromacs/legacyheaders/types/enums.h"
 +#include "gromacs/legacyheaders/types/forcerec.h"
 +#include "gromacs/legacyheaders/types/hw_info.h"
 +#include "gromacs/legacyheaders/types/ifunc.h"
 +#include "gromacs/legacyheaders/types/inputrec.h"
 +#include "gromacs/legacyheaders/types/mdatom.h"
 +#include "gromacs/legacyheaders/types/nrnb.h"
 +#include "gromacs/legacyheaders/types/ns.h"
 +#include "gromacs/legacyheaders/types/nsgrid.h"
 +#include "gromacs/legacyheaders/types/shellfc.h"
 +#include "gromacs/legacyheaders/types/simple.h"
 +#include "gromacs/legacyheaders/types/state.h"
 +#include "gromacs/math/vec.h"
 +#include "gromacs/math/vectypes.h"
 +#include "gromacs/mdlib/nb_verlet.h"
 +#include "gromacs/mdlib/nbnxn_search.h"
 +#include "gromacs/pbcutil/ishift.h"
 +#include "gromacs/pbcutil/pbc.h"
 +#include "gromacs/pulling/pull.h"
 +#include "gromacs/pulling/pull_rotation.h"
 +#include "gromacs/swap/swapcoords.h"
  #include "gromacs/timing/wallcycle.h"
 +#include "gromacs/topology/block.h"
 +#include "gromacs/topology/idef.h"
 +#include "gromacs/topology/mtop_util.h"
 +#include "gromacs/topology/topology.h"
 +#include "gromacs/utility/basedefinitions.h"
 +#include "gromacs/utility/basenetwork.h"
 +#include "gromacs/utility/cstringutil.h"
 +#include "gromacs/utility/fatalerror.h"
  #include "gromacs/utility/gmxmpi.h"
 -#include "gromacs/swap/swapcoords.h"
  #include "gromacs/utility/qsort_threadsafe.h"
 -#include "gromacs/pulling/pull.h"
 -#include "gromacs/pulling/pull_rotation.h"
 -#include "gromacs/imd/imd.h"
 +#include "gromacs/utility/real.h"
 +#include "gromacs/utility/smalloc.h"
  
  #define DDRANK(dd, rank)    (rank)
  #define DDMASTERRANK(dd)   (dd->masterrank)
@@@ -794,7 -768,7 +794,7 @@@ void dd_move_f(gmx_domdec_t *dd, rvec f
      rvec                  *buf, *sbuf;
      ivec                   vis;
      int                    is;
 -    gmx_bool               bPBC, bScrew;
 +    gmx_bool               bShiftForcesNeedPbc, bScrew;
  
      comm = dd->comm;
  
  
      buf = comm->vbuf.v;
  
 -    n       = 0;
      nzone   = comm->zones.n/2;
      nat_tot = dd->nat_tot;
      for (d = dd->ndim-1; d >= 0; d--)
      {
 -        bPBC   = (dd->ci[dd->dim[d]] == 0);
 -        bScrew = (bPBC && dd->bScrewPBC && dd->dim[d] == XX);
 +        /* Only forces in domains near the PBC boundaries need to
 +           consider PBC in the treatment of fshift */
 +        bShiftForcesNeedPbc   = (dd->ci[dd->dim[d]] == 0);
 +        bScrew                = (bShiftForcesNeedPbc && dd->bScrewPBC && dd->dim[d] == XX);
          if (fshift == NULL && !bScrew)
          {
 -            bPBC = FALSE;
 +            bShiftForcesNeedPbc = FALSE;
          }
          /* Determine which shift vector we need */
          clear_ivec(vis);
              index = ind->index;
              /* Add the received forces */
              n = 0;
 -            if (!bPBC)
 +            if (!bShiftForcesNeedPbc)
              {
                  for (i = 0; i < ind->nsend[nzone]; i++)
                  {
              }
              else if (!bScrew)
              {
 +                /* fshift should always be defined if this function is
 +                 * called when bShiftForcesNeedPbc is true */
 +                assert(NULL != fshift);
                  for (i = 0; i < ind->nsend[nzone]; i++)
                  {
                      at0 = cgindex[index[i]];
@@@ -986,6 -956,7 +986,6 @@@ void dd_atom_sum_real(gmx_domdec_t *dd
  
      buf = &comm->vbuf.v[0][0];
  
 -    n       = 0;
      nzone   = comm->zones.n/2;
      nat_tot = dd->nat_tot;
      for (d = dd->ndim-1; d >= 0; d--)
@@@ -1091,7 -1062,7 +1091,7 @@@ static void dd_sendrecv_ddzone(const gm
  static void dd_move_cellx(gmx_domdec_t *dd, gmx_ddbox_t *ddbox,
                            rvec cell_ns_x0, rvec cell_ns_x1)
  {
 -    int                d, d1, dim, dim1, pos, buf_size, i, j, k, p, npulse, npulse_min;
 +    int                d, d1, dim, pos, buf_size, i, j, p, npulse, npulse_min;
      gmx_ddzone_t      *zp;
      gmx_ddzone_t       buf_s[DDZONECOMM_MAXZONE];
      gmx_ddzone_t       buf_r[DDZONECOMM_MAXZONE];
          if (bPBC)
          {
              /* Take the minimum to avoid double communication */
 -            npulse_min = min(npulse, dd->nc[dim]-1-npulse);
 +            npulse_min = std::min(npulse, dd->nc[dim]-1-npulse);
          }
          else
          {
              {
                  for (d1 = d; d1 < dd->ndim-1; d1++)
                  {
 -                    extr_s[d1][0] = max(extr_s[d1][0], extr_r[d1][0]);
 -                    extr_s[d1][1] = min(extr_s[d1][1], extr_r[d1][1]);
 -                    extr_s[d1][2] = min(extr_s[d1][2], extr_r[d1][2]);
 +                    extr_s[d1][0] = std::max(extr_s[d1][0], extr_r[d1][0]);
 +                    extr_s[d1][1] = std::min(extr_s[d1][1], extr_r[d1][1]);
 +                    extr_s[d1][2] = std::min(extr_s[d1][2], extr_r[d1][2]);
                  }
              }
          }
                  {
                      if (bUse)
                      {
 -                        buf_e[i].min0 = min(buf_e[i].min0, buf_r[i].min0);
 -                        buf_e[i].max1 = max(buf_e[i].max1, buf_r[i].max1);
 -                        buf_e[i].min1 = min(buf_e[i].min1, buf_r[i].min1);
 +                        buf_e[i].min0 = std::min(buf_e[i].min0, buf_r[i].min0);
 +                        buf_e[i].max1 = std::max(buf_e[i].max1, buf_r[i].max1);
 +                        buf_e[i].min1 = std::min(buf_e[i].min1, buf_r[i].min1);
                      }
  
                      if (dd->ndim == 3 && d == 0 && i == buf_size - 1)
                      }
                      if (bUse && dh[d1] >= 0)
                      {
 -                        buf_e[i].mch0 = max(buf_e[i].mch0, buf_r[i].mch0-dh[d1]);
 -                        buf_e[i].mch1 = max(buf_e[i].mch1, buf_r[i].mch1-dh[d1]);
 +                        buf_e[i].mch0 = std::max(buf_e[i].mch0, buf_r[i].mch0-dh[d1]);
 +                        buf_e[i].mch1 = std::max(buf_e[i].mch1, buf_r[i].mch1-dh[d1]);
                      }
                  }
                  /* Copy the received buffer to the send buffer,
  
                  for (d1 = d; d1 < dd->ndim-1; d1++)
                  {
 -                    extr_s[d1][1] = min(extr_s[d1][1], buf_e[pos].min0);
 -                    extr_s[d1][0] = max(extr_s[d1][0], buf_e[pos].max1);
 -                    extr_s[d1][2] = min(extr_s[d1][2], buf_e[pos].min1);
 +                    extr_s[d1][1] = std::min(extr_s[d1][1], buf_e[pos].min0);
 +                    extr_s[d1][0] = std::max(extr_s[d1][0], buf_e[pos].max1);
 +                    extr_s[d1][2] = std::min(extr_s[d1][2], buf_e[pos].min1);
                      pos++;
                  }
  
              {
                  print_ddzone(debug, 1, i, 0, &comm->zone_d1[i]);
              }
 -            cell_ns_x0[dim] = min(cell_ns_x0[dim], comm->zone_d1[i].min0);
 -            cell_ns_x1[dim] = max(cell_ns_x1[dim], comm->zone_d1[i].max1);
 +            cell_ns_x0[dim] = std::min(cell_ns_x0[dim], comm->zone_d1[i].min0);
 +            cell_ns_x1[dim] = std::max(cell_ns_x1[dim], comm->zone_d1[i].max1);
          }
      }
      if (dd->ndim >= 3)
                  {
                      print_ddzone(debug, 2, i, j, &comm->zone_d2[i][j]);
                  }
 -                cell_ns_x0[dim] = min(cell_ns_x0[dim], comm->zone_d2[i][j].min0);
 -                cell_ns_x1[dim] = max(cell_ns_x1[dim], comm->zone_d2[i][j].max1);
 +                cell_ns_x0[dim] = std::min(cell_ns_x0[dim], comm->zone_d2[i][j].min0);
 +                cell_ns_x1[dim] = std::max(cell_ns_x1[dim], comm->zone_d2[i][j].max1);
              }
          }
      }
@@@ -1556,6 -1527,10 +1556,6 @@@ static void dd_collect_vec_gatherv(gmx_
  void dd_collect_vec(gmx_domdec_t *dd,
                      t_state *state_local, rvec *lv, rvec *v)
  {
 -    gmx_domdec_master_t *ma;
 -    int                  n, i, c, a, nalloc = 0;
 -    rvec                *buf = NULL;
 -
      dd_collect_cg(dd, state_local);
  
      if (dd->nnodes <= GMX_DD_NNODES_SENDRECV)
@@@ -1777,7 -1752,7 +1777,7 @@@ static void dd_distribute_vec_scatterv(
  {
      gmx_domdec_master_t *ma;
      int                 *scounts = NULL, *disps = NULL;
 -    int                  n, i, c, a, nalloc = 0;
 +    int                  n, i, c, a;
      rvec                *buf = NULL;
  
      if (DDMASTER(dd))
@@@ -1971,7 -1946,7 +1971,7 @@@ static void write_dd_grid_pdb(const cha
          snew(grid_r, 2*dd->nnodes);
      }
  
 -    dd_gather(dd, 2*sizeof(rvec), grid_s[0], DDMASTER(dd) ? grid_r[0] : NULL);
 +    dd_gather(dd, 2*sizeof(rvec), grid_s, DDMASTER(dd) ? grid_r : NULL);
  
      if (DDMASTER(dd))
      {
@@@ -2114,15 -2089,15 +2114,15 @@@ real dd_cutoff_mbody(gmx_domdec_t *dd
              r = comm->cellsize_min[dd->dim[0]];
              for (di = 1; di < dd->ndim; di++)
              {
 -                r = min(r, comm->cellsize_min[dd->dim[di]]);
 +                r = std::min(r, comm->cellsize_min[dd->dim[di]]);
              }
              if (comm->bBondComm)
              {
 -                r = max(r, comm->cutoff_mbody);
 +                r = std::max(r, comm->cutoff_mbody);
              }
              else
              {
 -                r = min(r, comm->cutoff);
 +                r = std::min(r, comm->cutoff);
              }
          }
      }
@@@ -2136,7 -2111,7 +2136,7 @@@ real dd_cutoff_twobody(gmx_domdec_t *dd
  
      r_mb = dd_cutoff_mbody(dd);
  
 -    return max(dd->comm->cutoff, r_mb);
 +    return std::max(dd->comm->cutoff, r_mb);
  }
  
  
@@@ -2198,7 -2173,7 +2198,7 @@@ static int *dd_pmenodes(t_commrec *cr
  static int gmx_ddcoord2pmeindex(t_commrec *cr, int x, int y, int z)
  {
      gmx_domdec_t *dd;
 -    ivec          coords, coords_pme, nc;
 +    ivec          coords;
      int           slab;
  
      dd = cr->dd;
@@@ -2267,6 -2242,7 +2267,6 @@@ static int dd_simnode2pmenode(t_commre
  {
      gmx_domdec_t      *dd;
      gmx_domdec_comm_t *comm;
 -    ivec               coord, coord_pme;
      int                i;
      int                pmenode = -1;
  
      if (comm->bCartesianPP_PME)
      {
  #ifdef GMX_MPI
 +        ivec coord, coord_pme;
          MPI_Cart_coords(cr->mpi_comm_mysim, sim_nodeid, DIM, coord);
          if (coord[comm->cartpmedim] < dd->nc[comm->cartpmedim])
          {
@@@ -2415,7 -2390,7 +2415,7 @@@ void get_pme_ddnodes(t_commrec *cr, in
  static gmx_bool receive_vir_ener(t_commrec *cr)
  {
      gmx_domdec_comm_t *comm;
 -    int                pmenode, coords[DIM], rank;
 +    int                pmenode;
      gmx_bool           bReceive;
  
      bReceive = TRUE;
          {
              pmenode = dd_simnode2pmenode(cr, cr->sim_nodeid);
  #ifdef GMX_MPI
 +            ivec coords;
              MPI_Cart_coords(cr->mpi_comm_mysim, cr->sim_nodeid, DIM, coords);
              coords[comm->cartpmedim]++;
              if (coords[comm->cartpmedim] < cr->dd->nc[comm->cartpmedim])
              {
 +                int rank;
                  MPI_Cart_rank(cr->mpi_comm_mysim, coords, &rank);
                  if (dd_simnode2pmenode(cr, rank) == pmenode)
                  {
@@@ -2539,8 -2512,12 +2539,8 @@@ static void make_dd_indices(gmx_domdec_
  {
      int          nzone, zone, zone1, cg0, cg1, cg1_p1, cg, cg_gl, a, a_gl;
      int         *zone2cg, *zone_ncg1, *index_gl, *gatindex;
 -    gmx_ga2la_t *ga2la;
 -    char        *bLocalCG;
      gmx_bool     bCGs;
  
 -    bLocalCG = dd->comm->bLocalCG;
 -
      if (dd->nat_tot > dd->gatindex_nalloc)
      {
          dd->gatindex_nalloc = over_alloc_dd(dd->nat_tot);
  static int check_bLocalCG(gmx_domdec_t *dd, int ncg_sys, const char *bLocalCG,
                            const char *where)
  {
 -    int ncg, i, ngl, nerr;
 +    int i, ngl, nerr;
  
      nerr = 0;
      if (bLocalCG == NULL)
@@@ -2765,12 -2742,12 +2765,12 @@@ static real cellsize_min_dlb(gmx_domdec
          /* The cut-off might have changed, e.g. by PME load balacning,
           * from the value used to set comm->cellsize_min, so check it.
           */
 -        cellsize_min = max(cellsize_min, comm->cutoff/comm->cd[dim_ind].np_dlb);
 +        cellsize_min = std::max(cellsize_min, comm->cutoff/comm->cd[dim_ind].np_dlb);
  
          if (comm->bPMELoadBalDLBLimits)
          {
              /* Check for the cut-off limit set by the PME load balancing */
 -            cellsize_min = max(cellsize_min, comm->PMELoadBal_max_cutoff/comm->cd[dim_ind].np_dlb);
 +            cellsize_min = std::max(cellsize_min, comm->PMELoadBal_max_cutoff/comm->cd[dim_ind].np_dlb);
          }
      }
  
@@@ -2793,10 -2770,10 +2793,10 @@@ static real grid_jump_limit(gmx_domdec_
      {
          if (comm->bPMELoadBalDLBLimits)
          {
 -            cutoff = max(cutoff, comm->PMELoadBal_max_cutoff);
 +            cutoff = std::max(cutoff, comm->PMELoadBal_max_cutoff);
          }
 -        grid_jump_limit = max(grid_jump_limit,
 -                              cutoff/comm->cd[dim_ind].np);
 +        grid_jump_limit = std::max(grid_jump_limit,
 +                                   cutoff/comm->cd[dim_ind].np);
      }
  
      return grid_jump_limit;
@@@ -2983,8 -2960,8 +2983,8 @@@ static void init_ddpme(gmx_domdec_t *dd
              {
                  slab = pmeindex % ddpme->nslab;
              }
 -            ddpme->pp_min[slab] = min(ddpme->pp_min[slab], xyz[dimind]);
 -            ddpme->pp_max[slab] = max(ddpme->pp_max[slab], xyz[dimind]);
 +            ddpme->pp_min[slab] = std::min(ddpme->pp_min[slab], xyz[dimind]);
 +            ddpme->pp_max[slab] = std::max(ddpme->pp_max[slab], xyz[dimind]);
          }
      }
  
@@@ -3185,7 -3162,7 +3185,7 @@@ static void set_dd_cell_sizes_slb(gmx_d
                  {
                      npulse[d]++;
                  }
 -                cellsize_min[d] = min(cellsize_min[d], cellsize);
 +                cellsize_min[d] = std::min(cellsize_min[d], cellsize);
              }
              if (setmode == setcellsizeslbLOCAL)
              {
@@@ -3455,7 -3432,7 +3455,7 @@@ static void set_dd_cell_sizes_dlb_root(
                                         gmx_bool bUniform, gmx_int64_t step)
  {
      gmx_domdec_comm_t *comm;
 -    int                ncd, d1, i, j, pos;
 +    int                ncd, d1, i, pos;
      real              *cell_size;
      real               load_aver, load_i, imbalance, change, change_max, sc;
      real               cellsize_limit_f, dist_min_f, dist_min_f_hard, space;
              imbalance = (load_i - load_aver)/(load_aver > 0 ? load_aver : 1);
              /* Determine the change of the cell size using underrelaxation */
              change     = -relax*imbalance;
 -            change_max = max(change_max, max(change, -change));
 +            change_max = std::max(change_max, std::max(change, -change));
          }
          /* Limit the amount of scaling.
           * We need to use the same rescaling for all cells in one row,
@@@ -3641,7 -3618,7 +3641,7 @@@ static void distribute_dd_cell_sizes_dl
                                           gmx_ddbox_t *ddbox)
  {
      gmx_domdec_comm_t *comm;
 -    int                d1, dim1, pos;
 +    int                d1, pos;
  
      comm = dd->comm;
  
@@@ -3929,13 -3906,13 +3929,13 @@@ static void distribute_cg(FILE *fplog, 
  {
      gmx_domdec_master_t *ma;
      int                **tmp_ind = NULL, *tmp_nalloc = NULL;
 -    int                  i, icg, j, k, k0, k1, d, npbcdim;
 +    int                  i, icg, j, k, k0, k1, d;
      matrix               tcm;
 -    rvec                 box_size, cg_cm;
 +    rvec                 cg_cm;
      ivec                 ind;
      real                 nrcg, inv_ncg, pos_d;
      atom_id             *cgindex;
 -    gmx_bool             bUnbounded, bScrew;
 +    gmx_bool             bScrew;
  
      ma = dd->ma;
  
@@@ -4143,9 -4120,9 +4143,9 @@@ static void get_cg_distribution(FILE *f
      }
  
      dd_scatterv(dd,
 -                DDMASTER(dd) ? ma->ibuf : NULL,
 -                DDMASTER(dd) ? ma->ibuf+dd->nnodes : NULL,
 -                DDMASTER(dd) ? ma->cg : NULL,
 +                bMaster ? ma->ibuf : NULL,
 +                bMaster ? ma->ibuf+dd->nnodes : NULL,
 +                bMaster ? ma->cg : NULL,
                  dd->ncg_home*sizeof(int), dd->index_gl);
  
      /* Determine the home charge group sizes */
@@@ -4480,8 -4457,8 +4480,8 @@@ static void calc_cg_move(FILE *fplog, g
                           int *move)
  {
      int      npbcdim;
 -    int      c, i, cg, k, k0, k1, d, dim, dim2, dir, d2, d3, d4, cell_d;
 -    int      mc, cdd, nrcg, ncg_recv, nat_recv, nvs, nvr, nvec, vec;
 +    int      cg, k, k0, k1, d, dim, d2;
 +    int      mc, nrcg;
      int      flag;
      gmx_bool bScrew;
      ivec     dev;
@@@ -4657,15 -4634,17 +4657,15 @@@ static void dd_redistribute_cg(FILE *fp
      int               *move;
      int                npbcdim;
      int                ncg[DIM*2], nat[DIM*2];
 -    int                c, i, cg, k, k0, k1, d, dim, dim2, dir, d2, d3, d4, cell_d;
 -    int                mc, cdd, nrcg, ncg_recv, nat_recv, nvs, nvr, nvec, vec;
 +    int                c, i, cg, k, d, dim, dim2, dir, d2, d3;
 +    int                mc, cdd, nrcg, ncg_recv, nvs, nvr, nvec, vec;
      int                sbuf[2], rbuf[2];
      int                home_pos_cg, home_pos_at, buf_pos;
      int                flag;
      gmx_bool           bV = FALSE, bSDX = FALSE, bCGP = FALSE;
 -    gmx_bool           bScrew;
 -    ivec               dev;
 -    real               inv_ncg, pos_d;
 +    real               pos_d;
      matrix             tcm;
 -    rvec              *cg_cm = NULL, cell_x0, cell_x1, limitd, limit0, limit1, cm_new;
 +    rvec              *cg_cm = NULL, cell_x0, cell_x1, limitd, limit0, limit1;
      atom_id           *cgindex;
      cginfo_mb_t       *cginfo_mb;
      gmx_domdec_comm_t *comm;
      {
          dim      = dd->dim[d];
          ncg_recv = 0;
 -        nat_recv = 0;
          nvr      = 0;
          for (dir = 0; dir < (dd->nc[dim] == 2 ? 1 : 2); dir++)
          {
                               comm->cgcm_state[cdd], nvs,
                               comm->vbuf.v+nvr, i);
              ncg_recv += rbuf[0];
 -            nat_recv += rbuf[1];
              nvr      += i;
          }
  
@@@ -5259,7 -5240,7 +5259,7 @@@ static void get_load_distribution(gmx_d
      gmx_domdec_comm_t *comm;
      gmx_domdec_load_t *load;
      gmx_domdec_root_t *root = NULL;
 -    int                d, dim, cid, i, pos;
 +    int                d, dim, i, pos;
      float              cell_frac = 0, sbuf[DD_NLOAD_MAX];
      gmx_bool           bSepPME;
  
                  for (i = 0; i < dd->nc[dim]; i++)
                  {
                      load->sum += load->load[pos++];
 -                    load->max  = max(load->max, load->load[pos]);
 +                    load->max  = std::max(load->max, load->load[pos]);
                      pos++;
                      if (dd->bGridJump)
                      {
                              /* This direction could not be load balanced properly,
                               * therefore we need to use the maximum iso the average load.
                               */
 -                            load->sum_m = max(load->sum_m, load->load[pos]);
 +                            load->sum_m = std::max(load->sum_m, load->load[pos]);
                          }
                          else
                          {
                              load->sum_m += load->load[pos];
                          }
                          pos++;
 -                        load->cvol_min = min(load->cvol_min, load->load[pos]);
 +                        load->cvol_min = std::min(load->cvol_min, load->load[pos]);
                          pos++;
                          if (d < dd->ndim-1)
                          {
                      }
                      if (bSepPME)
                      {
 -                        load->mdf = max(load->mdf, load->load[pos]);
 +                        load->mdf = std::max(load->mdf, load->load[pos]);
                          pos++;
 -                        load->pme = max(load->pme, load->load[pos]);
 +                        load->pme = std::max(load->pme, load->load[pos]);
                          pos++;
                      }
                  }
@@@ -5781,6 -5762,7 +5781,6 @@@ static void make_load_communicators(gmx
  
  void setup_dd_grid(FILE *fplog, gmx_domdec_t *dd)
  {
 -    gmx_bool                bZYX;
      int                     d, dim, i, j, m;
      ivec                    tmp, s;
      int                     nzone, nzonep;
  static void make_pp_communicator(FILE *fplog, t_commrec *cr, int gmx_unused reorder)
  {
      gmx_domdec_t      *dd;
 +    dd   = cr->dd;
 +
 +#ifdef GMX_MPI
      gmx_domdec_comm_t *comm;
 -    int                i, rank, *buf;
 +    int                rank, *buf;
      ivec               periods;
 -#ifdef GMX_MPI
      MPI_Comm           comm_cart;
 -#endif
  
 -    dd   = cr->dd;
      comm = dd->comm;
  
 -#ifdef GMX_MPI
      if (comm->bCartesianPP)
      {
          /* Set up cartesian communication for the particle-particle part */
                      dd->nc[XX], dd->nc[YY], dd->nc[ZZ]);
          }
  
 -        for (i = 0; i < DIM; i++)
 +        for (int i = 0; i < DIM; i++)
          {
              periods[i] = TRUE;
          }
          /* Determine the master coordinates and rank.
           * The DD master should be the same node as the master of this sim.
           */
 -        for (i = 0; i < dd->nnodes; i++)
 +        for (int i = 0; i < dd->nnodes; i++)
          {
              if (comm->ddindex2simnodeid[i] == 0)
              {
      }
  }
  
 -static void receive_ddindex2simnodeid(t_commrec *cr)
 +static void receive_ddindex2simnodeid(t_commrec gmx_unused *cr)
  {
 +#ifdef GMX_MPI
      gmx_domdec_t      *dd;
 -
      gmx_domdec_comm_t *comm;
 -    int               *buf;
  
      dd   = cr->dd;
      comm = dd->comm;
  
 -#ifdef GMX_MPI
      if (!comm->bCartesianPP_PME && comm->bCartesianPP)
      {
 +        int *buf;
          snew(comm->ddindex2simnodeid, dd->nnodes);
          snew(buf, dd->nnodes);
          if (cr->duty & DUTY_PP)
          {
              buf[dd_index(dd->nc, dd->ci)] = cr->sim_nodeid;
          }
 -#ifdef GMX_MPI
          /* Communicate the ddindex to simulation nodeid index */
          MPI_Allreduce(buf, comm->ddindex2simnodeid, dd->nnodes, MPI_INT, MPI_SUM,
                        cr->mpi_comm_mysim);
 -#endif
          sfree(buf);
      }
  #endif
@@@ -6122,8 -6108,9 +6122,8 @@@ static void split_communicator(FILE *fp
  {
      gmx_domdec_t      *dd;
      gmx_domdec_comm_t *comm;
 -    int                i, rank;
 +    int                i;
      gmx_bool           bDiv[DIM];
 -    ivec               periods;
  #ifdef GMX_MPI
      MPI_Comm           comm_cart;
  #endif
  #ifdef GMX_MPI
      if (comm->bCartesianPP_PME)
      {
 +        int  rank;
 +        ivec periods;
 +
          if (fplog)
          {
              fprintf(fplog, "Will use a Cartesian communicator for PP <-> PME: %d x %d x %d\n", comm->ntot[XX], comm->ntot[YY], comm->ntot[ZZ]);
          }
          MPI_Cart_create(cr->mpi_comm_mysim, DIM, comm->ntot, periods, reorder,
                          &comm_cart);
 -
          MPI_Comm_rank(comm_cart, &rank);
          if (MASTERNODE(cr) && rank != 0)
          {
@@@ -6368,7 -6353,7 +6368,7 @@@ static real *get_slb_frac(FILE *fplog, 
          for (i = 0; i < nc; i++)
          {
              dbl = 0;
 -            sscanf(size_string, "%lf%n", &dbl, &n);
 +            sscanf(size_string, "%20lf%n", &dbl, &n);
              if (dbl == 0)
              {
                  gmx_fatal(FARGS, "Incorrect or not enough DD cell size entries for direction %s: '%s'", dir, size_string);
@@@ -6431,7 -6416,7 +6431,7 @@@ static int dd_getenv(FILE *fplog, cons
      val = getenv(env_var);
      if (val)
      {
 -        if (sscanf(val, "%d", &nst) <= 0)
 +        if (sscanf(val, "%20d", &nst) <= 0)
          {
              nst = 1;
          }
@@@ -6492,7 -6477,7 +6492,7 @@@ static real average_cellsize_min(gmx_do
      {
          d = dd->dim[di];
          /* Check using the initial average cell size */
 -        r = min(r, ddbox->box_size[d]*ddbox->skew_fac[d]/dd->nc[d]);
 +        r = std::min(r, ddbox->box_size[d]*ddbox->skew_fac[d]/dd->nc[d]);
      }
  
      return r;
@@@ -6502,6 -6487,7 +6502,6 @@@ static int check_dlb_support(FILE *fplo
                               const char *dlb_opt, gmx_bool bRecordLoad,
                               unsigned long Flags, t_inputrec *ir)
  {
 -    gmx_domdec_t *dd;
      int           eDLB = -1;
      char          buf[STRLEN];
  
@@@ -6643,10 -6629,10 +6643,10 @@@ gmx_domdec_t *init_domain_decomposition
      gmx_domdec_t      *dd;
      gmx_domdec_comm_t *comm;
      int                recload;
 -    int                d, i, j;
      real               r_2b, r_mb, r_bonded = -1, r_bonded_limit = -1, limit, acs;
      gmx_bool           bC;
      char               buf[STRLEN];
 +    const real         tenPercentMargin = 1.1;
  
      if (fplog)
      {
       * within nstlist steps. Since boundaries are allowed to displace by half
       * a cell size, DD cells should be at least the size of the list buffer.
       */
 -    comm->cellsize_limit = max(comm->cellsize_limit,
 -                               ir->rlistlong - max(ir->rvdw, ir->rcoulomb));
 +    comm->cellsize_limit = std::max(comm->cellsize_limit,
 +                                    ir->rlistlong - std::max(ir->rvdw, ir->rcoulomb));
  
      if (comm->bInterCGBondeds)
      {
              }
              else
              {
 -                comm->cutoff = max(comm->cutoff, comm->cutoff_mbody);
 +                comm->cutoff = std::max(comm->cutoff, comm->cutoff_mbody);
              }
              r_bonded_limit = comm->cutoff_mbody;
          }
               */
              if (Flags & MD_DDBONDCOMM)
              {
 -                if (max(r_2b, r_mb) > comm->cutoff)
 +                if (std::max(r_2b, r_mb) > comm->cutoff)
                  {
 -                    r_bonded        = max(r_2b, r_mb);
 -                    r_bonded_limit  = 1.1*r_bonded;
 +                    r_bonded        = std::max(r_2b, r_mb);
 +                    r_bonded_limit  = tenPercentMargin*r_bonded;
                      comm->bBondComm = TRUE;
                  }
                  else
                  {
                      r_bonded       = r_mb;
 -                    r_bonded_limit = min(1.1*r_bonded, comm->cutoff);
 +                    r_bonded_limit = std::min(tenPercentMargin*r_bonded, comm->cutoff);
                  }
                  /* We determine cutoff_mbody later */
              }
                  /* No special bonded communication,
                   * simply increase the DD cut-off.
                   */
 -                r_bonded_limit     = 1.1*max(r_2b, r_mb);
 +                r_bonded_limit     = tenPercentMargin*std::max(r_2b, r_mb);
                  comm->cutoff_mbody = r_bonded_limit;
 -                comm->cutoff       = max(comm->cutoff, comm->cutoff_mbody);
 +                comm->cutoff       = std::max(comm->cutoff, comm->cutoff_mbody);
              }
          }
-         comm->cellsize_limit = std::max(comm->cellsize_limit, r_bonded_limit);
          if (fplog)
          {
              fprintf(fplog,
                      "Minimum cell size due to bonded interactions: %.3f nm\n",
-                     comm->cellsize_limit);
+                     r_bonded_limit);
          }
 -        comm->cellsize_limit = max(comm->cellsize_limit, r_bonded_limit);
++        comm->cellsize_limit = std::max(comm->cellsize_limit, r_bonded_limit);
      }
  
      if (dd->bInterCGcons && rconstr <= 0)
                  "User supplied maximum distance required for P-LINCS: %.3f nm\n",
                  rconstr);
      }
 -    comm->cellsize_limit = max(comm->cellsize_limit, rconstr);
 +    comm->cellsize_limit = std::max(comm->cellsize_limit, rconstr);
  
      comm->cgs_gl = gmx_mtop_global_cgs(mtop);
  
              if (comm->eDLB != edlbNO)
              {
                  /* Check if this does not limit the scaling */
 -                comm->cutoff_mbody = min(comm->cutoff_mbody, dlb_scale*acs);
 +                comm->cutoff_mbody = std::min(comm->cutoff_mbody, dlb_scale*acs);
              }
              if (!comm->bBondComm)
              {
                  /* Without bBondComm do not go beyond the n.b. cut-off */
 -                comm->cutoff_mbody = min(comm->cutoff_mbody, comm->cutoff);
 +                comm->cutoff_mbody = std::min(comm->cutoff_mbody, comm->cutoff);
                  if (comm->cellsize_limit >= comm->cutoff)
                  {
                      /* We don't loose a lot of efficieny
                  }
              }
              /* Check if we did not end up below our original limit */
 -            comm->cutoff_mbody = max(comm->cutoff_mbody, r_bonded_limit);
 +            comm->cutoff_mbody = std::max(comm->cutoff_mbody, r_bonded_limit);
  
              if (comm->cutoff_mbody > comm->cellsize_limit)
              {
@@@ -7105,7 -7091,7 +7105,7 @@@ static void turn_on_dlb(FILE *fplog, t_
      cellsize_min = comm->cellsize_min[dd->dim[0]];
      for (d = 1; d < dd->ndim; d++)
      {
 -        cellsize_min = min(cellsize_min, comm->cellsize_min[dd->dim[d]]);
 +        cellsize_min = std::min(cellsize_min, comm->cellsize_min[dd->dim[d]]);
      }
  
      if (cellsize_min < comm->cellsize_limit*1.05)
@@@ -7170,6 -7156,8 +7170,6 @@@ void dd_init_bondeds(FILE *fplog
                       t_inputrec *ir, gmx_bool bBCheck, cginfo_mb_t *cginfo_mb)
  {
      gmx_domdec_comm_t *comm;
 -    gmx_bool           bBondComm;
 -    int                d;
  
      dd_make_reverse_top(fplog, dd, mtop, vsite, ir, bBCheck);
  
@@@ -7280,7 -7268,7 +7280,7 @@@ static void print_dd_settings(FILE *fpl
              limit = dd->comm->cellsize_min[XX];
              for (d = 1; d < DIM; d++)
              {
 -                limit = min(limit, dd->comm->cellsize_min[d]);
 +                limit = std::min(limit, dd->comm->cellsize_min[d]);
              }
          }
  
          {
              fprintf(fplog, "%40s  %-7s %6.3f nm\n",
                      "two-body bonded interactions", "(-rdd)",
 -                    max(comm->cutoff, comm->cutoff_mbody));
 +                    std::max(comm->cutoff, comm->cutoff_mbody));
              fprintf(fplog, "%40s  %-7s %6.3f nm\n",
                      "multi-body bonded interactions", "(-rdd)",
 -                    (comm->bBondComm || dd->bGridJump) ? comm->cutoff_mbody : min(comm->cutoff, limit));
 +                    (comm->bBondComm || dd->bGridJump) ? comm->cutoff_mbody : std::min(comm->cutoff, limit));
          }
          if (dd->vsite_comm)
          {
@@@ -7326,7 -7314,7 +7326,7 @@@ static void set_cell_limits_dlb(gmx_dom
  
      /* Determine the maximum number of comm. pulses in one dimension */
  
 -    comm->cellsize_limit = max(comm->cellsize_limit, comm->cutoff_mbody);
 +    comm->cellsize_limit = std::max(comm->cellsize_limit, comm->cutoff_mbody);
  
      /* Determine the maximum required number of grid pulses */
      if (comm->cellsize_limit >= comm->cutoff)
      else
      {
          /* There is no cell size limit */
 -        npulse = max(dd->nc[XX]-1, max(dd->nc[YY]-1, dd->nc[ZZ]-1));
 +        npulse = std::max(dd->nc[XX]-1, std::max(dd->nc[YY]-1, dd->nc[ZZ]-1));
      }
  
      if (!bNoCutOff && npulse > 1)
              dim      = dd->dim[d];
              npulse_d = (int)(1 + dd->nc[dim]*comm->cutoff
                               /(ddbox->box_size[dim]*ddbox->skew_fac[dim]*dlb_scale));
 -            npulse_d_max = max(npulse_d_max, npulse_d);
 +            npulse_d_max = std::max(npulse_d_max, npulse_d);
          }
 -        npulse = min(npulse, npulse_d_max);
 +        npulse = std::min(npulse, npulse_d_max);
      }
  
      /* This env var can override npulse */
      comm->bVacDLBNoLimit = (ir->ePBC == epbcNONE);
      for (d = 0; d < dd->ndim; d++)
      {
 -        comm->cd[d].np_dlb    = min(npulse, dd->nc[dd->dim[d]]-1);
 +        comm->cd[d].np_dlb    = std::min(npulse, dd->nc[dd->dim[d]]-1);
          comm->cd[d].np_nalloc = comm->cd[d].np_dlb;
          snew(comm->cd[d].ind, comm->cd[d].np_nalloc);
 -        comm->maxpulse = max(comm->maxpulse, comm->cd[d].np_dlb);
 +        comm->maxpulse = std::max(comm->maxpulse, comm->cd[d].np_dlb);
          if (comm->cd[d].np_dlb < dd->nc[dd->dim[d]]-1)
          {
              comm->bVacDLBNoLimit = FALSE;
      /* cellsize_limit is set for LINCS in init_domain_decomposition */
      if (!comm->bVacDLBNoLimit)
      {
 -        comm->cellsize_limit = max(comm->cellsize_limit,
 -                                   comm->cutoff/comm->maxpulse);
 +        comm->cellsize_limit = std::max(comm->cellsize_limit,
 +                                        comm->cutoff/comm->maxpulse);
      }
 -    comm->cellsize_limit = max(comm->cellsize_limit, comm->cutoff_mbody);
 +    comm->cellsize_limit = std::max(comm->cellsize_limit, comm->cutoff_mbody);
      /* Set the minimum cell size for each DD dimension */
      for (d = 0; d < dd->ndim; d++)
      {
      }
      if (comm->cutoff_mbody <= 0)
      {
 -        comm->cutoff_mbody = min(comm->cutoff, comm->cellsize_limit);
 +        comm->cutoff_mbody = std::min(comm->cutoff, comm->cellsize_limit);
      }
      if (comm->bDynLoadBal)
      {
@@@ -7499,7 -7487,7 +7499,7 @@@ void set_dd_parameters(FILE *fplog, gmx
      }
      natoms_tot = comm->cgs_gl.index[comm->cgs_gl.nr];
  
 -    dd->ga2la = ga2la_init(natoms_tot, vol_frac*natoms_tot);
 +    dd->ga2la = ga2la_init(natoms_tot, static_cast<int>(vol_frac*natoms_tot));
  }
  
  static gmx_bool test_dd_cutoff(t_commrec *cr,
@@@ -7777,11 -7765,11 +7777,11 @@@ set_dd_corners(const gmx_domdec_t *dd
          c->c[1][1] = comm->cell_x0[dim1];
          if (dd->bGridJump)
          {
 -            c->c[1][1] = max(comm->cell_x0[dim1], comm->zone_d1[1].mch0);
 +            c->c[1][1] = std::max(comm->cell_x0[dim1], comm->zone_d1[1].mch0);
              if (bDistMB)
              {
                  /* For the multi-body distance we need the maximum */
 -                c->bc[1] = max(comm->cell_x0[dim1], comm->zone_d1[1].p1_0);
 +                c->bc[1] = std::max(comm->cell_x0[dim1], comm->zone_d1[1].p1_0);
              }
          }
          /* Set the upper-right corner for rounding */
                          if (j >= 4)
                          {
                              c->c[2][j-4] =
 -                                max(c->c[2][j-4],
 -                                    comm->zone_d2[zones->shift[i][dim0]][zones->shift[i][dim1]].mch0);
 +                                std::max(c->c[2][j-4],
 +                                         comm->zone_d2[zones->shift[i][dim0]][zones->shift[i][dim1]].mch0);
                          }
                      }
                  }
                      {
                          for (j = 0; j < 2; j++)
                          {
 -                            c->bc[2] = max(c->bc[2], comm->zone_d2[i][j].p1_0);
 +                            c->bc[2] = std::max(c->bc[2], comm->zone_d2[i][j].p1_0);
                          }
                      }
                  }
              c->cr1[3] = comm->cell_x1[dim1];
              if (dd->bGridJump)
              {
 -                c->cr1[0] = max(comm->cell_x1[dim1], comm->zone_d1[1].mch1);
 +                c->cr1[0] = std::max(comm->cell_x1[dim1], comm->zone_d1[1].mch1);
                  if (bDistMB)
                  {
                      /* For the multi-body distance we need the maximum */
 -                    c->bcr1 = max(comm->cell_x1[dim1], comm->zone_d1[1].p1_1);
 +                    c->bcr1 = std::max(comm->cell_x1[dim1], comm->zone_d1[1].p1_1);
                  }
              }
          }
@@@ -8115,7 -8103,7 +8115,7 @@@ static void setup_dd_communication(gmx_
  {
      int                    dim_ind, dim, dim0, dim1, dim2, dimd, p, nat_tot;
      int                    nzone, nzone_send, zone, zonei, cg0, cg1;
 -    int                    c, i, j, cg, cg_gl, nrcg;
 +    int                    c, i, cg, cg_gl, nrcg;
      int                   *zone_cg_range, pos_cg, *index_gl, *cgindex, *recv_i;
      gmx_domdec_comm_t     *comm;
      gmx_domdec_zones_t    *zones;
      gmx_domdec_ind_t      *ind;
      cginfo_mb_t           *cginfo_mb;
      gmx_bool               bBondComm, bDist2B, bDistMB, bDistBonded;
 -    real                   r_mb, r_comm2, r_scomm2, r_bcomm2, r_0, r_1, r2inc, inv_ncg;
 +    real                   r_comm2, r_bcomm2;
      dd_corners_t           corners;
      ivec                   tric_dist;
      rvec                  *cg_cm, *normal, *v_d, *v_0 = NULL, *v_1 = NULL, *recv_vr;
  
      for (dim_ind = 0; dim_ind < dd->ndim; dim_ind++)
      {
 -        dim = dd->dim[dim_ind];
 -
          /* Check if we need to use triclinic distances */
          tric_dist[dim_ind] = 0;
          for (i = 0; i <= dim_ind; i++)
                      /* The rvec buffer is also required for atom buffers
                       * of size nrecv in dd_move_x and dd_move_f.
                       */
 -                    i = max(cd->ind[0].nrecv[nzone+1], ind->nrecv[nzone+1]);
 +                    i = std::max(cd->ind[0].nrecv[nzone+1], ind->nrecv[nzone+1]);
                      vec_rvec_check_alloc(&comm->vbuf2, i);
                  }
              }
@@@ -8590,9 -8580,10 +8590,9 @@@ static void set_zones_size(gmx_domdec_
      gmx_domdec_comm_t  *comm;
      gmx_domdec_zones_t *zones;
      gmx_bool            bDistMB;
 -    int                 z, zi, zj0, zj1, d, dim;
 +    int                 z, zi, d, dim;
      real                rcs, rcmbs;
      int                 i, j;
 -    real                size_j, add_tric;
      real                vol;
  
      comm = dd->comm;
                                   * With multiple pulses this will lead
                                   * to a larger zone then strictly necessary.
                                   */
 -                                zones->size[z].x1[dim] = max(zones->size[z].x1[dim],
 -                                                             zones->size[zi].x1[dim]+rcmbs);
 +                                zones->size[z].x1[dim] = std::max(zones->size[z].x1[dim],
 +                                                                  zones->size[zi].x1[dim]+rcmbs);
                              }
                          }
                      }
                  {
                      if (zones->shift[z][dim] > 0)
                      {
 -                        zones->size[z].x1[dim] = max(zones->size[z].x1[dim],
 -                                                     zones->size[zi].x1[dim]+rcs);
 +                        zones->size[z].x1[dim] = std::max(zones->size[z].x1[dim],
 +                                                          zones->size[zi].x1[dim]+rcs);
                      }
                  }
              }
              {
                  for (i = 0; i < DIM; i++)
                  {
 -                    corner_min[i] = min(corner_min[i], corner[i]);
 -                    corner_max[i] = max(corner_max[i], corner[i]);
 +                    corner_min[i] = std::min(corner_min[i], corner[i]);
 +                    corner_max[i] = std::max(corner_max[i], corner[i]);
                  }
              }
          }
@@@ -8946,7 -8937,8 +8946,7 @@@ static int dd_sort_order(gmx_domdec_t *
  {
      gmx_domdec_sort_t *sort;
      gmx_cgsort_t      *cgsort, *sort_i;
 -    int                ncg_new, nsort2, nsort_new, i, *a, moved, *ibuf;
 -    int                sort_last, sort_skip;
 +    int                ncg_new, nsort2, nsort_new, i, *a, moved;
  
      sort = dd->comm->sort;
  
@@@ -9057,7 -9049,7 +9057,7 @@@ static void dd_sort_state(gmx_domdec_t 
                            int ncg_home_old)
  {
      gmx_domdec_sort_t *sort;
 -    gmx_cgsort_t      *cgsort, *sort_i;
 +    gmx_cgsort_t      *cgsort;
      int               *cgindex;
      int                ncg_new, i, *ibuf, cgsize;
      rvec              *vbuf;
@@@ -9314,7 -9306,7 +9314,7 @@@ void dd_partition_system(FIL
      t_block           *cgs_gl;
      gmx_int64_t        step_pcoupl;
      rvec               cell_ns_x0, cell_ns_x1;
 -    int                i, j, n, ncgindex_set, ncg_home_old = -1, ncg_moved, nat_f_novirsum;
 +    int                i, n, ncgindex_set, ncg_home_old = -1, ncg_moved, nat_f_novirsum;
      gmx_bool           bBoxChanged, bNStGlobalComm, bDoDLB, bCheckDLB, bTurnOnDLB, bLogLoad;
      gmx_bool           bRedist, bSortCG, bResortAll;
      ivec               ncells_old = {0, 0, 0}, ncells_new = {0, 0, 0}, np;
               * so the extra communication cost is negligible.
               */
              const int nddp_chk_dlb = 100;
 -
              bCheckDLB = (comm->n_load_collect == 0 ||
                           comm->n_load_have % nddp_chk_dlb == nddp_chk_dlb - 1);
          }
index e13434df5543e92f9d6659320b0a1836c13aaabc,0e905041ce57f5456e80830e4bb506dcd1b7b29d..4a47131b235f1339e4968e2aaa0c4b6f81e58ce4
   * To help us fund GROMACS development, we humbly ask that you cite
   * the research papers on the package. Check out http://www.gromacs.org.
   */
 -#ifdef HAVE_CONFIG_H
 -#include <config.h>
 -#endif
 +#include "gmxpre.h"
  
 -#include <string.h>
  #include <math.h>
 -#include "macros.h"
 -#include "sysstuff.h"
 -#include "gromacs/utility/smalloc.h"
 -#include "typedefs.h"
 +#include <stdlib.h>
 +#include <string.h>
 +
 +#include "gromacs/commandline/pargs.h"
 +#include "gromacs/fileio/confio.h"
  #include "gromacs/fileio/gmxfio.h"
 -#include "gromacs/fileio/tpxio.h"
 -#include "gromacs/fileio/trxio.h"
 -#include "gromacs/fileio/trnio.h"
 +#include "gromacs/fileio/pdbio.h"
  #include "gromacs/fileio/tngio.h"
  #include "gromacs/fileio/tngio_for_tools.h"
 -#include "gromacs/commandline/pargs.h"
 -#include "gromacs/fileio/futil.h"
 -#include "gromacs/fileio/pdbio.h"
 -#include "gromacs/fileio/confio.h"
 -#include "names.h"
 -#include "index.h"
 -#include "vec.h"
 +#include "gromacs/fileio/tpxio.h"
 +#include "gromacs/fileio/trnio.h"
 +#include "gromacs/fileio/trxio.h"
  #include "gromacs/fileio/xtcio.h"
 -#include "rmpbc.h"
 -#include "pbc.h"
 -#include "xvgr.h"
 -#include "gmx_ana.h"
 -
 -#include "gmx_fatal.h"
 +#include "gromacs/fileio/xvgr.h"
 +#include "gromacs/gmxana/gmx_ana.h"
 +#include "gromacs/legacyheaders/macros.h"
 +#include "gromacs/legacyheaders/names.h"
 +#include "gromacs/legacyheaders/typedefs.h"
 +#include "gromacs/math/vec.h"
 +#include "gromacs/topology/index.h"
 +#include "gromacs/utility/fatalerror.h"
 +#include "gromacs/utility/futil.h"
 +#include "gromacs/utility/smalloc.h"
  
  #define TIME_EXPLICIT 0
  #define TIME_CONTINUE 1
@@@ -502,7 -506,7 +502,7 @@@ int gmx_trjcat(int argc, char *argv[]
  
  #define NFILE asize(fnm)
  
 -    if (!parse_common_args(&argc, argv, PCA_BE_NICE | PCA_TIME_UNIT, NFILE, fnm,
 +    if (!parse_common_args(&argc, argv, PCA_TIME_UNIT, NFILE, fnm,
                             asize(pa), pa, asize(desc), desc, 0, NULL, &oenv))
      {
          return 0;
      {
          if (nfile_out != nset)
          {
 -            char *buf = strdup(fnms_out[0]);
 +            char *buf = gmx_strdup(fnms_out[0]);
              snew(fnms_out, nset);
              for (i = 0; (i < nset); i++)
              {
          {
              if (ftpout == efTNG)
              {
+                 if (ftpout != ftpin)
+                 {
+                     gmx_fatal(FARGS, "When writing TNG the input file format must also be TNG");
+                 }
                  if (bIndex)
                  {
                      trjtools_gmx_prepare_tng_writing(out_file, 'w', NULL, &trxout,
index a70c1c005c5921f050fbca682919b295ce704b7b,9743fc16a80a507da78fd471234a4a389d612195..c22bdb311da2f4078c7aa4d03992b64750a2cf02
   * To help us fund GROMACS development, we humbly ask that you cite
   * the research papers on the package. Check out http://www.gromacs.org.
   */
 -#ifdef HAVE_CONFIG_H
 -#include <config.h>
 -#endif
 +#include "gmxpre.h"
  
  /* This file is completely threadsafe - please keep it that way! */
  
 +#include "gromacs/legacyheaders/txtdump.h"
 +
  #include <stdio.h>
 -#include "typedefs.h"
 -#include "types/commrec.h"
 -#include "names.h"
 -#include "txtdump.h"
 -#include "vec.h"
 -#include "macros.h"
 -
 -#include "gmx_fatal.h"
 +#include <stdlib.h>
 +
 +#include "gromacs/legacyheaders/macros.h"
 +#include "gromacs/legacyheaders/names.h"
 +#include "gromacs/legacyheaders/typedefs.h"
 +#include "gromacs/legacyheaders/types/commrec.h"
 +#include "gromacs/math/vec.h"
 +#include "gromacs/utility/fatalerror.h"
  #include "gromacs/utility/smalloc.h"
  
  int pr_indent(FILE *fp, int n)
@@@ -792,8 -792,8 +792,8 @@@ static void pr_rotgrp(FILE *fp, int ind
      PR("rot-min-gauss", rotg->min_gaussian);
      PR("rot-eps", rotg->eps);
      PS("rot-fit-method", EROTFIT(rotg->eFittype));
-     PI("rot_potfit_nstep", rotg->PotAngle_nstep);
-     PR("rot_potfit_step", rotg->PotAngle_step);
+     PI("rot-potfit-nstep", rotg->PotAngle_nstep);
+     PR("rot-potfit-step", rotg->PotAngle_step);
  }
  
  static void pr_rot(FILE *fp, int indent, t_rot *rot)
@@@ -1411,7 -1411,7 +1411,7 @@@ void pr_ilist(FILE *fp, int indent, con
                  j++;
                  for (k = 0; k < interaction_function[ftype].nratoms; k++)
                  {
 -                    (void) fprintf(fp, " %u", *(iatoms++));
 +                    (void) fprintf(fp, " %d", *(iatoms++));
                  }
                  (void) fprintf(fp, "\n");
                  i += 1+interaction_function[ftype].nratoms;
@@@ -1554,13 -1554,13 +1554,13 @@@ static void low_pr_blocka(FILE *fp, in
          for (i = 0; i <= block->nr; i++)
          {
              (void) pr_indent(fp, indent+INDENT);
 -            (void) fprintf(fp, "%s->index[%d]=%u\n",
 +            (void) fprintf(fp, "%s->index[%d]=%d\n",
                             title, bShowNumbers ? i : -1, block->index[i]);
          }
          for (i = 0; i < block->nra; i++)
          {
              (void) pr_indent(fp, indent+INDENT);
 -            (void) fprintf(fp, "%s->a[%d]=%u\n",
 +            (void) fprintf(fp, "%s->a[%d]=%d\n",
                             title, bShowNumbers ? i : -1, block->a[i]);
          }
      }
@@@ -1641,7 -1641,7 +1641,7 @@@ void pr_blocka(FILE *fp, int indent, co
                          (void) fprintf(fp, "\n");
                          size = pr_indent(fp, indent+INDENT);
                      }
 -                    size += fprintf(fp, "%u", block->a[j]);
 +                    size += fprintf(fp, "%d", block->a[j]);
                  }
                  (void) fprintf(fp, "}\n");
                  start = end;
index 5c8329344c8434cb468cda962731219e7d12259e,0387dd1483353b22b6f26afddde41037e57d32f8..a2763fe1faae93d9bbb12fe98693d47c9bb26bf8
   * To help us fund GROMACS development, we humbly ask that you cite
   * the research papers on the package. Check out http://www.gromacs.org.
   */
 -#ifdef HAVE_CONFIG_H
 -#include <config.h>
 -#endif
 +#include "gmxpre.h"
 +
 +#include "calc_verletbuf.h"
  
  #include <assert.h>
 +#include <math.h>
 +#include <stdlib.h>
  
  #include <sys/types.h>
 -#include <math.h>
 -#include "typedefs.h"
 -#include "physics.h"
 +
 +#include "gromacs/ewald/ewald-util.h"
 +#include "gromacs/legacyheaders/macros.h"
 +#include "gromacs/legacyheaders/typedefs.h"
 +#include "gromacs/math/units.h"
 +#include "gromacs/math/vec.h"
 +#include "gromacs/mdlib/nbnxn_consts.h"
 +#include "gromacs/utility/fatalerror.h"
  #include "gromacs/utility/smalloc.h"
 -#include "gmx_fatal.h"
 -#include "macros.h"
 -#include "vec.h"
 -#include "coulomb.h"
 -#include "calc_verletbuf.h"
 -#include "../mdlib/nbnxn_consts.h"
  
  #ifdef GMX_NBNXN_SIMD
  /* The include below sets the SIMD instruction type (precision+width)
@@@ -208,7 -207,7 +208,7 @@@ static void get_vsite_masses(const gmx_
              for (i = 0; i < il->nr; i += 1+NRAL(ft))
              {
                  const t_iparams *ip;
 -                real             cam[5], inv_mass, coeff, m_aj;
 +                real             cam[5] = {0}, inv_mass, coeff, m_aj;
                  int              a1, j, aj;
  
                  ip = &ffparams->iparams[il->iatoms[i]];
@@@ -414,7 -413,6 +414,7 @@@ static void get_verlet_buffer_atomtypes
              add_at(&att, &natt, &prop[a], nmol);
          }
  
 +        /* cppcheck-suppress uninitvar Fixed in cppcheck 1.65 */
          sfree(vsite_m);
          sfree(prop);
      }
@@@ -556,13 -554,18 +556,18 @@@ static real ener_drift(const verletbuf_
                         real r_buffer,
                         real rlist, real boxvol)
  {
-     double drift_tot, pot1, pot2, pot3, pot;
-     int    i, j;
-     real   s2i_2d, s2i_3d, s2j_2d, s2j_3d, s2, s;
-     int    ti, tj;
-     real   md1, d2, md3;
-     real   sc_fac, rsh, rsh2;
-     double c_exp, c_erfc;
+     /* Erfc(8)=1e-29, use this limit so we have some space for arithmetic
+      * on the result when using float precision.
+      */
+     const real erfc_arg_max = 8.0;
+     double     drift_tot, pot1, pot2, pot3, pot;
+     int        i, j;
+     real       s2i_2d, s2i_3d, s2j_2d, s2j_3d, s2, s;
+     int        ti, tj;
+     real       md1, d2, md3;
+     real       sc_fac, rsh, rsh2;
+     double     c_exp, c_erfc;
  
      drift_tot = 0;
  
  
              rsh    = r_buffer;
              sc_fac = 1.0;
-             /* For constraints: adapt r and scaling for the Gaussian */
-             if (att[i].prop.bConstr)
-             {
-                 real sh, sc;
  
-                 approx_2dof(s2i_2d, r_buffer*s2i_2d/s2, &sh, &sc);
-                 rsh    += sh;
-                 sc_fac *= sc;
+             if (rsh*rsh > 2*s2*erfc_arg_max*erfc_arg_max)
+             {
+                 /* Erfc might run out of float and become 0, somewhat before
+                  * c_exp becomes 0. To avoid this and to avoid NaN in
+                  * approx_2dof, we set both c_expc and c_erfc to zero.
+                  * In any relevant case this has no effect on the results,
+                  * since c_exp < 6e-29, so the displacement is completely
+                  * negligible for such atom pairs (and an overestimate).
+                  * In nearly all use cases, there will be other atom
+                  * pairs that contribute much more to the total, so zeroing
+                  * this particular contribution has no effect at all.
+                  */
+                 c_exp  = 0;
+                 c_erfc = 0;
              }
-             if (att[j].prop.bConstr)
+             else
              {
-                 real sh, sc;
+                 /* For constraints: adapt r and scaling for the Gaussian */
+                 if (att[i].prop.bConstr)
+                 {
+                     real sh, sc;
  
-                 approx_2dof(s2j_2d, r_buffer*s2j_2d/s2, &sh, &sc);
-                 rsh    += sh;
-                 sc_fac *= sc;
-             }
+                     approx_2dof(s2i_2d, r_buffer*s2i_2d/s2, &sh, &sc);
+                     rsh    += sh;
+                     sc_fac *= sc;
+                 }
+                 if (att[j].prop.bConstr)
+                 {
+                     real sh, sc;
  
-             /* Exact contribution of an atom pair with Gaussian displacement
-              * with sigma s to the energy drift for a potential with
-              * derivative -md and second derivative dd at the cut-off.
-              * The only catch is that for potentials that change sign
-              * near the cut-off there could be an unlucky compensation
-              * of positive and negative energy drift.
-              * Such potentials are extremely rare though.
-              *
-              * Note that pot has unit energy*length, as the linear
-              * atom density still needs to be put in.
-              */
-             c_exp  = exp(-rsh*rsh/(2*s2))/sqrt(2*M_PI);
-             c_erfc = 0.5*gmx_erfc(rsh/(sqrt(2*s2)));
+                     approx_2dof(s2j_2d, r_buffer*s2j_2d/s2, &sh, &sc);
+                     rsh    += sh;
+                     sc_fac *= sc;
+                 }
+                 /* Exact contribution of an atom pair with Gaussian displacement
+                  * with sigma s to the energy drift for a potential with
+                  * derivative -md and second derivative dd at the cut-off.
+                  * The only catch is that for potentials that change sign
+                  * near the cut-off there could be an unlucky compensation
+                  * of positive and negative energy drift.
+                  * Such potentials are extremely rare though.
+                  *
+                  * Note that pot has unit energy*length, as the linear
+                  * atom density still needs to be put in.
+                  */
+                 c_exp  = exp(-rsh*rsh/(2*s2))/sqrt(2*M_PI);
+                 c_erfc = 0.5*gmx_erfc(rsh/(sqrt(2*s2)));
+             }
              s      = sqrt(s2);
              rsh2   = rsh*rsh;
  
                  md1/2*((rsh2 + s2)*c_erfc - rsh*s*c_exp);
              pot2 = sc_fac*
                  d2/6*(s*(rsh2 + 2*s2)*c_exp - rsh*(rsh2 + 3*s2)*c_erfc);
-             pot3 =
+             pot3 = sc_fac*
                  md3/24*((rsh2*rsh2 + 6*rsh2*s2 + 3*s2*s2)*c_erfc - rsh*s*(rsh2 + 5*s2)*c_exp);
              pot = pot1 + pot2 + pot3;
  
@@@ -1048,7 -1070,7 +1072,7 @@@ void calc_verlet_buffer_size(const gmx_
  
          if (debug)
          {
-             fprintf(debug, "ib %3d %3d %3d rb %.3f %dx%d fac %.3f drift %f\n",
+             fprintf(debug, "ib %3d %3d %3d rb %.3f %dx%d fac %.3f drift %.1e\n",
                      ib0, ib, ib1, rb,
                      list_setup->cluster_size_i, list_setup->cluster_size_j,
                      nb_clust_frac_pairs_not_in_list_at_cutoff,
index f6289af803e0604b26d8a4a18f208dcd87f777d6,237ea878dbafb8bcaac6eda4fdc5b4790b75ea30..5d9715f60deebf0169695e9aee0b0b177e874c9f
   * To help us fund GROMACS development, we humbly ask that you cite
   * the research papers on the package. Check out http://www.gromacs.org.
   */
 -#ifdef HAVE_CONFIG_H
 -#include <config.h>
 -#endif
 +#include "gmxpre.h"
  
 -#include "vec.h"
 -#include "gromacs/utility/smalloc.h"
 -#include "readir.h"
 -#include "names.h"
 -#include "gromacs/fileio/futil.h"
  #include "gromacs/fileio/trnio.h"
 -#include "txtdump.h"
 -
 -#include "gmx_fatal.h"
 +#include "gromacs/gmxpreprocess/readir.h"
 +#include "gromacs/legacyheaders/names.h"
 +#include "gromacs/legacyheaders/txtdump.h"
 +#include "gromacs/math/vec.h"
 +#include "gromacs/utility/fatalerror.h"
 +#include "gromacs/utility/futil.h"
 +#include "gromacs/utility/smalloc.h"
  
  static char *RotStr = {"Enforced rotation:"};
  
@@@ -77,15 -80,15 +77,15 @@@ extern char **read_rotparams(int *ninp_
  
      /* read rotation parameters */
      CTYPE("Output frequency for angle, torque and rotation potential energy for the whole group");
-     ITYPE("rot_nstrout",     rot->nstrout, 100);
+     ITYPE("rot-nstrout",     rot->nstrout, 100);
      CTYPE("Output frequency for per-slab data (angles, torques and slab centers)");
-     ITYPE("rot_nstsout",     rot->nstsout, 1000);
+     ITYPE("rot-nstsout",     rot->nstsout, 1000);
      CTYPE("Number of rotation groups");
-     ITYPE("rot_ngroups",     rot->ngrp, 1);
+     ITYPE("rot-ngroups",     rot->ngrp, 1);
  
      if (rot->ngrp < 1)
      {
-         gmx_fatal(FARGS, "rot_ngroups should be >= 1");
+         gmx_fatal(FARGS, "rot-ngroups should be >= 1");
      }
  
      snew(rot->grp, rot->ngrp);
          rotg = &rot->grp[g];
          snew(grpbuf[g], STRLEN);
          CTYPE("Rotation group name");
-         sprintf(buf, "rot_group%d", g);
+         sprintf(buf, "rot-group%d", g);
          STYPE(buf, grpbuf[g], "");
  
          CTYPE("Rotation potential. Can be iso, iso-pf, pm, pm-pf, rm, rm-pf, rm2, rm2-pf, flex, flex-t, flex2, flex2-t");
-         sprintf(buf, "rot_type%d", g);
+         sprintf(buf, "rot-type%d", g);
          ETYPE(buf, rotg->eType, erotg_names);
  
          CTYPE("Use mass-weighting of the rotation group positions");
-         sprintf(buf, "rot_massw%d", g);
+         sprintf(buf, "rot-massw%d", g);
          ETYPE(buf, rotg->bMassW, yesno_names);
  
          CTYPE("Rotation vector, will get normalized");
-         sprintf(buf, "rot_vec%d", g);
+         sprintf(buf, "rot-vec%d", g);
          STYPE(buf, s_vec, "1.0 0.0 0.0");
          string2dvec(s_vec, vec);
          /* Normalize the rotation vector */
          }
          else
          {
-             sprintf(warn_buf, "rot_vec%d = 0", g);
+             sprintf(warn_buf, "rot-vec%d = 0", g);
              warning_error(wi, warn_buf);
          }
          fprintf(stderr, "%s Group %d (%s) normalized rot. vector: %f %f %f\n",
          }
  
          CTYPE("Pivot point for the potentials iso, pm, rm, and rm2 (nm)");
-         sprintf(buf, "rot_pivot%d", g);
+         sprintf(buf, "rot-pivot%d", g);
          STYPE(buf, s_vec, "0.0 0.0 0.0");
          clear_dvec(vec);
          if ( (rotg->eType == erotgISO) || (rotg->eType == erotgPM) || (rotg->eType == erotgRM) || (rotg->eType == erotgRM2) )
          }
  
          CTYPE("Rotation rate (degree/ps) and force constant (kJ/(mol*nm^2))");
-         sprintf(buf, "rot_rate%d", g);
+         sprintf(buf, "rot-rate%d", g);
          RTYPE(buf, rotg->rate, 0.0);
  
-         sprintf(buf, "rot_k%d", g);
+         sprintf(buf, "rot-k%d", g);
          RTYPE(buf, rotg->k, 0.0);
          if (rotg->k <= 0.0)
          {
-             sprintf(warn_buf, "rot_k%d <= 0", g);
+             sprintf(warn_buf, "rot-k%d <= 0", g);
              warning_note(wi, warn_buf);
          }
  
          CTYPE("Slab distance for flexible axis rotation (nm)");
-         sprintf(buf, "rot_slab_dist%d", g);
+         sprintf(buf, "rot-slab-dist%d", g);
          RTYPE(buf, rotg->slab_dist, 1.5);
          if (rotg->slab_dist <= 0.0)
          {
-             sprintf(warn_buf, "rot_slab_dist%d <= 0", g);
+             sprintf(warn_buf, "rot-slab-dist%d <= 0", g);
              warning_error(wi, warn_buf);
          }
  
          CTYPE("Minimum value of Gaussian function for the force to be evaluated (for flex* potentials)");
-         sprintf(buf, "rot_min_gauss%d", g);
+         sprintf(buf, "rot-min-gauss%d", g);
          RTYPE(buf, rotg->min_gaussian, 1e-3);
          if (rotg->min_gaussian <= 0.0)
          {
-             sprintf(warn_buf, "rot_min_gauss%d <= 0", g);
+             sprintf(warn_buf, "rot-min-gauss%d <= 0", g);
              warning_error(wi, warn_buf);
          }
  
          CTYPE("Value of additive constant epsilon' (nm^2) for rm2* and flex2* potentials");
-         sprintf(buf, "rot_eps%d", g);
+         sprintf(buf, "rot-eps%d", g);
          RTYPE(buf, rotg->eps, 1e-4);
          if ( (rotg->eps <= 0.0) && (rotg->eType == erotgRM2 || rotg->eType == erotgFLEX2) )
          {
-             sprintf(warn_buf, "rot_eps%d <= 0", g);
+             sprintf(warn_buf, "rot-eps%d <= 0", g);
              warning_error(wi, warn_buf);
          }
  
          CTYPE("Fitting method to determine angle of rotation group (rmsd, norm, or potential)");
-         sprintf(buf, "rot_fit_method%d", g);
+         sprintf(buf, "rot-fit-method%d", g);
          ETYPE(buf, rotg->eFittype, erotg_fitnames);
          CTYPE("For fit type 'potential', nr. of angles around the reference for which the pot. is evaluated");
-         sprintf(buf, "rot_potfit_nsteps%d", g);
+         sprintf(buf, "rot-potfit-nsteps%d", g);
          ITYPE(buf, rotg->PotAngle_nstep, 21);
          if ( (rotg->eFittype == erotgFitPOT) && (rotg->PotAngle_nstep < 1) )
          {
-             sprintf(warn_buf, "rot_potfit_nsteps%d < 1", g);
+             sprintf(warn_buf, "rot-potfit-nsteps%d < 1", g);
              warning_error(wi, warn_buf);
          }
          CTYPE("For fit type 'potential', distance in degrees between two consecutive angles");
-         sprintf(buf, "rot_potfit_step%d", g);
+         sprintf(buf, "rot-potfit-step%d", g);
          RTYPE(buf, rotg->PotAngle_step, 0.25);
      }
  
index 7338b892ff9ef0df8524da693c80b80609c4e3c0,e8e94fdfdedde8797668401b9a68b3062d17063f..db548c115fe8ae54790f07832793f8a5ca7f62a0
  #ifndef _nbnxn_internal_h
  #define _nbnxn_internal_h
  
 -#include "typedefs.h"
 -#include "nbnxn_simd.h"
 -#include "domdec.h"
 +#include "gromacs/domdec/domdec.h"
 +#include "gromacs/legacyheaders/typedefs.h"
 +#include "gromacs/mdlib/nbnxn_pairlist.h"
 +#include "gromacs/mdlib/nbnxn_simd.h"
  #include "gromacs/timing/cyclecounter.h"
  
  
@@@ -98,6 -97,7 +98,7 @@@ typedef struct 
  typedef struct {
      rvec          c0;               /* The lower corner of the (local) grid        */
      rvec          c1;               /* The upper corner of the (local) grid        */
+     rvec          size;             /* c1 - c0                                     */
      real          atom_density;     /* The atom number density for the local grid  */
  
      gmx_bool      bSimple;          /* Is this grid simple or super/sub            */
index 5c368b2ddf3b0ffc532f0d1e0e94acbc48e63d26,5d182ee11c15a1202243e9bc3f3e4428cffbf2df..856ae32d0f0d8c00a97f736ab14d5d66b64026d1
   * the research papers on the package. Check out http://www.gromacs.org.
   */
  
 -#ifdef HAVE_CONFIG_H
 -#include <config.h>
 -#endif
 +#include "gmxpre.h"
 +
 +#include "nbnxn_search.h"
 +
 +#include "config.h"
  
 +#include <assert.h>
  #include <math.h>
  #include <string.h>
 -#include <assert.h>
  
 -#include "sysstuff.h"
 -#include "gromacs/utility/smalloc.h"
 -#include "types/commrec.h"
 -#include "macros.h"
 +#include "gromacs/legacyheaders/gmx_omp_nthreads.h"
 +#include "gromacs/legacyheaders/macros.h"
 +#include "gromacs/legacyheaders/nrnb.h"
 +#include "gromacs/legacyheaders/ns.h"
 +#include "gromacs/legacyheaders/types/commrec.h"
  #include "gromacs/math/utilities.h"
 -#include "vec.h"
 -#include "pbc.h"
 -#include "nbnxn_consts.h"
 +#include "gromacs/math/vec.h"
 +#include "gromacs/mdlib/nb_verlet.h"
 +#include "gromacs/mdlib/nbnxn_atomdata.h"
 +#include "gromacs/mdlib/nbnxn_consts.h"
 +#include "gromacs/pbcutil/ishift.h"
 +#include "gromacs/pbcutil/pbc.h"
 +#include "gromacs/utility/smalloc.h"
 +
  /* nbnxn_internal.h included gromacs/simd/macros.h */
 -#include "nbnxn_internal.h"
 -#ifdef GMX_NBNXN_SIMD
 +#include "gromacs/mdlib/nbnxn_internal.h"
 +#ifdef GMX_SIMD
  #include "gromacs/simd/vector_operations.h"
  #endif
 -#include "nbnxn_atomdata.h"
 -#include "nbnxn_search.h"
 -#include "gmx_omp_nthreads.h"
 -#include "nrnb.h"
 -#include "ns.h"
 -
 -#include "gromacs/fileio/gmxfio.h"
  
  #ifdef NBNXN_SEARCH_BB_SIMD4
  /* Always use 4-wide SIMD for bounding box calculations */
@@@ -572,6 -571,7 +572,7 @@@ static int set_grid_size_xy(const nbnxn
  
      copy_rvec(corner0, grid->c0);
      copy_rvec(corner1, grid->c1);
+     copy_rvec(size,    grid->size);
  
      return nc_max;
  }
@@@ -1008,7 -1008,6 +1009,6 @@@ static void combine_bounding_box_pairs(
  
  /* Prints the average bb size, used for debug output */
  static void print_bbsizes_simple(FILE                *fp,
-                                  const nbnxn_search_t nbs,
                                   const nbnxn_grid_t  *grid)
  {
      int  c, d;
      }
      dsvmul(1.0/grid->nc, ba, ba);
  
-     fprintf(fp, "ns bb: %4.2f %4.2f %4.2f  %4.2f %4.2f %4.2f rel %4.2f %4.2f %4.2f\n",
-             nbs->box[XX][XX]/grid->ncx,
-             nbs->box[YY][YY]/grid->ncy,
-             nbs->box[ZZ][ZZ]*grid->ncx*grid->ncy/grid->nc,
+     fprintf(fp, "ns bb: grid %4.2f %4.2f %4.2f abs %4.2f %4.2f %4.2f rel %4.2f %4.2f %4.2f\n",
+             grid->sx,
+             grid->sy,
+             grid->na_c/(grid->atom_density*grid->sx*grid->sy),
              ba[XX], ba[YY], ba[ZZ],
-             ba[XX]*grid->ncx/nbs->box[XX][XX],
-             ba[YY]*grid->ncy/nbs->box[YY][YY],
-             ba[ZZ]*grid->nc/(grid->ncx*grid->ncy*nbs->box[ZZ][ZZ]));
+             ba[XX]/grid->sx,
+             ba[YY]/grid->sy,
+             ba[ZZ]/(grid->na_c/(grid->atom_density*grid->sx*grid->sy)));
  }
  
  /* Prints the average bb size, used for debug output */
  static void print_bbsizes_supersub(FILE                *fp,
-                                    const nbnxn_search_t nbs,
                                     const nbnxn_grid_t  *grid)
  {
      int  ns, c, s;
      }
      dsvmul(1.0/ns, ba, ba);
  
-     fprintf(fp, "ns bb: %4.2f %4.2f %4.2f  %4.2f %4.2f %4.2f rel %4.2f %4.2f %4.2f\n",
-             nbs->box[XX][XX]/(grid->ncx*GPU_NSUBCELL_X),
-             nbs->box[YY][YY]/(grid->ncy*GPU_NSUBCELL_Y),
-             nbs->box[ZZ][ZZ]*grid->ncx*grid->ncy/(grid->nc*GPU_NSUBCELL_Z),
+     fprintf(fp, "ns bb: grid %4.2f %4.2f %4.2f abs %4.2f %4.2f %4.2f rel %4.2f %4.2f %4.2f\n",
+             grid->sx/GPU_NSUBCELL_X,
+             grid->sy/GPU_NSUBCELL_Y,
+             grid->na_sc/(grid->atom_density*grid->sx*grid->sy*GPU_NSUBCELL_Z),
              ba[XX], ba[YY], ba[ZZ],
-             ba[XX]*grid->ncx*GPU_NSUBCELL_X/nbs->box[XX][XX],
-             ba[YY]*grid->ncy*GPU_NSUBCELL_Y/nbs->box[YY][YY],
-             ba[ZZ]*grid->nc*GPU_NSUBCELL_Z/(grid->ncx*grid->ncy*nbs->box[ZZ][ZZ]));
+             ba[XX]*GPU_NSUBCELL_X/grid->sx,
+             ba[YY]*GPU_NSUBCELL_Y/grid->sy,
+             ba[ZZ]/(grid->na_sc/(grid->atom_density*grid->sx*grid->sy*GPU_NSUBCELL_Z)));
  }
  
  /* Potentially sorts atoms on LJ coefficients !=0 and ==0.
@@@ -1333,7 -1331,7 +1332,7 @@@ static void sort_columns_simple(const n
          sort_atoms(ZZ, FALSE, dd_zone,
                     nbs->a+ash, na, x,
                     grid->c0[ZZ],
-                    1.0/nbs->box[ZZ][ZZ], ncz*grid->na_sc,
+                    1.0/grid->size[ZZ], ncz*grid->na_sc,
                     sort_work);
  
          /* Fill the ncz cells in this column */
@@@ -1383,13 -1381,14 +1382,13 @@@ static void sort_columns_supersub(cons
                                    int cxy_start, int cxy_end,
                                    int *sort_work)
  {
 -    int  cxy;
 -    int  cx, cy, cz = -1, c = -1, ncz;
 -    int  na, ash, na_c, ind, a;
 -    int  subdiv_z, sub_z, na_z, ash_z;
 -    int  subdiv_y, sub_y, na_y, ash_y;
 -    int  subdiv_x, sub_x, na_x, ash_x;
 +    int        cxy;
 +    int        cx, cy, cz = -1, c = -1, ncz;
 +    int        na, ash, na_c, ind, a;
 +    int        subdiv_z, sub_z, na_z, ash_z;
 +    int        subdiv_y, sub_y, na_y, ash_y;
 +    int        subdiv_x, sub_x, na_x, ash_x;
  
 -    /* cppcheck-suppress unassignedVariable */
      nbnxn_bb_t bb_work_array[2], *bb_work_aligned;
  
      bb_work_aligned = (nbnxn_bb_t *)(((size_t)(bb_work_array+1)) & (~((size_t)15)));
          sort_atoms(ZZ, FALSE, dd_zone,
                     nbs->a+ash, na, x,
                     grid->c0[ZZ],
-                    1.0/nbs->box[ZZ][ZZ], ncz*grid->na_sc,
+                    1.0/grid->size[ZZ], ncz*grid->na_sc,
                     sort_work);
  
          /* This loop goes over the supercells and subcells along z at once */
@@@ -1742,14 -1741,14 +1741,14 @@@ static void calc_cell_indices(const nbn
      {
          if (grid->bSimple)
          {
-             print_bbsizes_simple(debug, nbs, grid);
+             print_bbsizes_simple(debug, grid);
          }
          else
          {
              fprintf(debug, "ns non-zero sub-cells: %d average atoms %.2f\n",
                      grid->nsubc_tot, (a1-a0)/(double)grid->nsubc_tot);
  
-             print_bbsizes_supersub(debug, nbs, grid);
+             print_bbsizes_supersub(debug, grid);
          }
      }
  }
@@@ -1767,7 -1766,7 +1766,7 @@@ static void init_buffer_flags(nbnxn_buf
      }
      for (b = 0; b < flags->nflag; b++)
      {
 -        flags->flag[b] = 0;
 +        bitmask_clear(&(flags->flag[b]));
      }
  }
  
@@@ -2623,7 -2622,7 +2622,7 @@@ static void print_nblist_statistics_sim
      fprintf(fp, "nbl na_sc %d rl %g ncp %d per cell %.1f atoms %.1f ratio %.2f\n",
              nbl->na_sc, rl, nbl->ncj, nbl->ncj/(double)grid->nc,
              nbl->ncj/(double)grid->nc*grid->na_sc,
-             nbl->ncj/(double)grid->nc*grid->na_sc/(0.5*4.0/3.0*M_PI*rl*rl*rl*grid->nc*grid->na_sc/det(nbs->box)));
+             nbl->ncj/(double)grid->nc*grid->na_sc/(0.5*4.0/3.0*M_PI*rl*rl*rl*grid->nc*grid->na_sc/(grid->size[XX]*grid->size[YY]*grid->size[ZZ])));
  
      fprintf(fp, "nbl average j cell list length %.1f\n",
              0.25*nbl->ncj/(double)nbl->nci);
@@@ -2673,7 -2672,7 +2672,7 @@@ static void print_nblist_statistics_sup
      fprintf(fp, "nbl na_c %d rl %g ncp %d per cell %.1f atoms %.1f ratio %.2f\n",
              nbl->na_ci, rl, nbl->nci_tot, nbl->nci_tot/(double)grid->nsubc_tot,
              nbl->nci_tot/(double)grid->nsubc_tot*grid->na_c,
-             nbl->nci_tot/(double)grid->nsubc_tot*grid->na_c/(0.5*4.0/3.0*M_PI*rl*rl*rl*grid->nsubc_tot*grid->na_c/det(nbs->box)));
+             nbl->nci_tot/(double)grid->nsubc_tot*grid->na_c/(0.5*4.0/3.0*M_PI*rl*rl*rl*grid->nsubc_tot*grid->na_c/(grid->size[XX]*grid->size[YY]*grid->size[ZZ])));
  
      fprintf(fp, "nbl average j super cell list length %.1f\n",
              0.25*nbl->ncj4/(double)nbl->nsci);
@@@ -2944,10 -2943,10 +2943,10 @@@ static void make_cluster_list_simple(co
  }
  
  #ifdef GMX_NBNXN_SIMD_4XN
 -#include "nbnxn_search_simd_4xn.h"
 +#include "gromacs/mdlib/nbnxn_search_simd_4xn.h"
  #endif
  #ifdef GMX_NBNXN_SIMD_2XNN
 -#include "nbnxn_search_simd_2xnn.h"
 +#include "gromacs/mdlib/nbnxn_search_simd_2xnn.h"
  #endif
  
  /* Plain C or SIMD4 code for making a pair list of super-cell sci vs scj.
@@@ -4846,7 -4845,7 +4845,7 @@@ static void nbnxn_make_pairlist_part(co
      int               ndistc;
      int               ncpcheck;
      int               gridi_flag_shift = 0, gridj_flag_shift = 0;
 -    unsigned int     *gridj_flag       = NULL;
 +    gmx_bitmask_t    *gridj_flag       = NULL;
      int               ncj_old_i, ncj_old_j;
  
      nbs_cycle_start(&work->cc[enbsCCsearch]);
                                          cbl = nbl->cj[nbl->ncj-1].cj >> gridj_flag_shift;
                                          for (cb = cbf; cb <= cbl; cb++)
                                          {
 -                                            gridj_flag[cb] = 1U<<th;
 +                                            bitmask_init_bit(&gridj_flag[cb], th);
                                          }
                                      }
                                  }
  
          if (bFBufferFlag && nbl->ncj > ncj_old_i)
          {
 -            work->buffer_flags.flag[(gridi->cell0+ci)>>gridi_flag_shift] = 1U<<th;
 +            bitmask_init_bit(&(work->buffer_flags.flag[(gridi->cell0+ci)>>gridi_flag_shift]), th);
          }
      }
  
@@@ -5451,8 -5450,8 +5450,8 @@@ static void reduce_buffer_flags(const n
                                  int                         nsrc,
                                  const nbnxn_buffer_flags_t *dest)
  {
 -    int                 s, b;
 -    const unsigned int *flag;
 +    int            s, b;
 +    gmx_bitmask_t *flag;
  
      for (s = 0; s < nsrc; s++)
      {
  
          for (b = 0; b < dest->nflag; b++)
          {
 -            dest->flag[b] |= flag[b];
 +            bitmask_union(&(dest->flag[b]), flag[b]);
          }
      }
  }
  
  static void print_reduction_cost(const nbnxn_buffer_flags_t *flags, int nout)
  {
 -    int nelem, nkeep, ncopy, nred, b, c, out;
 +    int           nelem, nkeep, ncopy, nred, b, c, out;
 +    gmx_bitmask_t mask_0;
  
      nelem = 0;
      nkeep = 0;
      ncopy = 0;
      nred  = 0;
 +    bitmask_init_bit(&mask_0, 0);
      for (b = 0; b < flags->nflag; b++)
      {
 -        if (flags->flag[b] == 1)
 +        if (bitmask_is_equal(flags->flag[b], mask_0))
          {
              /* Only flag 0 is set, no copy of reduction required */
              nelem++;
              nkeep++;
          }
 -        else if (flags->flag[b] > 0)
 +        else if (!bitmask_is_zero(flags->flag[b]))
          {
              c = 0;
              for (out = 0; out < nout; out++)
              {
 -                if (flags->flag[b] & (1U<<out))
 +                if (bitmask_is_set(flags->flag[b], out))
                  {
                      c++;
                  }
index 6fb10ad92d13a972627c4a1d2c4afed71a0b912d,c4bdfe3892908ec4b191ae80967367bb469df2a7..97712b74282ddb2120257a3777c995218963f615
   * To help us fund GROMACS development, we humbly ask that you cite
   * the research papers on the package. Check out http://www.gromacs.org.
   */
 -#ifdef HAVE_CONFIG_H
 -#include <config.h>
 -#endif
 +#include "gmxpre.h"
  
 +#include "gromacs/legacyheaders/update.h"
  
 -#include <stdio.h>
  #include <math.h>
 +#include <stdio.h>
  
 -#include "types/commrec.h"
 -#include "sysstuff.h"
 -#include "gromacs/utility/smalloc.h"
 -#include "typedefs.h"
 -#include "nrnb.h"
 -#include "physics.h"
 -#include "macros.h"
 -#include "vec.h"
 -#include "main.h"
 -#include "update.h"
 -#include "gromacs/random/random.h"
 -#include "mshift.h"
 -#include "tgroup.h"
 -#include "force.h"
 -#include "names.h"
 -#include "txtdump.h"
 -#include "mdrun.h"
 -#include "constr.h"
 -#include "disre.h"
 -#include "orires.h"
 -#include "gmx_omp_nthreads.h"
 +#include <algorithm>
  
  #include "gromacs/fileio/confio.h"
 -#include "gromacs/fileio/futil.h"
 +#include "gromacs/legacyheaders/constr.h"
 +#include "gromacs/legacyheaders/disre.h"
 +#include "gromacs/legacyheaders/force.h"
 +#include "gromacs/legacyheaders/gmx_omp_nthreads.h"
 +#include "gromacs/legacyheaders/macros.h"
 +#include "gromacs/legacyheaders/mdrun.h"
 +#include "gromacs/legacyheaders/names.h"
 +#include "gromacs/legacyheaders/nrnb.h"
 +#include "gromacs/legacyheaders/orires.h"
 +#include "gromacs/legacyheaders/tgroup.h"
 +#include "gromacs/legacyheaders/txtdump.h"
 +#include "gromacs/legacyheaders/typedefs.h"
 +#include "gromacs/legacyheaders/types/commrec.h"
 +#include "gromacs/math/units.h"
 +#include "gromacs/math/vec.h"
 +#include "gromacs/pbcutil/mshift.h"
 +#include "gromacs/pbcutil/pbc.h"
 +#include "gromacs/pulling/pull.h"
 +#include "gromacs/random/random.h"
  #include "gromacs/timing/wallcycle.h"
 +#include "gromacs/utility/futil.h"
  #include "gromacs/utility/gmxomp.h"
 -#include "gromacs/pulling/pull.h"
 +#include "gromacs/utility/smalloc.h"
  
  /*For debugging, start at v(-dt/2) for velolcity verlet -- uncomment next line */
  /*#define STARTFROMDT2*/
@@@ -264,8 -266,10 +264,8 @@@ static void do_update_vv_vel(int start
                               unsigned short cACC[], rvec v[], rvec f[],
                               gmx_bool bExtended, real veta, real alpha)
  {
 -    double imass, w_dt;
 +    double w_dt;
      int    gf = 0, ga = 0;
 -    rvec   vrel;
 -    real   u, vn, vv, va, vb, vnrel;
      int    n, d;
      double g, mv1, mv2;
  
@@@ -312,6 -316,7 +312,6 @@@ static void do_update_vv_pos(int start
                               rvec x[], rvec xprime[], rvec v[],
                               gmx_bool bExtended, real veta)
  {
 -    double imass, w_dt;
      int    gf = 0;
      int    n, d;
      double g, mr1, mr2;
@@@ -396,6 -401,8 +396,6 @@@ static void do_update_visc(int start, i
              }
              for (d = 0; d < DIM; d++)
              {
 -                vn             = v[n][d];
 -
                  if ((ptype[n] != eptVSite) && (ptype[n] != eptShell))
                  {
                      vn  = (lg*vrel[d] + dt*(imass*f[n][d] - 0.5*vxi*vrel[d]
@@@ -462,7 -469,7 +462,7 @@@ static gmx_stochd_t *init_stochd(t_inpu
  {
      gmx_stochd_t   *sd;
      gmx_sd_const_t *sdc;
 -    int             ngtc, n, th;
 +    int             ngtc, n;
      real            y;
  
      snew(sd, 1);
  
          for (n = 0; n < ngtc; n++)
          {
 -            reft = max(0.0, opts->ref_t[n]);
 +            reft = std::max<real>(0, opts->ref_t[n]);
              if ((opts->tau_t[n] > 0) && (reft > 0))  /* tau_t or ref_t = 0 means that no randomization is done */
              {
                  sd->randomize_group[n] = TRUE;
@@@ -661,6 -668,10 +661,6 @@@ static void do_update_sd1(gmx_stochd_t 
                  {
                      ga  = cACC[n];
                  }
 -                if (cTC)
 -                {
 -                    gt  = cTC[n];
 -                }
  
                  for (d = 0; d < DIM; d++)
                  {
                  {
                      gf  = cFREEZE[n];
                  }
 -                if (cACC)
 -                {
 -                    ga  = cACC[n];
 -                }
                  if (cTC)
                  {
                      gt  = cTC[n];
@@@ -768,6 -783,7 +768,6 @@@ static void do_update_sd2(gmx_stochd_t 
       * half of the update, needs to be remembered for the second half.
       */
      rvec  *sd_V;
 -    real   kT;
      int    gf = 0, ga = 0, gt = 0;
      real   vn = 0, Vmh, Xmh;
      real   ism;
@@@ -1321,7 -1337,7 +1321,7 @@@ void update_tcouple(gmx_int64_t       s
  {
      gmx_bool   bTCouple = FALSE;
      real       dttc;
 -    int        i, start, end, homenr, offset;
 +    int        i, offset;
  
      /* if using vv with trotter decomposition methods, we do this elsewhere in the code */
      if (inputrec->etc != etcNO &&
@@@ -1499,6 -1515,10 +1499,10 @@@ static void combine_forces(gmx_update_
                  {
                      xp[i][d] = state->x[i][d] + fac*f_lr[i][d]*md->invmass[i];
                  }
+                 else
+                 {
+                     xp[i][d] = state->x[i][d];
+                 }
              }
          }
          constrain(NULL, FALSE, FALSE, constr, idef, ir, NULL, cr, step, 0, 1.0, md,
@@@ -1539,11 -1559,12 +1543,11 @@@ void update_constraints(FIL
                          gmx_bool          bCalcVir,
                          real              vetanew)
  {
 -    gmx_bool             bExtended, bLastStep, bLog = FALSE, bEner = FALSE, bDoConstr = FALSE;
 +    gmx_bool             bLastStep, bLog = FALSE, bEner = FALSE, bDoConstr = FALSE;
      double               dt;
 -    real                 dt_1;
 -    int                  start, homenr, nrend, i, n, m, g, d;
 +    int                  start, homenr, nrend, i, m;
      tensor               vir_con;
 -    rvec                *vbuf, *xprime = NULL;
 +    rvec                *xprime = NULL;
      int                  nth, th;
  
      if (constr)
      nrend  = start+homenr;
  
      dt   = inputrec->delta_t;
 -    dt_1 = 1.0/dt;
  
      /*
       *  Steps (7C, 8C)
          }
          else
          {
 +            // cppcheck-suppress unreadVariable
              nth = gmx_omp_nthreads_get(emntUpdate);
  
  #pragma omp parallel for num_threads(nth) schedule(static)
@@@ -1780,11 -1801,20 +1784,11 @@@ void update_box(FILE             *fplog
                  t_nrnb           *nrnb,
                  gmx_update_t      upd)
  {
 -    gmx_bool             bExtended, bLastStep, bLog = FALSE, bEner = FALSE;
      double               dt;
 -    real                 dt_1;
 -    int                  start, homenr, nrend, i, n, m, g;
 -    tensor               vir_con;
 +    int                  start, homenr, i, n, m;
  
      start  = 0;
      homenr = md->homenr;
 -    nrend  = start+homenr;
 -
 -    bExtended =
 -        (inputrec->etc == etcNOSEHOOVER) ||
 -        (inputrec->epc == epcPARRINELLORAHMAN) ||
 -        (inputrec->epc == epcMTTK);
  
      dt = inputrec->delta_t;
  
@@@ -1886,11 -1916,16 +1890,11 @@@ void update_coords(FILE             *fp
                     gmx_constr_t      constr,
                     t_idef           *idef)
  {
 -    gmx_bool          bNH, bPR, bLastStep, bLog = FALSE, bEner = FALSE, bDoConstr = FALSE;
 +    gmx_bool          bNH, bPR, bDoConstr = FALSE;
      double            dt, alpha;
 -    real             *imass, *imassin;
      rvec             *force;
 -    real              dt_1;
 -    int               start, homenr, nrend, i, j, d, n, m, g;
 -    int               blen0, blen1, iatom, jatom, nshake, nsettle, nconstr, nexpand;
 -    int              *icom = NULL;
 -    tensor            vir_con;
 -    rvec             *vcom, *xcom, *vall, *xall, *xin, *vin, *forcein, *fall, *xpall, *xprimein, *xprime;
 +    int               start, homenr, nrend;
 +    rvec             *xprime;
      int               nth, th;
  
      bDoConstr = (NULL != constr);
      xprime = get_xprime(state, upd);
  
      dt   = inputrec->delta_t;
 -    dt_1 = 1.0/dt;
  
      /* We need to update the NMR restraint history when time averaging is used */
      if (state->flags & (1<<estDISRE_RM3TAV))
@@@ -2124,6 -2160,7 +2128,6 @@@ extern gmx_bool update_randomize_veloci
                                              t_mdatoms *md, t_state *state, gmx_update_t upd, gmx_constr_t constr)
  {
  
 -    int  i;
      real rate = (ir->delta_t)/ir->opts.tau_t[0];
  
      if (ir->etc == etcANDERSEN && constr != NULL)
index faa5912527cca004da5178ef6a8046c031bc2408,4ab6b30cb8c7d73b3ec2429a60e5addaa23538d9..5b7358a761d8f770dd0c611431a1b05818175c0f
   * To help us fund GROMACS development, we humbly ask that you cite
   * the research papers on the package. Check out http://www.gromacs.org.
   */
 -#ifdef HAVE_CONFIG_H
 -#include <config.h>
 -#endif
 +#include "gmxpre.h"
 +
 +#include "pull_rotation.h"
 +
 +#include "config.h"
  
  #include <stdio.h>
  #include <stdlib.h>
  #include <string.h>
  
 -#include "domdec.h"
 -#include "gromacs/utility/smalloc.h"
 -#include "network.h"
 -#include "pbc.h"
 -#include "mdrun.h"
 -#include "txtdump.h"
 -#include "names.h"
 -#include "mtop_util.h"
 -#include "names.h"
 -#include "vec.h"
 -#include "gmx_ga2la.h"
 -#include "xvgr.h"
 -#include "copyrite.h"
 -#include "macros.h"
 -
 -#include "gromacs/fileio/futil.h"
 +#include "gromacs/domdec/domdec.h"
  #include "gromacs/fileio/gmxfio.h"
  #include "gromacs/fileio/trnio.h"
 +#include "gromacs/fileio/xvgr.h"
 +#include "gromacs/legacyheaders/copyrite.h"
 +#include "gromacs/legacyheaders/gmx_ga2la.h"
 +#include "gromacs/legacyheaders/macros.h"
 +#include "gromacs/legacyheaders/mdrun.h"
 +#include "gromacs/legacyheaders/names.h"
 +#include "gromacs/legacyheaders/network.h"
 +#include "gromacs/legacyheaders/txtdump.h"
  #include "gromacs/linearalgebra/nrjac.h"
 +#include "gromacs/math/utilities.h"
 +#include "gromacs/math/vec.h"
 +#include "gromacs/mdlib/groupcoord.h"
 +#include "gromacs/pbcutil/pbc.h"
  #include "gromacs/timing/cyclecounter.h"
  #include "gromacs/timing/wallcycle.h"
 +#include "gromacs/topology/mtop_util.h"
 +#include "gromacs/utility/futil.h"
  #include "gromacs/utility/qsort_threadsafe.h"
 -#include "gromacs/pulling/pull_rotation.h"
 -#include "gromacs/mdlib/groupcoord.h"
 -#include "gromacs/math/utilities.h"
 +#include "gromacs/utility/smalloc.h"
  
  static char *RotStr = {"Enforced rotation:"};
  
@@@ -873,7 -874,7 +873,7 @@@ static FILE *open_rot_out(const char *f
          fp = xvgropen(fn, "Rotation angles and energy", "Time (ps)", "angles (degrees) and energies (kJ/mol)", oenv);
          fprintf(fp, "# Output of enforced rotation data is written in intervals of %d time step%s.\n#\n", rot->nstrout, rot->nstrout > 1 ? "s" : "");
          fprintf(fp, "# The scalar tau is the torque (kJ/mol) in the direction of the rotation vector v.\n");
-         fprintf(fp, "# To obtain the vectorial torque, multiply tau with the group's rot_vec.\n");
+         fprintf(fp, "# To obtain the vectorial torque, multiply tau with the group's rot-vec.\n");
          fprintf(fp, "# For flexible groups, tau(t,n) from all slabs n have been summed in a single value tau(t) here.\n");
          fprintf(fp, "# The torques tau(t,n) are found in the rottorque.log (-rt) output file\n");
  
  
              fprintf(fp, "#\n");
              fprintf(fp, "# ROTATION GROUP %d, potential type '%s':\n", g, erotg_names[rotg->eType]);
-             fprintf(fp, "# rot_massw%d          %s\n", g, yesno_names[rotg->bMassW]);
-             fprintf(fp, "# rot_vec%d            %12.5e %12.5e %12.5e\n", g, rotg->vec[XX], rotg->vec[YY], rotg->vec[ZZ]);
-             fprintf(fp, "# rot_rate%d           %12.5e degrees/ps\n", g, rotg->rate);
-             fprintf(fp, "# rot_k%d              %12.5e kJ/(mol*nm^2)\n", g, rotg->k);
+             fprintf(fp, "# rot-massw%d          %s\n", g, yesno_names[rotg->bMassW]);
+             fprintf(fp, "# rot-vec%d            %12.5e %12.5e %12.5e\n", g, rotg->vec[XX], rotg->vec[YY], rotg->vec[ZZ]);
+             fprintf(fp, "# rot-rate%d           %12.5e degrees/ps\n", g, rotg->rate);
+             fprintf(fp, "# rot-k%d              %12.5e kJ/(mol*nm^2)\n", g, rotg->k);
              if (rotg->eType == erotgISO || rotg->eType == erotgPM || rotg->eType == erotgRM || rotg->eType == erotgRM2)
              {
-                 fprintf(fp, "# rot_pivot%d          %12.5e %12.5e %12.5e  nm\n", g, rotg->pivot[XX], rotg->pivot[YY], rotg->pivot[ZZ]);
+                 fprintf(fp, "# rot-pivot%d          %12.5e %12.5e %12.5e  nm\n", g, rotg->pivot[XX], rotg->pivot[YY], rotg->pivot[ZZ]);
              }
  
              if (bFlex)
              {
-                 fprintf(fp, "# rot_slab_distance%d   %f nm\n", g, rotg->slab_dist);
-                 fprintf(fp, "# rot_min_gaussian%d   %12.5e\n", g, rotg->min_gaussian);
+                 fprintf(fp, "# rot-slab-distance%d   %f nm\n", g, rotg->slab_dist);
+                 fprintf(fp, "# rot-min-gaussian%d   %12.5e\n", g, rotg->min_gaussian);
              }
  
              /* Output the centers of the rotation groups for the pivot-free potentials */
  
              if ( (rotg->eType == erotgRM2) || (rotg->eType == erotgFLEX2) || (rotg->eType == erotgFLEX2T) )
              {
-                 fprintf(fp, "# rot_eps%d            %12.5e nm^2\n", g, rotg->eps);
+                 fprintf(fp, "# rot-eps%d            %12.5e nm^2\n", g, rotg->eps);
              }
              if (erotgFitPOT == rotg->eFittype)
              {
              add_to_string_aligned(&LegendStr, buf);
  
              sprintf(buf2, "%s (degrees)", buf);
 -            setname[nsets] = strdup(buf2);
 +            setname[nsets] = gmx_strdup(buf2);
              nsets++;
          }
          for (g = 0; g < rot->ngrp; g++)
              }
              add_to_string_aligned(&LegendStr, buf);
              sprintf(buf2, "%s (degrees)", buf);
 -            setname[nsets] = strdup(buf2);
 +            setname[nsets] = gmx_strdup(buf2);
              nsets++;
  
              sprintf(buf, "tau%d", g);
              add_to_string_aligned(&LegendStr, buf);
              sprintf(buf2, "%s (kJ/mol)", buf);
 -            setname[nsets] = strdup(buf2);
 +            setname[nsets] = gmx_strdup(buf2);
              nsets++;
  
              sprintf(buf, "energy%d", g);
              add_to_string_aligned(&LegendStr, buf);
              sprintf(buf2, "%s (kJ/mol)", buf);
 -            setname[nsets] = strdup(buf2);
 +            setname[nsets] = gmx_strdup(buf2);
              nsets++;
          }
          fprintf(fp, "#\n");
@@@ -1106,7 -1107,7 +1106,7 @@@ static FILE *open_torque_out(const cha
                  fprintf(fp, "# Rotation group %d (%s), slab distance %f nm.\n", g, erotg_names[rotg->eType], rotg->slab_dist);
                  fprintf(fp, "# The scalar tau is the torque (kJ/mol) in the direction of the rotation vector.\n");
                  fprintf(fp, "# To obtain the vectorial torque, multiply tau with\n");
-                 fprintf(fp, "# rot_vec%d            %10.3e %10.3e %10.3e\n", g, rotg->vec[XX], rotg->vec[YY], rotg->vec[ZZ]);
+                 fprintf(fp, "# rot-vec%d            %10.3e %10.3e %10.3e\n", g, rotg->vec[XX], rotg->vec[YY], rotg->vec[ZZ]);
                  fprintf(fp, "#\n");
              }
          }