Merge branch 'origin/release-2020' into merge-2020-into-2021
authorPaul Bauer <paul.bauer.q@gmail.com>
Tue, 24 Nov 2020 13:38:39 +0000 (14:38 +0100)
committerPaul Bauer <paul.bauer.q@gmail.com>
Tue, 24 Nov 2020 14:15:08 +0000 (15:15 +0100)
Resolved Conflicts:
admin/containers/scripted_gmx_docker_builds.py
docs/conf.cmakein.py
docs/dev-manual/change-management.rst
docs/dev-manual/contribute.rst
docs/dev-manual/releng/index.rst
docs/dev-manual/tools.rst
docs/install-guide/index.rst
docs/reference-manual/preface.rst
docs/release-notes/index.rst
scripts/GMXRC.bash.cmakein
src/gromacs/gmxana/gmx_covar.cpp
src/gromacs/gmxana/gmx_do_dssp.cpp
src/gromacs/gmxana/gmx_h2order.cpp
src/gromacs/gmxana/nsfactor.cpp
src/gromacs/hardware/detecthardware.cpp
src/gromacs/mdlib/coupling.cpp
src/gromacs/mdrun/md.cpp
src/gromacs/mdtypes/state.h
src/gromacs/trajectoryanalysis/modules/pairdist.cpp
src/gromacs/utility/basedefinitions.h

20 files changed:
1  2 
admin/containers/scripted_gmx_docker_builds.py
docs/conf.cmakein.py
docs/dev-manual/contribute.rst
docs/dev-manual/releng/index.rst
docs/dev-manual/reportstyle.rst
docs/dev-manual/tools.rst
docs/install-guide/index.rst
docs/reference-manual/preface.rst
docs/release-notes/2016/major/bugs-fixed.rst
docs/release-notes/index.rst
src/gromacs/domdec/partition.cpp
src/gromacs/gmxana/gmx_covar.cpp
src/gromacs/gmxana/gmx_do_dssp.cpp
src/gromacs/gmxana/gmx_h2order.cpp
src/gromacs/hardware/detecthardware.cpp
src/gromacs/hardware/prepare_detection.cpp
src/gromacs/mdlib/coupling.cpp
src/gromacs/mdlib/mdebin_bar.cpp
src/gromacs/mdtypes/state.h
src/gromacs/trajectoryanalysis/modules/pairdist.cpp

index 5464a7f17c5b62a9ac3519cb5c66e5cd1a4af9a8,26a20c720a16eb6eb74a0f9aac51945e98ed9672..1b820805bdce2683f03e6fb535d1705fa76a237e
mode 100755,100644..100755
@@@ -33,8 -33,7 +33,8 @@@
  # To help us fund GROMACS development, we humbly ask that you cite
  # the research papers on the package. Check out http://www.gromacs.org.
  
 -"""
 +"""Building block based Dockerfile generation for CI testing images.
 +
  Generates a set of docker images used for running GROMACS CI on Gitlab.
  The images are prepared according to a selection of build configuration targets
  that hope to cover a broad enough scope of different possible systems,
@@@ -45,9 -44,6 +45,9 @@@ analysing the logic and adding build st
  
  Based on the example script provided by the NVidia HPCCM repository.
  
 +Reference:
 +    `NVidia HPC Container Maker <https://github.com/NVIDIA/hpc-container-maker>`__
 +
  Authors:
      * Paul Bauer <paul.bauer.q@gmail.com>
      * Eric Irrgang <ericirrgang@gmail.com>
@@@ -60,9 -56,6 +60,9 @@@ Usage:
      $ python3 scripted_gmx_docker_builds.py --format docker > Dockerfile && docker build .
      $ python3 scripted_gmx_docker_builds.py | docker build -
  
 +See Also:
 +    :file:`buildall.sh`
 +
  """
  
  import argparse
@@@ -87,7 -80,6 +87,7 @@@ _common_packages = ['build-essential'
                      'ccache',
                      'git',
                      'gnupg',
 +                    'gpg-agent',
                      'libfftw3-dev',
                      'libhwloc-dev',
                      'liblapack-dev',
                      'wget',
                      'xsltproc']
  
 +_opencl_extra_packages = [
 +    'nvidia-opencl-dev',
 +    # The following require apt_ppas=['ppa:intel-opencl/intel-opencl']
 +    'intel-opencl-icd',
 +    'ocl-icd-libopencl1',
 +    'ocl-icd-opencl-dev',
 +    'opencl-headers',
 +    # The following require
 +    #             apt_keys=['http://repo.radeon.com/rocm/apt/debian/rocm.gpg.key'],
 +    #             apt_repositories=['deb [arch=amd64] http://repo.radeon.com/rocm/apt/debian/ xenial main']
 +    'libelf1',
 +    'rocm-opencl',
 +    'rocm-dev',
 +    'clinfo'
 +]
 +
  # Extra packages needed to build Python installations from source.
  _python_extra_packages = ['build-essential',
                            'ca-certificates',
@@@ -151,20 -127,25 +151,20 @@@ _docs_extra_packages = ['autoconf'
                          'linkchecker',
                          'mscgen',
                          'm4',
 +                        'openssh-client',
                          'texinfo',
                          'texlive-latex-base',
                          'texlive-latex-extra',
                          'texlive-fonts-recommended',
                          'texlive-fonts-extra']
  
 -# Supported Python versions for maintained branches.
 -# TODO: Remove '3.5.9' from defaults in master once script in release-2020 diverges.
 -_python_versions = ['3.5.9', '3.6.10', '3.7.7', '3.8.2']
 -
  # Parse command line arguments
 -parser = argparse.ArgumentParser(description='GROMACS CI image creation script', parents=[utility.parser])
 +parser = argparse.ArgumentParser(description='GROMACS CI image creation script',
 +                                 parents=[utility.parser])
  
  parser.add_argument('--format', type=str, default='docker',
                      choices=['docker', 'singularity'],
                      help='Container specification format (default: docker)')
 -parser.add_argument('--venvs', nargs='*', type=str, default=_python_versions,
 -                    help='List of Python versions ("major.minor.patch") for which to install venvs. '
 -                         'Default: {}'.format(' '.join(_python_versions)))
  
  
  def base_image_tag(args) -> str:
@@@ -193,20 -174,12 +193,18 @@@ def get_llvm_packages(args) -> typing.I
      # If we use the package version of LLVM, we need to install extra packages for it.
      if (args.llvm is not None) and (args.tsan is None):
          return ['libomp-dev',
 +                'libomp5',
                  'clang-format-' + str(args.llvm),
                  'clang-tidy-' + str(args.llvm)]
      else:
          return []
  
 +def get_opencl_packages(args) -> typing.Iterable[str]:
 +    if (args.doxygen is None) and (args.oneapi is None):
 +        return _opencl_extra_packages
 +    else:
 +        return []
  
  def get_compiler(args, compiler_build_stage: hpccm.Stage = None) -> bb_base:
      # Compiler
      if args.icc is not None:
          else:
              compiler = hpccm.building_blocks.llvm(extra_repository=True, version=args.llvm)
  
 +    elif args.oneapi is not None:
 +        if compiler_build_stage is not None:
 +            compiler = compiler_build_stage.runtime(_from='oneapi')
 +            # Prepare the toolchain (needed only for builds done within the Dockerfile, e.g.
 +            # OpenMPI builds, which don't currently work for other reasons)
 +            oneapi_toolchain = hpccm.toolchain(CC='/opt/intel/oneapi/compiler/latest/linux/bin/intel64/icc',
 +                                               CXX='/opt/intel/oneapi/compiler/latest/linux/bin/intel64/icpc')
 +            setattr(compiler, 'toolchain', oneapi_toolchain)
 +
 +        else:
 +            raise RuntimeError('No oneAPI compiler build stage!')
 +
      elif args.gcc is not None:
          compiler = hpccm.building_blocks.gnu(extra_repository=True,
                                               version=args.gcc,
@@@ -253,17 -214,11 +251,17 @@@ def get_mpi(args, compiler)
                  use_cuda = True
  
              if hasattr(compiler, 'toolchain'):
 +                if args.oneapi is not None:
 +                    raise RuntimeError('oneAPI building OpenMPI is not supported')
                  return hpccm.building_blocks.openmpi(toolchain=compiler.toolchain, cuda=use_cuda, infiniband=False)
              else:
                  raise RuntimeError('compiler is not an HPCCM compiler building block!')
  
          elif args.mpi == 'impi':
 +            # TODO Intel MPI from the oneAPI repo is not working reliably,
 +            # reasons are unclear. When solved, add packagages called:
 +            # 'intel-oneapi-mpi', 'intel-oneapi-mpi-devel'
 +            # during the compiler stage.
              # TODO also consider hpccm's intel_mpi package if that doesn't need
              # a license to run.
              raise RuntimeError('Intel MPI recipe not implemented yet.')
          return None
  
  
 -def get_opencl(args):
 -    # Add OpenCL environment if needed
 -    if (args.opencl is not None):
 -        if args.opencl == 'nvidia':
 -            if (args.cuda is None):
 -                raise RuntimeError('Need Nvidia environment for Nvidia OpenCL image')
 -
 -            return hpccm.building_blocks.packages(ospackages=['nvidia-opencl-dev'])
 -
 -        elif args.opencl == 'intel':
 -            return hpccm.building_blocks.packages(
 -                    apt_ppas=['ppa:intel-opencl/intel-opencl'],
 -                    ospackages=['opencl-headers', 'ocl-icd-libopencl1',
 -                                'ocl-icd-opencl-dev', 'intel-opencl-icd'])
 -
 -        elif args.opencl == 'amd':
 -            # libelf1 is a necessary dependency for something in the ROCm stack,
 -            # which they should set up, but seem to have omitted.
 -            return hpccm.building_blocks.packages(
 -                    apt_keys=['http://repo.radeon.com/rocm/apt/debian/rocm.gpg.key'],
 -                    apt_repositories=['deb [arch=amd64] http://repo.radeon.com/rocm/apt/debian/ xenial main'],
 -                    ospackages=['ocl-icd-libopencl1', 'ocl-icd-opencl-dev', 'opencl-headers', 'libelf1', 'rocm-opencl'])
 -    else:
 -        return None
 -
 -
  def get_clfft(args):
      if (args.clfft is not None):
          return hpccm.building_blocks.generic_cmake(
@@@ -319,53 -300,10 +317,53 @@@ def add_tsan_compiler_build_stage(input
                       'ln -s /usr/local/libexec/c++-analyzer /usr/local/bin/c++-analyzer-' + str(input_args.llvm)])
      output_stages['compiler_build'] = tsan_stage
  
 +def oneapi_runtime(_from='0'):
 +    oneapi_runtime_stage = hpccm.Stage()
 +    oneapi_runtime_stage += hpccm.primitives.copy(_from='oneapi-build',
 +                                                  files={"/opt/intel": "/opt/intel",
 +                                                         "/etc/bash.bashrc": "/etc/bash.bashrc"})
 +    return oneapi_runtime_stage
 +
 +def add_oneapi_compiler_build_stage(input_args, output_stages: typing.Mapping[str, hpccm.Stage]):
 +    """Isolate the oneAPI preparation stage.
 +
 +    This stage is isolated so that its installed components are minimized in the
 +    final image (chiefly /opt/intel) and its environment setup script can be
 +    sourced. This also helps with rebuild time and final image size.
 +
 +    Note that the ICC compiler inside oneAPI on linux also needs
 +    gcc to build other components and provide libstdc++.
 +    """
 +    if not isinstance(output_stages, collections.abc.MutableMapping):
 +        raise RuntimeError('Need output_stages container.')
 +    oneapi_stage = hpccm.Stage()
 +    oneapi_stage += hpccm.primitives.baseimage(image=base_image_tag(input_args), _as='oneapi-build')
 +
 +    version = str(input_args.oneapi)
 +
 +    # Add required components for the next stage (both for hpccm and Intel's setvars.sh script)
 +    oneapi_stage += hpccm.building_blocks.packages(ospackages=['wget', 'gnupg2', 'ca-certificates', 'lsb-release'])
 +    oneapi_stage += hpccm.building_blocks.packages(
 +        apt_keys=['https://apt.repos.intel.com/intel-gpg-keys/GPG-PUB-KEY-INTEL-SW-PRODUCTS-2023.PUB'],
 +        apt_repositories=['deb https://apt.repos.intel.com/oneapi all main'],
 +        # Add minimal packages (not the whole HPC toolkit!)
 +        ospackages=['intel-oneapi-dpcpp-cpp-compiler-pro-{}'.format(version),
 +            'intel-oneapi-openmp-{}'.format(version),
 +            'intel-oneapi-mkl-{}'.format(version),
 +            'intel-oneapi-mkl-devel-{}'.format(version)]
 +    )
 +    # Ensure that all bash shells on the final container will have access to oneAPI
 +    oneapi_stage += hpccm.primitives.shell(
 +            commands=['echo "source /opt/intel/oneapi/setvars.sh" >> /etc/bash.bashrc']
 +            )
 +    setattr(oneapi_stage, 'runtime', oneapi_runtime)
 +
 +    output_stages['compiler_build'] = oneapi_stage
 +
  def prepare_venv(version: StrictVersion) -> typing.Sequence[str]:
      """Get shell commands to set up the venv for the requested Python version."""
      major = version.version[0]
 -    minor = version.version[1]
 +    minor = version.version[1]  # type: int
  
      pyenv = '$HOME/.pyenv/bin/pyenv'
  
      # WARNING: Please keep this list synchronized with python_packaging/requirements-test.txt
      # TODO: Get requirements.txt from an input argument.
      commands.append("""{path}/bin/python -m pip install --upgrade \
 -            'cmake>=3.9.6' \
 +            'cmake>=3.13' \
              'flake8>=3.7.7' \
 -            'mpi4py>=2' \
 +            'mpi4py>=3.0.3' \
              'networkx>=2.0' \
              'numpy>=1' \
              'pip>=10.1' \
              'pytest>=3.9' \
 -            'setuptools>=28.0.0' \
 -            'scikit-build>=0.7'""".format(path=venv_path))
 +            'setuptools>=42' \
 +            'scikit-build>=0.10'""".format(path=venv_path))
 +
 +    # TODO: Remove 'importlib_resources' dependency when Python >=3.7 is required.
 +    if minor == 6:
 +        commands.append("""{path}/bin/python -m pip install --upgrade \
 +                'importlib_resources'""".format(path=venv_path))
  
      return commands
  
@@@ -472,47 -405,6 +470,47 @@@ def add_python_stages(building_blocks: 
      output_stages['pyenv'] = pyenv_stage
  
  
 +def add_documentation_dependencies(input_args,
 +                                   output_stages: typing.MutableMapping[str, hpccm.Stage]):
 +    """Add appropriate layers according to doxygen input arguments."""
 +    if input_args.doxygen is None:
 +        return
 +    output_stages['main'] += hpccm.primitives.shell(
 +        commands=['sed -i \'/\"XPS\"/d;/\"PDF\"/d;/\"PS\"/d;/\"EPS\"/d;/disable ghostscript format types/d\' /etc/ImageMagick-6/policy.xml'])
 +    output_stages['main'] += hpccm.building_blocks.pip(pip='pip3', packages=['sphinx==1.6.1', 'gcovr'])
 +    if input_args.doxygen == '1.8.5':
 +        doxygen_commit = 'ed4ed873ab0e7f15116e2052119a6729d4589f7a'
 +        output_stages['main'] += hpccm.building_blocks.generic_autotools(
 +            repository='https://github.com/westes/flex.git',
 +            commit='f7788a9a0ecccdc953ed12043ccb59ca25714018',
 +            prefix='/tmp/install-of-flex',
 +            configure_opts=['--disable-shared'],
 +            preconfigure=['./autogen.sh'])
 +        output_stages['main'] += hpccm.building_blocks.generic_autotools(
 +            repository='https://github.com/doxygen/doxygen.git',
 +            commit=doxygen_commit,
 +            prefix='',
 +            configure_opts=[
 +                '--flex /tmp/install-of-flex/bin/flex',
 +                '--static'])
 +    else:
 +        version = input_args.doxygen
 +        archive_name = 'doxygen-{}.linux.bin.tar.gz'.format(version)
 +        archive_url = 'https://sourceforge.net/projects/doxygen/files/rel-{}/{}'.format(
 +            version,
 +            archive_name
 +        )
 +        binary_path = 'doxygen-{}/bin/doxygen'.format(version)
 +        commands = [
 +            'mkdir doxygen && cd doxygen',
 +            'wget {}'.format(archive_url),
 +            'tar xf {} {}'.format(archive_name, binary_path),
 +            'cp {} /usr/local/bin/'.format(binary_path),
 +            'cd .. && rm -rf doxygen'
 +        ]
 +        output_stages['main'] += hpccm.primitives.shell(commands=commands)
 +
 +
  def build_stages(args) -> typing.Iterable[hpccm.Stage]:
      """Define and sequence the stages for the recipe corresponding to *args*."""
  
      # installed.
      if args.llvm is not None and args.tsan is not None:
          add_tsan_compiler_build_stage(input_args=args, output_stages=stages)
 +    if args.oneapi is not None:
 +        add_oneapi_compiler_build_stage(input_args=args, output_stages=stages)
  
      # Building blocks are chunks of container-builder instructions that can be
      # copied to any build stage with the addition operator.
      building_blocks = collections.OrderedDict()
 +    building_blocks['base_packages'] = hpccm.building_blocks.packages(
 +        ospackages=_common_packages)
  
      # These are the most expensive and most reusable layers, so we put them first.
      building_blocks['compiler'] = get_compiler(args, compiler_build_stage=stages.get('compiler_build'))
      building_blocks['mpi'] = get_mpi(args, building_blocks['compiler'])
 +    for i, cmake in enumerate(args.cmake):
 +        building_blocks['cmake' + str(i)] = hpccm.building_blocks.cmake(
 +            eula=True,
 +            prefix='/usr/local/cmake-{}'.format(cmake),
 +            version=cmake)
  
      # Install additional packages early in the build to optimize Docker build layer cache.
 -    os_packages = _common_packages + get_llvm_packages(args)
 +    os_packages = list(get_llvm_packages(args)) + get_opencl_packages(args)
      if args.doxygen is not None:
          os_packages += _docs_extra_packages
 -    building_blocks['ospackages'] = hpccm.building_blocks.packages(ospackages=os_packages)
 +    if args.oneapi is not None:
 +        os_packages += ['lsb-release']
 +    building_blocks['extra_packages'] = hpccm.building_blocks.packages(
 +        ospackages=os_packages,
 +        apt_ppas=['ppa:intel-opencl/intel-opencl'],
 +        apt_keys=['http://repo.radeon.com/rocm/apt/debian/rocm.gpg.key'],
 +        apt_repositories=['deb [arch=amd64] http://repo.radeon.com/rocm/apt/debian/ xenial main']
 +    )
  
 -    building_blocks['cmake'] = hpccm.building_blocks.cmake(eula=True, version=args.cmake)
 -    building_blocks['opencl'] = get_opencl(args)
      building_blocks['clfft'] = get_clfft(args)
  
      # Add Python environments to MPI images, only, so we don't have to worry
                                                  packages=['pytest', 'networkx', 'numpy'])
  
      # Add documentation requirements (doxygen and sphinx + misc).
 -    if (args.doxygen is not None):
 -        if (args.doxygen == '1.8.5'):
 -            doxygen_commit = 'ed4ed873ab0e7f15116e2052119a6729d4589f7a'
 -        else:
 -            doxygen_commit = 'a6d4f4df45febe588c38de37641513fd576b998f'
 -        stages['main'] += hpccm.building_blocks.generic_autotools(
 -            repository='https://github.com/westes/flex.git',
 -            commit='f7788a9a0ecccdc953ed12043ccb59ca25714018',
 -            prefix='/tmp/install-of-flex',
 -            configure_opts=['--disable-shared'],
 -            preconfigure=['./autogen.sh'])
 -        stages['main'] += hpccm.building_blocks.generic_autotools(
 -            repository='https://github.com/doxygen/doxygen.git',
 -            commit=doxygen_commit,
 -            prefix='',
 -            configure_opts=[
 -                '--flex /tmp/install-of-flex/bin/flex',
 -                '--static'],
 -            postinstall=[
 -                'sed -i \'/\"XPS\"/d;/\"PDF\"/d;/\"PS\"/d;/\"EPS\"/d;/disable ghostscript format types/d\' /etc/ImageMagick-6/policy.xml'])
 -        stages['main'] += hpccm.building_blocks.pip(pip='pip3', packages=['sphinx==1.6.1'])
 +    if args.doxygen is not None:
 +        add_documentation_dependencies(args, stages)
  
      if 'pyenv' in stages and stages['pyenv'] is not None:
          stages['main'] += hpccm.primitives.copy(_from='pyenv', _mkdir=True, src=['/root/.pyenv/'],
          # stages['main'] += hpccm.primitives.copy(_from='pyenv', src=['/root/.bashrc'],
          #                                         dest='/root/')
  
 +    # Make sure that `python` resolves to something.
 +    stages['main'] += hpccm.primitives.shell(commands=['test -x /usr/bin/python || '
 +                                                       'update-alternatives --install /usr/bin/python python /usr/bin/python3 1 && '
 +                                                       '/usr/bin/python --version'])
 +
      # Note that the list of stages should be sorted in dependency order.
      for build_stage in stages.values():
          if build_stage is not None:
diff --combined docs/conf.cmakein.py
index 5ca284b39426b725831f8f9851511c45b05ef6f6,0000000000000000000000000000000000000000..c5993c6a09e4f8f864fd7c47c93499f2c61fd25c
mode 100644,000000..100644
--- /dev/null
@@@ -1,434 -1,0 +1,429 @@@
- .. _continuous integration server used by GROMACS: http://jenkins.gromacs.org
- .. _Jenkins: http://jenkins-ci.org
 +#
 +# This file is part of the GROMACS molecular simulation package.
 +#
 +# Copyright (c) 2015,2016,2017,2018,2019 by the GROMACS development team.
 +# Copyright (c) 2020, by the GROMACS development team, led by
 +# Mark Abraham, David van der Spoel, Berk Hess, and Erik Lindahl,
 +# and including many others, as listed in the AUTHORS file in the
 +# top-level source directory and at http://www.gromacs.org.
 +#
 +# GROMACS is free software; you can redistribute it and/or
 +# modify it under the terms of the GNU Lesser General Public License
 +# as published by the Free Software Foundation; either version 2.1
 +# of the License, or (at your option) any later version.
 +#
 +# GROMACS is distributed in the hope that it will be useful,
 +# but WITHOUT ANY WARRANTY; without even the implied warranty of
 +# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
 +# Lesser General Public License for more details.
 +#
 +# You should have received a copy of the GNU Lesser General Public
 +# License along with GROMACS; if not, see
 +# http://www.gnu.org/licenses, or write to the Free Software Foundation,
 +# Inc., 51 Franklin Street, Fifth Floor, Boston, MA  02110-1301  USA.
 +#
 +# If you want to redistribute modifications to GROMACS, please
 +# consider that scientific software is very special. Version
 +# control is crucial - bugs must be traceable. We will be happy to
 +# consider code for inclusion in the official distribution, but
 +# derived work must not be called official GROMACS. Details are found
 +# in the README & COPYING files - if they are missing, get the
 +# official version at http://www.gromacs.org.
 +#
 +# To help us fund GROMACS development, we humbly ask that you cite
 +# the research papers on the package. Check out http://www.gromacs.org.
 +
 +# -*- coding: utf-8 -*-
 +#
 +# GROMACS documentation build configuration file, created by
 +# sphinx-quickstart on Tue Jan 13 14:28:44 2015.
 +#
 +# This file is execfile()d with the current directory set to its
 +# containing dir.
 +#
 +# Note that not all possible configuration values are present in this
 +# autogenerated file.
 +#
 +# All configuration values have a default; values that are commented out
 +# serve to show the default.
 +
 +import datetime
 +import os
 +import sys
 +
 +# The following definitions are completed via CMake machinery.
 +gmx_containers_path = '@GMX_ADMIN_DIR@/containers'
 +gmx_sphinx_extension_path = '@SPHINX_EXTENSION_PATH@'
 +gmxapi_staging_path = '@GMXAPI_PYTHON_STAGING_DIR@'
 +releng_path = '@RELENG_PATH@'
 +gmx_version_string = '@GMX_VERSION_STRING@'
 +gmx_version_string_full = '@GMX_VERSION_STRING_FULL@'
 +regressiontest_version = '@REGRESSIONTEST_VERSION@'
 +gmx_min_sphinx = '@EXPECTED_SPHINX_VERSION@'
 +gmx_image_convert = '@IMAGE_CONVERT_STRING@'
 +variables = [
 +    ('EXPECTED_DOXYGEN_VERSION', '@EXPECTED_DOXYGEN_VERSION@'),
 +    ('EXPECTED_SPHINX_VERSION', '@EXPECTED_SPHINX_VERSION@'),
 +    ('CMAKE_MINIMUM_REQUIRED_VERSION', '@CMAKE_MINIMUM_REQUIRED_VERSION@'),
 +    ('REQUIRED_CUDA_VERSION', '@REQUIRED_CUDA_VERSION@'),
 +    ('REQUIRED_CUDA_COMPUTE_CAPABILITY', '@REQUIRED_CUDA_COMPUTE_CAPABILITY@'),
 +    ('REQUIRED_OPENCL_MIN_VERSION', '@REQUIRED_OPENCL_MIN_VERSION@'),
 +    ('SOURCE_MD5SUM', '@SOURCE_MD5SUM@'),
 +    ('REGRESSIONTEST_MD5SUM', '@REGRESSIONTEST_MD5SUM_STRING@'),
 +    ('GMX_TNG_MINIMUM_REQUIRED_VERSION', '@GMX_TNG_MINIMUM_REQUIRED_VERSION@'),
 +    ('GMX_LMFIT_REQUIRED_VERSION', '@GMX_LMFIT_REQUIRED_VERSION@'),
 +    ('GMX_MANUAL_DOI_STRING', '@GMX_MANUAL_DOI_STRING@'),
 +    ('GMX_SOURCE_DOI_STRING', '@GMX_SOURCE_DOI_STRING@')
 +]
 +# End of build-time substitutions.
 +
 +sys.path.append(gmx_containers_path)
 +sys.path.append(gmx_sphinx_extension_path)
 +if releng_path and os.path.isdir(releng_path):
 +    sys.path.append(releng_path)
 +if gmxapi_staging_path and os.path.isdir(gmxapi_staging_path):
 +    sys.path.append(gmxapi_staging_path)
 +
 +# If extensions (or modules to document with autodoc) are in another directory,
 +# add these directories to sys.path here. If the directory is relative to the
 +# documentation root, use os.path.abspath to make it absolute, like shown here.
 +#sys.path.insert(0, os.path.abspath('.'))
 +
 +# -- General configuration ------------------------------------------------
 +
 +# If your documentation needs a minimal Sphinx version, state it here.
 +# gmx_min_sphinx is set from the expected minimum version of Sphinx
 +# in CMakeLists.txt
 +needs_sphinx = gmx_min_sphinx
 +
 +# Add any Sphinx extension module names here, as strings. They can be
 +# extensions coming with Sphinx (named 'sphinx.ext.*') or your custom
 +# ones.
 +extensions = [
 +    'sphinx.ext.autodoc',
 +    'sphinx.ext.graphviz',
 +    'sphinx.ext.extlinks',
 +    'sphinx.ext.napoleon',
 +    'sphinx.ext.doctest',
 +    'sphinx.ext.ifconfig',
 +    'sphinx.ext.intersphinx',
 +    'sphinx.ext.mathjax',
 +    'sphinx.ext.napoleon',
 +    'sphinx.ext.todo',
 +    # The viewcode extension requires extra configuration or consideration to
 +    # avoid collecting sources for external packages (through intersphninx) or
 +    # generating broken links for compiled extension modules that do not have
 +    # inspectable 'code' member data.
 +    # Ref: https://www.sphinx-doc.org/en/master/usage/extensions/viewcode.html
 +    # Please run linkchecker and inspect the generated
 +    # docs/html/_modules/index.html page before committing a change that enables
 +    # 'sphinx.ext.viewcode',
 +    'gmxsphinx'
 +]
 +extlinks = {'issue': ('https://gitlab.com/gromacs/gromacs/-/issues/%s',
 +                      'Issue ')}
 +
 +# Add any paths that contain templates here, relative to this directory.
 +# templates_path = ['_templates']
 +
 +# The suffix of source filenames.
 +source_suffix = '.rst'
 +
 +# The encoding of source files.
 +#source_encoding = 'utf-8-sig'
 +
 +# The master toctree document.
 +master_doc = 'index'
 +
 +# General information about the project.
 +project = u'GROMACS'
 +copyright = str(datetime.datetime.now().year) + u', GROMACS development team'
 +thisyear_string = str(datetime.datetime.now().year)
 +# The version info for the project you're documenting, acts as replacement for
 +# |version| and |release|, also used in various other places throughout the
 +# built documents.
 +#
 +# The short X.Y version.
 +version = gmx_version_string
 +# The full version, including alpha/beta/rc tags.
 +release = gmx_version_string_full
 +# default file extension for plots
 +plotext = u'.eps'
 +
 +# The language for content autogenerated by Sphinx. Refer to documentation
 +# for a list of supported languages.
 +#language = None
 +
 +# There are two options for replacing |today|: either, you set today to some
 +# non-false value, then it is used:
 +#today = ''
 +# Else, today_fmt is used as the format for a strftime call.
 +#today_fmt = '%B %d, %Y'
 +
 +# List of patterns, relative to source directory, that match files and
 +# directories to ignore when looking for source files.
 +exclude_patterns = ['fragments']
 +if not tags.has('do_man'):
 +    exclude_patterns += ['man']
 +
 +# Set variable if documentation build can convert images or not
 +# to selectively include files in reference manual
 +def setup(app):
 +    app.add_config_value('gmx_image_convert', 'impossible', 'env')
 +
 +# The reST default role (used for this markup: `text`) to use for all
 +# documents.
 +default_role = 'any'
 +
 +# If true, '()' will be appended to :func: etc. cross-reference text.
 +#add_function_parentheses = True
 +
 +# If true, the current module name will be prepended to all description
 +# unit titles (such as .. function::).
 +#add_module_names = True
 +
 +# If true, sectionauthor and moduleauthor directives will be shown in the
 +# output. They are ignored by default.
 +#show_authors = False
 +
 +# The name of the Pygments (syntax highlighting) style to use.
 +pygments_style = 'sphinx'
 +
 +# A list of ignored prefixes for module index sorting.
 +#modindex_common_prefix = []
 +
 +# If true, keep warnings as "system message" paragraphs in the built documents.
 +#keep_warnings = False
 +
 +# Configure the values for all the variables that might later configure any of the .rst files.
 +substitutions = ['.. |{0}| replace:: {1}'.format(*x) for x in variables if x[1]]
 +substitutions.extend(['.. |{0}| replace:: unknown'.format(x[0]) for x in variables if x[1] == ''])
 +rst_epilog = "\n".join(substitutions)
 +rst_epilog += """
 +.. |Gromacs| replace:: GROMACS
 +.. _gmx-manual: manual-{gmx_version_string}.pdf
 +.. _gmx-manual-parent-dir: ../manual-{gmx_version_string}.pdf
 +.. |gmx-source-package-ftp| replace:: As ftp ftp://ftp.gromacs.org/pub/gromacs/gromacs-{gmx_version_string}.tar.gz
 +.. |gmx-source-package-http| replace:: As http http://ftp.gromacs.org/pub/gromacs/gromacs-{gmx_version_string}.tar.gz
 +.. |gmx-regressiontests-package| replace:: http://ftp.gromacs.org/pub/regressiontests/regressiontests-{regressiontest_version}.tar.gz
 +.. _up-to-date installation instructions: http://manual.gromacs.org/documentation/current/install-guide/index.html
 +.. _CUDA: http://www.nvidia.com/object/cuda_home_new.html
 +.. _OpenCL: https://www.khronos.org/opencl/
 +.. _OpenMPI: http://www.open-mpi.org
 +.. _MPICH: http://www.mpich.org
 +.. _LAM-MPI: http://www.lam-mpi.org
 +.. _OpenMP: http://en.wikipedia.org/wiki/OpenMP
 +.. _CMake installation page: http://www.cmake.org/install/
 +.. _Ubuntu toolchain ppa page: https://launchpad.net/~ubuntu-toolchain-r/+archive/ubuntu/test
 +.. _EPEL page: https://fedoraproject.org/wiki/EPEL
 +.. _running CMake: http://www.cmake.org/runningcmake/
 +.. _CMake environment variables: http://cmake.org/Wiki/CMake_Useful_Variables#Environment_Variables
 +.. _FFTW: http://www.fftw.org
 +.. _FFTW installation guide: http://www.fftw.org/doc/Installation-and-Customization.html#Installation-and-Customization
 +.. _MKL: https://software.intel.com/en-us/intel-mkl
 +.. _VMD: http://www.ks.uiuc.edu/Research/vmd/
 +.. _PyMOL: http://www.pymol.org
- .. _issues: https://gitlab.com/gromacs/gromacs/-/issues/
 +.. _webpage: http://www.gromacs.org
 +.. _ftp site: ftp://ftp.gromacs.org/pub/gromacs/
 +.. _tutorials: http://www.gromacs.org/Documentation/Tutorials
 +.. _issue tracker: https://gitlab.com/gromacs/gromacs/-/issues/
- .. _redmine: http://redmine.gromacs.org
- .. _gerrit: http://gerrit.gromacs.org
 +.. _gitlab: https://gitlab.com/gromacs/gromacs/
 +.. _download: ../download.html
 +.. |thisyear| replace:: {thisyear_string}
 +""".format(gmx_version_string=gmx_version_string, regressiontest_version=regressiontest_version, thisyear_string=thisyear_string)
 +
 +# -- Options for HTML output ----------------------------------------------
 +
 +# The theme to use for HTML and HTML Help pages.  See the documentation for
 +# a list of builtin themes.
 +html_theme = 'classic'
 +
 +# Theme options are theme-specific and customize the look and feel of a theme
 +# further.  For a list of options available for each theme, see the
 +# documentation.
 +#html_theme_options = {}
 +
 +# Add any paths that contain custom themes here, relative to this directory.
 +html_theme_path = ['']
 +
 +# The name for this set of Sphinx documents.  If None, it defaults to
 +# "<project> v<release> documentation".
 +html_title = u'GROMACS ' + release + ' documentation'
 +
 +# A shorter title for the navigation bar.  Default is the same as html_title.
 +html_short_title = u'GROMACS ' + version
 +
 +# The name of an image file (relative to this directory) to place at the top
 +# of the sidebar.
 +#html_logo = None
 +
 +# The name of an image file (within the static path) to use as favicon of the
 +# docs.  This file should be a Windows icon file (.ico) being 16x16 or 32x32
 +# pixels large.
 +#html_favicon = None
 +
 +# Add any paths that contain custom static files (such as style sheets) here,
 +# relative to this directory. They are copied after the builtin static files,
 +# so a file named "default.css" will overwrite the builtin "default.css".
 +html_static_path = ['_static']
 +
 +# Add any extra paths that contain custom files (such as robots.txt or
 +# .htaccess) here, relative to this directory. These files are copied
 +# directly to the root of the documentation.
 +#html_extra_path = []
 +
 +# If not '', a 'Last updated on:' timestamp is inserted at every page bottom,
 +# using the given strftime format.
 +#html_last_updated_fmt = '%b %d, %Y'
 +
 +# If true, SmartyPants will be used to convert quotes and dashes to
 +# typographically correct entities.
 +#html_use_smartypants = True
 +
 +# Custom sidebar templates, maps document names to template names.
 +#html_sidebars = {}
 +
 +# Additional templates that should be rendered to pages, maps page names to
 +# template names.
 +#html_additional_pages = {}
 +
 +# If false, no module index is generated.
 +#html_domain_indices = True
 +
 +# If false, no index is generated.
 +#html_use_index = True
 +
 +# If true, the index is split into individual pages for each letter.
 +#html_split_index = False
 +
 +# If true, links to the reST sources are added to the pages.
 +#html_show_sourcelink = True
 +
 +# If true, "Created using Sphinx" is shown in the HTML footer. Default is True.
 +#html_show_sphinx = True
 +
 +# If true, "(C) Copyright ..." is shown in the HTML footer. Default is True.
 +#html_show_copyright = True
 +
 +# If true, an OpenSearch description file will be output, and all pages will
 +# contain a <link> tag referring to it.  The value of this option must be the
 +# base URL from which the finished HTML is served.
 +#html_use_opensearch = ''
 +
 +# This is the file name suffix for HTML files (e.g. ".xhtml").
 +#html_file_suffix = None
 +
 +# Output file base name for HTML help builder.
 +htmlhelp_basename = 'Gromacsdoc'
 +
 +
 +# -- Options for LaTeX output ---------------------------------------------
 +
 +latex_elements = {
 +# The paper size ('letterpaper' or 'a4paper').
 +'papersize': 'a4paper',
 +
 +
 +# The font size ('10pt', '11pt' or '12pt').
 +'pointsize': '11',
 +
 +# Additional stuff for the LaTeX preamble.
 +# The tocdepth setting is needed to overwrite the default value set by Sphinx
 +# to get a more detailed toctree in the pdf version of the manual.
 +    'preamble': r'''
 +    \usepackage{here}
 +    \usepackage{picins}
 +    \usepackage{underscore}
 +    \usepackage{tabularx}
 +    \usepackage{multicol}
 +    \usepackage{dcolumn}
 +    \usepackage{makeidx}
 +    \usepackage{times}
 +    \usepackage{ifthen}
 +    \usepackage{enumitem}
 +    \usepackage{longtable}
 +    \usepackage{pdflscape}
 +    \pagenumbering{roman}
 +    \usepackage{array}
 +    \setcounter{tocdepth}{2}
 +    ''',
 +# Format output to avoid empty pages
 +  'classoptions': ',openany,oneside'
 +}
 +
 +# Grouping the document tree into LaTeX files. List of tuples
 +# (source start file, target name, title,
 +#  author, documentclass [howto, manual, or own class]).
 +latex_documents = [
 +  ('index', 'gromacs.tex', u'GROMACS Documentation',
 +   u'GROMACS development team', 'manual'),
 +]
 +
 +# The name of an image file (relative to this directory) to place at the top of
 +# the title page.
 +latex_logo = 'reference-manual/plots/peregrine.png'
 +
 +# For "manual" documents, if this is true, then toplevel headings are parts,
 +# not chapters.
 +latex_use_parts = True
 +
 +# If true, show page references after internal links.
 +latex_show_pagerefs = True
 +
 +# If true, show URL addresses after external links.
 +#latex_show_urls = False
 +
 +# Documents to append as an appendix to all manuals.
 +#latex_appendices = []
 +
 +# If false, no module index is generated.
 +#latex_domain_indices = True
 +
 +
 +# -- Options for manual page output ---------------------------------------
 +
 +# One entry per manual page. List of tuples
 +# (source start file, name, description, authors, manual section).
 +if tags.has('do_man'):
 +    exec(open('conf-man.py').read())
 +
 +# If true, show URL addresses after external links.
 +#man_show_urls = False
 +
 +
 +# -- Options for Texinfo output -------------------------------------------
 +
 +# Grouping the document tree into Texinfo files. List of tuples
 +# (source start file, target name, title, author,
 +#  dir menu entry, description, category)
 +texinfo_documents = [
 +  ('index', 'GROMACS', u'GROMACS Documentation',
 +   u'GROMACS development team', 'GROMACS', 'One line description of project.',
 +   'Miscellaneous'),
 +]
 +
 +# Documents to append as an appendix to all manuals.
 +#texinfo_appendices = []
 +
 +# If false, no module index is generated.
 +#texinfo_domain_indices = True
 +
 +# How to display URL addresses: 'footnote', 'no', or 'inline'.
 +#texinfo_show_urls = 'footnote'
 +
 +# If true, do not generate a @detailmenu in the "Top" node's menu.
 +#texinfo_no_detailmenu = False
 +
 +# Make it possible to use numbered labels for figures and tables
 +numfig = True
 +
 +# -- Options for autodoc extension ----------------------------------------
 +
 +# http://www.sphinx-doc.org/en/master/usage/extensions/autodoc.html#confval-autodoc-mock-imports
 +autodoc_mock_imports = ['hpccm', 'hpccm.config']
 +
 +# -- Options for intersphinx extension ------------------------------------
 +
 +intersphinx_mapping = {'python': ('https://docs.python.org/3', None)}
 +intersphinx_cache_limit = -1
 +intersphinx_timeout = 10
index 22db6df9f58f0e51a7aa78316cf2282281891be5,53db811dd7625e76ad0d030a627f6e986ec6f441..dfc8862a531117e5b1bce8f00ae4cc6eecc791af
@@@ -40,7 -40,8 +40,8 @@@ please make sure that you have checked 
    widespread adoption of the method.
  
  * *Advance discussion*: Please communicate with the other developers,
-   e.g.  on the `developer mailing list`_ mailing list, or `issue tracker`_ to let them know of the general
 -  e.g.  on the `developer mailing list`_ mailing list, or `GitLab`_
 -  to let them know of the general
++  e.g.  on the `developer mailing list`_ mailing list, or
++  `issue tracker`_ to let them know of the general
    nature of your plans. This will prevent duplicate or wasted
    effort. It is also a good idea to search those resources as well as
    the literature and WWW for other projects that may be relevant.
  Preparing code for submission
  -----------------------------
  
 -|Gromacs| revision control uses a ``git`` repository managed by :ref:`Gerrit <gmx-gerrit>`.
 +|Gromacs| uses ``git`` for :doc:`change-management`.
  Instead of accepting "pull requests", |Gromacs| changes are submitted as individual
 -commits on the tip of the ``master`` branch hosted at https://gerrit.gromacs.org.
 +commits on the tip of the ``master`` branch hosted at `gitlab`_.
  Preparing, submitting, and managing patches for a change requires a little bit
  of set-up. Refer to :doc:`change-management` for information about
  
 -* accessing the |Gromacs| Gerrit server
 +* accessing the |Gromacs| *git* repository
  * structure of the repository
  * source control without merge commits
  * ``git`` usage that may be less common in other development work flows
index 74272b9ab7414af2538cf3ebf7d88a1b9bd243a2,8b7e8065bd633219ec4a02a7b95123aec73dae6a..5b9ccdade4d23b1b6c83e7810dc6d10cbbfeaf53
@@@ -3,23 -3,25 +3,25 @@@
  
  .. _releng-workflow-release:
  
- Release engineering with Gitlab
+ Release engineering with GitLab
  ===============================
  
  .. toctree::
     :hidden:
  
- We are currently switching our build and testing system to use Gitlab
- CI pipelines run on GitLab Runner. This section is going to be extended
- with individual build information as it comes available. For now we are
- using a combination of building with the previous system on Jenkins
- and post-submit verification on Gitlab.
+ We are currently switching our build and testing system to use GitLab
+ and the integrated CI system, with information for the general system found
+ in the official `GitLab documentation <https://docs.gitlab.com/ee/ci/yaml/>`_.
+ The new configuration for the builds and tests can be found in the file
+ ``.gitlab-ci.yml``, with the templates for configuring is found in the files in the
+ ``admin/ci-templates/`` directory. This section is going to be extended
 -with individual build information as it comes available. For now we are
 -using a combination of building with the previous system on Jenkins
 -and post-submit verification on GitLab.
++with individual build information as it comes available.
 +
 +.. seealso:: :doc:`../infrastructure`
  
  .. _releng-triggering-builds:
  
- Triggering builds on Gitlab
+ Triggering builds on GitLab
  ---------------------------
  
  Pipelines can be triggered through the web interface, with different
index 93e58492fe599fd6a66100fd257da1a7315b9ede,2d9109b99347bebb9e24f7f35d3b06c9154cc41c..07f11f72041a82fbfd9e5182bc66d9cb6f6f0818
@@@ -1,9 -1,11 +1,9 @@@
 -Guidelines for creating meaningful redmine issue reports
 -========================================================
 +Guidelines for creating meaningful issue reports
 +================================================
  This section gives some started on how to generate useful issues on the
 -|Gromacs| `redmine issue tracker`_. The information here comes to a large extent
 +|Gromacs| `issue tracker`_. The information here comes to a large extent
  directly from there, to help you in preparing your reports.
  
 -.. _redmine issue tracker: https://redmine.gromacs.org
 -
  What to report
  ^^^^^^^^^^^^^^
  Please only report issues you have confirmed to be caused by |Gromacs| behaving in an
@@@ -34,8 -36,8 +34,8 @@@ does not result in any warnings or erro
  not be considered as *real*, or at the minimum it will be much harder to analyse to find the actual issue.
  
  
- If your inputs are sensitive, then it is possible to create private `issues`_ so that the
 -If your inputs are sensitive, then it is possible to create private Redmine issues so that the
--developer team can have access to solve the problem, while preventing widespread
++If your inputs are sensitive, then it is possible to create private `issues <issue tracker>`_
++so that the developer team can have access to solve the problem, while preventing widespread
  visibility on the internet.
  
  
@@@ -59,7 -61,7 +59,7 @@@ General issue workflo
  The general issue workflow is shown in the figure below:
  
  .. image:: redmine-states.png
 -   :alt:  Sample procedure pathway for issues reported in redmine.
 +   :alt:  Sample procedure pathway for reported issues.
  
  
  .. Text below is stolen from the old Gromacs web page
index fd90412c285334743a55380d9308631d3bc7deae,a4037096c1e2d602edaca0f962f1694a4cb743a4..c6f7d6f3c58a00f734453ab0dc8fb536a0d0c2c0
@@@ -8,28 -8,24 +8,28 @@@ Several tools have their own individua
  
     doxygen
     change-management
 -   jenkins
 +   infrastructure
     releng/index
     gmxtree
     code-formatting
     testutils
     physical_validation
  
 -.. TODO: Consider what is the most reasonable structure; currently, this list
 +.. todo:: :issue:`3032`
 +
 +   Consider what is the most reasonable structure; currently, this list
     here does not make much sense in the overall organization and creates a
     confusing TOC for the developer guide.
  
 -.. TODO: Add details for most of the tools, either in the form of links to wiki,
 +.. todo:: :issue:`3267`
 +
 +   Add details for most of the tools, either in the form of links to wiki,
     or to a separate page that explains more details.
  
  Change management
  -----------------
  
- |Gromacs| change management is supported by the following tools.
+ |Gromacs| change management uses git and `GitLab`_ for code uploading and testing as well as issues tracking.
  (For change submission guidelines, refer to :doc:`contribute`.)
  
  git
    Other basic tutorial material for ``git`` can be found on the `web <https://git-scm.com/doc/ext>`__.
  
  GitLab
 -  All code changes go through a code review system at
 +  Bugs and issues, as well as some random features and discussions,
 +  are tracked, and all code changes go through a code review system at
    https://gitlab.com/gromacs/gromacs.
  
 -Continuous Integration
 +Build testing
    All changes pushed to GitLab are automatically compiled and otherwise
 -  checked on various platforms using a continuous integration system.
 +  checked on various platforms.
 +  :doc:`infrastructure` documents how builds are automated,
 +  providing information on how to replicate the builds (e.g., to
 +  diagnose issues).
    :doc:`releng/index` provides more information on the technical implementation
    of the builds.
  
 -Issue tracking
 -  Bugs and issues, as well as some random features and discussions,
 -  are tracked at https://gitlab.com/gromacs/gromacs/issues.
 -
  .. _Git Tips & Tricks: http://www.gromacs.org/index.php?title=Developer_Zone/Git/Git_Tips_%26_Tricks
  
  Build system
  ------------
  
 -.. TODO: details, ASAN, others?
 +.. todo:: details, ASAN, others?
  
  CMake
    Main tool used in the build system.
@@@ -86,46 -82,49 +86,46 @@@ The tools and scripts listed below are 
  formatting that follows |Gromacs| style guidelines described on a separate page:
  :doc:`style`.
  
 -uncrustify
 -  `uncrustify <http://uncrustify.sourceforge.net>`_ is used for automatic
 -  indentation and other formatting of the source code to follow
 -  :doc:`formatting`.  All code must remain invariant under uncrustify
 -  with the config at ``admin/uncrustify.cfg``.  A patched version of uncrustify is
 -  used.  See :ref:`gmx-uncrustify` for details.
 -
  clang-format
    We use clang-format to enforce a consistent coding style, with the
    settings recorded in ``.clang-format`` in the main tree.
    See :ref:`gmx-clang-format` for details.
  
 +clang-tidy
 +  The source code linter clang-tidy is used to enforce common restrictions to the
 +  code, with the checks collected under ``.clang-tidy`` at the top of the main tree.
 +  See :ref:`gmx-clang-tidy` for details.
 +
  ``admin/copyright.py``
    This Python script adds and formats copyright headers in source files.
    ``copyright.sh`` (see below) uses the script to check/update copyright years on
    changed files automatically.
  
 -``admin/uncrustify.sh``
 -  This ``bash`` script runs uncrustify for all
 -  files that have local changes and checks that they conform to the prescribed
 -  style.  Optionally, the script can also apply changes to make the files
 -  conform. It is included only for historical reasons.
 -  See :doc:`formatting` for details.
 -
  ``admin/copyright.sh``
    This ``bash`` script runs the ``copyright.py`` python script to enforce
    correct copyright information in all files that have local changes
    and checks that they conform to the prescribed
    style.  Optionally, the script can also apply changes to make the files
    conform.
 -  This script is automatically run by Jenkins to ensure that all commits adhere
 +  This script is automatically run by the CI to ensure that all commits adhere
    to :doc:`formatting`.  If the copyright job does not succeed, it
    means that this script has something to complain.
    See :doc:`code-formatting` for details.
  
  ``admin/clang-format.sh``
    This script enforces coding style using clang-format.
 -  This script is automatically run by Jenkins to ensure that all commits adhere
 +  This script is automatically run by our CI to ensure that all commits adhere
    to :doc:`formatting`.
  
 +``admin/clang-tidy.sh``
 +  The clang-tidy code correctness restrictions are enforced by this script.
 +  The script is also used by the CI to verify the code, in addition to nightly
 +  compilations using clang-tidy on the whole tree.
 +
  ``admin/git-pre-commit``
    This sample git pre-commit hook can be used if one wants to apply
 -  ``uncrustify.sh`` and ``clang-format.sh`` automatically before every commit to check for formatting
 +  ``clang-tidy.sh``, ``copyright.sh`` and ``clang-format.sh`` automatically
 +  before every commit to check for formatting
    issues.  See :doc:`code-formatting` for details.
  
  ``docs/doxygen/includesorter.py``
@@@ -139,12 -138,12 +139,12 @@@ include directive checke
    applied in the formatting script.  To check for issues, it is instead integrated into
    a ``check-source`` build target.  When this target is built, it also checks for
    include formatting issues.  Internally, it uses the sorter script.  This check
 -  is run in Jenkins as part of the Documentation job.
 +  is run in the CI as part of the Documentation job.
    Details for the checking mechanism are on a separate page (common for several
    checkers): :doc:`gmxtree`.
  
  ``admin/reformat_all.sh``
 -  This ``bash`` script runs uncrustify/clang-format/``copyright.py``/include sorter
 +  This ``bash`` script runs clang-format/``copyright.py``/include sorter
    on all relevant files in the source tree (or in a particular directory).
    The script can also produce the list of files where these scripts are applied,
    for use with other scripts.  See :doc:`code-formatting` for details.
@@@ -157,3 -156,4 +157,3 @@@ git attribute
    checking/formatting to apply.  Custom attributes are used for specifying some
    build system dependencies for easier processing in CMake.
  
 -include-what-you-use
index 69daf6f2266f975e25dcd2fa8f090c9a4a2b4272,6f332dff78f33cc2670539335a7a2f9ed7c18d2b..02c05b950e655234e7e942d34af48afad55b5c5c
@@@ -72,8 -72,8 +72,8 @@@ appropriate value instead of ``xxx`` 
  * ``-DCMAKE_C_COMPILER=xxx`` equal to the name of the C99 `Compiler`_ you wish to use (or the environment variable ``CC``)
  * ``-DCMAKE_CXX_COMPILER=xxx`` equal to the name of the C++98 `compiler`_ you wish to use (or the environment variable ``CXX``)
  * ``-DGMX_MPI=on`` to build using `MPI support`_ (generally good to combine with `building only mdrun`_)
 -* ``-DGMX_GPU=on`` to build using nvcc to run using NVIDIA `CUDA GPU acceleration`_ or an OpenCL_ GPU
 -* ``-DGMX_USE_OPENCL=on`` to build with OpenCL_ support enabled. ``GMX_GPU`` must also be set.
 +* ``-DGMX_GPU=CUDA`` to build with NVIDIA CUDA support enabled.
 +* ``-DGMX_GPU=OpenCL`` to build with OpenCL_ support enabled.
  * ``-DGMX_SIMD=xxx`` to specify the level of `SIMD support`_ of the node on which |Gromacs| will run
  * ``-DGMX_BUILD_MDRUN_ONLY=on`` for `building only mdrun`_, e.g. for compute cluster back-end nodes
  * ``-DGMX_DOUBLE=on`` to build |Gromacs| in double precision (slower, and not normally useful)
@@@ -99,25 -99,25 +99,25 @@@ Platfor
  |Gromacs| can be compiled for many operating systems and
  architectures.  These include any distribution of Linux, Mac OS X or
  Windows, and architectures including x86, AMD64/x86-64, several
 -PowerPC including POWER8, ARM v7, ARM v8, and SPARC VIII.
 +PowerPC including POWER8, ARM v8, and SPARC VIII.
  
  Compiler
  ^^^^^^^^
  
 -|Gromacs| can be compiled on any platform with ANSI C99 and C++14
 +|Gromacs| can be compiled on any platform with ANSI C99 and C++17
  compilers, and their respective standard C/C++ libraries. Good
  performance on an OS and architecture requires choosing a good
  compiler. We recommend gcc, because it is free, widely available and
  frequently provides the best performance.
  
  You should strive to use the most recent version of your
 -compiler. Since we require full C++14 support the minimum supported
 +compiler. Since we require full C++17 support the minimum supported
  compiler versions are
  
 -* GNU (gcc) 5.1
 -* Intel (icc) 17.0.1
 -* LLVM (clang) 3.6
 -* Microsoft (MSVC) 2017
 +* GNU (gcc/libstdc++) 7
 +* Intel (icc) 19.1
 +* LLVM (clang/libc++) 5
 +* Microsoft (MSVC) 2017 15.7
  
  Other compilers may work (Cray, Pathscale, older clang) but do
  not offer competitive performance. We recommend against PGI because
@@@ -131,7 -131,7 +131,7 @@@ You may also need the most recent versi
  components beside the compiler itself (e.g. assembler or linker);
  these are often shipped by your OS distribution's binutils package.
  
 -C++14 support requires adequate support in both the compiler and the
 +C++17 support requires adequate support in both the compiler and the
  C++ library. The gcc and MSVC compilers include their own standard
  libraries and require no further configuration. If your vendor's
  compiler also manages the standard library library via compiler flags,
@@@ -139,12 -139,12 +139,12 @@@ these will be honored. For configuratio
  
  On Linux, both the Intel and clang compiler use the libstdc++ which
  comes with gcc as the default C++ library. For |Gromacs|, we require
 -the compiler to support libstc++ version 5.1 or higher. To select a
 +the compiler to support libstc++ version 7.1 or higher. To select a
  particular libstdc++ library, provide the path to g++ with
  ``-DGMX_GPLUSPLUS_PATH=/path/to/g++``.
  
  On Windows with the Intel compiler, the MSVC standard library is used,
 -and at least MSVC 2017 is required. Load the enviroment variables with
 +and at least MSVC 2017 15.7 is required. Load the enviroment variables with
  vcvarsall.bat.
  
  To build with clang and llvm's libcxx standard library, use
@@@ -237,7 -237,7 +237,7 @@@ library. LAM-MPI_ might work, but sinc
  been deprecated for years, it is not supported.
  
  For example, depending on your actual MPI library, use ``cmake
 --DCMAKE_C_COMPILER=mpicc -DCMAKE_CXX_COMPILER=mpicxx -DGMX_MPI=on``.
 +-DMPI_C_COMPILER=mpicc -DGMX_MPI=on``.
  
  
  CMake
@@@ -301,7 -301,7 +301,7 @@@ should also add ``--enable-avx2`` also
  512-wide AVX, including KNL, add ``--enable-avx512`` also.
  FFTW will create a fat library with codelets for all different instruction sets,
  and pick the fastest supported one at runtime.
 -On ARM architectures with NEON SIMD support and IBM Power8 and later, you
 +On ARM architectures with SIMD support and IBM Power8 and later, you
  definitely want version 3.3.5 or later,
  and to compile it with ``--enable-neon`` and ``--enable-vsx``, respectively, for
  SIMD support. If you are using a Cray, there is a special modified
@@@ -348,9 -348,7 +348,9 @@@ Other optional build component
  ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
  
  * Run-time detection of hardware capabilities can be improved by
 -  linking with hwloc, which is automatically enabled if detected.
 +  linking with hwloc. By default this is turned off since it might
 +  not be supported everywhere, but if you have hwloc installed it
 +  should work by just setting ``-DGMX_HWLOC=ON``
  * Hardware-optimized BLAS and LAPACK libraries are useful
    for a few of the |Gromacs| utilities focused on normal modes and
    matrix manipulation, but they do not provide any benefits for normal
    ``-DGMX_USE_LMFIT=none``.
  * zlib is used by TNG for compressing some kinds of trajectory data
  * Building the |Gromacs| documentation is optional, and requires
 -  ImageMagick, pdflatex, bibtex, doxygen, python 3.5, sphinx
 +  ImageMagick, pdflatex, bibtex, doxygen, python 3.6, sphinx
    |EXPECTED_SPHINX_VERSION|, and pygments.
  * The |Gromacs| utility programs often write data files in formats
    suitable for the Grace plotting tool, but it is straightforward to
@@@ -502,7 -500,7 +502,7 @@@ For example, the following command lin
  
  ::
  
 -    cmake .. -DGMX_GPU=ON -DGMX_MPI=ON -DCMAKE_INSTALL_PREFIX=/home/marydoe/programs
 +    cmake .. -DGMX_GPU=CUDA -DGMX_MPI=ON -DCMAKE_INSTALL_PREFIX=/home/marydoe/programs
  
  can be used to build with CUDA GPUs, MPI and install in a custom
  location. You can even save that in a shell script to make it even
@@@ -575,10 -573,6 +575,10 @@@ lead to performance loss, e.g. on Inte
  12. ``IBM_VSX`` Power7, Power8, Power9 and later have this.
  13. ``ARM_NEON`` 32-bit ARMv7 with NEON support.
  14. ``ARM_NEON_ASIMD`` 64-bit ARMv8 and later.
 +15. ``ARM_SVE`` 64-bit ARMv8 and later with the Scalable Vector Extensions (SVE).
 +    The SVE vector length is fixed at CMake configure time. The default vector
 +    length is automatically detected, and this can be changed via the
 +    ``GMX_SIMD_ARM_SVE_LENGTH`` CMake variable.
  
  The CMake configure system will check that the compiler you have
  chosen can target the architecture you have chosen. mdrun will check
@@@ -654,7 -648,7 +654,7 @@@ If you have the CUDA_ Toolkit installed
  
  ::
  
 -    cmake .. -DGMX_GPU=ON -DCUDA_TOOLKIT_ROOT_DIR=/usr/local/cuda
 +    cmake .. -DGMX_GPU=CUDA -DCUDA_TOOLKIT_ROOT_DIR=/usr/local/cuda
  
  (or whichever path has your installation). In some cases, you might
  need to specify manually which of your C++ compilers should be used,
@@@ -675,7 -669,7 +675,7 @@@ manual
  
  The GPU acceleration has been tested on AMD64/x86-64 platforms with
  Linux, Mac OS X and Windows operating systems, but Linux is the
 -best-tested and supported of these. Linux running on POWER 8, ARM v7 and v8
 +best-tested and supported of these. Linux running on POWER 8 and ARM v8
  CPUs also works well.
  
  Experimental support is available for compiling CUDA code, both for host and
@@@ -729,7 -723,7 +729,7 @@@ To trigger an OpenCL_ build the followi
  
  ::
  
 -    cmake .. -DGMX_GPU=ON -DGMX_USE_OPENCL=ON
 +    cmake .. -DGMX_GPU=OpenCL
  
  To build with support for Intel integrated GPUs, it is required
  to add ``-DGMX_OPENCL_NB_CLUSTER_SIZE=4`` to the cmake command line,
@@@ -747,7 -741,7 +747,7 @@@ external library, us
  
  ::
  
 -    cmake .. -DGMX_GPU=ON -DGMX_USE_OPENCL=ON -DclFFT_ROOT_DIR=/path/to/your/clFFT -DGMX_EXTERNAL_CLFFT=TRUE
 +    cmake .. -DGMX_GPU=OpenCL -DclFFT_ROOT_DIR=/path/to/your/clFFT -DGMX_EXTERNAL_CLFFT=TRUE
  
  Static linking
  ~~~~~~~~~~~~~~
@@@ -808,9 -802,8 +808,9 @@@ of the build host machine or otherwise 
  configuration.
  
  Often it is possible to ensure portability by choosing the least
 -common denominator of SIMD support, e.g. SSE2 for x86, and ensuring
 -the you use ``cmake -DGMX_USE_RDTSCP=off`` if any of the target CPU
 +common denominator of SIMD support, e.g. SSE2 for x86. In rare cases
 +of very old x86 machines, ensure that
 +you use ``cmake -DGMX_USE_RDTSCP=off`` if any of the target CPU
  architectures does not support the ``RDTSCP`` instruction.  However, we
  discourage attempts to use a single |Gromacs| installation when the
  execution environment is heterogeneous, such as a mix of AVX and
@@@ -1054,15 -1047,11 +1054,15 @@@ are individual failed tests it could b
  that a tolerance is just a tiny bit too tight. Check the output files
  the script directs you too, and try a different or newer compiler if
  the errors appear to be real. If you cannot get it to pass the
 -regression tests, you might try dropping a line to the gmx-users
 -mailing list, but then you should include a detailed description of
 +regression tests, you might try dropping a line to the
 +`|Gromacs| users forum <https://gromacs.bioexcel.eu/c/gromacs-user-forum>`__,
 +but then you should include a detailed description of
  your hardware, and the output of ``gmx mdrun -version`` (which contains
  valuable diagnostic information in the header).
  
 +Testing for MDRUN_ONLY executables
 +~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
 +
  A build with ``-DGMX_BUILD_MDRUN_ONLY`` cannot be tested with
  ``make check`` from the build tree, because most of the tests
  require a full build to run things like ``grompp``. To test such an
@@@ -1076,24 -1065,19 +1076,24 @@@ directory
  
      mkdir build-normal
      cd build-normal
 -    cmake .. -DCMAKE_INSTALL_PREFIX=/your/installation/prefix/here
 +    # First, build and install normally to allow full testing of the standalone simulator.
 +    cmake .. -DGMX_MPI=ON -DCMAKE_INSTALL_PREFIX=/your/installation/prefix/here
      make -j 4
      make install
      cd ..
      mkdir build-mdrun-only
      cd build-mdrun-only
 -    cmake .. -DGMX_MPI=ON -DGMX_GPU=ON -DGMX_BUILD_MDRUN_ONLY=ON -DCMAKE_INSTALL_PREFIX=/your/installation/prefix/here
 +    # Next, build and install the GMX_BUILD_MDRUN_ONLY version (optional).
 +    cmake .. -DGMX_MPI=ON -DGMX_BUILD_MDRUN_ONLY=ON -DCMAKE_INSTALL_PREFIX=/your/installation/prefix/here
      make -j 4
      make install
      cd /to/your/unpacked/regressiontests
      source /your/installation/prefix/here/bin/GMXRC
      ./gmxtest.pl all -np 2
  
 +Non-standard suffix
 +~~~~~~~~~~~~~~~~~~~
 +
  If your mdrun program has been suffixed in a non-standard way, then
  the ``./gmxtest.pl -mdrun`` option will let you specify that name to the
  test machinery. You can use ``./gmxtest.pl -double`` to test the
@@@ -1102,9 -1086,6 +1102,9 @@@ to stop the test harness attempting to 
  be run. You can use ``./gmxtest.pl -mpirun srun`` if your command to
  run an MPI program is called ``srun``.
  
 +Running MPI-enabled tests
 +~~~~~~~~~~~~~~~~~~~~~~~~~
 +
  The ``make check`` target also runs integration-style tests that may run
  with MPI if ``GMX_MPI=ON`` was set. To make these work with various possible
  MPI libraries, you may need to
@@@ -1311,17 -1292,17 +1311,17 @@@ much everywhere, it is important that w
  it works because we have tested it.
  Every commit in our git source code repository
  is currently tested with a range of configuration options on x86 with
 -gcc versions 6 and 7,
 -clang versions 3.6 and 8,
 +gcc versions 7 and 8,
 +clang versions 8 and 9,
  and
 -For this testing, we use Ubuntu 16.04 or 18.04 operating system.
 +a beta version of oneAPI containing Intel's compiler.
 +For this testing, we use Ubuntu 18.04 or 20.04 operating system.
  Other compiler, library, and OS versions are tested less frequently.
--For details, you can
--have a look at the `continuous integration server used by GROMACS`_,
++For details, you can have a look at the
++`continuous integration server used by GROMACS <gitlab>`_,
  which uses GitLab runner on a local k8s x86 cluster with NVIDIA and
  AMD GPU support.
  
 -We test irregularly on ARM v7, ARM v8, Cray, Fujitsu
 -PRIMEHPC, Power8, Power9,
 +We test irregularly on ARM v8, Cray, Power8, Power9,
  Google Native Client and other environments, and
  with other compilers and compiler versions, too.
index 1fc4bcc6e126d92086585ef35f95670c46f27c9c,795b11ad1fa2c7b9590719ac455159db17250f68..b675c224318eed4524e32606d0cb1038dc1f8a60
@@@ -38,8 -38,8 +38,8 @@@ which in some cases might mean the info
  Comments on form and content are welcome, please send them to one of
  the mailing lists (see our `webpage`_ or this section on
  how to :ref:`contribute <gmx-contribute>`), or open an issue
 -on `GitLab`_. Corrections can also be made in the |Gromacs| git
 +on our `issue tracker`_. Corrections can also be made in the |Gromacs| git
- source repository and uploaded to the |Gromacs| `gerrit`_.
+ source repository and uploaded to the |Gromacs| `GitLab`_.
  
  We release an updated version of the manual whenever
  we release a new version of the software, so in general 
@@@ -49,7 -49,7 +49,7 @@@ minor release number as your |Gromacs| 
  Citation information
  --------------------
  
 -.. TODO needs link to ref list
 +.. todo:: needs link to ref list
  
  |GMX_MANUAL_DOI_STRING|
  
index d79603c6708459dc16da65e7cecbc1d851ca7399,64e252520968ddb67b7d2e1ba48a4a5ed4885351..62108c17874a3bd41a767f8c22c534deca2fb44f
@@@ -160,7 -160,7 +160,7 @@@ These document fixes for issues that we
  introduced into the release-2016 branch since it diverged from
  release-5-1. These will not appear in the final release notes, because
  no formal release is thought to have had the problem. Of course, the
- tracked `issues`_ remain available should further discussion arise.
 -Redmine issues remain available should further discussion arise.
++tracked `issues <issue tracker>`_ remain available should further discussion arise.
  
  Fixed bug in v-rescale thermostat & replica exchange
  """"""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""
index dfb12c7e5a2ae8c581cdd8be8d075811cf75af58,e7de4da185648909d612a79c0d7a165cd78e9b2c..6b41402eb98eb2ad7f5b0c7f8c9bd22fc1d504f0
@@@ -8,40 -8,18 +8,39 @@@ releases of |Gromacs|. Major releases c
  functionality supported, whereas patch releases contain only fixes for
  issues identified in the corresponding major releases.
  
 -Two versions of |Gromacs| are under active maintenance, the 2020
 -series and the 2019 series. In the latter, only highly conservative
 +Two versions of |Gromacs| are under active maintenance, the 2021
 +series and the 2020 series. In the latter, only highly conservative
  fixes will be made, and only to address issues that affect scientific
  correctness. Naturally, some of those releases will be made after the
 -year 2019 ends, but we keep 2019 in the name so users understand how
 +year 2020 ends, but we keep 2019 in the name so users understand how
  up to date their version is. Such fixes will also be incorporated into
 -the 2020 release series, as appropriate. Around the time the 2021
 -release is made, the 2019 series will no longer be maintained.
 +the 2021 release series, as appropriate. Around the time the 2022
 +release is made, the 2020 series will no longer be maintained.
  
  Where issue numbers are reported in these release notes, more details
- can be found at https://gitlab.com/gromacs/gromacs/-/issues at that issue number.
 -can be found at https://gitlab.com/gromacs/gromacs/issues at that issue number.
++can be found on the `issue tracker`_ at that issue number.
 +
 +|Gromacs| 2021 series
 +---------------------
 +
 +.. todolist::
 +
 +Major release
 +^^^^^^^^^^^^^
 +
 +.. toctree::
 +   :maxdepth: 1
 +
 +   2021/major/highlights
 +   2021/major/features
 +   2021/major/performance
 +   2021/major/tools
 +   2021/major/bugs-fixed
 +   2021/major/deprecated-functionality
 +   2021/major/removed-functionality
 +   2021/major/portability
 +   2021/major/miscellaneous
  
  |Gromacs| 2020 series
  ---------------------
  
index d1488c855dec91163532b7ddfe33e0a8130c05bb,19ab95707d1503271d1f2f913c41f31b68ee14b1..b295040db745d94abe1231988ed7041458fc4360
@@@ -60,7 -60,7 +60,7 @@@
  #include "gromacs/domdec/ga2la.h"
  #include "gromacs/domdec/localatomsetmanager.h"
  #include "gromacs/domdec/mdsetup.h"
 -#include "gromacs/ewald/pme.h"
 +#include "gromacs/ewald/pme_pp.h"
  #include "gromacs/gmxlib/network.h"
  #include "gromacs/gmxlib/nrnb.h"
  #include "gromacs/imd/imd.h"
@@@ -75,7 -75,6 +75,7 @@@
  #include "gromacs/mdtypes/forcerec.h"
  #include "gromacs/mdtypes/inputrec.h"
  #include "gromacs/mdtypes/md_enums.h"
 +#include "gromacs/mdtypes/mdatom.h"
  #include "gromacs/mdtypes/nblist.h"
  #include "gromacs/mdtypes/state.h"
  #include "gromacs/nbnxm/nbnxm.h"
@@@ -210,7 -209,7 +210,7 @@@ static void dd_move_cellx(gmx_domdec_t
          if (applyPbc)
          {
              /* Take the minimum to avoid double communication */
 -            numPulsesMin = std::min(numPulses, dd->nc[dim] - 1 - numPulses);
 +            numPulsesMin = std::min(numPulses, dd->numCells[dim] - 1 - numPulses);
          }
          else
          {
          for (int pulse = 0; pulse < numPulses; pulse++)
          {
              /* Communicate all the zone information backward */
 -            bool receiveValidData = (applyPbc || dd->ci[dim] < dd->nc[dim] - 1);
 +            bool receiveValidData = (applyPbc || dd->ci[dim] < dd->numCells[dim] - 1);
  
              static_assert(
                      sizeof(gmx_ddzone_t) == c_ddzoneNumReals * sizeof(real),
                   */
                  buf_s[i] = buf_r[i];
              }
 -            if (((applyPbc || dd->ci[dim] + numPulses < dd->nc[dim]) && pulse == numPulses - 1)
 -                || (!applyPbc && dd->ci[dim] + 1 + pulse == dd->nc[dim] - 1))
 +            if (((applyPbc || dd->ci[dim] + numPulses < dd->numCells[dim]) && pulse == numPulses - 1)
 +                || (!applyPbc && dd->ci[dim] + 1 + pulse == dd->numCells[dim] - 1))
              {
                  /* Store the extremes */
                  int pos = 0;
@@@ -463,8 -462,8 +463,8 @@@ static void dd_set_cginfo(gmx::ArrayRef
  {
      if (fr != nullptr)
      {
 -        const cginfo_mb_t* cginfo_mb = fr->cginfo_mb;
 -        gmx::ArrayRef<int> cginfo    = fr->cginfo;
 +        gmx::ArrayRef<cginfo_mb_t> cginfo_mb = fr->cginfo_mb;
 +        gmx::ArrayRef<int>         cginfo    = fr->cginfo;
  
          for (int cg = cg0; cg < cg1; cg++)
          {
@@@ -733,7 -732,7 +733,7 @@@ static void comm_dd_ns_cell_sizes(gmx_d
          dim = dd->dim[dim_ind];
  
          /* Without PBC we don't have restrictions on the outer cells */
 -        if (!(dim >= ddbox->npbcdim && (dd->ci[dim] == 0 || dd->ci[dim] == dd->nc[dim] - 1))
 +        if (!(dim >= ddbox->npbcdim && (dd->ci[dim] == 0 || dd->ci[dim] == dd->numCells[dim] - 1))
              && isDlbOn(comm)
              && (comm->cell_x1[dim] - comm->cell_x0[dim]) * ddbox->skew_fac[dim] < comm->cellsize_min[dim])
          {
@@@ -864,7 -863,7 +864,7 @@@ static void get_load_distribution(gmx_d
                  load->mdf      = 0;
                  load->pme      = 0;
                  int pos        = 0;
 -                for (int i = 0; i < dd->nc[dim]; i++)
 +                for (int i = 0; i < dd->numCells[dim]; i++)
                  {
                      load->sum += load->load[pos++];
                      load->max = std::max(load->max, load->load[pos]);
                  }
                  if (isDlbOn(comm) && rowMaster->dlbIsLimited)
                  {
 -                    load->sum_m *= dd->nc[dim];
 +                    load->sum_m *= dd->numCells[dim];
                      load->flags |= (1 << d);
                  }
              }
@@@ -1259,7 -1258,7 +1259,7 @@@ static void turn_on_dlb(const gmx::MDLo
          {
              comm->load[d].sum_m = comm->load[d].sum;
  
 -            int nc = dd->nc[dd->dim[d]];
 +            int nc = dd->numCells[dd->dim[d]];
              for (int i = 0; i < nc; i++)
              {
                  rowMaster->cellFrac[i] = i / static_cast<real>(nc);
@@@ -1328,7 -1327,7 +1328,7 @@@ static void merge_cg_buffers(in
                               const int*                     recv_i,
                               gmx::ArrayRef<gmx::RVec>       x,
                               gmx::ArrayRef<const gmx::RVec> recv_vr,
 -                             cginfo_mb_t*                   cginfo_mb,
 +                             gmx::ArrayRef<cginfo_mb_t>     cginfo_mb,
                               gmx::ArrayRef<int>             cginfo)
  {
      gmx_domdec_ind_t *ind, *ind_p;
@@@ -1807,7 -1806,12 +1807,7 @@@ static void clearCommSetupData(dd_comm_
  }
  
  //! Prepare DD communication.
 -static void setup_dd_communication(gmx_domdec_t*                dd,
 -                                   matrix                       box,
 -                                   gmx_ddbox_t*                 ddbox,
 -                                   t_forcerec*                  fr,
 -                                   t_state*                     state,
 -                                   PaddedHostVector<gmx::RVec>* f)
 +static void setup_dd_communication(gmx_domdec_t* dd, matrix box, gmx_ddbox_t* ddbox, t_forcerec* fr, t_state* state)
  {
      int                    dim_ind, dim, dim0, dim1, dim2, dimd, nat_tot;
      int                    nzone, nzone_send, zone, zonei, cg0, cg1;
      gmx_domdec_comm_t*     comm;
      gmx_domdec_zones_t*    zones;
      gmx_domdec_comm_dim_t* cd;
 -    cginfo_mb_t*           cginfo_mb;
      gmx_bool               bBondComm, bDist2B, bDistMB, bDistBonded;
      dd_corners_t           corners;
      rvec *                 normal, *v_d, *v_0 = nullptr, *v_1 = nullptr;
          v_1 = ddbox->v[dim1];
      }
  
 -    zone_cg_range = zones->cg_range;
 -    cginfo_mb     = fr->cginfo_mb;
 +    zone_cg_range                        = zones->cg_range;
 +    gmx::ArrayRef<cginfo_mb_t> cginfo_mb = fr->cginfo_mb;
  
      zone_cg_range[0]   = 0;
      zone_cg_range[1]   = dd->ncg_home;
              ddSendrecv<int>(dd, dim_ind, dddirBackward, work.atomGroupBuffer, integerBufferRef);
  
              /* Make space for cg_cm */
 -            dd_check_alloc_ncg(fr, state, f, pos_cg + ind->nrecv[nzone]);
 +            dd_resize_atominfo_and_state(fr, state, pos_cg + ind->nrecv[nzone]);
  
              /* Communicate the coordinates */
              gmx::ArrayRef<gmx::RVec> rvecBufferRef;
@@@ -2624,27 -2629,27 +2624,27 @@@ void print_dd_statistics(const t_commre
  }
  
  //!\brief TODO Remove fplog when group scheme and charge groups are gone
 -void dd_partition_system(FILE*                        fplog,
 -                         const gmx::MDLogger&         mdlog,
 -                         int64_t                      step,
 -                         const t_commrec*             cr,
 -                         gmx_bool                     bMasterState,
 -                         int                          nstglobalcomm,
 -                         t_state*                     state_global,
 -                         const gmx_mtop_t&            top_global,
 -                         const t_inputrec*            ir,
 -                         gmx::ImdSession*             imdSession,
 -                         pull_t*                      pull_work,
 -                         t_state*                     state_local,
 -                         PaddedHostVector<gmx::RVec>* f,
 -                         gmx::MDAtoms*                mdAtoms,
 -                         gmx_localtop_t*              top_local,
 -                         t_forcerec*                  fr,
 -                         gmx_vsite_t*                 vsite,
 -                         gmx::Constraints*            constr,
 -                         t_nrnb*                      nrnb,
 -                         gmx_wallcycle*               wcycle,
 -                         gmx_bool                     bVerbose)
 +void dd_partition_system(FILE*                     fplog,
 +                         const gmx::MDLogger&      mdlog,
 +                         int64_t                   step,
 +                         const t_commrec*          cr,
 +                         gmx_bool                  bMasterState,
 +                         int                       nstglobalcomm,
 +                         t_state*                  state_global,
 +                         const gmx_mtop_t&         top_global,
 +                         const t_inputrec*         ir,
 +                         gmx::ImdSession*          imdSession,
 +                         pull_t*                   pull_work,
 +                         t_state*                  state_local,
 +                         gmx::ForceBuffers*        f,
 +                         gmx::MDAtoms*             mdAtoms,
 +                         gmx_localtop_t*           top_local,
 +                         t_forcerec*               fr,
 +                         gmx::VirtualSitesHandler* vsite,
 +                         gmx::Constraints*         constr,
 +                         t_nrnb*                   nrnb,
 +                         gmx_wallcycle*            wcycle,
 +                         gmx_bool                  bVerbose)
  {
      gmx_domdec_t*      dd;
      gmx_domdec_comm_t* comm;
  
          set_ddbox(*dd, true, DDMASTER(dd) ? state_global->box : nullptr, true, xGlobal, &ddbox);
  
 -        distributeState(mdlog, dd, top_global, state_global, ddbox, state_local, f);
 +        distributeState(mdlog, dd, top_global, state_global, ddbox, state_local);
  
          /* Ensure that we have space for the new distribution */
 -        dd_check_alloc_ncg(fr, state_local, f, dd->ncg_home);
 +        dd_resize_atominfo_and_state(fr, state_local, dd->ncg_home);
  
          inc_nrnb(nrnb, eNR_CGCM, comm->atomRanges.numHomeAtoms());
  
          wallcycle_sub_start(wcycle, ewcsDD_REDIST);
  
          ncgindex_set = dd->ncg_home;
 -        dd_redistribute_cg(fplog, step, dd, ddbox.tric_dir, state_local, f, fr, nrnb, &ncg_moved);
 +        dd_redistribute_cg(fplog, step, dd, ddbox.tric_dir, state_local, fr, nrnb, &ncg_moved);
  
          GMX_RELEASE_ASSERT(bSortCG, "Sorting is required after redistribution");
  
          dd_sort_state(dd, fr, state_local);
  
          /* After sorting and compacting we set the correct size */
 -        dd_resize_state(state_local, f, comm->atomRanges.numHomeAtoms());
 +        state_change_natoms(state_local, comm->atomRanges.numHomeAtoms());
  
          /* Rebuild all the indices */
          dd->ga2la->clear();
      make_dd_indices(dd, ncgindex_set);
  
      /* Setup up the communication and communicate the coordinates */
 -    setup_dd_communication(dd, state_local->box, &ddbox, fr, state_local, f);
 +    setup_dd_communication(dd, state_local->box, &ddbox, fr, state_local);
  
      /* Set the indices for the halo atoms */
      make_dd_indices(dd, dd->ncg_home);
      /* Set the charge group boundaries for neighbor searching */
      set_cg_boundaries(&comm->zones);
  
 -    if (fr->cutoff_scheme == ecutsVERLET)
 -    {
 -        /* When bSortCG=true, we have already set the size for zone 0 */
 -        set_zones_size(dd, state_local->box, &ddbox, bSortCG ? 1 : 0, comm->zones.n, 0);
 -    }
 +    /* When bSortCG=true, we have already set the size for zone 0 */
 +    set_zones_size(dd, state_local->box, &ddbox, bSortCG ? 1 : 0, comm->zones.n, 0);
  
      wallcycle_sub_stop(wcycle, ewcsDD_SETUPCOMM);
  
          switch (range)
          {
              case DDAtomRanges::Type::Vsites:
 -                if (vsite && vsite->numInterUpdategroupVsites)
 +                if (vsite && vsite->numInterUpdategroupVirtualSites())
                  {
                      n = dd_make_local_vsites(dd, n, top_local->idef.il);
                  }
       */
      state_local->natoms = comm->atomRanges.numAtomsTotal();
  
 -    dd_resize_state(state_local, f, state_local->natoms);
 +    state_change_natoms(state_local, state_local->natoms);
  
 -    if (fr->haveDirectVirialContributions)
 +    if (vsite && vsite->numInterUpdategroupVirtualSites())
      {
 -        if (vsite && vsite->numInterUpdategroupVsites)
 +        nat_f_novirsum = comm->atomRanges.end(DDAtomRanges::Type::Vsites);
 +    }
 +    else
 +    {
 +        if (EEL_FULL(ir->coulombtype) && dd->haveExclusions)
          {
 -            nat_f_novirsum = comm->atomRanges.end(DDAtomRanges::Type::Vsites);
 +            nat_f_novirsum = comm->atomRanges.end(DDAtomRanges::Type::Zones);
          }
          else
          {
 -            if (EEL_FULL(ir->coulombtype) && dd->haveExclusions)
 -            {
 -                nat_f_novirsum = comm->atomRanges.end(DDAtomRanges::Type::Zones);
 -            }
 -            else
 -            {
 -                nat_f_novirsum = comm->atomRanges.numHomeAtoms();
 -            }
 +            nat_f_novirsum = comm->atomRanges.numHomeAtoms();
          }
      }
 -    else
 -    {
 -        nat_f_novirsum = 0;
 -    }
  
      /* Set the number of atoms required for the force calculation.
       * Forces need to be constrained when doing energy
                          comm->atomRanges.end(DDAtomRanges::Type::Constraints), nat_f_novirsum);
  
      /* Update atom data for mdatoms and several algorithms */
 -    mdAlgorithmsSetupAtomData(cr, ir, top_global, top_local, fr, nullptr, mdAtoms, constr, vsite, nullptr);
 +    mdAlgorithmsSetupAtomData(cr, ir, top_global, top_local, fr, f, mdAtoms, constr, vsite, nullptr);
  
      auto mdatoms = mdAtoms->mdatoms();
      if (!thisRankHasDuty(cr, DUTY_PME))
                                  mdatoms->sigmaB, dd_pme_maxshift_x(dd), dd_pme_maxshift_y(dd));
      }
  
-     if (ir->bPull)
-     {
-         /* Update the local pull groups */
-         dd_make_local_pull_groups(cr, pull_work);
-     }
      if (dd->atomSets != nullptr)
      {
          /* Update the local atom sets */
          dd->atomSets->setIndicesInDomainDecomposition(*(dd->ga2la));
      }
  
+     // The pull group construction can need the atom sets updated above
+     if (ir->bPull)
+     {
+         /* Update the local pull groups */
+         dd_make_local_pull_groups(cr, pull_work);
+     }
      /* Update the local atoms to be communicated via the IMD protocol if bIMD is TRUE. */
      imdSession->dd_make_local_IMD_atoms(dd);
  
       * the last vsite construction, we need to communicate the constructing
       * atom coordinates again (for spreading the forces this MD step).
       */
 -    dd_move_x_vsites(dd, state_local->box, state_local->x.rvec_array());
 +    dd_move_x_vsites(*dd, state_local->box, state_local->x.rvec_array());
  
      wallcycle_sub_stop(wcycle, ewcsDD_TOPOTHER);
  
  }
  
  /*! \brief Check whether bonded interactions are missing, if appropriate */
 -void checkNumberOfBondedInteractions(const gmx::MDLogger&  mdlog,
 -                                     t_commrec*            cr,
 -                                     int                   totalNumberOfBondedInteractions,
 -                                     const gmx_mtop_t*     top_global,
 -                                     const gmx_localtop_t* top_local,
 -                                     const rvec*           x,
 -                                     const matrix          box,
 -                                     bool*                 shouldCheckNumberOfBondedInteractions)
 +void checkNumberOfBondedInteractions(const gmx::MDLogger&           mdlog,
 +                                     t_commrec*                     cr,
 +                                     int                            totalNumberOfBondedInteractions,
 +                                     const gmx_mtop_t*              top_global,
 +                                     const gmx_localtop_t*          top_local,
 +                                     gmx::ArrayRef<const gmx::RVec> x,
 +                                     const matrix                   box,
 +                                     bool* shouldCheckNumberOfBondedInteractions)
  {
      if (*shouldCheckNumberOfBondedInteractions)
      {
index 5fa50372bf70a467c0642d32076df33cc5b1f655,828edb4168faa635058e19c32dc34a22fa116c74..559c1f96372f7686f65fc3fe7ccc00bb92586c9a
@@@ -3,8 -3,7 +3,8 @@@
   *
   * Copyright (c) 1991-2000, University of Groningen, The Netherlands.
   * Copyright (c) 2001-2004, The GROMACS development team.
 - * Copyright (c) 2013,2014,2015,2016,2017,2018,2019,2020, by the GROMACS development team, led by
 + * Copyright (c) 2013,2014,2015,2016,2017 by the GROMACS development team.
 + * Copyright (c) 2018,2019,2020, by the GROMACS development team, led by
   * Mark Abraham, David van der Spoel, Berk Hess, and Erik Lindahl,
   * and including many others, as listed in the AUTHORS file in the
   * top-level source directory and at http://www.gromacs.org.
  #include "gromacs/pbcutil/rmpbc.h"
  #include "gromacs/topology/index.h"
  #include "gromacs/topology/topology.h"
+ #include "gromacs/utility/arrayref.h"
  #include "gromacs/utility/arraysize.h"
  #include "gromacs/utility/cstringutil.h"
+ #include "gromacs/utility/exceptions.h"
  #include "gromacs/utility/fatalerror.h"
  #include "gromacs/utility/futil.h"
  #include "gromacs/utility/smalloc.h"
  #include "gromacs/utility/sysinfo.h"
  
+ namespace gmx
+ {
+ namespace
+ {
+ /*! \brief Throw an error if any element in index exceeds a given number.
+  *
+  * \param[in] indices to be acessed
+  * \param[in] largestOkayIndex to be accessed
+  * \param[in] indexUsagePurpose to be more explicit in the error message
+  *
+  * \throws RangeError if largestOkayIndex is larger than any element in indices
+  *
+  */
+ void throwErrorIfIndexOutOfBounds(ArrayRef<const int> indices,
+                                   const int           largestOkayIndex,
+                                   const std::string&  indexUsagePurpose)
+ {
+     // do nothing if index is empty
+     if (indices.ssize() == 0)
+     {
+         return;
+     }
+     const int largestIndex = *std::max_element(indices.begin(), indices.end());
+     if (largestIndex < largestOkayIndex)
+     {
+         GMX_THROW(RangeError("The provided structure file only contains "
+                              + std::to_string(largestOkayIndex) + " coordinates, but coordinate index "
+                              + std::to_string(largestIndex) + " was requested for " + indexUsagePurpose
+                              + ". Make sure to update structure files "
+                                "and index files if you store only a part of your system."));
+     }
+ };
+ } // namespace
+ } // namespace gmx
  int gmx_covar(int argc, char* argv[])
  {
      const char* desc[] = {
      FILE*             out = nullptr; /* initialization makes all compilers happy */
      t_trxstatus*      status;
      t_topology        top;
 -    int               ePBC;
 +    PbcType           pbcType;
      t_atoms*          atoms;
      rvec *            x, *xread, *xref, *xav, *xproj;
      matrix            box, zerobox;
      xpmfile    = opt2fn_null("-xpm", NFILE, fnm);
      xpmafile   = opt2fn_null("-xpma", NFILE, fnm);
  
 -    read_tps_conf(fitfile, &top, &ePBC, &xref, nullptr, box, TRUE);
 +    read_tps_conf(fitfile, &top, &pbcType, &xref, nullptr, box, TRUE);
      atoms = &top.atoms;
  
      if (bFit)
      {
          printf("\nChoose a group for the least squares fit\n");
          get_index(atoms, ndxfile, 1, &nfit, &ifit, &fitname);
+         // Make sure that we never attempt to access a coordinate out of range
+         gmx::throwErrorIfIndexOutOfBounds({ ifit, ifit + nfit }, atoms->nr, "fitting");
          if (nfit < 3)
          {
              gmx_fatal(FARGS, "Need >= 3 points to fit!\n");
      }
      printf("\nChoose a group for the covariance analysis\n");
      get_index(atoms, ndxfile, 1, &natoms, &index, &ananame);
+     gmx::throwErrorIfIndexOutOfBounds({ index, index + natoms }, atoms->nr, "analysis");
  
      bDiffMass1 = FALSE;
      if (bFit)
      /* Prepare reference frame */
      if (bPBC)
      {
 -        gpbc = gmx_rmpbc_init(&top.idef, ePBC, atoms->nr);
 +        gpbc = gmx_rmpbc_init(&top.idef, pbcType, atoms->nr);
          gmx_rmpbc(gpbc, atoms->nr, box, xref);
      }
      if (bFit)
      nat      = read_first_x(oenv, &status, trxfile, &t, &xread, box);
      if (nat != atoms->nr)
      {
-         fprintf(stderr, "\nWARNING: number of atoms in tpx (%d) and trajectory (%d) do not match\n",
+         fprintf(stderr,
+                 "\nWARNING: number of atoms in structure file (%d) and trajectory (%d) do not "
+                 "match\n",
                  natoms, nat);
      }
+     gmx::throwErrorIfIndexOutOfBounds({ ifit, ifit + nfit }, nat, "fitting");
+     gmx::throwErrorIfIndexOutOfBounds({ index, index + natoms }, nat, "analysis");
      do
      {
          nframes0++;
          }
      }
      write_sto_conf_indexed(opt2fn("-av", NFILE, fnm), "Average structure", atoms, xread, nullptr,
 -                           epbcNONE, zerobox, natoms, index);
 +                           PbcType::No, zerobox, natoms, index);
      sfree(xread);
  
      fprintf(stderr, "Constructing covariance matrix (%dx%d) ...\n", static_cast<int>(ndim),
index b1c9aa11c52be1727d7ac054babecdc8d1235ff0,0d31a64b3aa671d72bf0d126e8cc12f9f3e72da2..1cfcde6f3876650d12fe9f2c2f8951bbd8cacf8f
@@@ -3,8 -3,7 +3,8 @@@
   *
   * Copyright (c) 1991-2000, University of Groningen, The Netherlands.
   * Copyright (c) 2001-2004, The GROMACS development team.
 - * Copyright (c) 2012,2013,2014,2015,2017,2018,2019,2020, by the GROMACS development team, led by
 + * Copyright (c) 2012,2013,2014,2015,2017 by the GROMACS development team.
 + * Copyright (c) 2018,2019,2020, by the GROMACS development team, led by
   * Mark Abraham, David van der Spoel, Berk Hess, and Erik Lindahl,
   * and including many others, as listed in the AUTHORS file in the
   * top-level source directory and at http://www.gromacs.org.
@@@ -259,17 -258,14 +259,14 @@@ static gmx_bool* bPhobics(t_atoms* atom
  
  static void check_oo(t_atoms* atoms)
  {
-     char* OOO;
+     char* OOO = gmx_strdup("O");
  
-     int i;
-     OOO = gmx_strdup("O");
-     for (i = 0; (i < atoms->nr); i++)
+     for (int i = 0; (i < atoms->nr); i++)
      {
          if ((std::strcmp(*(atoms->atomname[i]), "OXT") == 0)
              || (std::strcmp(*(atoms->atomname[i]), "O1") == 0)
-             || (std::strcmp(*(atoms->atomname[i]), "OC1") == 0))
+             || (std::strcmp(*(atoms->atomname[i]), "OC1") == 0)
+             || (std::strcmp(*(atoms->atomname[i]), "OT1") == 0))
          {
              *atoms->atomname[i] = OOO;
          }
@@@ -512,7 -508,7 +509,7 @@@ int gmx_do_dssp(int argc, char* argv[]
      const char *      fnSCount, *fnArea, *fnTArea, *fnAArea;
      const char*       leg[] = { "Phobic", "Phylic" };
      t_topology        top;
 -    int               ePBC;
 +    PbcType           pbcType;
      t_atoms*          atoms;
      t_matrix          mat;
      int               nres, nr0, naccr, nres_plus_separators;
      fnAArea    = opt2fn_null("-aa", NFILE, fnm);
      bDoAccSurf = ((fnArea != nullptr) || (fnTArea != nullptr) || (fnAArea != nullptr));
  
 -    read_tps_conf(ftp2fn(efTPS, NFILE, fnm), &top, &ePBC, &xp, nullptr, box, FALSE);
 +    read_tps_conf(ftp2fn(efTPS, NFILE, fnm), &top, &pbcType, &xp, nullptr, box, FALSE);
      atoms = &(top.atoms);
      check_oo(atoms);
      bPhbres = bPhobics(atoms);
      accr  = nullptr;
      naccr = 0;
  
 -    gpbc = gmx_rmpbc_init(&top.idef, ePBC, natoms);
 +    gpbc = gmx_rmpbc_init(&top.idef, pbcType, natoms);
      do
      {
          t = output_env_conv_time(oenv, t);
          }
          gmx_rmpbc(gpbc, natoms, box, x);
          tapein = gmx_ffopen(pdbfile, "w");
 -        write_pdbfile_indexed(tapein, nullptr, atoms, x, ePBC, box, ' ', -1, gnx, index, nullptr, FALSE);
 +        write_pdbfile_indexed(tapein, nullptr, atoms, x, pbcType, box, ' ', -1, gnx, index, nullptr, FALSE);
          gmx_ffclose(tapein);
          /* strip_dssp returns the number of lines found in the dssp file, i.e.
           * the number of residues plus the separator lines */
index 1e52fe0b491ac09cc4d5024cae6931318e39455b,889016981933f32132612ffb994899016ab1774f..509801a57292a95c46647a7da56057469c0ffa9d
@@@ -3,8 -3,7 +3,8 @@@
   *
   * Copyright (c) 1991-2000, University of Groningen, The Netherlands.
   * Copyright (c) 2001-2004, The GROMACS development team.
 - * Copyright (c) 2013,2014,2015,2016,2017,2018,2019,2020, by the GROMACS development team, led by
 + * Copyright (c) 2013,2014,2015,2016,2017 by the GROMACS development team.
 + * Copyright (c) 2018,2019,2020, by the GROMACS development team, led by
   * Mark Abraham, David van der Spoel, Berk Hess, and Erik Lindahl,
   * and including many others, as listed in the AUTHORS file in the
   * top-level source directory and at http://www.gromacs.org.
@@@ -72,7 -71,7 +72,7 @@@ static void calc_h2order(const char
                           real*                   slWidth,
                           int*                    nslices,
                           const t_topology*       top,
 -                         int                     ePBC,
 +                         PbcType                 pbcType,
                           int                     axis,
                           gmx_bool                bMicel,
                           int                     micel[],
  
      teller = 0;
  
 -    gpbc = gmx_rmpbc_init(&top->idef, ePBC, natoms);
 +    gpbc = gmx_rmpbc_init(&top->idef, pbcType, natoms);
      /*********** Start processing trajectory ***********/
      do
      {
@@@ -259,6 -258,15 +259,15 @@@ static void h2order_plot(rvec dipole[]
      xvgrclose(ord);
  }
  
+ enum
+ {
+     axisSEL,
+     axisZ,
+     axisY,
+     axisX,
+     axisNR
+ };
  int gmx_h2order(int argc, char* argv[])
  {
      const char* desc[] = {
          "dipole and the axis from the center of mass to the oxygen is calculated",
          "instead of the angle between the dipole and a box axis."
      };
-     static int         axis    = 2; /* normal to memb. default z  */
-     static const char* axtitle = "Z";
-     static int         nslices = 0; /* nr of slices defined       */
-     t_pargs            pa[]    = { { "-d",
+     static const char* axisOption[axisNR + 1] = { nullptr, "Z", "Y", "X", nullptr };
+     static int         nslices                = 0; /* nr of slices defined       */
+     // The struct that will hold the parsed user input
+     t_pargs     pa[]   = { { "-d",
                         FALSE,
-                        etSTR,
-                        { &axtitle },
+                        etENUM,
+                        { axisOption },
                         "Take the normal on the membrane in direction X, Y or Z." },
                       { "-sl",
                         FALSE,
                         { &nslices },
                         "Calculate order parameter as function of boxlength, dividing the box"
                         " in this number of slices." } };
-     const char*        bugs[]  = {
+     const char* bugs[] = {
          "The program assigns whole water molecules to a slice, based on the first "
          "atom of three in the index file group. It assumes an order O,H,H. "
          "Name is not important, but the order is. If this demand is not met, "
      int ngx,          /* nr. of atomsin sol group   */
              nmic = 0; /* nr. of atoms in micelle    */
      t_topology* top;  /* topology           */
 -    int         ePBC;
 +    PbcType     pbcType;
      int *       index, /* indices for solvent group  */
              *micelle = nullptr;
      gmx_bool bMicel  = FALSE; /* think we're a micel        */
  
  #define NFILE asize(fnm)
  
+     // Parse the user input in argv into pa
      if (!parse_common_args(&argc, argv, PCA_CAN_VIEW | PCA_CAN_TIME, NFILE, fnm, asize(pa), pa,
                             asize(desc), desc, asize(bugs), bugs, &oenv))
      {
          return 0;
      }
+     // Process the axis option chosen by the user to set the
+     // axis used for the computation. The useful choice is an
+     // axis normal to the membrane. Default is z-axis.
+     int axis = ZZ;
+     switch (nenum(axisOption))
+     {
+         case axisX: axis = XX; break;
+         case axisY: axis = YY; break;
+         case axisZ: axis = ZZ; break;
+         default: axis = ZZ;
+     }
      bMicel = opt2bSet("-nm", NFILE, fnm);
  
 -    top = read_top(ftp2fn(efTPR, NFILE, fnm), &ePBC); /* read topology file */
 +    top = read_top(ftp2fn(efTPR, NFILE, fnm), &pbcType); /* read topology file */
  
      rd_index(ftp2fn(efNDX, NFILE, fnm), 1, &ngx, &index, &grpname);
  
      }
  
      calc_h2order(ftp2fn(efTRX, NFILE, fnm), index, ngx, &slDipole, &slOrder, &slWidth, &nslices,
 -                 top, ePBC, axis, bMicel, micelle, nmic, oenv);
 +                 top, pbcType, axis, bMicel, micelle, nmic, oenv);
  
      h2order_plot(slDipole, slOrder, opt2fn("-o", NFILE, fnm), nslices, slWidth, oenv);
  
index 54e2eefb3e7ff4d7d55e1079800d4d776114acf1,e5c754e78f20b7d69bd7aa58e861d69246661fd0..d808249678abdb704be7cdbb13562ce846b41ae7
@@@ -1,7 -1,7 +1,7 @@@
  /*
   * This file is part of the GROMACS molecular simulation package.
   *
-  * Copyright (c) 2012,2013,2014,2015,2016 by the GROMACS development team.
+  * Copyright (c) 2012,2013,2014,2015,2016, The GROMACS development team.
   * Copyright (c) 2017,2018,2019,2020, by the GROMACS development team, led by
   * Mark Abraham, David van der Spoel, Berk Hess, and Erik Lindahl,
   * and including many others, as listed in the AUTHORS file in the
  
  #include <algorithm>
  #include <array>
 -#include <chrono>
  #include <memory>
  #include <string>
 -#include <thread>
  #include <vector>
  
 -#include "gromacs/compat/pointers.h"
 -#include "gromacs/gpu_utils/gpu_utils.h"
  #include "gromacs/hardware/cpuinfo.h"
 +#include "gromacs/hardware/device_management.h"
  #include "gromacs/hardware/hardwaretopology.h"
  #include "gromacs/hardware/hw_info.h"
  #include "gromacs/simd/support.h"
  #include "gromacs/utility/exceptions.h"
  #include "gromacs/utility/fatalerror.h"
  #include "gromacs/utility/gmxassert.h"
 +#include "gromacs/utility/inmemoryserializer.h"
  #include "gromacs/utility/logger.h"
 -#include "gromacs/utility/mutex.h"
  #include "gromacs/utility/physicalnodecommunicator.h"
  
  #include "architecture.h"
 +#include "device_information.h"
 +#include "prepare_detection.h"
  
  #ifdef HAVE_UNISTD_H
  #    include <unistd.h> // sysconf()
@@@ -76,7 -77,10 +76,7 @@@ gmx_hw_info_t::gmx_hw_info_t(std::uniqu
  {
  }
  
 -gmx_hw_info_t::~gmx_hw_info_t()
 -{
 -    free_gpu_info(&gpu_info);
 -}
 +gmx_hw_info_t::~gmx_hw_info_t() = default;
  
  namespace gmx
  {
  #    define _SC_NPROCESSORS_CONF _SC_NPROC_CONF
  #endif
  
 -/*! \brief Information about the hardware of all nodes (common to all threads in this process).
 +/*! \brief The result of device detection
   *
 - * This information is constructed only when required, but thereafter
 - * its lifetime is that of the whole process, potentially across
 - * multiple successive simulation parts. It's wise to ensure that only
 - * one thread can create the information, but thereafter they can all
 - * read it without e.g. needing a std::shared_ptr to ensure its
 - * lifetime exceeds that of the thread. */
 -static std::unique_ptr<gmx_hw_info_t> g_hardwareInfo;
 -//! A mutex to protect the hwinfo structure
 -static Mutex g_hardwareInfoMutex;
 -
 -//! Detect GPUs, if that makes sense to attempt.
 -static void gmx_detect_gpus(const gmx::MDLogger&             mdlog,
 -                            const PhysicalNodeCommunicator&  physicalNodeComm,
 -                            compat::not_null<gmx_hw_info_t*> hardwareInfo)
 + * Note that non-functional device detection still produces
 + * a detection result, ie. of no devices. This might not be
 + * what the user wanted, so it makes sense to log later when
 + * that is possible. */
 +struct DeviceDetectionResult
  {
 -    hardwareInfo->gpu_info.bDetectGPUs = canPerformGpuDetection();
 +    //! The device information detected
 +    std::vector<std::unique_ptr<DeviceInformation>> deviceInfoList_;
 +    //! Container of possible warnings to issue when that is possible
 +    std::vector<std::string> deviceDetectionWarnings_;
 +};
  
 -    if (!hardwareInfo->gpu_info.bDetectGPUs)
 +/*! \brief Detect GPUs when that makes sense to attempt.
 + *
 + * \param[in]  physicalNodeComm  The communicator across this physical node
 + * \return The result of the detection, perhaps including diagnostic messages
 + *         to issue later.
 + *
 + * \todo Coordinating the efficient detection of devices across
 + * multiple ranks per node should be separated from the lower-level
 + * hardware detection. See
 + * https://gitlab.com/gromacs/gromacs/-/issues/3650.
 + */
 +static DeviceDetectionResult detectAllDeviceInformation(const PhysicalNodeCommunicator& physicalNodeComm)
 +{
 +    DeviceDetectionResult deviceDetectionResult;
 +
 +    if (!isDeviceDetectionEnabled())
      {
 -        return;
 +        return deviceDetectionResult;
      }
  
 +    std::string errorMessage;
 +
      bool isMasterRankOfPhysicalNode = true;
  #if GMX_LIB_MPI
      isMasterRankOfPhysicalNode = (physicalNodeComm.rank_ == 0);
  #else
 -    // We choose to run the detection only once with thread-MPI and
 -    // use a mutex to enforce it.
 +    // Without an MPI library, this process is trivially the only one
 +    // on the physical node. This code runs before e.g. thread-MPI
 +    // ranks are spawned, so detection is race-free by construction.
 +    // Read-only access is enforced with providing those ranks with a
 +    // handle to a const object, so usage is also free of races.
      GMX_UNUSED_VALUE(physicalNodeComm);
 -    isMasterRankOfPhysicalNode = true;
 +    isMasterRankOfPhysicalNode        = true;
  #endif
  
 -    /* The OpenCL support requires us to run detection on all ranks.
 +    /* The SYCL and OpenCL support requires us to run detection on all
 +     * ranks.
 +     *
       * With CUDA we don't need to, and prefer to detect on one rank
 -     * and send the information to the other ranks over MPI. */
 -    bool allRanksMustDetectGpus = (GMX_GPU == GMX_GPU_OPENCL);
 -    bool gpusCanBeDetected      = false;
 +     * and send the information to the other ranks over MPI. This
 +     * avoids creating a start-up bottleneck with each MPI rank on a
 +     * node making the same GPU API calls. */
 +    constexpr bool allRanksMustDetectGpus = (GMX_GPU_OPENCL != 0 || GMX_GPU_SYCL != 0);
 +    bool           gpusCanBeDetected      = false;
      if (isMasterRankOfPhysicalNode || allRanksMustDetectGpus)
      {
          std::string errorMessage;
 -        gpusCanBeDetected = isGpuDetectionFunctional(&errorMessage);
 +        gpusCanBeDetected = isDeviceDetectionFunctional(&errorMessage);
          if (!gpusCanBeDetected)
          {
 -            GMX_LOG(mdlog.info)
 -                    .asParagraph()
 -                    .appendTextFormatted(
 -                            "NOTE: Detection of GPUs failed. The API reported:\n"
 -                            "      %s\n"
 -                            "      GROMACS cannot run tasks on a GPU.",
 -                            errorMessage.c_str());
 +            deviceDetectionResult.deviceDetectionWarnings_.emplace_back(
 +                    "Detection of GPUs failed. The API reported:\n" + errorMessage);
          }
      }
  
      if (gpusCanBeDetected)
      {
 -        findGpus(&hardwareInfo->gpu_info);
 +        deviceDetectionResult.deviceInfoList_ = findDevices();
          // No need to tell the user anything at this point, they get a
          // hardware report later.
      }
  
  #if GMX_LIB_MPI
 -    if (!allRanksMustDetectGpus)
 +    if (!allRanksMustDetectGpus && (physicalNodeComm.size_ > 1))
      {
 -        /* Broadcast the GPU info to the other ranks within this node */
 -        MPI_Bcast(&hardwareInfo->gpu_info.n_dev, 1, MPI_INT, 0, physicalNodeComm.comm_);
 -
 -        if (hardwareInfo->gpu_info.n_dev > 0)
 +        // Master rank must serialize the device information list and
 +        // send it to the other ranks on this node.
 +        std::vector<char> buffer;
 +        int               sizeOfBuffer;
 +        if (isMasterRankOfPhysicalNode)
          {
 -            int dev_size;
 -
 -            dev_size = hardwareInfo->gpu_info.n_dev * sizeof_gpu_dev_info();
 -
 +            gmx::InMemorySerializer writer;
 +            serializeDeviceInformations(deviceDetectionResult.deviceInfoList_, &writer);
 +            buffer       = writer.finishAndGetBuffer();
 +            sizeOfBuffer = buffer.size();
 +        }
 +        // Ensure all ranks agree on the size of list to be sent
 +        MPI_Bcast(&sizeOfBuffer, 1, MPI_INT, 0, physicalNodeComm.comm_);
 +        buffer.resize(sizeOfBuffer);
 +        if (!buffer.empty())
 +        {
 +            // Send the list and deserialize it
 +            MPI_Bcast(buffer.data(), buffer.size(), MPI_BYTE, 0, physicalNodeComm.comm_);
              if (!isMasterRankOfPhysicalNode)
              {
 -                hardwareInfo->gpu_info.gpu_dev = (struct gmx_device_info_t*)malloc(dev_size);
 +                gmx::InMemoryDeserializer reader(buffer, false);
 +                deviceDetectionResult.deviceInfoList_ = deserializeDeviceInformations(&reader);
              }
 -            MPI_Bcast(hardwareInfo->gpu_info.gpu_dev, dev_size, MPI_BYTE, 0, physicalNodeComm.comm_);
 -            MPI_Bcast(&hardwareInfo->gpu_info.n_dev_compatible, 1, MPI_INT, 0, physicalNodeComm.comm_);
          }
      }
  #endif
 +    return deviceDetectionResult;
  }
  
  //! Reduce the locally collected \p hardwareInfo over MPI ranks
 -static void gmx_collect_hardware_mpi(const gmx::CpuInfo&              cpuInfo,
 -                                     const PhysicalNodeCommunicator&  physicalNodeComm,
 -                                     compat::not_null<gmx_hw_info_t*> hardwareInfo)
 +static void gmx_collect_hardware_mpi(const gmx::CpuInfo&             cpuInfo,
 +                                     const PhysicalNodeCommunicator& physicalNodeComm,
 +                                     gmx_hw_info_t*                  hardwareInfo)
  {
      const int ncore = hardwareInfo->hardwareTopology->numberOfCores();
      /* Zen1 is assumed for:
                                  && (cpuInfo.model() == 1 || cpuInfo.model() == 17
                                      || cpuInfo.model() == 8 || cpuInfo.model() == 24))
                                 || cpuInfo.vendor() == CpuInfo::Vendor::Hygon);
 +
 +    int numCompatibleDevices = getCompatibleDevices(hardwareInfo->deviceInfoList).size();
  #if GMX_LIB_MPI
 -    int nhwthread, ngpu, i;
 +    int nhwthread;
      int gpu_hash;
  
      nhwthread = hardwareInfo->nthreads_hw_avail;
 -    ngpu      = hardwareInfo->gpu_info.n_dev_compatible;
      /* Create a unique hash of the GPU type(s) in this node */
      gpu_hash = 0;
      /* Here it might be better to only loop over the compatible GPU, but we
       * don't have that information available and it would also require
       * removing the device ID from the device info string.
       */
 -    for (i = 0; i < hardwareInfo->gpu_info.n_dev; i++)
 +    for (const auto& deviceInfo : hardwareInfo->deviceInfoList)
      {
 -        char stmp[STRLEN];
 -
          /* Since the device ID is incorporated in the hash, the order of
           * the GPUs affects the hash. Also two identical GPUs won't give
           * a gpu_hash of zero after XORing.
           */
 -        get_gpu_device_info_string(stmp, hardwareInfo->gpu_info, i);
 -        gpu_hash ^= gmx_string_fullhash_func(stmp, gmx_string_hash_init);
 +        std::string deviceInfoString = getDeviceInformationString(*deviceInfo);
 +        gpu_hash ^= gmx_string_fullhash_func(deviceInfoString.c_str(), gmx_string_hash_init);
      }
  
      constexpr int                      numElementsCounts = 4;
              countsLocal[0] = 1;
              countsLocal[1] = ncore;
              countsLocal[2] = nhwthread;
 -            countsLocal[3] = ngpu;
 +            countsLocal[3] = numCompatibleDevices;
          }
  
          MPI_Allreduce(countsLocal.data(), countsReduced.data(), countsLocal.size(), MPI_INT,
           */
          maxMinLocal[0]  = ncore;
          maxMinLocal[1]  = nhwthread;
 -        maxMinLocal[2]  = ngpu;
 +        maxMinLocal[2]  = numCompatibleDevices;
          maxMinLocal[3]  = static_cast<int>(gmx::simdSuggested(cpuInfo));
          maxMinLocal[4]  = gpu_hash;
          maxMinLocal[5]  = -maxMinLocal[0];
      hardwareInfo->bIdenticalGPUs      = (maxMinReduced[4] == -maxMinReduced[9]);
      hardwareInfo->haveAmdZen1Cpu      = (maxMinReduced[10] > 0);
  #else
 -    /* All ranks use the same pointer, protected by a mutex in the caller */
      hardwareInfo->nphysicalnode       = 1;
      hardwareInfo->ncore_tot           = ncore;
      hardwareInfo->ncore_min           = ncore;
      hardwareInfo->nhwthread_tot       = hardwareInfo->nthreads_hw_avail;
      hardwareInfo->nhwthread_min       = hardwareInfo->nthreads_hw_avail;
      hardwareInfo->nhwthread_max       = hardwareInfo->nthreads_hw_avail;
 -    hardwareInfo->ngpu_compatible_tot = hardwareInfo->gpu_info.n_dev_compatible;
 -    hardwareInfo->ngpu_compatible_min = hardwareInfo->gpu_info.n_dev_compatible;
 -    hardwareInfo->ngpu_compatible_max = hardwareInfo->gpu_info.n_dev_compatible;
 +    hardwareInfo->ngpu_compatible_tot = numCompatibleDevices;
 +    hardwareInfo->ngpu_compatible_min = numCompatibleDevices;
 +    hardwareInfo->ngpu_compatible_max = numCompatibleDevices;
      hardwareInfo->simd_suggest_min    = static_cast<int>(simdSuggested(cpuInfo));
      hardwareInfo->simd_suggest_max    = static_cast<int>(simdSuggested(cpuInfo));
      hardwareInfo->bIdenticalGPUs      = TRUE;
  #endif
  }
  
 -/*! \brief Utility that does dummy computing for max 2 seconds to spin up cores
 - *
 - *  This routine will check the number of cores configured and online
 - *  (using sysconf), and the spins doing dummy compute operations for up to
 - *  2 seconds, or until all cores have come online. This can be used prior to
 - *  hardware detection for platforms that take unused processors offline.
 - *
 - *  This routine will not throw exceptions. In principle it should be
 - *  declared noexcept, but at least icc 19.1 and 21-beta08 with the
 - *  libstdc++-7.5 has difficulty implementing a std::vector of
 - *  std::thread started with this function when declared noexcept. It
 - *  is a known compiler bug that should be fixed after 19.1.
 - *  Fortunately, this function is not performance sensitive,
 - *  and only runs on platforms other than x86 and POWER (ie ARM),
 - *  so the possible overhead introduced by omitting noexcept is not
 - *  important.
 - */
 -static void spinUpCore()
 -{
 -#if defined(HAVE_SYSCONF) && defined(_SC_NPROCESSORS_CONF) && defined(_SC_NPROCESSORS_ONLN)
 -    float dummy           = 0.1;
 -    int   countConfigured = sysconf(_SC_NPROCESSORS_CONF);    // noexcept
 -    auto  start           = std::chrono::steady_clock::now(); // noexcept
 -
 -    while (sysconf(_SC_NPROCESSORS_ONLN) < countConfigured
 -           && std::chrono::steady_clock::now() - start < std::chrono::seconds(2))
 -    {
 -        for (int i = 1; i < 10000; i++)
 -        {
 -            dummy /= i;
 -        }
 -    }
 -
 -    if (dummy < 0)
 -    {
 -        printf("This cannot happen, but prevents loop from being optimized away.");
 -    }
 -#endif
 -}
 -
 -/*! \brief Prepare the system before hardware topology detection
 - *
 - * This routine should perform any actions we want to put the system in a state
 - * where we want it to be before detecting the hardware topology. For most
 - * processors there is nothing to do, but some architectures (in particular ARM)
 - * have support for taking configured cores offline, which will make them disappear
 - * from the online processor count.
 - *
 - * This routine checks if there is a mismatch between the number of cores
 - * configured and online, and in that case we issue a small workload that
 - * attempts to wake sleeping cores before doing the actual detection.
 - *
 - * This type of mismatch can also occur for x86 or PowerPC on Linux, if SMT has only
 - * been disabled in the kernel (rather than bios). Since those cores will never
 - * come online automatically, we currently skip this test for x86 & PowerPC to
 - * avoid wasting 2 seconds. We also skip the test if there is no thread support.
 - *
 - * \note Cores will sleep relatively quickly again, so it's important to issue
 - *       the real detection code directly after this routine.
 - */
 -static void hardwareTopologyPrepareDetection()
 -{
 -#if defined(HAVE_SYSCONF) && defined(_SC_NPROCESSORS_CONF) \
 -        && (defined(THREAD_PTHREADS) || defined(THREAD_WINDOWS))
 -
 -    // Modify this conditional when/if x86 or PowerPC starts to sleep some cores
 -    if (c_architecture != Architecture::X86 && c_architecture != Architecture::PowerPC)
 -    {
 -        int                      countConfigured = sysconf(_SC_NPROCESSORS_CONF);
 -        std::vector<std::thread> workThreads(countConfigured);
 -
 -        for (auto& t : workThreads)
 -        {
 -            t = std::thread(spinUpCore);
 -        }
 -
 -        for (auto& t : workThreads)
 -        {
 -            t.join();
 -        }
 -    }
 -#endif
 -}
 -
 -/*! \brief Sanity check hardware topology and print some notes to log
 - *
 - *  \param mdlog            Logger.
 - *  \param hardwareTopology Reference to hardwareTopology object.
 - */
 -static void hardwareTopologyDoubleCheckDetection(const gmx::MDLogger gmx_unused& mdlog,
 -                                                 const gmx::HardwareTopology gmx_unused& hardwareTopology)
 +void hardwareTopologyDoubleCheckDetection(const gmx::MDLogger gmx_unused& mdlog,
 +                                          const gmx::HardwareTopology gmx_unused& hardwareTopology)
  {
  #if defined HAVE_SYSCONF && defined(_SC_NPROCESSORS_CONF)
      if (hardwareTopology.supportLevel() < gmx::HardwareTopology::SupportLevel::LogicalProcessorCount)
                              "performance.");
          }
      }
 +#else
 +    GMX_UNUSED_VALUE(mdlog);
 +    GMX_UNUSED_VALUE(hardwareTopology);
  #endif
  }
  
 -gmx_hw_info_t* gmx_detect_hardware(const gmx::MDLogger& mdlog, const PhysicalNodeCommunicator& physicalNodeComm)
 +std::unique_ptr<gmx_hw_info_t> gmx_detect_hardware(const PhysicalNodeCommunicator& physicalNodeComm)
  {
 -    // By construction, only one thread ever runs hardware detection,
 -    // but we may as well prevent issues arising if that would change.
 -    // Taking the lock early ensures that exactly one thread can
 -    // attempt to construct g_hardwareInfo.
 -    lock_guard<Mutex> lock(g_hardwareInfoMutex);
 -
 -    // If we already have the information, just return a handle to it.
 -    if (g_hardwareInfo != nullptr)
 -    {
 -        return g_hardwareInfo.get();
 -    }
 -
 -    // Make the new hardwareInfo in a temporary.
 +    // Ensure all cores have spun up, where applicable.
      hardwareTopologyPrepareDetection();
  
      // TODO: We should also do CPU hardware detection only once on each
              std::make_unique<CpuInfo>(CpuInfo::detect()),
              std::make_unique<HardwareTopology>(HardwareTopology::detect()));
  
 -    // If we detected the topology on this system, double-check that it makes sense
 -    if (hardwareInfo->hardwareTopology->isThisSystem())
 -    {
 -        hardwareTopologyDoubleCheckDetection(mdlog, *hardwareInfo->hardwareTopology);
 -    }
 -
      // TODO: Get rid of this altogether.
      hardwareInfo->nthreads_hw_avail = hardwareInfo->hardwareTopology->machine().logicalProcessorCount;
  
      // Detect GPUs
 -    hardwareInfo->gpu_info.n_dev            = 0;
 -    hardwareInfo->gpu_info.n_dev_compatible = 0;
 -    hardwareInfo->gpu_info.gpu_dev          = nullptr;
 -
 -    gmx_detect_gpus(mdlog, physicalNodeComm, compat::make_not_null(hardwareInfo));
 -    gmx_collect_hardware_mpi(*hardwareInfo->cpuInfo, physicalNodeComm, compat::make_not_null(hardwareInfo));
 +    // Open a nested scope so no temporary variables can
 +    // be mis-used later.
 +    {
 +        DeviceDetectionResult deviceDetectionResult = detectAllDeviceInformation(physicalNodeComm);
 +        hardwareInfo->deviceInfoList.swap(deviceDetectionResult.deviceInfoList_);
 +        std::swap(hardwareInfo->hardwareDetectionWarnings_, deviceDetectionResult.deviceDetectionWarnings_);
 +    }
  
 -    // Now that the temporary is fully constructed, swap it to become
 -    // the real thing.
 -    g_hardwareInfo.swap(hardwareInfo);
 +    gmx_collect_hardware_mpi(*hardwareInfo->cpuInfo, physicalNodeComm, hardwareInfo.get());
  
 -    return g_hardwareInfo.get();
 +    return hardwareInfo;
  }
  
 -bool compatibleGpusFound(const gmx_gpu_info_t& gpu_info)
 +void logHardwareDetectionWarnings(const gmx::MDLogger& mdlog, const gmx_hw_info_t& hardwareInformation)
  {
 -    return gpu_info.n_dev_compatible > 0;
 +    for (const std::string& warningString : hardwareInformation.hardwareDetectionWarnings_)
 +    {
 +        GMX_LOG(mdlog.warning).asParagraph().appendText(warningString);
 +    }
  }
  
  } // namespace gmx
index 885afdc3b2ae64c0c2dad6461f4cbc795ddc057b,0000000000000000000000000000000000000000..853192f52a70c810c54d8f5ee6ab502fedc727b1
mode 100644,000000..100644
--- /dev/null
@@@ -1,134 -1,0 +1,134 @@@
-  *  is not clear whether the problem is the compiler or the standard
-  *  library. Fortunately, this function is not performance sensitive,
 +/*
 + * This file is part of the GROMACS molecular simulation package.
 + *
 + * Copyright (c) 2012,2013,2014,2015,2016 by the GROMACS development team.
 + * Copyright (c) 2017,2018,2019,2020, by the GROMACS development team, led by
 + * Mark Abraham, David van der Spoel, Berk Hess, and Erik Lindahl,
 + * and including many others, as listed in the AUTHORS file in the
 + * top-level source directory and at http://www.gromacs.org.
 + *
 + * GROMACS is free software; you can redistribute it and/or
 + * modify it under the terms of the GNU Lesser General Public License
 + * as published by the Free Software Foundation; either version 2.1
 + * of the License, or (at your option) any later version.
 + *
 + * GROMACS is distributed in the hope that it will be useful,
 + * but WITHOUT ANY WARRANTY; without even the implied warranty of
 + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
 + * Lesser General Public License for more details.
 + *
 + * You should have received a copy of the GNU Lesser General Public
 + * License along with GROMACS; if not, see
 + * http://www.gnu.org/licenses, or write to the Free Software Foundation,
 + * Inc., 51 Franklin Street, Fifth Floor, Boston, MA  02110-1301  USA.
 + *
 + * If you want to redistribute modifications to GROMACS, please
 + * consider that scientific software is very special. Version
 + * control is crucial - bugs must be traceable. We will be happy to
 + * consider code for inclusion in the official distribution, but
 + * derived work must not be called official GROMACS. Details are found
 + * in the README & COPYING files - if they are missing, get the
 + * official version at http://www.gromacs.org.
 + *
 + * To help us fund GROMACS development, we humbly ask that you cite
 + * the research papers on the package. Check out http://www.gromacs.org.
 + */
 +/*! \internal \file
 + * \brief Defines routine for activating potentially deactivated cores
 + * so they can be detected.
 + *
 + * The use of std::thread makes for brittle interaction with std
 + * library headers. Its caller also handles GPU detection and
 + * allocation of device-specific data structures. This is more
 + * manageable when separated into two distinct translation units.
 + *
 + * \author Erik Lindahl <erik.lindahl@scilifelab.se>
 + * \author Mark Abraham <mark.j.abraham@gmail.com>
 + * \ingroup module_hardware
 + */
 +#include "gmxpre.h"
 +
 +#include "prepare_detection.h"
 +
 +#include "config.h"
 +
 +#include <cstdio>
 +
 +#include <chrono>
 +#include <thread>
 +#include <vector>
 +
 +#include "architecture.h"
 +
 +#ifdef HAVE_UNISTD_H
 +#    include <unistd.h> // sysconf()
 +#endif
 +
 +namespace gmx
 +{
 +
 +/*! \brief Utility that does dummy computing for max 2 seconds to spin up cores
 + *
 + *  This routine will check the number of cores configured and online
 + *  (using sysconf), and the spins doing dummy compute operations for up to
 + *  2 seconds, or until all cores have come online. This can be used prior to
 + *  hardware detection for platforms that take unused processors offline.
 + *
 + *  This routine will not throw exceptions. In principle it should be
 + *  declared noexcept, but at least icc 19.1 and 21-beta08 with the
 + *  libstdc++-7.5 has difficulty implementing a std::vector of
 + *  std::thread started with this function when declared noexcept. It
++ *  is a known compiler bug that should be fixed after 19.1.
++ *  Fortunately, this function is not performance sensitive,
 + *  and only runs on platforms other than x86 and POWER (ie ARM),
 + *  so the possible overhead introduced by omitting noexcept is not
 + *  important.
 + */
 +static void spinUpCore()
 +{
 +#if defined(HAVE_SYSCONF) && defined(_SC_NPROCESSORS_CONF) && defined(_SC_NPROCESSORS_ONLN)
 +    float dummy           = 0.1;
 +    int   countConfigured = sysconf(_SC_NPROCESSORS_CONF);    // noexcept
 +    auto  start           = std::chrono::steady_clock::now(); // noexcept
 +
 +    while (sysconf(_SC_NPROCESSORS_ONLN) < countConfigured
 +           && std::chrono::steady_clock::now() - start < std::chrono::seconds(2))
 +    {
 +        for (int i = 1; i < 10000; i++)
 +        {
 +            dummy /= i;
 +        }
 +    }
 +
 +    if (dummy < 0)
 +    {
 +        printf("This cannot happen, but prevents loop from being optimized away.");
 +    }
 +#endif
 +}
 +
 +void hardwareTopologyPrepareDetection()
 +{
 +#if defined(HAVE_SYSCONF) && defined(_SC_NPROCESSORS_CONF) \
 +        && (defined(THREAD_PTHREADS) || defined(THREAD_WINDOWS))
 +
 +    // Modify this conditional when/if x86 or PowerPC starts to sleep some cores
 +    if (c_architecture != Architecture::X86 && c_architecture != Architecture::PowerPC)
 +    {
 +        int                      countConfigured = sysconf(_SC_NPROCESSORS_CONF);
 +        std::vector<std::thread> workThreads(countConfigured);
 +
 +        for (auto& t : workThreads)
 +        {
 +            t = std::thread(spinUpCore);
 +        }
 +
 +        for (auto& t : workThreads)
 +        {
 +            t.join();
 +        }
 +    }
 +#endif
 +}
 +
 +} // namespace gmx
index 43ef9d2459be5a85aa439296a3edd83fc2ae9a00,bdf47ab3e84324b059e8a586b94e5f2ae9bf3ac4..e56bcada93ffee70985837cd004719c09761a05a
@@@ -3,8 -3,7 +3,8 @@@
   *
   * Copyright (c) 1991-2000, University of Groningen, The Netherlands.
   * Copyright (c) 2001-2004, The GROMACS development team.
 - * Copyright (c) 2013,2014,2015,2016,2017,2018,2019,2020, by the GROMACS development team, led by
 + * Copyright (c) 2013,2014,2015,2016,2017 by the GROMACS development team.
 + * Copyright (c) 2018,2019,2020, by the GROMACS development team, led by
   * Mark Abraham, David van der Spoel, Berk Hess, and Erik Lindahl,
   * and including many others, as listed in the AUTHORS file in the
   * top-level source directory and at http://www.gromacs.org.
@@@ -37,8 -36,6 +37,8 @@@
   */
  #include "gmxpre.h"
  
 +#include "coupling.h"
 +
  #include <cassert>
  #include <cmath>
  
@@@ -51,7 -48,6 +51,7 @@@
  #include "gromacs/math/units.h"
  #include "gromacs/math/vec.h"
  #include "gromacs/math/vecdump.h"
 +#include "gromacs/mdlib/boxdeformation.h"
  #include "gromacs/mdlib/expanded.h"
  #include "gromacs/mdlib/gmx_omp_nthreads.h"
  #include "gromacs/mdlib/stat.h"
@@@ -92,239 -88,6 +92,239 @@@ static const double sy_const_5[] = { 0.
  
  static const double* sy_const[] = { nullptr, sy_const_1, nullptr, sy_const_3, nullptr, sy_const_5 };
  
 +
 +void update_tcouple(int64_t           step,
 +                    const t_inputrec* inputrec,
 +                    t_state*          state,
 +                    gmx_ekindata_t*   ekind,
 +                    const t_extmass*  MassQ,
 +                    const t_mdatoms*  md)
 +
 +{
 +    // This condition was explicitly checked in previous version, but should have never been satisfied
 +    GMX_ASSERT(!(EI_VV(inputrec->eI)
 +                 && (inputrecNvtTrotter(inputrec) || inputrecNptTrotter(inputrec)
 +                     || inputrecNphTrotter(inputrec))),
 +               "Temperature coupling was requested with velocity verlet and trotter");
 +
 +    bool doTemperatureCoupling = false;
 +
 +    // For VV temperature coupling parameters are updated on the current
 +    // step, for the others - one step before.
 +    if (inputrec->etc == etcNO)
 +    {
 +        doTemperatureCoupling = false;
 +    }
 +    else if (EI_VV(inputrec->eI))
 +    {
 +        doTemperatureCoupling = do_per_step(step, inputrec->nsttcouple);
 +    }
 +    else
 +    {
 +        doTemperatureCoupling = do_per_step(step + inputrec->nsttcouple - 1, inputrec->nsttcouple);
 +    }
 +
 +    if (doTemperatureCoupling)
 +    {
 +        real dttc = inputrec->nsttcouple * inputrec->delta_t;
 +
 +        // TODO: berendsen_tcoupl(...), nosehoover_tcoupl(...) and vrescale_tcoupl(...) update
 +        //      temperature coupling parameters, which should be reflected in the name of these
 +        //      subroutines
 +        switch (inputrec->etc)
 +        {
 +            case etcNO: break;
 +            case etcBERENDSEN:
 +                berendsen_tcoupl(inputrec, ekind, dttc, state->therm_integral);
 +                break;
 +            case etcNOSEHOOVER:
 +                nosehoover_tcoupl(&(inputrec->opts), ekind, dttc, state->nosehoover_xi.data(),
 +                                  state->nosehoover_vxi.data(), MassQ);
 +                break;
 +            case etcVRESCALE:
 +                vrescale_tcoupl(inputrec, step, ekind, dttc, state->therm_integral.data());
 +                break;
 +        }
 +        /* rescale in place here */
 +        if (EI_VV(inputrec->eI))
 +        {
 +            rescale_velocities(ekind, md, 0, md->homenr, state->v.rvec_array());
 +        }
 +    }
 +    else
 +    {
 +        // Set the T scaling lambda to 1 to have no scaling
 +        // TODO: Do we have to do it on every non-t-couple step?
 +        for (int i = 0; (i < inputrec->opts.ngtc); i++)
 +        {
 +            ekind->tcstat[i].lambda = 1.0;
 +        }
 +    }
 +}
 +
 +void update_pcouple_before_coordinates(FILE*             fplog,
 +                                       int64_t           step,
 +                                       const t_inputrec* inputrec,
 +                                       t_state*          state,
 +                                       matrix            parrinellorahmanMu,
 +                                       matrix            M,
 +                                       gmx_bool          bInitStep)
 +{
 +    /* Berendsen P-coupling is completely handled after the coordinate update.
 +     * Trotter P-coupling is handled by separate calls to trotter_update().
 +     */
 +    if (inputrec->epc == epcPARRINELLORAHMAN
 +        && do_per_step(step + inputrec->nstpcouple - 1, inputrec->nstpcouple))
 +    {
 +        real dtpc = inputrec->nstpcouple * inputrec->delta_t;
 +
 +        parrinellorahman_pcoupl(fplog, step, inputrec, dtpc, state->pres_prev, state->box,
 +                                state->box_rel, state->boxv, M, parrinellorahmanMu, bInitStep);
 +    }
 +}
 +
 +void update_pcouple_after_coordinates(FILE*                fplog,
 +                                      int64_t              step,
 +                                      const t_inputrec*    inputrec,
 +                                      const t_mdatoms*     md,
 +                                      const matrix         pressure,
 +                                      const matrix         forceVirial,
 +                                      const matrix         constraintVirial,
 +                                      matrix               pressureCouplingMu,
 +                                      t_state*             state,
 +                                      t_nrnb*              nrnb,
 +                                      gmx::BoxDeformation* boxDeformation,
 +                                      const bool           scaleCoordinates)
 +{
 +    int start  = 0;
 +    int homenr = md->homenr;
 +
 +    /* Cast to real for faster code, no loss in precision (see comment above) */
 +    real dt = inputrec->delta_t;
 +
 +
 +    /* now update boxes */
 +    switch (inputrec->epc)
 +    {
 +        case (epcNO): break;
 +        case (epcBERENDSEN):
 +            if (do_per_step(step, inputrec->nstpcouple))
 +            {
 +                real dtpc = inputrec->nstpcouple * dt;
 +                berendsen_pcoupl(fplog, step, inputrec, dtpc, pressure, state->box, forceVirial,
 +                                 constraintVirial, pressureCouplingMu, &state->baros_integral);
 +                berendsen_pscale(inputrec, pressureCouplingMu, state->box, state->box_rel, start,
 +                                 homenr, state->x.rvec_array(), md->cFREEZE, nrnb, scaleCoordinates);
 +            }
 +            break;
 +        case (epcCRESCALE):
 +            if (do_per_step(step, inputrec->nstpcouple))
 +            {
 +                real dtpc = inputrec->nstpcouple * dt;
 +                crescale_pcoupl(fplog, step, inputrec, dtpc, pressure, state->box, forceVirial,
 +                                constraintVirial, pressureCouplingMu, &state->baros_integral);
 +                crescale_pscale(inputrec, pressureCouplingMu, state->box, state->box_rel, start,
 +                                homenr, state->x.rvec_array(), state->v.rvec_array(), md->cFREEZE,
 +                                nrnb, scaleCoordinates);
 +            }
 +            break;
 +        case (epcPARRINELLORAHMAN):
 +            if (do_per_step(step + inputrec->nstpcouple - 1, inputrec->nstpcouple))
 +            {
 +                /* The box velocities were updated in do_pr_pcoupl,
 +                 * but we dont change the box vectors until we get here
 +                 * since we need to be able to shift/unshift above.
 +                 */
 +                real dtpc = inputrec->nstpcouple * dt;
 +                for (int i = 0; i < DIM; i++)
 +                {
 +                    for (int m = 0; m <= i; m++)
 +                    {
 +                        state->box[i][m] += dtpc * state->boxv[i][m];
 +                    }
 +                }
 +                preserve_box_shape(inputrec, state->box_rel, state->box);
 +
 +                /* Scale the coordinates */
 +                if (scaleCoordinates)
 +                {
 +                    auto x = state->x.rvec_array();
 +                    for (int n = start; n < start + homenr; n++)
 +                    {
 +                        tmvmul_ur0(pressureCouplingMu, x[n], x[n]);
 +                    }
 +                }
 +            }
 +            break;
 +        case (epcMTTK):
 +            switch (inputrec->epct)
 +            {
 +                case (epctISOTROPIC):
 +                    /* DIM * eta = ln V.  so DIM*eta_new = DIM*eta_old + DIM*dt*veta =>
 +                       ln V_new = ln V_old + 3*dt*veta => V_new = V_old*exp(3*dt*veta) =>
 +                       Side length scales as exp(veta*dt) */
 +
 +                    msmul(state->box, std::exp(state->veta * dt), state->box);
 +
 +                    /* Relate veta to boxv.  veta = d(eta)/dT = (1/DIM)*1/V dV/dT.
 +                       o               If we assume isotropic scaling, and box length scaling
 +                       factor L, then V = L^DIM (det(M)).  So dV/dt = DIM
 +                       L^(DIM-1) dL/dt det(M), and veta = (1/L) dL/dt.  The
 +                       determinant of B is L^DIM det(M), and the determinant
 +                       of dB/dt is (dL/dT)^DIM det (M).  veta will be
 +                       (det(dB/dT)/det(B))^(1/3).  Then since M =
 +                       B_new*(vol_new)^(1/3), dB/dT_new = (veta_new)*B(new). */
 +
 +                    msmul(state->box, state->veta, state->boxv);
 +                    break;
 +                default: break;
 +            }
 +            break;
 +        default: break;
 +    }
 +
 +    if (boxDeformation)
 +    {
 +        auto localX = makeArrayRef(state->x).subArray(start, homenr);
 +        boxDeformation->apply(localX, state->box, step);
 +    }
 +}
 +
 +extern gmx_bool update_randomize_velocities(const t_inputrec*        ir,
 +                                            int64_t                  step,
 +                                            const t_commrec*         cr,
 +                                            const t_mdatoms*         md,
 +                                            gmx::ArrayRef<gmx::RVec> v,
 +                                            const gmx::Update*       upd,
 +                                            const gmx::Constraints*  constr)
 +{
 +
 +    real rate = (ir->delta_t) / ir->opts.tau_t[0];
 +
 +    if (ir->etc == etcANDERSEN && constr != nullptr)
 +    {
 +        /* Currently, Andersen thermostat does not support constrained
 +           systems. Functionality exists in the andersen_tcoupl
 +           function in GROMACS 4.5.7 to allow this combination. That
 +           code could be ported to the current random-number
 +           generation approach, but has not yet been done because of
 +           lack of time and resources. */
 +        gmx_fatal(FARGS,
 +                  "Normal Andersen is currently not supported with constraints, use massive "
 +                  "Andersen instead");
 +    }
 +
 +    /* proceed with andersen if 1) it's fixed probability per
 +       particle andersen or 2) it's massive andersen and it's tau_t/dt */
 +    if ((ir->etc == etcANDERSEN) || do_per_step(step, gmx::roundToInt(1.0 / rate)))
 +    {
 +        andersen_tcoupl(ir, step, cr, md, v, rate, upd->getAndersenRandomizeGroup(),
 +                        upd->getBoltzmanFactor());
 +        return TRUE;
 +    }
 +    return FALSE;
 +}
 +
  /*
     static const double sy_const[MAX_SUZUKI_YOSHIDA_NUM+1][MAX_SUZUKI_YOSHIDA_NUM+1] = {
      {},
@@@ -525,7 -288,7 +525,7 @@@ static void boxv_trotter(const t_inputr
      /* for now, we use Elr = 0, because if you want to get it right, you
         really should be using PME. Maybe print a warning? */
  
 -    pscal = calc_pres(ir->ePBC, nwall, box, ekinmod, vir, localpres) + pcorr;
 +    pscal = calc_pres(ir->pbcType, nwall, box, ekinmod, vir, localpres) + pcorr;
  
      vol = det(box);
      GW  = (vol * (MassQ->Winv / PRESFAC)) * (DIM * pscal - trace(ir->ref_p)); /* W is in ps^2 * bar * nm^3 */
   *
   */
  
 -real calc_pres(int ePBC, int nwall, const matrix box, const tensor ekin, const tensor vir, tensor pres)
 +real calc_pres(PbcType pbcType, int nwall, const matrix box, const tensor ekin, const tensor vir, tensor pres)
  {
      int  n, m;
      real fac;
  
 -    if (ePBC == epbcNONE || (ePBC == epbcXY && nwall != 2))
 +    if (pbcType == PbcType::No || (pbcType == PbcType::XY && nwall != 2))
      {
          clear_mat(pres);
      }
@@@ -811,6 -574,7 +811,6 @@@ void berendsen_pcoupl(FILE
                        matrix            mu,
                        double*           baros_integral)
  {
 -    int  d, n;
      real scalar_pressure, xy_pressure, p_corr_z;
      char buf[STRLEN];
  
       */
      scalar_pressure = 0;
      xy_pressure     = 0;
 -    for (d = 0; d < DIM; d++)
 +    for (int d = 0; d < DIM; d++)
      {
          scalar_pressure += pres[d][d] / DIM;
          if (d != ZZ)
      switch (ir->epct)
      {
          case epctISOTROPIC:
 -            for (d = 0; d < DIM; d++)
 +            for (int d = 0; d < DIM; d++)
              {
                  mu[d][d] = 1.0 - factor(d, d) * (ir->ref_p[d][d] - scalar_pressure) / DIM;
              }
              break;
          case epctSEMIISOTROPIC:
 -            for (d = 0; d < ZZ; d++)
 +            for (int d = 0; d < ZZ; d++)
              {
                  mu[d][d] = 1.0 - factor(d, d) * (ir->ref_p[d][d] - xy_pressure) / DIM;
              }
              mu[ZZ][ZZ] = 1.0 - factor(ZZ, ZZ) * (ir->ref_p[ZZ][ZZ] - pres[ZZ][ZZ]) / DIM;
              break;
          case epctANISOTROPIC:
 -            for (d = 0; d < DIM; d++)
 +            for (int d = 0; d < DIM; d++)
              {
 -                for (n = 0; n < DIM; n++)
 +                for (int n = 0; n < DIM; n++)
                  {
                      mu[d][n] = (d == n ? 1.0 : 0.0) - factor(d, n) * (ir->ref_p[d][n] - pres[d][n]) / DIM;
                  }
                  p_corr_z = 0;
              }
              mu[ZZ][ZZ] = 1.0 - ir->compress[ZZ][ZZ] * p_corr_z;
 -            for (d = 0; d < DIM - 1; d++)
 +            for (int d = 0; d < DIM - 1; d++)
              {
                  mu[d][d] = 1.0
                             + factor(d, d)
      }
  }
  
 +void crescale_pcoupl(FILE*             fplog,
 +                     int64_t           step,
 +                     const t_inputrec* ir,
 +                     real              dt,
 +                     const tensor      pres,
 +                     const matrix      box,
 +                     const matrix      force_vir,
 +                     const matrix      constraint_vir,
 +                     matrix            mu,
 +                     double*           baros_integral)
 +{
 +    /*
 +     *  Calculate the scaling matrix mu
 +     */
 +    real scalar_pressure = 0;
 +    real xy_pressure     = 0;
 +    for (int d = 0; d < DIM; d++)
 +    {
 +        scalar_pressure += pres[d][d] / DIM;
 +        if (d != ZZ)
 +        {
 +            xy_pressure += pres[d][d] / (DIM - 1);
 +        }
 +    }
 +    clear_mat(mu);
 +
 +    gmx::ThreeFry2x64<64>         rng(ir->ld_seed, gmx::RandomDomain::Barostat);
 +    gmx::NormalDistribution<real> normalDist;
 +    rng.restart(step, 0);
 +    real vol = 1.0;
 +    for (int d = 0; d < DIM; d++)
 +    {
 +        vol *= box[d][d];
 +    }
 +    real gauss;
 +    real gauss2;
 +    real kt = ir->opts.ref_t[0] * BOLTZ;
 +    if (kt < 0.0)
 +    {
 +        kt = 0.0;
 +    }
 +
 +    switch (ir->epct)
 +    {
 +        case epctISOTROPIC:
 +            gauss = normalDist(rng);
 +            for (int d = 0; d < DIM; d++)
 +            {
 +                const real compressibilityFactor = ir->compress[d][d] * dt / ir->tau_p;
 +                mu[d][d] = std::exp(-compressibilityFactor * (ir->ref_p[d][d] - scalar_pressure) / DIM
 +                                    + std::sqrt(2.0 * kt * compressibilityFactor * PRESFAC / vol)
 +                                              * gauss / DIM);
 +            }
 +            break;
 +        case epctSEMIISOTROPIC:
 +            gauss  = normalDist(rng);
 +            gauss2 = normalDist(rng);
 +            for (int d = 0; d < ZZ; d++)
 +            {
 +                const real compressibilityFactor = ir->compress[d][d] * dt / ir->tau_p;
 +                mu[d][d]                         = std::exp(
 +                        -compressibilityFactor * (ir->ref_p[d][d] - xy_pressure) / DIM
 +                        + std::sqrt((DIM - 1) * 2.0 * kt * compressibilityFactor * PRESFAC / vol / DIM)
 +                                  / (DIM - 1) * gauss);
 +            }
 +            {
 +                const real compressibilityFactor = ir->compress[ZZ][ZZ] * dt / ir->tau_p;
 +                mu[ZZ][ZZ]                       = std::exp(
 +                        -compressibilityFactor * (ir->ref_p[ZZ][ZZ] - pres[ZZ][ZZ]) / DIM
 +                        + std::sqrt(2.0 * kt * compressibilityFactor * PRESFAC / vol / DIM) * gauss2);
 +            }
 +            break;
 +        case epctSURFACETENSION:
 +            gauss  = normalDist(rng);
 +            gauss2 = normalDist(rng);
 +            for (int d = 0; d < ZZ; d++)
 +            {
 +                const real compressibilityFactor = ir->compress[d][d] * dt / ir->tau_p;
 +                /* Notice: we here use ref_p[ZZ][ZZ] as isotropic pressure and ir->ref_p[d][d] as surface tension */
 +                mu[d][d] = std::exp(
 +                        -compressibilityFactor
 +                                * (ir->ref_p[ZZ][ZZ] - ir->ref_p[d][d] / box[ZZ][ZZ] - xy_pressure) / DIM
 +                        + std::sqrt(4.0 / 3.0 * kt * compressibilityFactor * PRESFAC / vol)
 +                                  / (DIM - 1) * gauss);
 +            }
 +            {
 +                const real compressibilityFactor = ir->compress[ZZ][ZZ] * dt / ir->tau_p;
 +                mu[ZZ][ZZ]                       = std::exp(
 +                        -compressibilityFactor * (ir->ref_p[ZZ][ZZ] - pres[ZZ][ZZ]) / DIM
 +                        + std::sqrt(2.0 / 3.0 * kt * compressibilityFactor * PRESFAC / vol) * gauss2);
 +            }
 +            break;
 +        default:
 +            gmx_fatal(FARGS, "C-rescale pressure coupling type %s not supported yet\n",
 +                      EPCOUPLTYPETYPE(ir->epct));
 +    }
 +    /* To fullfill the orientation restrictions on triclinic boxes
 +     * we will set mu_yx, mu_zx and mu_zy to 0 and correct
 +     * the other elements of mu to first order.
 +     */
 +    mu[YY][XX] += mu[XX][YY];
 +    mu[ZZ][XX] += mu[XX][ZZ];
 +    mu[ZZ][YY] += mu[YY][ZZ];
 +    mu[XX][YY] = 0;
 +    mu[XX][ZZ] = 0;
 +    mu[YY][ZZ] = 0;
 +
 +    /* Keep track of the work the barostat applies on the system.
 +     * Without constraints force_vir tells us how Epot changes when scaling.
 +     * With constraints constraint_vir gives us the constraint contribution
 +     * to both Epot and Ekin. Although we are not scaling velocities, scaling
 +     * the coordinates leads to scaling of distances involved in constraints.
 +     * This in turn changes the angular momentum (even if the constrained
 +     * distances are corrected at the next step). The kinetic component
 +     * of the constraint virial captures the angular momentum change.
 +     */
 +    for (int d = 0; d < DIM; d++)
 +    {
 +        for (int n = 0; n <= d; n++)
 +        {
 +            *baros_integral -=
 +                    2 * (mu[d][n] - (n == d ? 1 : 0)) * (force_vir[d][n] + constraint_vir[d][n]);
 +        }
 +    }
 +
 +    if (debug)
 +    {
 +        pr_rvecs(debug, 0, "PC: pres ", pres, 3);
 +        pr_rvecs(debug, 0, "PC: mu   ", mu, 3);
 +    }
 +
 +    if (mu[XX][XX] < 0.99 || mu[XX][XX] > 1.01 || mu[YY][YY] < 0.99 || mu[YY][YY] > 1.01
 +        || mu[ZZ][ZZ] < 0.99 || mu[ZZ][ZZ] > 1.01)
 +    {
 +        char buf[STRLEN];
 +        char buf2[22];
 +        sprintf(buf,
 +                "\nStep %s  Warning: pressure scaling more than 1%%, "
 +                "mu: %g %g %g\n",
 +                gmx_step_str(step, buf2), mu[XX][XX], mu[YY][YY], mu[ZZ][ZZ]);
 +        if (fplog)
 +        {
 +            fprintf(fplog, "%s", buf);
 +        }
 +        fprintf(stderr, "%s", buf);
 +    }
 +}
 +
 +void crescale_pscale(const t_inputrec*    ir,
 +                     const matrix         mu,
 +                     matrix               box,
 +                     matrix               box_rel,
 +                     int                  start,
 +                     int                  nr_atoms,
 +                     rvec                 x[],
 +                     rvec                 v[],
 +                     const unsigned short cFREEZE[],
 +                     t_nrnb*              nrnb,
 +                     const bool           scaleCoordinates)
 +{
 +    ivec* nFreeze = ir->opts.nFreeze;
 +    int nthreads gmx_unused;
 +    matrix       inv_mu;
 +
 +#ifndef __clang_analyzer__
 +    nthreads = gmx_omp_nthreads_get(emntUpdate);
 +#endif
 +
 +    gmx::invertBoxMatrix(mu, inv_mu);
 +
 +    /* Scale the positions and the velocities */
 +    if (scaleCoordinates)
 +    {
 +#pragma omp parallel for num_threads(nthreads) schedule(static)
 +        for (int n = start; n < start + nr_atoms; n++)
 +        {
 +            // Trivial OpenMP region that does not throw
 +            int g;
 +
 +            if (cFREEZE == nullptr)
 +            {
 +                g = 0;
 +            }
 +            else
 +            {
 +                g = cFREEZE[n];
 +            }
 +
 +            if (!nFreeze[g][XX])
 +            {
 +                x[n][XX] = mu[XX][XX] * x[n][XX] + mu[YY][XX] * x[n][YY] + mu[ZZ][XX] * x[n][ZZ];
 +                v[n][XX] = inv_mu[XX][XX] * v[n][XX] + inv_mu[YY][XX] * v[n][YY]
 +                           + inv_mu[ZZ][XX] * v[n][ZZ];
 +            }
 +            if (!nFreeze[g][YY])
 +            {
 +                x[n][YY] = mu[YY][YY] * x[n][YY] + mu[ZZ][YY] * x[n][ZZ];
 +                v[n][YY] = inv_mu[YY][YY] * v[n][YY] + inv_mu[ZZ][YY] * v[n][ZZ];
 +            }
 +            if (!nFreeze[g][ZZ])
 +            {
 +                x[n][ZZ] = mu[ZZ][ZZ] * x[n][ZZ];
 +                v[n][ZZ] = inv_mu[ZZ][ZZ] * v[n][ZZ];
 +            }
 +        }
 +    }
 +    /* compute final boxlengths */
 +    for (int d = 0; d < DIM; d++)
 +    {
 +        box[d][XX] = mu[XX][XX] * box[d][XX] + mu[YY][XX] * box[d][YY] + mu[ZZ][XX] * box[d][ZZ];
 +        box[d][YY] = mu[YY][YY] * box[d][YY] + mu[ZZ][YY] * box[d][ZZ];
 +        box[d][ZZ] = mu[ZZ][ZZ] * box[d][ZZ];
 +    }
 +
 +    preserve_box_shape(ir, box_rel, box);
 +
 +    /* (un)shifting should NOT be done after this,
 +     * since the box vectors might have changed
 +     */
 +    inc_nrnb(nrnb, eNR_PCOUPL, nr_atoms);
 +}
 +
  void berendsen_pscale(const t_inputrec*    ir,
                        const matrix         mu,
                        matrix               box,
@@@ -1770,7 -1312,7 +1770,7 @@@ static real energyNoseHoover(const t_in
  
          if (nd > 0.0)
          {
-             if (inputrecNvtTrotter(ir))
+             if (inputrecNvtTrotter(ir) || inputrecNptTrotter(ir))
              {
                  /* contribution from the thermal momenta of the NH chain */
                  for (int j = 0; j < nh; j++)
@@@ -1821,7 -1363,7 +1821,7 @@@ static real energyPressureMTTK(const t_
              double iQinv = MassQ->QPinv[i * nh + j];
              if (iQinv > 0)
              {
-                 energy += 0.5 * gmx::square(state->nhpres_vxi[i * nh + j] / iQinv);
+                 energy += 0.5 * gmx::square(state->nhpres_vxi[i * nh + j]) / iQinv;
                  /* contribution from the thermal variable of the NH chain */
                  energy += state->nhpres_xi[i * nh + j] * kT;
              }
@@@ -1899,8 -1441,7 +1899,8 @@@ real NPT_energy(const t_inputrec* ir, c
                      energyNPT += energyPressureMTTK(ir, state, MassQ);
                  }
                  break;
 -            case epcBERENDSEN: energyNPT += state->baros_integral; break;
 +            case epcBERENDSEN:
 +            case epcCRESCALE: energyNPT += state->baros_integral; break;
              default:
                  GMX_RELEASE_ASSERT(
                          false,
@@@ -2201,7 -1742,7 +2201,7 @@@ void update_annealing_target_temp(t_inp
          }
      }
  
 -    update_temperature_constants(upd->sd(), ir);
 +    upd->update_temperature_constants(*ir);
  }
  
  void pleaseCiteCouplingAlgorithms(FILE* fplog, const t_inputrec& ir)
          {
              please_cite(fplog, "Bussi2007a");
          }
 +        if (ir.epc == epcCRESCALE)
 +        {
 +            please_cite(fplog, "Bernetti2020");
 +        }
          // TODO this is actually an integrator, not a coupling algorithm
          if (ir.eI == eiSD1)
          {
index 394c225ddeb819aa704156142ed1c2f76e1a8d2d,4589a85683fa4177f7c215ad53196b6ee0379128..bed78d90c1699df9169929577556099483e8e246
@@@ -730,7 -730,7 +730,7 @@@ void mde_delta_h_coll_update_energyhist
      {
          std::vector<real>& dh = deltaH->dh[i];
          dh.resize(dhc->dh[i].ndh);
-         std::copy(dh.begin(), dh.end(), dhc->dh[i].dh);
+         std::copy(dhc->dh[i].dh, dhc->dh[i].dh + dhc->dh[i].ndh, dh.begin());
      }
      deltaH->start_time   = dhc->start_time;
      deltaH->start_lambda = dhc->start_lambda;
@@@ -759,5 -759,12 +759,5 @@@ void mde_delta_h_coll_restore_energyhis
      {
          dhc->start_lambda = deltaH->start_lambda;
      }
 -    if (dhc->dh[0].ndh > 0)
 -    {
 -        dhc->start_time_set = TRUE;
 -    }
 -    else
 -    {
 -        dhc->start_time_set = FALSE;
 -    }
 +    dhc->start_time_set = (dhc->dh[0].ndh > 0);
  }
index 9b2db3b779d1c17592746a333d58f396dfc13693,a54bff29bbed91811711587af795d3195ea34fcc..e38f3f7dbc7941fb8c283b75a18ee71251dffff0
@@@ -3,8 -3,7 +3,8 @@@
   *
   * Copyright (c) 1991-2000, University of Groningen, The Netherlands.
   * Copyright (c) 2001-2004, The GROMACS development team.
 - * Copyright (c) 2013,2014,2015,2016,2017,2018,2019,2020, by the GROMACS development team, led by
 + * Copyright (c) 2013,2014,2015,2016,2017 by the GROMACS development team.
 + * Copyright (c) 2018,2019,2020, by the GROMACS development team, led by
   * Mark Abraham, David van der Spoel, Berk Hess, and Erik Lindahl,
   * and including many others, as listed in the AUTHORS file in the
   * top-level source directory and at http://www.gromacs.org.
@@@ -72,10 -71,7 +72,10 @@@ struct t_inputrec
  namespace gmx
  {
  struct AwhHistory;
 -}
 +enum class CheckpointDataOperation;
 +template<CheckpointDataOperation operation>
 +class CheckpointData;
 +} // namespace gmx
  
  //! Convenience alias for until all is moved in the gmx namespace
  template<class T>
@@@ -179,14 -175,6 +179,14 @@@ public
       * before we enter the MD loop should compute these quantities
       * fresh, or not. */
      bool hasReadEkinState;
 +
 +    /*!
 +     * \brief Allows to read and write checkpoint within modular simulator
 +     * \tparam operation  Whether we're reading or writing
 +     * \param checkpointData  The CheckpointData object
 +     */
 +    template<gmx::CheckpointDataOperation operation>
 +    void doCheckpoint(gmx::CheckpointData<operation> checkpointData);
  };
  
  /*! \brief Free-energy sampling history struct
@@@ -234,8 -222,8 +234,8 @@@ public
      // All things public
      int natoms; //!< Number of atoms, local + non-local; this is the size of \p x, \p v and \p cg_p, when used
      int ngtc;          //!< The number of temperature coupling groups
-     int nnhpres;       //!< The NH-chain length for the MTTK barostat
-     int nhchainlength; //!< The NH-chain length for temperature coupling
+     int nnhpres;       //!< The number of NH-chains for the MTTK barostat (always 1 or 0)
+     int nhchainlength; //!< The NH-chain length for temperature coupling and MTTK barostat
      int flags; //!< Set of bit-flags telling which entries are present, see enum at the top of the file
      int                      fep_state;      //!< indicates which of the alchemical states we are in
      std::array<real, efptNR> lambda;         //!< Free-energy lambda vector
@@@ -352,21 -340,19 +352,21 @@@ static inline gmx::ArrayRef<const gmx::
      }
  };
  
 -/*! \brief Fills fep_state, lambda, and lam0 if needed
 +/*! \brief Prints the current lambda state to the log file.
   *
 - * If FEP or simulated tempering is in use:
 + * \param[in] fplog  The log file. If fplog == nullptr there will be no output.
 + * \param[in] lambda The array of lambda values.
 + * \param[in] isInitialOutput Whether this output is the initial lambda state or not.
 + */
 +void printLambdaStateToLog(FILE* fplog, gmx::ArrayRef<real> lambda, bool isInitialOutput);
 +
 +
 +/*! \brief Fills fep_state and lambda if needed
   *
 - *    fills non-null lam0 with the initial lambda values, and
 - *    on master rank fills fep_state and lambda.
 + * If FEP or simulated tempering is in use,  fills fep_state
 + * and lambda on master rank.
   *
   * Reports the initial lambda state to the log file. */
 -void initialize_lambdas(FILE*               fplog,
 -                        const t_inputrec&   ir,
 -                        bool                isMaster,
 -                        int*                fep_state,
 -                        gmx::ArrayRef<real> lambda,
 -                        double*             lam0);
 +void initialize_lambdas(FILE* fplog, const t_inputrec& ir, bool isMaster, int* fep_state, gmx::ArrayRef<real> lambda);
  
  #endif
index f1d470779a7d743d5a5f15617fbb7bba056cbf71,0f28b29a56ac6c177a62cf13b5b3e174ebd38851..e5e27165d34119a584ffb699dbc9576e492b467a
@@@ -78,28 -78,25 +78,28 @@@ namespac
  //! \{
  
  //! Enum value to store the selected value for `-type`.
 -enum DistanceType
 +enum class DistanceType : int
  {
 -    eDistanceType_Min,
 -    eDistanceType_Max
 +    Min,
 +    Max,
 +    Count
  };
  
  //! Enum value to store the selected value for `-refgrouping`/`-selgrouping`.
 -enum GroupType
 +enum class GroupType : int
  {
 -    eGroupType_All,
 -    eGroupType_Residue,
 -    eGroupType_Molecule,
 -    eGroupType_None
 +    All,
 +    Residue,
 +    Molecule,
 +    None,
 +    Count
  };
  
  //! Strings corresponding to DistanceType.
 -const char* const c_distanceTypes[] = { "min", "max" };
 +const EnumerationArray<DistanceType, const char*> c_distanceTypeNames = { { "min", "max" } };
  //! Strings corresponding to GroupType.
 -const char* const c_groupTypes[] = { "all", "res", "mol", "none" };
 +const EnumerationArray<GroupType, const char*> c_groupTypeNames = { { "all", "res", "mol",
 +                                                                      "none" } };
  
  /*! \brief
   * Implements `gmx pairdist` trajectory analysis module.
@@@ -168,9 -165,9 +168,9 @@@ private
  
  PairDistance::PairDistance() :
      cutoff_(0.0),
 -    distanceType_(eDistanceType_Min),
 -    refGroupType_(eGroupType_All),
 -    selGroupType_(eGroupType_All),
 +    distanceType_(DistanceType::Min),
 +    refGroupType_(GroupType::All),
 +    selGroupType_(GroupType::All),
      refGroupCount_(0),
      maxGroupCount_(0),
      initialDist2_(0.0),
@@@ -231,17 -228,17 +231,17 @@@ void PairDistance::initOptions(IOptions
              DoubleOption("cutoff").store(&cutoff_).description("Maximum distance to consider"));
      options->addOption(EnumOption<DistanceType>("type")
                                 .store(&distanceType_)
 -                               .enumValue(c_distanceTypes)
 +                               .enumValue(c_distanceTypeNames)
                                 .description("Type of distances to calculate"));
      options->addOption(
              EnumOption<GroupType>("refgrouping")
                      .store(&refGroupType_)
 -                    .enumValue(c_groupTypes)
 +                    .enumValue(c_groupTypeNames)
                      .description("Grouping of -ref positions to compute the min/max over"));
      options->addOption(
              EnumOption<GroupType>("selgrouping")
                      .store(&selGroupType_)
 -                    .enumValue(c_groupTypes)
 +                    .enumValue(c_groupTypeNames)
                      .description("Grouping of -sel positions to compute the min/max over"));
  
      options->addOption(SelectionOption("ref").store(&refSel_).required().description(
  }
  
  //! Helper function to initialize the grouping for a selection.
 -int initSelectionGroups(Selection* sel, const gmx_mtop_t* top, int type)
 +int initSelectionGroups(Selection* sel, const gmx_mtop_t* top, GroupType type)
  {
      e_index_t indexType = INDEX_UNKNOWN;
-     switch (type)
+     // If the selection type is INDEX_UNKNOWN (e.g. a position not associated
+     // with a set of particles), don't overwrite the selection type.
+     if (sel->type() != INDEX_UNKNOWN)
      {
-         case GroupType::All: indexType = INDEX_ALL; break;
-         case GroupType::Residue: indexType = INDEX_RES; break;
-         case GroupType::Molecule: indexType = INDEX_MOL; break;
-         case GroupType::None: indexType = INDEX_ATOM; break;
-         case GroupType::Count: GMX_THROW(InternalError("Invalid GroupType"));
+         switch (type)
+         {
 -            case eGroupType_All: indexType = INDEX_ALL; break;
 -            case eGroupType_Residue: indexType = INDEX_RES; break;
 -            case eGroupType_Molecule: indexType = INDEX_MOL; break;
 -            case eGroupType_None: indexType = INDEX_ATOM; break;
++            case GroupType::All: indexType = INDEX_ALL; break;
++            case GroupType::Residue: indexType = INDEX_RES; break;
++            case GroupType::Molecule: indexType = INDEX_MOL; break;
++            case GroupType::None: indexType = INDEX_ATOM; break;
++            case GroupType::Count: GMX_THROW(InternalError("Invalid GroupType"));
+         }
      }
      return sel->initOriginalIdsToGroup(top, indexType);
  }
@@@ -284,7 -285,7 +289,7 @@@ void PairDistance::initAnalysis(const T
      {
          AnalysisDataPlotModulePointer plotm(new AnalysisDataPlotModule(settings.plotSettings()));
          plotm->setFileName(fnDist_);
 -        if (distanceType_ == eDistanceType_Max)
 +        if (distanceType_ == DistanceType::Max)
          {
              plotm->setTitle("Maximum distance");
          }
      {
          initialDist2_ = cutoff_ * cutoff_;
      }
 -    if (distanceType_ == eDistanceType_Max)
 +    if (distanceType_ == DistanceType::Max)
      {
          initialDist2_ = 0.0;
      }
@@@ -410,8 -411,8 +415,8 @@@ TrajectoryAnalysisModuleDataPointer Pai
  void PairDistance::analyzeFrame(int frnr, const t_trxframe& fr, t_pbc* pbc, TrajectoryAnalysisModuleData* pdata)
  {
      AnalysisDataHandle      dh         = pdata->dataHandle(distances_);
 -    const Selection&        refSel     = pdata->parallelSelection(refSel_);
 -    const SelectionList&    sel        = pdata->parallelSelections(sel_);
 +    const Selection&        refSel     = TrajectoryAnalysisModuleData::parallelSelection(refSel_);
 +    const SelectionList&    sel        = TrajectoryAnalysisModuleData::parallelSelections(sel_);
      PairDistanceModuleData& frameData  = *static_cast<PairDistanceModuleData*>(pdata);
      std::vector<real>&      distArray  = frameData.distArray_;
      std::vector<int>&       countArray = frameData.countArray_;
              const int                selIndex = selPos.mappedId();
              const int                index    = selIndex * refGroupCount_ + refIndex;
              const real               r2       = pair.distance2();
 -            if (distanceType_ == eDistanceType_Min)
 +            if (distanceType_ == DistanceType::Min)
              {
                  if (distArray[index] > r2)
                  {
                      // update the distance if necessary and the count.
                      if (countArray[index] < totalCount)
                      {
 -                        if (distanceType_ == eDistanceType_Max)
 +                        if (distanceType_ == DistanceType::Max)
                          {
                              distArray[index] = cutoff2_;
                          }