Merge branch release-2016
authorMark Abraham <mark.j.abraham@gmail.com>
Tue, 14 Mar 2017 14:40:05 +0000 (15:40 +0100)
committerMark Abraham <mark.j.abraham@gmail.com>
Tue, 14 Mar 2017 14:40:05 +0000 (15:40 +0100)
Change-Id: I14f7cd22447b8a76110d8b41487666eb49835835

14 files changed:
1  2 
cmake/gmxManageSimd.cmake
docs/install-guide/index.rst
src/gromacs/gmxana/gmx_do_dssp.cpp
src/gromacs/gmxana/gmx_mindist.cpp
src/gromacs/gmxana/gmx_traj.cpp
src/gromacs/gmxana/gmx_trjconv.cpp
src/gromacs/gmxana/gmx_tune_pme.cpp
src/gromacs/gmxpreprocess/readir.cpp
src/gromacs/gmxpreprocess/solvate.cpp
src/gromacs/listed-forces/pairs.cpp
src/gromacs/mdlib/sim_util.cpp
src/gromacs/utility/binaryinformation.cpp
src/programs/mdrun/md.cpp
src/testutils/CMakeLists.txt

index b862a71cd61344571c3dc699903248cc19164bb6,01150d19be31b069112ac4170629227438564883..85249369cf682805ce64ed2faec38d8b87a76937
@@@ -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, led by
+ # Copyright (c) 2012,2013,2014,2015,2016,2017, 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.
@@@ -32,6 -32,8 +32,6 @@@
  # To help us fund GROMACS development, we humbly ask that you cite
  # the research papers on the package. Check out http://www.gromacs.org.
  
 -# include avx test source, used if the AVX flags are set below
 -include(gmxTestAVXMaskload)
  include(gmxFindFlagsForSource)
  
  # Macro that manages setting the respective C and C++ toolchain
@@@ -72,6 -74,11 +72,11 @@@ macro(prepare_power_vsx_toolchain TOOLC
              message(FATAL_ERROR "Using VSX SIMD in double precision with GCC requires GCC-4.9 or later.")
          endif()
      endif()
+     if(${CMAKE_CXX_COMPILER_ID} MATCHES "XL" OR ${CMAKE_C_COMPILER_ID} MATCHES "XL")
+         if(CMAKE_CXX_COMPILER_VERSION VERSION_LESS "13.1.5" OR CMAKE_C_COMPILER_VERSION VERSION_LESS "13.1.5")
+             message(FATAL_ERROR "Using VSX SIMD requires XL compiler version 13.1.5 or later.")
+         endif()
+     endif()
  endmacro()
  
  # Issue a fatal error with an appropriate message, when the toolchain
@@@ -245,6 -252,8 +250,6 @@@ elseif(GMX_SIMD STREQUAL "AVX_128_FMA"
      set(GMX_SIMD_X86_${GMX_SIMD} 1)
      set(SIMD_STATUS_MESSAGE "Enabling 128-bit AVX SIMD GROMACS SIMD (with fused-multiply add)")
  
 -    gmx_test_avx_gcc_maskload_bug(GMX_SIMD_X86_AVX_GCC_MASKLOAD_BUG "${SIMD_C_FLAGS}")
 -
  elseif(GMX_SIMD STREQUAL "AVX_256")
  
      prepare_x86_toolchain(TOOLCHAIN_C_FLAGS TOOLCHAIN_CXX_FLAGS)
      set(GMX_SIMD_X86_${GMX_SIMD} 1)
      set(SIMD_STATUS_MESSAGE "Enabling 256-bit AVX SIMD instructions")
  
 -    gmx_test_avx_gcc_maskload_bug(GMX_SIMD_X86_AVX_GCC_MASKLOAD_BUG "${SIMD_C_FLAGS}")
 -
  elseif(GMX_SIMD STREQUAL "AVX2_256")
  
      prepare_x86_toolchain(TOOLCHAIN_C_FLAGS TOOLCHAIN_CXX_FLAGS)
      set(GMX_SIMD_X86_${GMX_SIMD} 1)
      set(SIMD_STATUS_MESSAGE "Enabling 256-bit AVX2 SIMD instructions")
  
 -    # No need to test for Maskload bug - it was fixed before gcc added AVX2 support
 -
  elseif(GMX_SIMD STREQUAL "MIC")
  
      # No flags needed. Not testing.
@@@ -378,7 -391,7 +383,7 @@@ elseif(GMX_SIMD STREQUAL "IBM_QPX"
          set(SIMD_STATUS_MESSAGE "Enabling IBM QPX SIMD instructions")
  
      else()
 -        gmx_give_fatal_error_when_simd_support_not_found("IBM QPX" "or 'cmake .. -DCMAKE_TOOLCHAIN_FILE=Platform/BlueGeneQ-static-XL-CXX' to set up the tool chain" "${SUGGEST_BINUTILS_UPDATE}")
 +        gmx_give_fatal_error_when_simd_support_not_found("IBM QPX" "or 'cmake .. -DCMAKE_TOOLCHAIN_FILE=Platform/BlueGeneQ-static-bgclang-CXX' to set up the tool chain" "${SUGGEST_BINUTILS_UPDATE}")
      endif()
  
  elseif(GMX_SIMD STREQUAL "IBM_VMX")
@@@ -410,7 -423,12 +415,12 @@@ elseif(GMX_SIMD STREQUAL "IBM_VSX"
          SIMD_${GMX_SIMD}_C_FLAGS SIMD_${GMX_SIMD}_CXX_FLAGS
          "-mvsx" "-maltivec -mabi=altivec" "-qarch=auto -qaltivec")
  
-     if(NOT SIMD_${GMX_SIMD}_C_FLAGS OR NOT SIMD_${GMX_SIMD}_CXX_FLAGS)
+     # Usually we check also for the C compiler here, but a C compiler
+     # is not required for SIMD support on this platform. cmake through
+     # at least version 3.7 cannot pass this check with the C compiler
+     # in the latest xlc 13.1.5, but the C++ compiler has different
+     # behaviour and is OK. See Redmine #2102.
+     if(NOT SIMD_${GMX_SIMD}_CXX_FLAGS)
          gmx_give_fatal_error_when_simd_support_not_found("IBM VSX" "disable SIMD support (slower)" "${SUGGEST_BINUTILS_UPDATE}")
      endif()
  
@@@ -468,11 -486,16 +478,16 @@@ endif(
  # rather than actual vector data. For now we disable __vectorcall with clang
  # when using the reference build.
  # 
+ # xlc 13.1.5 does not seem recognize any attribute, and warns about invalid ones
+ # so we avoid searching for any.
+ #
  if(NOT DEFINED GMX_SIMD_CALLING_CONVENTION)
      if(GMX_TARGET_BGQ)
          set(CALLCONV_LIST " ")
      elseif(CMAKE_CXX_COMPILER_ID MATCHES "Clang" AND GMX_SIMD STREQUAL "REFERENCE")
          set(CALLCONV_LIST __regcall " ")
+    elseif(CMAKE_CXX_COMPILER_ID MATCHES "XL")
+         set(CALLCONV_LIST " ")
      else()
          set(CALLCONV_LIST __vectorcall __regcall " ")
      endif()
index 9fc01b432e03b38570d5663cf36777408a9f0f88,eab5e872849fc70d256a7ccd21fe1288eec69675..09950835793bd507346058a0194a5fd7e608f8eb
@@@ -15,7 -15,7 +15,7 @@@ These instructions pertain to building 
  Quick and dirty installation
  ----------------------------
  1. Get the latest version of your C and C++ compilers.
 -2. Check that you have CMake version |GMX_CMAKE_MINIMUM_REQUIRED_VERSION| or later.
 +2. Check that you have CMake version |CMAKE_MINIMUM_REQUIRED_VERSION| or later.
  3. Get and unpack the latest version of the |Gromacs| tarball.
  4. Make a separate build directory and change to it. 
  5. Run ``cmake`` with the path to the source as an argument
@@@ -100,12 -100,13 +100,14 @@@ compiler. We recommend gcc, because it 
  frequently provides the best performance.
  
  You should strive to use the most recent version of your
 -compiler. Minimum supported compiler versions are
 +compiler. Since we require full C++11 support the minimum supported
 +compiler versions are
 -* GNU (gcc) 4.6
 -* Intel (icc) 14
 -* LLVM (clang) 3.4
 +* GNU (gcc) 4.8.1
 +* Intel (icc) 15.0
 +* LLVM (clang) 3.3
  * Microsoft (MSVC) 2015
  Other compilers may work (Cray, Pathscale, older clang) but do
  not offer competitive performance. We recommend against PGI because
  the performance with C++ is very bad.
@@@ -121,7 -122,7 +123,7 @@@ other compilers, read on
  
  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 4.6.1 or higher. To select a
 +the compiler to support libstc++ version 4.8.1 or higher. To select a
  particular libstdc++ library, use:
  
  * For Intel: ``-DGMX_STDLIB_CXX_FLAGS=-gcc-name=/path/to/gcc/binary``
@@@ -150,11 -151,6 +152,11 @@@ For all non-x86 platforms, your best op
  the vendor's default or recommended compiler, and check for
  specialized information below.
  
 +For updated versions of gcc to add to your Linux OS, see
 +
 +* Ubuntu: `Ubuntu toolchain ppa page`_
 +* RHEL/CentOS: `EPEL page`_ or the RedHat Developer Toolset
 +
  Compiling with parallelization options
  --------------------------------------
  
@@@ -166,10 -162,8 +168,10 @@@ generally built into your compiler and 
  GPU support
  ^^^^^^^^^^^
  |Gromacs| has excellent support for NVIDIA GPUs supported via CUDA.
 -NVIDIA's CUDA_ version |REQUIRED_CUDA_VERSION| software development kit is required,
 -and the latest version is strongly encouraged. NVIDIA GPUs with at
 +On Linux with gcc, NVIDIA's CUDA_ version |REQUIRED_CUDA_VERSION|
 +software development kit is required, and the latest
 +version is strongly encouraged. Using Intel or Microsoft compilers
 +requires version 7.0 and 8.0, respectively. NVIDIA GPUs with at
  least NVIDIA compute capability |REQUIRED_CUDA_COMPUTE_CAPABILITY| are
  required, e.g. Fermi, Kepler, Maxwell or Pascal cards. You are strongly recommended to
  get the latest CUDA version and driver supported by your hardware, but
@@@ -226,7 -220,7 +228,7 @@@ CMak
  -----
  
  |Gromacs| builds with the CMake build system, requiring at least
 -version |GMX_CMAKE_MINIMUM_REQUIRED_VERSION|. You can check whether
 +version |CMAKE_MINIMUM_REQUIRED_VERSION|. You can check whether
  CMake is installed, and what version it is, with ``cmake
  --version``. If you need to install CMake, then first check whether
  your platform's package management system provides a suitable version,
@@@ -1147,9 -1141,9 +1149,9 @@@ much everywhere, it is important that w
  it works because we have tested it. We do test on Linux, Windows, and
  Mac with a range of compilers and libraries for a range of our
  configuration options. Every commit in our git source code repository
 -is currently tested on x86 with gcc versions ranging from 4.6 through
 -5.2, and versions 16 of the Intel compiler as well as Clang
 -version 3.4 through 3.8. For this, we use a variety of GNU/Linux
 +is currently tested on x86 with a number of gcc versions ranging from 4.8.1
 +through 6.1, versions 16 of the Intel compiler, and Clang
 +versions 3.4 through 3.8. For this, we use a variety of GNU/Linux
  flavors and versions as well as recent versions of Windows. Under
  Windows, we test both MSVC 2015 and version 16 of the Intel compiler.
  For details, you can
index f10286b1421f4aca09f8271330d981829943ab06,60619c54995e36c54d7fdafa2763d6b4bbcdc11c..c7e1210b36fc0c4043aa1560d8d2ef4bb018351a
@@@ -85,7 -85,7 +85,7 @@@ static int strip_dssp(char *dsspfile, i
      {
          fgets2(buf, STRLEN, tapeout);
      }
 -    while (std::strstr(buf, "KAPPA") == NULL);
 +    while (std::strstr(buf, "KAPPA") == nullptr);
      if (bFirst)
      {
          /* Since we also have empty lines in the dssp output (temp) file,
      nresidues = 0;
      naccf     = 0;
      naccb     = 0;
 -    for (nr = 0; (fgets2(buf, STRLEN, tapeout) != NULL); nr++)
 +    for (nr = 0; (fgets2(buf, STRLEN, tapeout) != nullptr); nr++)
      {
          if (buf[13] == '!') /* Chain separator line has '!' at pos. 13 */
          {
          buf[39] = '\0';
  
          /* Only calculate solvent accessible area if needed */
 -        if ((NULL != acc) && (buf[13] != '!'))
 +        if ((nullptr != acc) && (buf[13] != '!'))
          {
              sscanf(&(buf[34]), "%d", &iacc);
              acc[nr] = iacc;
  
      if (bFirst)
      {
 -        if (0 != acc)
 +        if (nullptr != acc)
          {
              fprintf(stderr, "%d residues were classified as hydrophobic and %d as hydrophilic.\n", naccb, naccf);
          }
          {
              mat->axis_y[i] = i+1;
          }
 -        mat->axis_x = NULL;
 -        mat->matrix = NULL;
 +        mat->axis_x = nullptr;
 +        mat->matrix = nullptr;
          xsize       = 0;
          frame       = 0;
          bFirst      = FALSE;
@@@ -288,7 -288,7 +288,7 @@@ void prune_ss_legend(t_matrix *mat
      }
  
      newnmap = 0;
 -    newmap  = NULL;
 +    newmap  = nullptr;
      for (i = 0; i < mat->nmap; i++)
      {
          newnum[i] = -1;
@@@ -423,8 -423,8 +423,8 @@@ void analyse_ss(const char *outfile, t_
      }
      fprintf(fp, "\n");
  
-     /* now print percentages */
-     fprintf(fp, "%-8s %5.2f", "# SS %", total_count / static_cast<real>(mat->nx * mat->ny));
+     /* now print probabilities */
+     fprintf(fp, "%-8s %5.2f", "# SS pr.", total_count / static_cast<real>(mat->nx * mat->ny));
      for (s = 0; s < static_cast<size_t>(mat->nmap); s++)
      {
          fprintf(fp, " %5.2f", total[s] / static_cast<real>(mat->nx * mat->ny));
@@@ -505,17 -505,17 +505,17 @@@ int gmx_do_dssp(int argc, char *argv[]
      int               *index;
      rvec              *xp, *x;
      int               *average_area;
 -    real             **accr, *accr_ptr = NULL, *av_area, *norm_av_area;
 +    real             **accr, *accr_ptr = nullptr, *av_area, *norm_av_area;
      char               pdbfile[32], tmpfile[32];
      char               dssp[256];
      const char        *dptr;
      gmx_output_env_t  *oenv;
 -    gmx_rmpbc_t        gpbc = NULL;
 +    gmx_rmpbc_t        gpbc = nullptr;
  
      t_filenm           fnm[] = {
 -        { efTRX, "-f",   NULL,      ffREAD },
 -        { efTPS, NULL,   NULL,      ffREAD },
 -        { efNDX, NULL,   NULL,      ffOPTRD },
 +        { efTRX, "-f",   nullptr,      ffREAD },
 +        { efTPS, nullptr,   nullptr,      ffREAD },
 +        { efNDX, nullptr,   nullptr,      ffOPTRD },
          { efDAT, "-ssdump", "ssdump", ffOPTWR },
          { efMAP, "-map", "ss",      ffLIBRD },
          { efXPM, "-o",   "ss",      ffWRITE },
  
      if (!parse_common_args(&argc, argv,
                             PCA_CAN_TIME | PCA_CAN_VIEW | PCA_TIME_UNIT,
 -                           NFILE, fnm, asize(pa), pa, asize(desc), desc, 0, NULL, &oenv))
 +                           NFILE, fnm, asize(pa), pa, asize(desc), desc, 0, nullptr, &oenv))
      {
          return 0;
      }
      fnAArea    = opt2fn_null("-aa", NFILE, fnm);
      bDoAccSurf = (fnArea || fnTArea || fnAArea);
  
 -    read_tps_conf(ftp2fn(efTPS, NFILE, fnm), &top, &ePBC, &xp, NULL, box, FALSE);
 +    read_tps_conf(ftp2fn(efTPS, NFILE, fnm), &top, &ePBC, &xp, nullptr, box, FALSE);
      atoms = &(top.atoms);
      check_oo(atoms);
      bPhbres = bPhobics(atoms);
  
      std::strcpy(pdbfile, "ddXXXXXX");
      gmx_tmpnam(pdbfile);
 -    if ((tmpf = fopen(pdbfile, "w")) == NULL)
 +    if ((tmpf = fopen(pdbfile, "w")) == nullptr)
      {
          sprintf(pdbfile, "%ctmp%cfilterXXXXXX", DIR_SEPARATOR, DIR_SEPARATOR);
          gmx_tmpnam(pdbfile);
 -        if ((tmpf = fopen(pdbfile, "w")) == NULL)
 +        if ((tmpf = fopen(pdbfile, "w")) == nullptr)
          {
              gmx_fatal(FARGS, "Can not open tmp file %s", pdbfile);
          }
  
      std::strcpy(tmpfile, "ddXXXXXX");
      gmx_tmpnam(tmpfile);
 -    if ((tmpf = fopen(tmpfile, "w")) == NULL)
 +    if ((tmpf = fopen(tmpfile, "w")) == nullptr)
      {
          sprintf(tmpfile, "%ctmp%cfilterXXXXXX", DIR_SEPARATOR, DIR_SEPARATOR);
          gmx_tmpnam(tmpfile);
 -        if ((tmpf = fopen(tmpfile, "w")) == NULL)
 +        if ((tmpf = fopen(tmpfile, "w")) == nullptr)
          {
              gmx_fatal(FARGS, "Can not open tmp file %s", tmpfile);
          }
          fclose(tmpf);
      }
  
 -    if ((dptr = getenv("DSSP")) == NULL)
 +    if ((dptr = getenv("DSSP")) == nullptr)
      {
          dptr = "/usr/local/bin/dssp";
      }
      }
      else
      {
 -        fTArea = NULL;
 +        fTArea = nullptr;
      }
  
 -    mat.map  = NULL;
 +    mat.map  = nullptr;
      mat.nmap = readcmap(opt2fn("-map", NFILE, fnm), &(mat.map));
  
      natoms = read_first_x(oenv, &status, ftp2fn(efTRX, NFILE, fnm), &t, &x, box);
      snew(average_area, atoms->nres);
      snew(av_area, atoms->nres);
      snew(norm_av_area, atoms->nres);
 -    accr  = NULL;
 +    accr  = nullptr;
      naccr = 0;
  
      gpbc = gmx_rmpbc_init(&top.idef, ePBC, natoms);
          }
          gmx_rmpbc(gpbc, natoms, box, x);
          tapein = gmx_ffopen(pdbfile, "w");
 -        write_pdbfile_indexed(tapein, NULL, atoms, x, ePBC, box, ' ', -1, gnx, index, NULL, TRUE);
 +        write_pdbfile_indexed(tapein, nullptr, atoms, x, ePBC, box, ' ', -1, gnx, index, nullptr, TRUE);
          gmx_ffclose(tapein);
  
          if (0 != system(dssp))
index f19124b1f25d2ca07954d3e2977d551e602f6480,399256ba2bf528b0f9fb914bc14685a49f48ef9f..1229511dd733fc2feb22c3d1511828645efd3db6
@@@ -154,11 -154,11 +154,11 @@@ static void periodic_mindist_plot(cons
      int          natoms, ind_min[2] = {0, 0}, ind_mini = 0, ind_minj = 0;
      real         rmin, rmax, rmint, tmint;
      gmx_bool     bFirst;
 -    gmx_rmpbc_t  gpbc = NULL;
 +    gmx_rmpbc_t  gpbc = nullptr;
  
      natoms = read_first_x(oenv, &status, trxfn, &t, &x, box);
  
 -    check_index(NULL, n, index, NULL, natoms);
 +    check_index(nullptr, n, index, nullptr, natoms);
  
      out = xvgropen(outfn, "Minimum distance to periodic image",
                     output_env_get_time_label(oenv), "Distance (nm)", oenv);
      rmint = box[XX][XX];
      tmint = 0;
  
 -    if (NULL != top)
 +    if (nullptr != top)
      {
          gpbc = gmx_rmpbc_init(&top->idef, ePBC, natoms);
      }
      bFirst = TRUE;
      do
      {
 -        if (NULL != top)
 +        if (nullptr != top)
          {
              gmx_rmpbc(gpbc, natoms, box, x);
          }
      }
      while (read_next_x(oenv, status, &t, x, box));
  
 -    if (NULL != top)
 +    if (nullptr != top)
      {
          gmx_rmpbc_done(gpbc);
      }
@@@ -262,7 -262,7 +262,7 @@@ static void calc_dist(real rcut, gmx_bo
      for (j = 0; (j < j1); j++)
      {
          jx = index3[j];
 -        if (index2 == NULL)
 +        if (index2 == nullptr)
          {
              i0 = j + 1;
          }
                  {
                      nmin_j++;
                  }
 -                else if (r2 > rcut2)
 +                else
                  {
                      nmax_j++;
                  }
@@@ -337,7 -337,7 +337,7 @@@ void dist_plot(const char *fn, const ch
      t_trxstatus     *trxout;
      char             buf[256];
      char           **leg;
 -    real             t, dmin, dmax, **mindres = NULL, **maxdres = NULL;
 +    real             t, dmin, dmax, **mindres = nullptr, **maxdres = nullptr;
      int              nmin, nmax;
      t_trxstatus     *status;
      int              i = -1, j, k;
      rvec            *x0;
      matrix           box;
      gmx_bool         bFirst;
 -    FILE            *respertime = NULL;
 +    FILE            *respertime = nullptr;
  
      if (read_first_x(oenv, &status, fn, &t, &x0, box) == 0)
      {
      sprintf(buf, "%simum Distance", bMin ? "Min" : "Max");
      dist = xvgropen(dfile, buf, output_env_get_time_label(oenv), "Distance (nm)", oenv);
      sprintf(buf, "Number of Contacts %s %g nm", bMin ? "<" : ">", rcut);
 -    num    = nfile ? xvgropen(nfile, buf, output_env_get_time_label(oenv), "Number", oenv) : NULL;
 -    atm    = afile ? gmx_ffopen(afile, "w") : NULL;
 -    trxout = xfile ? open_trx(xfile, "w") : NULL;
 +    num    = nfile ? xvgropen(nfile, buf, output_env_get_time_label(oenv), "Number", oenv) : nullptr;
 +    atm    = afile ? gmx_ffopen(afile, "w") : nullptr;
 +    trxout = xfile ? open_trx(xfile, "w") : nullptr;
  
      if (bMat)
      {
          sprintf(buf, "%simum Distance", bMin ? "Min" : "Max");
          respertime = xvgropen(rfile, buf, output_env_get_time_label(oenv), "Distance (nm)", oenv);
          xvgr_legend(respertime, ng-1, (const char**)leg, oenv);
-         if (bPrintResName)
+         if (bPrintResName && output_env_get_print_xvgr_codes(oenv) )
          {
              fprintf(respertime, "# ");
+             for (j = 0; j < nres; j++)
+             {
+                 fprintf(respertime, "%s%d ", *(atoms->resinfo[atoms->atom[index[0][residue[j]]].resind].name), atoms->atom[index[0][residue[j]]].resind);
+             }
+             fprintf(respertime, "\n");
          }
-         for (j = 0; j < nres; j++)
-         {
-             fprintf(respertime, "%s%d ", *(atoms->resinfo[atoms->atom[index[0][residue[j]]].resind].name), atoms->atom[index[0][residue[j]]].resind);
-         }
-         fprintf(respertime, "\n");
  
      }
  
          {
              oindex[0] = bMin ? min1 : max1;
              oindex[1] = bMin ? min2 : max2;
 -            write_trx(trxout, 2, oindex, atoms, i, t, box, x0, NULL, NULL);
 +            write_trx(trxout, 2, oindex, atoms, i, t, box, x0, nullptr, nullptr);
          }
          bFirst = FALSE;
          /*dmin should be minimum distance for residue and group*/
@@@ -693,7 -694,7 +694,7 @@@ int gmx_mindist(int argc, char *argv[]
            "Write residue names" }
      };
      gmx_output_env_t *oenv;
 -    t_topology       *top  = NULL;
 +    t_topology       *top  = nullptr;
      int               ePBC = -1;
      rvec             *x;
      matrix            box;
      const char       *trxfnm, *tpsfnm, *ndxfnm, *distfnm, *numfnm, *atmfnm, *oxfnm, *resfnm;
      char            **grpname;
      int              *gnx;
 -    int             **index, *residues = NULL;
 +    int             **index, *residues = nullptr;
      t_filenm          fnm[] = {
 -        { efTRX, "-f",  NULL,      ffREAD },
 -        { efTPS,  NULL, NULL,      ffOPTRD },
 -        { efNDX,  NULL, NULL,      ffOPTRD },
 +        { efTRX, "-f",  nullptr,      ffREAD },
 +        { efTPS,  nullptr, nullptr,      ffOPTRD },
 +        { efNDX,  nullptr, nullptr,      ffOPTRD },
          { efXVG, "-od", "mindist",  ffWRITE },
          { efXVG, "-on", "numcont",  ffOPTWR },
          { efOUT, "-o", "atm-pair", ffOPTWR },
  
      if (!parse_common_args(&argc, argv,
                             PCA_CAN_VIEW | PCA_CAN_TIME | PCA_TIME_UNIT,
 -                           NFILE, fnm, asize(pa), pa, asize(desc), desc, 0, NULL, &oenv))
 +                           NFILE, fnm, asize(pa), pa, asize(desc), desc, 0, nullptr, &oenv))
      {
          return 0;
      }
      atmfnm  = ftp2fn_null(efOUT, NFILE, fnm);
      oxfnm   = opt2fn_null("-ox", NFILE, fnm);
      resfnm  = opt2fn_null("-or", NFILE, fnm);
 -    if (bPI || resfnm != NULL)
 +    if (bPI || resfnm != nullptr)
      {
          /* We need a tps file */
          tpsfnm = ftp2fn(efTPS, NFILE, fnm);
      if (tpsfnm || resfnm || !ndxfnm)
      {
          snew(top, 1);
 -        bTop = read_tps_conf(tpsfnm, top, &ePBC, &x, NULL, box, FALSE);
 +        bTop = read_tps_conf(tpsfnm, top, &ePBC, &x, nullptr, box, FALSE);
          if (bPI && !bTop)
          {
              printf("\nWARNING: Without a run input file a trajectory with broken molecules will not give the correct periodic image distance\n\n");
          }
      }
 -    get_index(top ? &(top->atoms) : NULL, ndxfnm, ng, gnx, index, grpname);
 +    get_index(top ? &(top->atoms) : nullptr, ndxfnm, ng, gnx, index, grpname);
  
      if (bMat && (ng == 1))
      {
  
      if (resfnm)
      {
 -        GMX_RELEASE_ASSERT(top != NULL, "top pointer cannot be NULL when finding residues");
 +        GMX_RELEASE_ASSERT(top != nullptr, "top pointer cannot be NULL when finding residues");
          nres = find_residues(&(top->atoms), gnx[0], index[0], &residues);
  
          if (debug)
      else
      {
          dist_plot(trxfnm, atmfnm, distfnm, numfnm, resfnm, oxfnm,
 -                  rcutoff, bMat, top ? &(top->atoms) : NULL,
 +                  rcutoff, bMat, top ? &(top->atoms) : nullptr,
                    ng, index, gnx, grpname, bSplit, !bMax, nres, residues, bPBC, ePBC,
                    bGroup, bEachResEachTime, bPrintResName, oenv);
      }
index 0d273adb5f6c731a9cb00086fe58c1f66e634121,e7025cdd2393e7aefb1063f1ef913228c6710313..445f8148be1f0f69e2f896adc17995c3b210b6b2
@@@ -72,7 -72,7 +72,7 @@@ static void low_print_data(FILE *fp, re
      fprintf(fp, " %g", time);
      for (i = 0; i < n; i++)
      {
 -        if (index != NULL)
 +        if (index != nullptr)
          {
              ii = index[i];
          }
@@@ -111,7 -111,7 +111,7 @@@ static void average_data(rvec x[], rve
          for (i = 0; i < isize[g]; i++)
          {
              ind = index[g][i];
 -            if (mass != NULL)
 +            if (mass != nullptr)
              {
                  m = mass[ind];
                  svmul(m, x[ind], tmp);
                  }
              }
          }
 -        if (mass != NULL)
 +        if (mass != nullptr)
          {
              for (d = 0; d < DIM; d++)
              {
@@@ -151,16 -151,16 +151,16 @@@ static void print_data(FILE *fp, real t
                         int ngrps, int isize[], int **index, gmx_bool bDim[],
                         const char *sffmt)
  {
 -    static rvec *xav = NULL;
 +    static rvec *xav = nullptr;
  
      if (bCom)
      {
 -        if (xav == NULL)
 +        if (xav == nullptr)
          {
              snew(xav, ngrps);
          }
          average_data(x, xav, mass, ngrps, isize, index);
 -        low_print_data(fp, time, xav, ngrps, NULL, bDim, sffmt);
 +        low_print_data(fp, time, xav, ngrps, nullptr, bDim, sffmt);
      }
      else
      {
  static void write_trx_x(t_trxstatus *status, const t_trxframe *fr, real *mass, gmx_bool bCom,
                          int ngrps, int isize[], int **index)
  {
 -    static rvec    *xav   = NULL;
 -    static t_atoms *atoms = NULL;
 +    static rvec    *xav   = nullptr;
 +    static t_atoms *atoms = nullptr;
      t_trxframe      fr_av;
      int             i;
  
      if (bCom)
      {
 -        if (xav == NULL)
 +        if (xav == nullptr)
          {
              snew(xav, ngrps);
              snew(atoms, 1);
          fr_av.natoms = ngrps;
          fr_av.atoms  = atoms;
          fr_av.x      = xav;
 -        write_trxframe(status, &fr_av, NULL);
 +        write_trxframe(status, &fr_av, nullptr);
      }
      else
      {
 -        write_trxframe_indexed(status, fr, isize[0], index[0], NULL);
 +        write_trxframe_indexed(status, fr, isize[0], index[0], nullptr);
      }
  }
  
@@@ -260,14 -260,14 +260,14 @@@ static void make_legend(FILE *fp, int n
  
  static real ekrot(rvec x[], rvec v[], real mass[], int isize, int index[])
  {
 -    static real **TCM = NULL, **L;
 +    static real **TCM = nullptr, **L;
      double        tm, m0, lxx, lxy, lxz, lyy, lyz, lzz, ekrot;
      rvec          a0, ocm;
      dvec          dx, b0;
      dvec          xcm, vcm, acm;
      int           i, j, m, n;
  
 -    if (TCM == NULL)
 +    if (TCM == nullptr)
      {
          snew(TCM, DIM);
          for (i = 0; i < DIM; i++)
@@@ -461,7 -461,7 +461,7 @@@ static void write_pdb_bfac(const char *
              svmul(scale, sum[index[i]], sum[index[i]]);
          }
  
-         fp = xvgropen(xname, title, "Atom", "", oenv);
+         fp = xvgropen(xname, title, "Atom", "Spatial component", oenv);
          for (i = 0; i < isize; i++)
          {
              fprintf(fp, "%-5d  %10.3f  %10.3f  %10.3f\n", 1+i,
                 *(atoms->resinfo[atoms->atom[maxi].resind].name),
                 atoms->resinfo[atoms->atom[maxi].resind].nr);
  
 -        if (atoms->pdbinfo == NULL)
 +        if (atoms->pdbinfo == nullptr)
          {
              snew(atoms->pdbinfo, atoms->nr);
          }
 +        atoms->havePdbInfo = TRUE;
 +
          if (onedim == -1)
          {
              for (i = 0; i < isize; i++)
                  atoms->pdbinfo[index[i]].bfac = sum[index[i]][onedim]*scale;
              }
          }
 -        write_sto_conf_indexed(fname, title, atoms, x, NULL, ePBC, box, isize, index);
 +        write_sto_conf_indexed(fname, title, atoms, x, nullptr, ePBC, box, isize, index);
      }
  }
  
@@@ -545,7 -543,7 +545,7 @@@ static void update_histo(int gnx, int i
      int  i, m, in, nnn;
      real vn, vnmax;
  
 -    if (*histo == NULL)
 +    if (*histo == nullptr)
      {
          vnmax = 0;
          for (i = 0; (i < gnx); i++)
@@@ -623,9 -621,7 +623,9 @@@ int gmx_traj(int argc, char *argv[]
          "(specified with [TT]-av[tt] or [TT]-af[tt]).[PAR]",
          "Option [TT]-vd[tt] computes a velocity distribution, i.e. the",
          "norm of the vector is plotted. In addition in the same graph",
 -        "the kinetic energy distribution is given."
 +        "the kinetic energy distribution is given.",
 +        "",
 +        "See [gmx-trajectory] for plotting similar data for selections."
      };
      static gmx_bool   bMol    = FALSE, bCom = FALSE, bPBC = TRUE, bNoJump = FALSE;
      static gmx_bool   bX      = TRUE, bY = TRUE, bZ = TRUE, bNorm = FALSE, bFP = FALSE;
          { "-scale", FALSE, etREAL, {&scale},
            "Scale factor for [REF].pdb[ref] output, 0 is autoscale" }
      };
 -    FILE             *outx   = NULL, *outv = NULL, *outf = NULL, *outb = NULL, *outt = NULL;
 -    FILE             *outekt = NULL, *outekr = NULL;
 +    FILE             *outx   = nullptr, *outv = nullptr, *outf = nullptr, *outb = nullptr, *outt = nullptr;
 +    FILE             *outekt = nullptr, *outekr = nullptr;
      t_topology        top;
      int               ePBC;
      real             *mass, time;
      const char       *indexfn;
      t_trxframe        fr;
 -    int               flags, nvhisto = 0, *vhisto = NULL;
 -    rvec             *xtop, *xp = NULL;
 -    rvec             *sumx = NULL, *sumv = NULL, *sumf = NULL;
 +    int               flags, nvhisto = 0, *vhisto = nullptr;
 +    rvec             *xtop, *xp = nullptr;
 +    rvec             *sumx = nullptr, *sumv = nullptr, *sumf = nullptr;
      matrix            topbox;
      t_trxstatus      *status;
 -    t_trxstatus      *status_out = NULL;
 -    gmx_rmpbc_t       gpbc       = NULL;
 +    t_trxstatus      *status_out = nullptr;
 +    gmx_rmpbc_t       gpbc       = nullptr;
      int               i, j;
      int               nr_xfr, nr_vfr, nr_ffr;
      char            **grpname;
      gmx_output_env_t *oenv;
  
      t_filenm          fnm[] = {
 -        { efTRX, "-f", NULL, ffREAD },
 -        { efTPS, NULL, NULL, ffREAD },
 -        { efNDX, NULL, NULL, ffOPTRD },
 +        { efTRX, "-f", nullptr, ffREAD },
 +        { efTPS, nullptr, nullptr, ffREAD },
 +        { efNDX, nullptr, nullptr, ffOPTRD },
          { efXVG, "-ox",  "coord",     ffOPTWR },
          { efTRX, "-oxt", "coord",     ffOPTWR },
          { efXVG, "-ov",  "veloc",     ffOPTWR },
  
      if (!parse_common_args(&argc, argv,
                             PCA_CAN_TIME | PCA_TIME_UNIT | PCA_CAN_VIEW,
 -                           NFILE, fnm, asize(pa), pa, asize(desc), desc, 0, NULL, &oenv))
 +                           NFILE, fnm, asize(pa), pa, asize(desc), desc, 0, nullptr, &oenv))
      {
          return 0;
      }
      sprintf(sffmt6, "%s%s%s%s%s%s", sffmt, sffmt, sffmt, sffmt, sffmt, sffmt);
  
      bTop = read_tps_conf(ftp2fn(efTPS, NFILE, fnm), &top, &ePBC,
 -                         &xtop, NULL, topbox,
 +                         &xtop, nullptr, topbox,
                           bCom && (bOX || bOXT || bOV || bOT || bEKT || bEKR));
      sfree(xtop);
      if ((bMol || bCV || bCF) && !bTop)
      }
      else
      {
 -        mass = NULL;
 +        mass = nullptr;
      }
  
      flags = 0;
          }
          if (bOF && fr.bF)
          {
 -            print_data(outf, time, fr.f, NULL, bCom, ngroups, isize, index, bDim, sffmt);
 +            print_data(outf, time, fr.f, nullptr, bCom, ngroups, isize, index, bDim, sffmt);
          }
          if (bOB && fr.bBox)
          {
      }
      while (read_next_frame(oenv, status, &fr));
  
 -    if (gpbc != NULL)
 +    if (gpbc != nullptr)
      {
          gmx_rmpbc_done(gpbc);
      }
index 92cae1eff5250383ca1415d2b0ae8013591ad5e9,0b9173470599d79cc61f89270d3cc031d777d69f..92c192830d97015a21247d12f83dc094dca83748
@@@ -49,7 -49,7 +49,7 @@@
  #include "gromacs/fileio/gmxfio.h"
  #include "gromacs/fileio/groio.h"
  #include "gromacs/fileio/pdbio.h"
 -#include "gromacs/fileio/tngio_for_tools.h"
 +#include "gromacs/fileio/tngio.h"
  #include "gromacs/fileio/tpxio.h"
  #include "gromacs/fileio/trrio.h"
  #include "gromacs/fileio/trxio.h"
@@@ -492,7 -492,7 +492,7 @@@ void do_trunc(const char *fn, real t0
  
      in   = gmx_trr_open(fn, "r");
      fp   = gmx_fio_getfp(in);
 -    if (fp == NULL)
 +    if (fp == nullptr)
      {
          fprintf(stderr, "Sorry, can not trunc %s, truncation of this filetype is not supported\n", fn);
          gmx_trr_close(in);
          bStop = FALSE;
          while (!bStop && gmx_trr_read_frame_header(in, &sh, &bOK))
          {
 -            gmx_trr_read_frame_data(in, &sh, NULL, NULL, NULL, NULL);
 +            gmx_trr_read_frame_data(in, &sh, nullptr, nullptr, nullptr, nullptr);
              fpos = gmx_ftell(fp);
              t    = sh.t;
              if (t >= t0)
@@@ -568,7 -568,7 +568,7 @@@ static gmx_mtop_t *read_mtop_for_tng(co
                                       const char *input_file,
                                       const char *output_file)
  {
 -    gmx_mtop_t *mtop = NULL;
 +    gmx_mtop_t *mtop = nullptr;
  
      if (fn2bTPX(tps_file) &&
          efTNG != fn2ftp(input_file) &&
      {
          int temp_natoms = -1;
          snew(mtop, 1);
 -        read_tpx(tps_file, NULL, NULL, &temp_natoms,
 -                 NULL, NULL, mtop);
 +        read_tpx(tps_file, nullptr, nullptr, &temp_natoms,
 +                 nullptr, nullptr, mtop);
      }
  
      return mtop;
@@@ -665,10 -665,6 +665,6 @@@ int gmx_trjconv(int argc, char *argv[]
          "   results if you in fact have a cluster. Luckily that can be checked",
          "   afterwards using a trajectory viewer. Note also that if your molecules",
          "   are broken this will not work either.",
-         "",
-         "   The separate option [TT]-clustercenter[tt] can be used to specify an",
-         "   approximate center for the cluster. This is useful e.g. if you have",
-         "   two big vesicles, and you want to maintain their relative positions.",
          " * [TT]whole[tt] only makes broken molecules whole.",
          "",
  
      };
      const char *pbc_opt[epNR + 1] =
      {
 -        NULL, "none", "mol", "res", "atom", "nojump", "cluster", "whole",
 -        NULL
 +        nullptr, "none", "mol", "res", "atom", "nojump", "cluster", "whole",
 +        nullptr
      };
  
      int         unitcell_enum;
      const char *unitcell_opt[euNR+1] =
 -    { NULL, "rect", "tric", "compact", NULL };
 +    { nullptr, "rect", "tric", "compact", nullptr };
  
      enum
      {
          ecSel, ecTric, ecRect, ecZero, ecNR
      };
      const char *center_opt[ecNR+1] =
 -    { NULL, "tric", "rect", "zero", NULL };
 +    { nullptr, "tric", "rect", "zero", nullptr };
      int         ecenter;
  
      int         fit_enum;
      };
      const char *fit[efNR + 1] =
      {
 -        NULL, "none", "rot+trans", "rotxy+transxy", "translation", "transxy",
 -        "progressive", NULL
 +        nullptr, "none", "rot+trans", "rotxy+transxy", "translation", "transxy",
 +        "progressive", nullptr
      };
  
      static gmx_bool  bSeparate     = FALSE, bVels = TRUE, bForce = FALSE, bCONECT = FALSE;
      static int       skip_nr       = 1, ndec = 3, nzero = 0;
      static real      tzero         = 0, delta_t = 0, timestep = 0, ttrunc = -1, tdump = -1, split_t = 0;
      static rvec      newbox        = {0, 0, 0}, shift = {0, 0, 0}, trans = {0, 0, 0};
 -    static char     *exec_command  = NULL;
 +    static char     *exec_command  = nullptr;
      static real      dropunder     = 0, dropover = 0;
      static gmx_bool  bRound        = FALSE;
  
      };
  #define NPA asize(pa)
  
 -    FILE             *out    = NULL;
 -    t_trxstatus      *trxout = NULL;
 +    FILE             *out    = nullptr;
 +    t_trxstatus      *trxout = nullptr;
      t_trxstatus      *trxin;
      int               file_nr;
      t_trxframe        fr, frout;
      int               flags;
 -    rvec             *xmem  = NULL, *vmem = NULL, *fmem = NULL;
 -    rvec             *xp    = NULL, x_shift, hbox;
 -    real             *w_rls = NULL;
 +    rvec             *xmem  = nullptr, *vmem = nullptr, *fmem = nullptr;
 +    rvec             *xp    = nullptr, x_shift, hbox;
 +    real             *w_rls = nullptr;
      int               m, i, d, frame, outframe, natoms, nout, ncent, newstep = 0, model_nr;
  #define SKIP 10
      t_topology        top;
 -    gmx_mtop_t       *mtop  = NULL;
 -    gmx_conect        gc    = NULL;
 +    gmx_mtop_t       *mtop  = nullptr;
 +    gmx_conect        gc    = nullptr;
      int               ePBC  = -1;
 -    t_atoms          *atoms = NULL, useatoms;
 +    t_atoms          *atoms = nullptr, useatoms;
      matrix            top_box;
      int              *index, *cindex;
      char             *grpnm;
      int               ifit, my_clust = -1;
      int              *ind_fit;
      char             *gn_fit;
 -    t_cluster_ndx    *clust           = NULL;
 -    t_trxstatus     **clust_status    = NULL;
 -    int              *clust_status_id = NULL;
 +    t_cluster_ndx    *clust           = nullptr;
 +    t_trxstatus     **clust_status    = nullptr;
 +    int              *clust_status_id = nullptr;
      int               ntrxopen        = 0;
 -    int              *nfwritten       = NULL;
 +    int              *nfwritten       = nullptr;
      int               ndrop           = 0, ncol, drop0 = 0, drop1 = 0, dropuse = 0;
      double          **dropval;
      real              tshift = 0, t0 = -1, dt = 0.001, prec;
      gmx_bool          bFit, bPFit, bReset;
      int               nfitdim;
 -    gmx_rmpbc_t       gpbc = NULL;
 +    gmx_rmpbc_t       gpbc = nullptr;
      gmx_bool          bRmPBC, bPBCWhole, bPBCcomRes, bPBCcomMol, bPBCcomAtom, bPBC, bNoJump, bCluster;
      gmx_bool          bCopy, bDoIt, bIndex, bTDump, bSetTime, bTPS = FALSE, bDTset = FALSE;
      gmx_bool          bExec, bTimeStep = FALSE, bDumpFrame = FALSE, bSetPrec, bNeedPrec;
      gmx_bool          bHaveFirstFrame, bHaveNextFrame, bSetBox, bSetUR, bSplit = FALSE;
      gmx_bool          bSubTraj = FALSE, bDropUnder = FALSE, bDropOver = FALSE, bTrans = FALSE;
      gmx_bool          bWriteFrame, bSplitHere;
 -    const char       *top_file, *in_file, *out_file = NULL;
 +    const char       *top_file, *in_file, *out_file = nullptr;
      char              out_file2[256], *charpt;
 -    char             *outf_base = NULL;
 -    const char       *outf_ext  = NULL;
 +    char             *outf_base = nullptr;
 +    const char       *outf_ext  = nullptr;
      char              top_title[256], title[256], filemode[5];
      gmx_output_env_t *oenv;
  
      t_filenm          fnm[] = {
 -        { efTRX, "-f",   NULL,      ffREAD  },
 -        { efTRO, "-o",   NULL,      ffWRITE },
 -        { efTPS, NULL,   NULL,      ffOPTRD },
 -        { efNDX, NULL,   NULL,      ffOPTRD },
 +        { efTRX, "-f",   nullptr,      ffREAD  },
 +        { efTRO, "-o",   nullptr,      ffWRITE },
 +        { efTPS, nullptr,   nullptr,      ffOPTRD },
 +        { efNDX, nullptr,   nullptr,      ffOPTRD },
          { efNDX, "-fr",  "frames",  ffOPTRD },
          { efNDX, "-sub", "cluster", ffOPTRD },
          { efXVG, "-drop", "drop",    ffOPTRD }
                             PCA_CAN_BEGIN | PCA_CAN_END | PCA_CAN_VIEW |
                             PCA_TIME_UNIT,
                             NFILE, fnm, NPA, pa, asize(desc), desc,
 -                           0, NULL, &oenv))
 +                           0, nullptr, &oenv))
      {
          return 0;
      }
          if (bSeparate || bSplit)
          {
              outf_ext = std::strrchr(out_file, '.');
 -            if (outf_ext == NULL)
 +            if (outf_ext == nullptr)
              {
                  gmx_fatal(FARGS, "Output file name '%s' does not contain a '.'", out_file);
              }
                  gmx_fatal(FARGS, "Can only use the sub option with output file types "
                            "xtc and trr");
              }
 -            clust = cluster_index(NULL, opt2fn("-sub", NFILE, fnm));
 +            clust = cluster_index(nullptr, opt2fn("-sub", NFILE, fnm));
  
              /* Check for number of files disabled, as FOPEN_MAX is not the correct
               * number to check for. In my linux box it is only 16.
              snew(nfwritten, clust->clust->nr);
              for (i = 0; (i < clust->clust->nr); i++)
              {
 -                clust_status[i]    = NULL;
 +                clust_status[i]    = nullptr;
                  clust_status_id[i] = -1;
              }
              bSeparate = bSplit = FALSE;
  
          if (bTPS)
          {
 -            read_tps_conf(top_file, &top, &ePBC, &xp, NULL, top_box,
 +            read_tps_conf(top_file, &top, &ePBC, &xp, nullptr, top_box,
                            bReset || bPBCcomRes);
              std::strncpy(top_title, *top.name, 255);
              top_title[255] = '\0';
          }
  
          /* get frame number index */
 -        frindex = NULL;
 +        frindex = nullptr;
          if (opt2bSet("-fr", NFILE, fnm))
          {
              printf("Select groups of frame number indices:\n");
                  gmx_rmpbc(gpbc, top.atoms.nr, top_box, xp);
              }
              copy_rvec(xp[index[0]], x_shift);
 -            reset_x_ndim(nfitdim, ifit, ind_fit, atoms->nr, NULL, xp, w_rls);
 +            reset_x_ndim(nfitdim, ifit, ind_fit, atoms->nr, nullptr, xp, w_rls);
              rvec_dec(x_shift, xp[index[0]]);
          }
          else
              /* get memory for stuff to go in .pdb file, and initialize
               * the pdbinfo structure part if the input has it.
               */
 -            init_t_atoms(&useatoms, atoms->nr, (atoms->pdbinfo != NULL));
 +            init_t_atoms(&useatoms, atoms->nr, atoms->havePdbInfo);
              sfree(useatoms.resinfo);
              useatoms.resinfo = atoms->resinfo;
              for (i = 0; (i < nout); i++)
              {
                  useatoms.atomname[i] = atoms->atomname[index[i]];
                  useatoms.atom[i]     = atoms->atom[index[i]];
 -                if (atoms->pdbinfo != NULL)
 +                if (atoms->havePdbInfo)
                  {
                      useatoms.pdbinfo[i]  = atoms->pdbinfo[index[i]];
                  }
                                                       filemode[0],
                                                       trxin,
                                                       &trxout,
 -                                                     NULL,
 +                                                     nullptr,
                                                       nout,
                                                       mtop,
                                                       index,
                      break;
                  case efXTC:
                  case efTRR:
 -                    out = NULL;
 +                    out = nullptr;
                      if (!bSplit && !bSubTraj)
                      {
                          trxout = open_trx(out_file, filemode);
                          gmx_rmpbc_trxfr(gpbc, &fr);
                      }
  
 -                    reset_x_ndim(nfitdim, ifit, ind_fit, natoms, NULL, fr.x, w_rls);
 +                    reset_x_ndim(nfitdim, ifit, ind_fit, natoms, nullptr, fr.x, w_rls);
                      do_fit(natoms, w_rls, xp, fr.x);
                  }
  
                  /* store this set of coordinates for future use */
                  if (bPFit || bNoJump)
                  {
 -                    if (xp == NULL)
 +                    if (xp == nullptr)
                      {
                          snew(xp, natoms);
                      }
  
                              if (bReset)
                              {
 -                                reset_x_ndim(nfitdim, ifit, ind_fit, natoms, NULL, fr.x, w_rls);
 +                                reset_x_ndim(nfitdim, ifit, ind_fit, natoms, nullptr, fr.x, w_rls);
                                  if (bFit)
                                  {
                                      do_fit_ndim(nfitdim, natoms, w_rls, xp, fr.x);
                                               clust->clust->index[my_clust]))
                                          {
                                              close_trx(clust_status[my_clust]);
 -                                            clust_status[my_clust]    = NULL;
 +                                            clust_status[my_clust]    = nullptr;
                                              clust_status_id[my_clust] = -2;
                                              ntrxopen--;
                                              if (ntrxopen < 0)
                                  {
                                      case efGRO:
                                          write_hconf_p(out, title, &useatoms,
 -                                                      frout.x, frout.bV ? frout.v : NULL, frout.box);
 +                                                      frout.x, frout.bV ? frout.v : nullptr, frout.box);
                                          break;
                                      case efPDB:
                                          fprintf(out, "REMARK    GENERATED BY TRJCONV\n");
                                              frout.bStep  = TRUE;
                                              frout.bTime  = TRUE;
                                          }
 -                                        write_g96_conf(out, &frout, -1, NULL);
 +                                        write_g96_conf(out, &frout, -1, nullptr);
                                  }
                                  if (bSeparate || bSplitHere)
                                  {
                                      gmx_ffclose(out);
 -                                    out = NULL;
 +                                    out = nullptr;
                                  }
                                  break;
                              default:
          {
              close_trx(trxout);
          }
 -        else if (out != NULL)
 +        else if (out != nullptr)
          {
              gmx_ffclose(out);
          }
  
      sfree(mtop);
  
 -    do_view(oenv, out_file, NULL);
 +    do_view(oenv, out_file, nullptr);
  
      return 0;
  }
index ed14fbc0c8f978ce2f154d524ecb401cc92a12a0,02cdca3fce04624de4085d98899454e49be20eee..e93b4298f126d99691910d60e3368feb9431b5bd
@@@ -48,9 -48,9 +48,9 @@@
  #endif
  
  #include "gromacs/commandline/pargs.h"
 +#include "gromacs/ewald/pme.h"
  #include "gromacs/fft/calcgrid.h"
  #include "gromacs/fileio/checkpoint.h"
 -#include "gromacs/fileio/readinp.h"
  #include "gromacs/fileio/tpxio.h"
  #include "gromacs/gmxana/gmx_ana.h"
  #include "gromacs/math/utilities.h"
@@@ -59,7 -59,6 +59,7 @@@
  #include "gromacs/mdtypes/commrec.h"
  #include "gromacs/mdtypes/inputrec.h"
  #include "gromacs/mdtypes/md_enums.h"
 +#include "gromacs/mdtypes/state.h"
  #include "gromacs/pbcutil/pbc.h"
  #include "gromacs/timing/walltime_accounting.h"
  #include "gromacs/topology/topology.h"
@@@ -169,7 -168,7 +169,7 @@@ static void finalize(const char *fn_out
      fp = fopen(fn_out, "r");
      fprintf(stdout, "\n\n");
  
 -    while (fgets(buf, STRLEN-1, fp) != NULL)
 +    while (fgets(buf, STRLEN-1, fp) != nullptr)
      {
          fprintf(stdout, "%s", buf);
      }
@@@ -219,13 -218,13 +219,13 @@@ static int parse_logfile(const char *lo
          iFound = eFoundDDStr; /* Skip some case statements */
      }
  
 -    while (fgets(line, STRLEN, fp) != NULL)
 +    while (fgets(line, STRLEN, fp) != nullptr)
      {
          /* Remove leading spaces */
          ltrim(line);
  
          /* Check for TERM and INT signals from user: */
 -        if (std::strstr(line, errSIG) != NULL)
 +        if (std::strstr(line, errSIG) != nullptr)
          {
              fclose(fp);
              cleandata(perfdata, test_nr);
          /* Check whether cycle resetting  worked */
          if (presteps > 0 && !bFoundResetStr)
          {
 -            if (std::strstr(line, matchstrcr) != NULL)
 +            if (std::strstr(line, matchstrcr) != nullptr)
              {
                  sprintf(dumstring, "step %s", "%" GMX_SCNd64);
                  sscanf(line, dumstring, &resetsteps);
      if (gmx_fexist(errfile))
      {
          fp = fopen(errfile, "r");
 -        while (fgets(line, STRLEN, fp) != NULL)
 +        while (fgets(line, STRLEN, fp) != nullptr)
          {
              if (str_starts(line, "Fatal error:") )
              {
 -                if (fgets(line, STRLEN, fp) != NULL)
 +                if (fgets(line, STRLEN, fp) != nullptr)
                  {
                      fprintf(stderr, "\nWARNING: An error occurred during this benchmark:\n"
                              "%s\n", line);
@@@ -608,7 -607,7 +608,7 @@@ static void get_program_paths(gmx_bool 
      /* Get the commands we need to set up the runs from environment variables */
      if (!bThreads)
      {
 -        if ( (cp = getenv("MPIRUN")) != NULL)
 +        if ( (cp = getenv("MPIRUN")) != nullptr)
          {
              *cmd_mpirun = gmx_strdup(cp);
          }
          *cmd_mpirun = gmx_strdup(empty_mpirun);
      }
  
 -    if (*cmd_mdrun == NULL)
 +    if (*cmd_mdrun == nullptr)
      {
          /* The use of MDRUN is deprecated, but made available in 5.1
             for backward compatibility. It may be removed in a future
             version. */
 -        if ( (cp = getenv("MDRUN" )) != NULL)
 +        if ( (cp = getenv("MDRUN" )) != nullptr)
          {
              *cmd_mdrun = gmx_strdup(cp);
          }
@@@ -650,7 -649,7 +650,7 @@@ static void check_mdrun_works(gmx_boo
                                const char *cmd_mdrun,
                                gmx_bool    bNeedGpuSupport)
  {
 -    char      *command = NULL;
 +    char      *command = nullptr;
      char      *cp;
      char       line[STRLEN];
      FILE      *fp;
       * gmx_print_version_info() in the GMX_MPI section */
      const char match_mpi[]     = "MPI library:        MPI";
      const char match_mdrun[]   = "Executable: ";
-     const char match_gpu[]     = "GPU support:        enabled";
+     const char match_nogpu[]   = "GPU support:        disabled";
      gmx_bool   bMdrun          = FALSE;
      gmx_bool   bMPI            = FALSE;
-     gmx_bool   bHaveGpuSupport = FALSE;
+     gmx_bool   bHaveGpuSupport = TRUE;
  
      /* Run a small test to see whether mpirun + mdrun work  */
      fprintf(stdout, "Making sure that mdrun can be executed. ");
      if (bThreads)
      {
          snew(command, std::strlen(cmd_mdrun) + std::strlen(cmd_np) + std::strlen(filename) + 50);
-         sprintf(command, "%s%s-version -maxh 0.001 1> %s 2>&1", cmd_mdrun, cmd_np, filename);
+         sprintf(command, "%s%s -version -maxh 0.001 1> %s 2>&1", cmd_mdrun, cmd_np, filename);
      }
      else
      {
      while (!feof(fp) )
      {
          cp = fgets(line, STRLEN, fp);
 -        if (cp != NULL)
 +        if (cp != nullptr)
          {
              if (str_starts(line, match_mdrun) )
              {
              {
                  bMPI = TRUE;
              }
-             if (str_starts(line, match_gpu) )
+             if (str_starts(line, match_nogpu) )
              {
-                 bHaveGpuSupport = TRUE;
+                 bHaveGpuSupport = FALSE;
              }
          }
      }
@@@ -870,12 -869,12 +870,12 @@@ static void modify_PMEsettings
          const char     *fn_best_tpr, /* tpr file with the best performance */
          const char     *fn_sim_tpr)  /* name of tpr file to be launched */
  {
 -    t_inputrec   *ir;
 -    t_state       state;
 -    gmx_mtop_t    mtop;
 -    char          buf[200];
 +    t_state        state;
 +    gmx_mtop_t     mtop;
 +    char           buf[200];
  
 -    snew(ir, 1);
 +    t_inputrec     irInstance;
 +    t_inputrec    *ir = &irInstance;
      read_tpx_state(fn_best_tpr, ir, &state, &mtop);
  
      /* Reset nsteps and init_step to the value of the input .tpr file */
      fprintf(stdout, buf, ir->nsteps);
      fflush(stdout);
      write_tpx_state(fn_sim_tpr, ir, &state, &mtop);
 -
 -    sfree(ir);
  }
  
  static gmx_bool can_scale_rvdw(int vdwtype)
@@@ -913,6 -914,7 +913,6 @@@ static void make_benchmark_tprs
          FILE           *fp)              /* Write the output here                         */
  {
      int           i, j, d;
 -    t_inputrec   *ir;
      t_state       state;
      gmx_mtop_t    mtop;
      real          nlist_buffer;     /* Thickness of the buffer regions for PME-switch potentials */
      }
      fprintf(stdout, ".\n");
  
 -
 -    snew(ir, 1);
 +    t_inputrec  irInstance;
 +    t_inputrec *ir = &irInstance;
      read_tpx_state(fn_sim_tpr, ir, &state, &mtop);
  
      /* Check if some kind of PME was chosen */
  
              /* Scale the Fourier grid spacing */
              ir->nkx = ir->nky = ir->nkz = 0;
 -            calc_grid(NULL, state.box, fourierspacing*fac, &ir->nkx, &ir->nky, &ir->nkz);
 +            calcFftGrid(nullptr, state.box, fourierspacing*fac, minimalPmeGridSize(ir->pme_order),
 +                        &ir->nkx, &ir->nky, &ir->nkz);
  
              /* Adjust other radii since various conditions need to be fulfilled */
              if (eelPME == ir->coulombtype)
      }
      fflush(stdout);
      fflush(fp);
 -    sfree(ir);
  }
  
  
@@@ -1178,7 -1180,7 +1178,7 @@@ static void cleanup(const t_filenm *fnm
  {
      char        numstring[STRLEN];
      char        newfilename[STRLEN];
 -    const char *fn = NULL;
 +    const char *fn = nullptr;
      int         i;
      const char *opt;
  
@@@ -1447,8 -1449,8 +1447,8 @@@ static void do_the_tests
                                                     * constructing mdrun command lines */
  {
      int      i, nr, k, ret, count = 0, totaltests;
 -    int     *nPMEnodes = NULL;
 -    t_perf  *pd        = NULL;
 +    int     *nPMEnodes = nullptr;
 +    t_perf  *pd        = nullptr;
      int      cmdline_length;
      char    *command, *cmd_stub;
      char     buf[STRLEN];
          /* Loop over various numbers of PME nodes: */
          for (i = 0; i < *pmeentries; i++)
          {
 -            char *cmd_gpu_ids = NULL;
 +            char *cmd_gpu_ids = nullptr;
  
              pd = &perfdata[k][i];
  
@@@ -2045,21 -2047,20 +2045,21 @@@ static float inspect_tpr(int nfile, t_f
      gmx_bool     bFree;     /* Is a free energy simulation requested?         */
      gmx_bool     bNM;       /* Is a normal mode analysis requested?           */
      gmx_bool     bSwap;     /* Is water/ion position swapping requested?      */
 -    t_inputrec   ir;
      t_state      state;
      gmx_mtop_t   mtop;
  
  
      /* Check tpr file for options that trigger extra output files */
 -    read_tpx_state(opt2fn("-s", nfile, fnm), &ir, &state, &mtop);
 -    bFree = (efepNO  != ir.efep );
 -    bNM   = (eiNM    == ir.eI   );
 -    bSwap = (eswapNO != ir.eSwapCoords);
 -    bTpi  = EI_TPI(ir.eI);
 +    t_inputrec  irInstance;
 +    t_inputrec *ir = &irInstance;
 +    read_tpx_state(opt2fn("-s", nfile, fnm), ir, &state, &mtop);
 +    bFree = (efepNO  != ir->efep );
 +    bNM   = (eiNM    == ir->eI   );
 +    bSwap = (eswapNO != ir->eSwapCoords);
 +    bTpi  = EI_TPI(ir->eI);
  
      /* Set these output files on the tuning command-line */
 -    if (ir.bPull)
 +    if (ir->bPull)
      {
          setopt("-pf", nfile, fnm);
          setopt("-px", nfile, fnm);
          setopt("-swap", nfile, fnm);
      }
  
 -    *rcoulomb = ir.rcoulomb;
 +    *rcoulomb = ir->rcoulomb;
  
      /* Return the estimate for the number of PME nodes */
 -    return pme_load_estimate(&mtop, &ir, state.box);
 +    float npme = pme_load_estimate(&mtop, ir, state.box);
 +    return npme;
  }
  
  
@@@ -2206,25 -2206,25 +2206,25 @@@ int gmx_tune_pme(int argc, char *argv[]
      int             presteps       = 1500; /* Do a full cycle reset after presteps steps */
      gmx_bool        bOverwrite     = FALSE, bKeepTPR;
      gmx_bool        bLaunch        = FALSE;
 -    char           *ExtraArgs      = NULL;
 -    char          **tpr_names      = NULL;
 -    const char     *simulation_tpr = NULL;
 -    char           *deffnm         = NULL;
 +    char           *ExtraArgs      = nullptr;
 +    char          **tpr_names      = nullptr;
 +    const char     *simulation_tpr = nullptr;
 +    char           *deffnm         = nullptr;
      int             best_npme, best_tpr;
      int             sim_part = 1; /* For benchmarks with checkpoint files */
      char            bbuf[STRLEN];
  
  
      /* Default program names if nothing else is found */
 -    char               *cmd_mpirun = NULL, *cmd_mdrun = NULL;
 +    char               *cmd_mpirun = nullptr, *cmd_mdrun = nullptr;
      char               *cmd_args_bench, *cmd_args_launch;
 -    char               *cmd_np           = NULL;
 +    char               *cmd_np           = nullptr;
  
      /* IDs of GPUs that are eligible for computation */
 -    char               *eligible_gpu_ids = NULL;
 -    t_eligible_gpu_ids *gpu_ids          = NULL;
 +    char               *eligible_gpu_ids = nullptr;
 +    t_eligible_gpu_ids *gpu_ids          = nullptr;
  
 -    t_perf            **perfdata = NULL;
 +    t_perf            **perfdata = nullptr;
      t_inputinfo        *info;
      int                 i;
      FILE               *fp;
          { efLOG, "-err",    "bencherr", ffWRITE },
          { efTPR, "-so",     "tuned",    ffWRITE },
          /* mdrun: */
 -        { efTPR, NULL,      NULL,       ffREAD },
 -        { efTRN, "-o",      NULL,       ffWRITE },
 -        { efCOMPRESSED, "-x", NULL,     ffOPTWR },
 -        { efCPT, "-cpi",    NULL,       ffOPTRD },
 -        { efCPT, "-cpo",    NULL,       ffOPTWR },
 +        { efTPR, nullptr,      nullptr,       ffREAD },
 +        { efTRN, "-o",      nullptr,       ffWRITE },
 +        { efCOMPRESSED, "-x", nullptr,     ffOPTWR },
 +        { efCPT, "-cpi",    nullptr,       ffOPTRD },
 +        { efCPT, "-cpo",    nullptr,       ffOPTWR },
          { efSTO, "-c",      "confout",  ffWRITE },
          { efEDR, "-e",      "ener",     ffWRITE },
          { efLOG, "-g",      "md",       ffWRITE },
      int             nthreads = 1;
  
      const char     *procstring[] =
 -    { NULL, "np", "n", "none", NULL };
 +    { nullptr, "np", "n", "none", nullptr };
      const char     *npmevalues_opt[] =
 -    { NULL, "auto", "all", "subset", NULL };
 +    { nullptr, "auto", "all", "subset", nullptr };
  
      gmx_bool          bAppendFiles          = TRUE;
      gmx_bool          bKeepAndNumCPT        = FALSE;
      gmx_bool          bBenchmark            = TRUE;
      gmx_bool          bCheck                = TRUE;
  
 -    gmx_output_env_t *oenv = NULL;
 +    gmx_output_env_t *oenv = nullptr;
  
      t_pargs           pa[] = {
          /***********************/
  
      if (!parse_common_args(&argc, argv, PCA_NOEXIT_ON_ARGS,
                             NFILE, fnm, asize(pa), pa, asize(desc), desc,
 -                           0, NULL, &oenv))
 +                           0, nullptr, &oenv))
      {
          return 0;
      }
  
      // procstring[0] is used inside two different conditionals further down
 -    GMX_RELEASE_ASSERT(procstring[0] != NULL, "Options inconsistency; procstring[0] is NULL");
 +    GMX_RELEASE_ASSERT(procstring[0] != nullptr, "Options inconsistency; procstring[0] is NULL");
  
      /* Store the remaining unparsed command line entries in a string which
       * is then attached to the mdrun command line */
      get_program_paths(bThreads, &cmd_mpirun, &cmd_mdrun);
      if (bBenchmark && repeats > 0)
      {
 -        check_mdrun_works(bThreads, cmd_mpirun, cmd_np, cmd_mdrun, NULL != eligible_gpu_ids);
 +        check_mdrun_works(bThreads, cmd_mpirun, cmd_np, cmd_mdrun, nullptr != eligible_gpu_ids);
      }
  
      /* Print some header info to file */
      snew(perfdata, ntprs);
      if (bBenchmark)
      {
 -        GMX_RELEASE_ASSERT(npmevalues_opt[0] != NULL, "Options inconsistency; npmevalues_opt[0] is NULL");
 +        GMX_RELEASE_ASSERT(npmevalues_opt[0] != nullptr, "Options inconsistency; npmevalues_opt[0] is NULL");
          do_the_tests(fp, tpr_names, maxPMEnodes, minPMEnodes, npme_fixed, npmevalues_opt[0], perfdata, &pmeentries,
                       repeats, nnodes, ntprs, bThreads, cmd_mpirun, cmd_np, cmd_mdrun,
                       cmd_args_bench, fnm, NFILE, presteps, cpt_steps, bCheck, gpu_ids);
index bb8e44bb2cb5ed001bb3d861100118eb0a8ec6ce,40585ad45e19a66ad360f6c003f0bbd445d11095..c90f0dfec62723e0ea3f961dc049e2447786b3e0
  #include "gromacs/math/units.h"
  #include "gromacs/math/vec.h"
  #include "gromacs/mdlib/calc_verletbuf.h"
 +#include "gromacs/mdrunutility/mdmodules.h"
  #include "gromacs/mdtypes/inputrec.h"
  #include "gromacs/mdtypes/md_enums.h"
  #include "gromacs/mdtypes/pull-params.h"
 +#include "gromacs/options/options.h"
 +#include "gromacs/options/treesupport.h"
  #include "gromacs/pbcutil/pbc.h"
  #include "gromacs/topology/block.h"
  #include "gromacs/topology/ifunc.h"
  #include "gromacs/topology/symtab.h"
  #include "gromacs/topology/topology.h"
  #include "gromacs/utility/cstringutil.h"
 +#include "gromacs/utility/exceptions.h"
  #include "gromacs/utility/fatalerror.h"
 +#include "gromacs/utility/filestream.h"
 +#include "gromacs/utility/gmxassert.h"
 +#include "gromacs/utility/ikeyvaluetreeerror.h"
 +#include "gromacs/utility/keyvaluetree.h"
 +#include "gromacs/utility/keyvaluetreetransform.h"
  #include "gromacs/utility/smalloc.h"
 +#include "gromacs/utility/stringcompare.h"
 +#include "gromacs/utility/stringutil.h"
  
  #define MAXPTR 254
  #define NOGID  255
@@@ -107,10 -96,12 +107,10 @@@ typedef struct t_inputrec_string
      char   QMmethod[STRLEN], QMbasis[STRLEN], QMcharge[STRLEN], QMmult[STRLEN],
             bSH[STRLEN], CASorbitals[STRLEN], CASelectrons[STRLEN], SAon[STRLEN],
             SAoff[STRLEN], SAsteps[STRLEN], bTS[STRLEN], bOPT[STRLEN];
 -    char efield_x[STRLEN], efield_xt[STRLEN], efield_y[STRLEN],
 -         efield_yt[STRLEN], efield_z[STRLEN], efield_zt[STRLEN];
  
  } gmx_inputrec_strings;
  
 -static gmx_inputrec_strings *is = NULL;
 +static gmx_inputrec_strings *is = nullptr;
  
  void init_inputrec_strings()
  {
  void done_inputrec_strings()
  {
      sfree(is);
 -    is = NULL;
 +    is = nullptr;
  }
  
  
@@@ -139,13 -130,22 +139,13 @@@ enum 
  };
  
  static const char *constraints[eshNR+1]    = {
 -    "none", "h-bonds", "all-bonds", "h-angles", "all-angles", NULL
 +    "none", "h-bonds", "all-bonds", "h-angles", "all-angles", nullptr
  };
  
  static const char *couple_lam[ecouplamNR+1]    = {
 -    "vdw-q", "vdw", "q", "none", NULL
 +    "vdw-q", "vdw", "q", "none", nullptr
  };
  
 -void init_ir(t_inputrec *ir, t_gromppopts *opts)
 -{
 -    snew(opts->include, STRLEN);
 -    snew(opts->define, STRLEN);
 -    snew(ir->fepvals, 1);
 -    snew(ir->expandedvals, 1);
 -    snew(ir->simtempvals, 1);
 -}
 -
  static void GetSimTemps(int ntemps, t_simtemp *simtemp, double *temperature_lambdas)
  {
  
@@@ -811,7 -811,7 +811,7 @@@ void check_ir(const char *mdparin, t_in
          if (expand->nstTij > 0)
          {
              sprintf(err_buf, "nstlog must be non-zero");
-             CHECK(ir->nstlog != 0);
+             CHECK(ir->nstlog == 0);
              sprintf(err_buf, "nst-transition-matrix (%d) must be an integer multiple of nstlog (%d)",
                      expand->nstTij, ir->nstlog);
              CHECK((expand->nstTij % ir->nstlog) != 0);
          CHECK(EEL_FULL(ir->coulombtype) || ir->implicit_solvent == eisGBSA);
      }
  
 -    if (getenv("GMX_DO_GALACTIC_DYNAMICS") == NULL)
 +    if (getenv("GMX_DO_GALACTIC_DYNAMICS") == nullptr)
      {
          sprintf(err_buf, "epsilon-r must be >= 0 instead of %g\n", ir->epsilon_r);
          CHECK(ir->epsilon_r < 0);
  
      if (EEL_PME(ir->coulombtype) || EVDW_PME(ir->vdwtype))
      {
 -        if (ir->pme_order < 3)
 +        // TODO: Move these checks into the ewald module with the options class
 +        int orderMin = 3;
 +        int orderMax = (ir->coulombtype == eelP3M_AD ? 8 : 12);
 +
 +        if (ir->pme_order < orderMin || ir->pme_order > orderMax)
          {
 -            warning_error(wi, "pme-order can not be smaller than 3");
 +            sprintf(warn_buf, "With coulombtype = %s, you should have %d <= pme-order <= %d", eel_names[ir->coulombtype], orderMin, orderMax);
 +            warning_error(wi, warn_buf);
          }
      }
  
@@@ -1430,7 -1425,7 +1430,7 @@@ int str_nelem(const char *str, int maxp
          }
          ltrim(copy);
      }
 -    if (ptr == NULL)
 +    if (ptr == nullptr)
      {
          sfree(copy0);
      }
@@@ -1647,8 -1642,8 +1647,8 @@@ static void do_wall_params(t_inputrec *
      char  *names[MAXPTR];
      double dbl;
  
 -    opts->wall_atomtype[0] = NULL;
 -    opts->wall_atomtype[1] = NULL;
 +    opts->wall_atomtype[0] = nullptr;
 +    opts->wall_atomtype[1] = nullptr;
  
      ir->wall_atomtype[0] = -1;
      ir->wall_atomtype[1] = -1;
@@@ -1764,53 -1759,9 +1764,53 @@@ static gmx_bool couple_lambda_has_vdw_o
              couple_lambda_value == ecouplamVDWQ);
  }
  
 +namespace
 +{
 +
 +class MdpErrorHandler : public gmx::IKeyValueTreeErrorHandler
 +{
 +    public:
 +        explicit MdpErrorHandler(warninp_t wi)
 +            : wi_(wi), mapping_(nullptr)
 +        {
 +        }
 +
 +        void setBackMapping(const gmx::IKeyValueTreeBackMapping &mapping)
 +        {
 +            mapping_ = &mapping;
 +        }
 +
 +        virtual bool onError(gmx::UserInputError *ex, const gmx::KeyValueTreePath &context)
 +        {
 +            ex->prependContext(gmx::formatString("Error in mdp option \"%s\":",
 +                                                 getOptionName(context).c_str()));
 +            std::string message = gmx::formatExceptionMessageToString(*ex);
 +            warning_error(wi_, message.c_str());
 +            return true;
 +        }
 +
 +    private:
 +        std::string getOptionName(const gmx::KeyValueTreePath &context)
 +        {
 +            if (mapping_ != nullptr)
 +            {
 +                gmx::KeyValueTreePath path = mapping_->originalPath(context);
 +                GMX_ASSERT(path.size() == 1, "Inconsistent mapping back to mdp options");
 +                return path[0];
 +            }
 +            GMX_ASSERT(context.size() == 1, "Inconsistent context for mdp option parsing");
 +            return context[0];
 +        }
 +
 +        warninp_t                            wi_;
 +        const gmx::IKeyValueTreeBackMapping *mapping_;
 +};
 +
 +} // namespace
 +
  void get_ir(const char *mdparin, const char *mdparout,
 -            t_inputrec *ir, t_gromppopts *opts,
 -            warninp_t wi)
 +            gmx::MDModules *mdModules, t_inputrec *ir, t_gromppopts *opts,
 +            WriteMdpHeader writeMdpHeader, warninp_t wi)
  {
      char       *dumstr[2];
      double      dumdub[2][6];
      t_expanded *expand = ir->expandedvals;
  
      init_inputrec_strings();
 -    inp = read_inpfile(mdparin, &ninp, wi);
 +    gmx::TextInputFile stream(mdparin);
 +    inp = read_inpfile(&stream, mdparin, &ninp, wi);
  
      snew(dumstr[0], STRLEN);
      snew(dumstr[1], STRLEN);
      CCTYPE ("VARIOUS PREPROCESSING OPTIONS");
      CTYPE ("Preprocessor information: use cpp syntax.");
      CTYPE ("e.g.: -I/home/joe/doe -I/home/mary/roe");
 -    STYPE ("include", opts->include,  NULL);
 +    STYPE ("include", opts->include,  nullptr);
      CTYPE ("e.g.: -DPOSRES -DFLEXIBLE (note these variable names are case sensitive)");
 -    STYPE ("define",  opts->define,   NULL);
 +    STYPE ("define",  opts->define,   nullptr);
  
      CCTYPE ("RUN CONTROL PARAMETERS");
      EETYPE("integrator",  ir->eI,         ei_names);
      CTYPE ("number of steps for center of mass motion removal");
      ITYPE ("nstcomm", ir->nstcomm,    100);
      CTYPE ("group(s) for center of mass motion removal");
 -    STYPE ("comm-grps",   is->vcm,            NULL);
 +    STYPE ("comm-grps",   is->vcm,            nullptr);
  
      CCTYPE ("LANGEVIN DYNAMICS OPTIONS");
      CTYPE ("Friction coefficient (amu/ps) and random seed");
      CTYPE ("This selects the subset of atoms for the compressed");
      CTYPE ("trajectory file. You can select multiple groups. By");
      CTYPE ("default, all atoms will be written.");
 -    STYPE ("compressed-x-grps", is->x_compressed_groups, NULL);
 +    STYPE ("compressed-x-grps", is->x_compressed_groups, nullptr);
      CTYPE ("Selection of energy groups");
 -    STYPE ("energygrps",  is->energy,         NULL);
 +    STYPE ("energygrps",  is->energy,         nullptr);
  
      /* Neighbor searching */
      CCTYPE ("NEIGHBORSEARCHING PARAMETERS");
      CTYPE ("Extension of the potential lookup tables beyond the cut-off");
      RTYPE ("table-extension", ir->tabext, 1.0);
      CTYPE ("Separate tables between energy group pairs");
 -    STYPE ("energygrp-table", is->egptable,   NULL);
 +    STYPE ("energygrp-table", is->egptable,   nullptr);
      CTYPE ("Spacing for the PME/PPPM FFT grid");
      RTYPE ("fourierspacing", ir->fourier_spacing, 0.12);
      CTYPE ("FFT grid size, when a value is 0 fourierspacing will be used");
      ITYPE("nh-chain-length",     ir->opts.nhchainlength, 10);
      EETYPE("print-nose-hoover-chain-variables", ir->bPrintNHChains, yesno_names);
      CTYPE ("Groups to couple separately");
 -    STYPE ("tc-grps",     is->tcgrps,         NULL);
 +    STYPE ("tc-grps",     is->tcgrps,         nullptr);
      CTYPE ("Time constant (ps) and reference temperature (K)");
 -    STYPE ("tau-t",   is->tau_t,      NULL);
 -    STYPE ("ref-t",   is->ref_t,      NULL);
 +    STYPE ("tau-t",   is->tau_t,      nullptr);
 +    STYPE ("ref-t",   is->ref_t,      nullptr);
      CTYPE ("pressure coupling");
      EETYPE("pcoupl",  ir->epc,        epcoupl_names);
      EETYPE("pcoupltype",  ir->epct,       epcoupltype_names);
      ITYPE ("nstpcouple", ir->nstpcouple,  -1);
      CTYPE ("Time constant (ps), compressibility (1/bar) and reference P (bar)");
      RTYPE ("tau-p",   ir->tau_p,  1.0);
 -    STYPE ("compressibility", dumstr[0],  NULL);
 -    STYPE ("ref-p",       dumstr[1],      NULL);
 +    STYPE ("compressibility", dumstr[0],  nullptr);
 +    STYPE ("ref-p",       dumstr[1],      nullptr);
      CTYPE ("Scaling of reference coordinates, No, All or COM");
      EETYPE ("refcoord-scaling", ir->refcoord_scaling, erefscaling_names);
  
      CCTYPE ("OPTIONS FOR QMMM calculations");
      EETYPE("QMMM", ir->bQMMM, yesno_names);
      CTYPE ("Groups treated Quantum Mechanically");
 -    STYPE ("QMMM-grps",  is->QMMM,          NULL);
 +    STYPE ("QMMM-grps",  is->QMMM,          nullptr);
      CTYPE ("QM method");
 -    STYPE("QMmethod",     is->QMmethod, NULL);
 +    STYPE("QMmethod",     is->QMmethod, nullptr);
      CTYPE ("QMMM scheme");
      EETYPE("QMMMscheme",  ir->QMMMscheme,    eQMMMscheme_names);
      CTYPE ("QM basisset");
 -    STYPE("QMbasis",      is->QMbasis, NULL);
 +    STYPE("QMbasis",      is->QMbasis, nullptr);
      CTYPE ("QM charge");
 -    STYPE ("QMcharge",    is->QMcharge, NULL);
 +    STYPE ("QMcharge",    is->QMcharge, nullptr);
      CTYPE ("QM multiplicity");
 -    STYPE ("QMmult",      is->QMmult, NULL);
 +    STYPE ("QMmult",      is->QMmult, nullptr);
      CTYPE ("Surface Hopping");
 -    STYPE ("SH",          is->bSH, NULL);
 +    STYPE ("SH",          is->bSH, nullptr);
      CTYPE ("CAS space options");
 -    STYPE ("CASorbitals",      is->CASorbitals,   NULL);
 -    STYPE ("CASelectrons",     is->CASelectrons,  NULL);
 -    STYPE ("SAon", is->SAon, NULL);
 -    STYPE ("SAoff", is->SAoff, NULL);
 -    STYPE ("SAsteps", is->SAsteps, NULL);
 +    STYPE ("CASorbitals",      is->CASorbitals,   nullptr);
 +    STYPE ("CASelectrons",     is->CASelectrons,  nullptr);
 +    STYPE ("SAon", is->SAon, nullptr);
 +    STYPE ("SAoff", is->SAoff, nullptr);
 +    STYPE ("SAsteps", is->SAsteps, nullptr);
      CTYPE ("Scale factor for MM charges");
      RTYPE ("MMChargeScaleFactor", ir->scalefactor, 1.0);
      CTYPE ("Optimization of QM subsystem");
 -    STYPE ("bOPT",          is->bOPT, NULL);
 -    STYPE ("bTS",          is->bTS, NULL);
 +    STYPE ("bOPT",          is->bOPT, nullptr);
 +    STYPE ("bTS",          is->bTS, nullptr);
  
      /* Simulated annealing */
      CCTYPE("SIMULATED ANNEALING");
      CTYPE ("Type of annealing for each temperature group (no/single/periodic)");
 -    STYPE ("annealing",   is->anneal,      NULL);
 +    STYPE ("annealing",   is->anneal,      nullptr);
      CTYPE ("Number of time points to use for specifying annealing in each group");
 -    STYPE ("annealing-npoints", is->anneal_npoints, NULL);
 +    STYPE ("annealing-npoints", is->anneal_npoints, nullptr);
      CTYPE ("List of times at the annealing points for each group");
 -    STYPE ("annealing-time",       is->anneal_time,       NULL);
 +    STYPE ("annealing-time",       is->anneal_time,       nullptr);
      CTYPE ("Temp. at each annealing point, for each group.");
 -    STYPE ("annealing-temp",  is->anneal_temp,  NULL);
 +    STYPE ("annealing-temp",  is->anneal_temp,  nullptr);
  
      /* Startup run */
      CCTYPE ("GENERATE VELOCITIES FOR STARTUP RUN");
      /* Energy group exclusions */
      CCTYPE ("ENERGY GROUP EXCLUSIONS");
      CTYPE ("Pairs of energy groups for which all non-bonded interactions are excluded");
 -    STYPE ("energygrp-excl", is->egpexcl,     NULL);
 +    STYPE ("energygrp-excl", is->egpexcl,     nullptr);
  
      /* Walls */
      CCTYPE ("WALLS");
      ITYPE ("nwall", ir->nwall, 0);
      EETYPE("wall-type",     ir->wall_type,   ewt_names);
      RTYPE ("wall-r-linpot", ir->wall_r_linpot, -1);
 -    STYPE ("wall-atomtype", is->wall_atomtype, NULL);
 -    STYPE ("wall-density",  is->wall_density,  NULL);
 +    STYPE ("wall-atomtype", is->wall_atomtype, nullptr);
 +    STYPE ("wall-density",  is->wall_density,  nullptr);
      RTYPE ("wall-ewald-zfac", ir->wall_ewald_zfac, 3);
  
      /* COM pulling */
      /* Interactive MD */
      ir->bIMD = FALSE;
      CCTYPE("Group to display and/or manipulate in interactive MD session");
 -    STYPE ("IMD-group", is->imd_grp, NULL);
 +    STYPE ("IMD-group", is->imd_grp, nullptr);
      if (is->imd_grp[0] != '\0')
      {
          snew(ir->imd, 1);
      CTYPE ("Orientation restraints force constant and tau for time averaging");
      RTYPE ("orire-fc",    ir->orires_fc,  0.0);
      RTYPE ("orire-tau",   ir->orires_tau, 0.0);
 -    STYPE ("orire-fitgrp", is->orirefitgrp,    NULL);
 +    STYPE ("orire-fitgrp", is->orirefitgrp,    nullptr);
      CTYPE ("Output frequency for trace(SD) and S to energy file");
      ITYPE ("nstorireout", ir->nstorireout, 100);
  
      /* free energy variables */
      CCTYPE ("Free energy variables");
      EETYPE("free-energy", ir->efep, efep_names);
 -    STYPE ("couple-moltype",  is->couple_moltype,  NULL);
 +    STYPE ("couple-moltype",  is->couple_moltype,  nullptr);
      EETYPE("couple-lambda0", opts->couple_lam0, couple_lam);
      EETYPE("couple-lambda1", opts->couple_lam1, couple_lam);
      EETYPE("couple-intramol", opts->bCoupleIntra, yesno_names);
      ITYPE ("init-lambda-state", fep->init_fep_state, -1);
      RTYPE ("delta-lambda", fep->delta_lambda, 0.0);
      ITYPE ("nstdhdl", fep->nstdhdl, 50);
 -    STYPE ("fep-lambdas", is->fep_lambda[efptFEP], NULL);
 -    STYPE ("mass-lambdas", is->fep_lambda[efptMASS], NULL);
 -    STYPE ("coul-lambdas", is->fep_lambda[efptCOUL], NULL);
 -    STYPE ("vdw-lambdas", is->fep_lambda[efptVDW], NULL);
 -    STYPE ("bonded-lambdas", is->fep_lambda[efptBONDED], NULL);
 -    STYPE ("restraint-lambdas", is->fep_lambda[efptRESTRAINT], NULL);
 -    STYPE ("temperature-lambdas", is->fep_lambda[efptTEMPERATURE], NULL);
 +    STYPE ("fep-lambdas", is->fep_lambda[efptFEP], nullptr);
 +    STYPE ("mass-lambdas", is->fep_lambda[efptMASS], nullptr);
 +    STYPE ("coul-lambdas", is->fep_lambda[efptCOUL], nullptr);
 +    STYPE ("vdw-lambdas", is->fep_lambda[efptVDW], nullptr);
 +    STYPE ("bonded-lambdas", is->fep_lambda[efptBONDED], nullptr);
 +    STYPE ("restraint-lambdas", is->fep_lambda[efptRESTRAINT], nullptr);
 +    STYPE ("temperature-lambdas", is->fep_lambda[efptTEMPERATURE], nullptr);
      ITYPE ("calc-lambda-neighbors", fep->lambda_neighbors, 1);
 -    STYPE ("init-lambda-weights", is->lambda_weights, NULL);
 +    STYPE ("init-lambda-weights", is->lambda_weights, nullptr);
      EETYPE("dhdl-print-energy", fep->edHdLPrintEnergy, edHdLPrintEnergy_names);
      RTYPE ("sc-alpha", fep->sc_alpha, 0.0);
      ITYPE ("sc-power", fep->sc_power, 1);
  
      /* Non-equilibrium MD stuff */
      CCTYPE("Non-equilibrium MD stuff");
 -    STYPE ("acc-grps",    is->accgrps,        NULL);
 -    STYPE ("accelerate",  is->acc,            NULL);
 -    STYPE ("freezegrps",  is->freeze,         NULL);
 -    STYPE ("freezedim",   is->frdim,          NULL);
 +    STYPE ("acc-grps",    is->accgrps,        nullptr);
 +    STYPE ("accelerate",  is->acc,            nullptr);
 +    STYPE ("freezegrps",  is->freeze,         nullptr);
 +    STYPE ("freezedim",   is->frdim,          nullptr);
      RTYPE ("cos-acceleration", ir->cos_accel, 0);
 -    STYPE ("deform",      is->deform,         NULL);
 +    STYPE ("deform",      is->deform,         nullptr);
  
      /* simulated tempering variables */
      CCTYPE("simulated tempering variables");
      }
  
      /* Electric fields */
 -    CCTYPE("Electric fields");
 -    CTYPE ("Format is number of terms (int) and for all terms an amplitude (real)");
 -    CTYPE ("and a phase angle (real)");
 -    STYPE ("E-x",     is->efield_x,   NULL);
 -    CTYPE ("Time dependent (pulsed) electric field. Format is omega, time for pulse");
 -    CTYPE ("peak, and sigma (width) for pulse. Sigma = 0 removes pulse, leaving");
 -    CTYPE ("the field to be a cosine function.");
 -    STYPE ("E-xt",    is->efield_xt,  NULL);
 -    STYPE ("E-y",     is->efield_y,   NULL);
 -    STYPE ("E-yt",    is->efield_yt,  NULL);
 -    STYPE ("E-z",     is->efield_z,   NULL);
 -    STYPE ("E-zt",    is->efield_zt,  NULL);
 +    {
 +        gmx::KeyValueTreeObject      convertedValues = flatKeyValueTreeFromInpFile(ninp, inp);
 +        gmx::KeyValueTreeTransformer transform;
 +        transform.rules()->addRule()
 +            .keyMatchType("/", gmx::StringCompareType::CaseAndDashInsensitive);
 +        mdModules->initMdpTransform(transform.rules());
 +        for (const auto &path : transform.mappedPaths())
 +        {
 +            GMX_ASSERT(path.size() == 1, "Inconsistent mapping back to mdp options");
 +            mark_einp_set(ninp, inp, path[0].c_str());
 +        }
 +        MdpErrorHandler              errorHandler(wi);
 +        auto                         result
 +                   = transform.transform(convertedValues, &errorHandler);
 +        ir->params = new gmx::KeyValueTreeObject(result.object());
 +        mdModules->adjustInputrecBasedOnModules(ir);
 +        errorHandler.setBackMapping(result.backMapping());
 +        mdModules->assignOptionsToModules(*ir->params, &errorHandler);
 +    }
  
      /* Ion/water position swapping ("computational electrophysiology") */
      CCTYPE("Ion/water position swapping for computational electrophysiology setups");
              snew(ir->swap->grp[i].molname, STRLEN);
          }
          CTYPE("Two index groups that contain the compartment-partitioning atoms");
 -        STYPE("split-group0", ir->swap->grp[eGrpSplit0].molname, NULL);
 -        STYPE("split-group1", ir->swap->grp[eGrpSplit1].molname, NULL);
 +        STYPE("split-group0", ir->swap->grp[eGrpSplit0].molname, nullptr);
 +        STYPE("split-group1", ir->swap->grp[eGrpSplit1].molname, nullptr);
          CTYPE("Use center of mass of split groups (yes/no), otherwise center of geometry is used");
          EETYPE("massw-split0", ir->swap->massw_split[0], yesno_names);
          EETYPE("massw-split1", ir->swap->massw_split[1], yesno_names);
  
          CTYPE("Name of solvent molecules");
 -        STYPE("solvent-group", ir->swap->grp[eGrpSolvent].molname, NULL);
 +        STYPE("solvent-group", ir->swap->grp[eGrpSolvent].molname, nullptr);
  
          CTYPE("Split cylinder: radius, upper and lower extension (nm) (this will define the channels)");
          CTYPE("Note that the split cylinder settings do not have an influence on the swapping protocol,");
              int ig = eSwapFixedGrpNR + i;
  
              sprintf(buf, "iontype%d-name", i);
 -            STYPE(buf, ir->swap->grp[ig].molname, NULL);
 +            STYPE(buf, ir->swap->grp[ig].molname, nullptr);
              sprintf(buf, "iontype%d-in-A", i);
              ITYPE(buf, ir->swap->grp[ig].nmolReq[0], -1);
              sprintf(buf, "iontype%d-in-B", i);
  
      /* User defined thingies */
      CCTYPE ("User defined thingies");
 -    STYPE ("user1-grps",  is->user1,          NULL);
 -    STYPE ("user2-grps",  is->user2,          NULL);
 +    STYPE ("user1-grps",  is->user1,          nullptr);
 +    STYPE ("user2-grps",  is->user2,          nullptr);
      ITYPE ("userint1",    ir->userint1,   0);
      ITYPE ("userint2",    ir->userint2,   0);
      ITYPE ("userint3",    ir->userint3,   0);
      RTYPE ("userreal4",   ir->userreal4,  0);
  #undef CTYPE
  
 -    write_inpfile(mdparout, ninp, inp, FALSE, wi);
 +    {
 +        gmx::TextOutputFile stream(mdparout);
 +        write_inpfile(&stream, mdparout, ninp, inp, FALSE, writeMdpHeader, wi);
 +        stream.close();
 +    }
 +
      for (i = 0; (i < ninp); i++)
      {
          sfree(inp[i].name);
          ir->nstcomm = 0;
      }
  
 -    opts->couple_moltype = NULL;
 +    opts->couple_moltype = nullptr;
      if (strlen(is->couple_moltype) > 0)
      {
          if (ir->efep != efepNO)
  
      /* ORIENTATION RESTRAINT PARAMETERS */
  
 -    if (opts->bOrire && str_nelem(is->orirefitgrp, MAXPTR, NULL) != 1)
 +    if (opts->bOrire && str_nelem(is->orirefitgrp, MAXPTR, nullptr) != 1)
      {
          warning_error(wi, "ERROR: Need one orientation restraint fit group\n");
      }
@@@ -2782,7 -2720,7 +2782,7 @@@ static gmx_bool do_numbering(int natoms
      {
          /* All atoms are part of one (or no) group, no index required */
          groups->ngrpnr[gtype] = 0;
 -        groups->grpnr[gtype]  = NULL;
 +        groups->grpnr[gtype]  = nullptr;
      }
      else
      {
@@@ -2810,6 -2748,7 +2810,6 @@@ static void calc_nrdf(gmx_mtop_t *mtop
      double                 *nrdf_tc, *nrdf_vcm, nrdf_uc, *nrdf_vcm_sub;
      ivec                   *dof_vcm;
      gmx_mtop_atomloop_all_t aloop;
 -    t_atom                 *atom;
      int                     mb, mol, ftype, as;
      gmx_molblock_t         *molb;
      gmx_moltype_t          *molt;
  
      snew(nrdf2, natoms);
      aloop = gmx_mtop_atomloop_all_init(mtop);
 +    const t_atom *atom;
      while (gmx_mtop_atomloop_all_next(aloop, &i, &atom))
      {
          nrdf2[i] = 0;
      sfree(nrdf_vcm_sub);
  }
  
 -static void decode_cos(char *s, t_cosines *cosine)
 -{
 -    char              *t;
 -    char               format[STRLEN], f1[STRLEN];
 -    double             a, phi;
 -    int                i;
 -
 -    t = gmx_strdup(s);
 -    trim(t);
 -
 -    cosine->n   = 0;
 -    cosine->a   = NULL;
 -    cosine->phi = NULL;
 -    if (strlen(t))
 -    {
 -        if (sscanf(t, "%d", &(cosine->n)) != 1)
 -        {
 -            gmx_fatal(FARGS, "Cannot parse cosine multiplicity from string '%s'", t);
 -        }
 -        if (cosine->n <= 0)
 -        {
 -            cosine->n = 0;
 -        }
 -        else
 -        {
 -            snew(cosine->a, cosine->n);
 -            snew(cosine->phi, cosine->n);
 -
 -            sprintf(format, "%%*d");
 -            for (i = 0; (i < cosine->n); i++)
 -            {
 -                double  gmx_unused canary;
 -
 -                strcpy(f1, format);
 -                strcat(f1, "%lf%lf%lf");
 -                if (sscanf(t, f1, &a, &phi, &canary) != 2)
 -                {
 -                    gmx_fatal(FARGS, "Invalid input for electric field shift: '%s'", t);
 -                }
 -                cosine->a[i]   = a;
 -                cosine->phi[i] = phi;
 -                strcat(format, "%*lf%*lf");
 -            }
 -        }
 -    }
 -    sfree(t);
 -}
 -
  static gmx_bool do_egp_flag(t_inputrec *ir, gmx_groups_t *groups,
                              const char *option, const char *val, int flag)
  {
@@@ -3232,7 -3218,7 +3232,7 @@@ void do_index(const char* mdparin, cons
      {
          fprintf(stderr, "processing index file...\n");
      }
 -    if (ndx == NULL)
 +    if (ndx == nullptr)
      {
          snew(grps, 1);
          snew(grps->index, 1);
          {
              ir->opts.annealing[i]      = eannNO;
              ir->opts.anneal_npoints[i] = 0;
 -            ir->opts.anneal_time[i]    = NULL;
 -            ir->opts.anneal_temp[i]    = NULL;
 +            ir->opts.anneal_time[i]    = nullptr;
 +            ir->opts.anneal_temp[i]    = nullptr;
          }
          if (nSA > 0)
          {
          gmx_fatal(FARGS, "Can only have energy group pair tables in combination with user tables for VdW and/or Coulomb");
      }
  
 -    decode_cos(is->efield_x, &(ir->ex[XX]));
 -    decode_cos(is->efield_xt, &(ir->et[XX]));
 -    decode_cos(is->efield_y, &(ir->ex[YY]));
 -    decode_cos(is->efield_yt, &(ir->et[YY]));
 -    decode_cos(is->efield_z, &(ir->ex[ZZ]));
 -    decode_cos(is->efield_zt, &(ir->et[ZZ]));
 -
      for (i = 0; (i < grps->nr); i++)
      {
          sfree(gnames[i]);
@@@ -3964,7 -3957,7 +3964,7 @@@ check_combination_rule_differences(cons
       */
      tol = 1e-5;
      ptr = getenv("GMX_LJCOMB_TOL");
 -    if (ptr != NULL)
 +    if (ptr != nullptr)
      {
          double            dbl;
          double gmx_unused canary;
@@@ -4082,6 -4075,7 +4082,6 @@@ void triple_check(const char *mdparin, 
      rvec                      acc;
      gmx_mtop_atomloop_block_t aloopb;
      gmx_mtop_atomloop_all_t   aloop;
 -    t_atom                   *atom;
      ivec                      AbsRef;
      char                      warn_buf[STRLEN];
  
  
      bCharge = FALSE;
      aloopb  = gmx_mtop_atomloop_block_init(sys);
 +    const t_atom *atom;
      while (gmx_mtop_atomloop_block_next(aloopb, &atom, &nmol))
      {
          if (atom->q != 0 || atom->qB != 0)
          clear_rvec(acc);
          snew(mgrp, sys->groups.grps[egcACC].nr);
          aloop = gmx_mtop_atomloop_all_init(sys);
 +        const t_atom *atom;
          while (gmx_mtop_atomloop_all_next(aloop, &i, &atom))
          {
              mgrp[ggrpnr(&sys->groups, egcACC, i)] += atom->m;
index 67e264e85457ba3efbab75d9ec86b9c5929698f0,e2c1d49ad9d28602826cfba70f393b783734536c..f8d979085bf2f18285f55f0c1d9c54e0feaf89ca
@@@ -86,7 -86,7 +86,7 @@@ static void sort_molecule(t_atoms **ato
      atoms = *atoms_solvt;
  
      /* copy each residue from *atoms to a molecule in *molecule */
 -    moltypes   = NULL;
 +    moltypes   = nullptr;
      nrmoltypes = 0;
      for (i = 0; i < atoms->nr; i++)
      {
@@@ -311,7 -311,7 +311,7 @@@ static void replicateSolventBox(t_atom
      // but not in all).
      t_atoms           newAtoms;
      init_t_atoms(&newAtoms, 0, FALSE);
 -    gmx::AtomsBuilder builder(&newAtoms, NULL);
 +    gmx::AtomsBuilder builder(&newAtoms, nullptr);
      builder.reserve(atoms->nr * nmol, atoms->nres * nmol);
      std::vector<RVec> newX(atoms->nr * nmol);
      std::vector<RVec> newV(!v->empty() ? atoms->nr * nmol : 0);
@@@ -491,50 -491,86 +491,86 @@@ static void removeSolventBoxOverlap(t_a
  }
  
  /*! \brief
-  * Removes solvent molecules that overlap with the solute, and optionally also
-  * those that are outside a given shell radius from the solute.
+  * Remove all solvent molecules outside a give radius from the solute.
   *
-  * \param[in,out] atoms      Solvent atoms.
-  * \param[in,out] x          Solvent positions.
-  * \param[in,out] v          Solvent velocities (can be empty).
-  * \param[in,out] r          Solvent exclusion radii.
-  * \param[in]     pbc        PBC information.
-  * \param[in]     x_solute   Solute positions.
-  * \param[in]     r_solute   Solute exclusion radii.
-  * \param[in]     rshell     If >0, only keep solvent atoms within a shell of
-  *     this size from the solute.
+  * \param[in,out] atoms     Solvent atoms.
+  * \param[in,out] x_solvent Solvent positions.
+  * \param[in,out] v_solvent Solvent velocities.
+  * \param[in,out] r         Atomic exclusion radii.
+  * \param[in]     pbc       PBC information.
+  * \param[in]     x_solute  Solute positions.
+  * \param[in]     rshell    The radius outside the solute molecule.
   */
- static void removeSoluteOverlap(t_atoms *atoms, std::vector<RVec> *x,
-                                 std::vector<RVec> *v, std::vector<real> *r,
-                                 const t_pbc &pbc,
-                                 const std::vector<RVec> &x_solute,
-                                 const std::vector<real> &r_solute,
-                                 real rshell)
+ static void removeSolventOutsideShell(t_atoms                 *atoms,
+                                       std::vector<RVec>       *x_solvent,
+                                       std::vector<RVec>       *v_solvent,
+                                       std::vector<real>       *r,
+                                       const t_pbc             &pbc,
+                                       const std::vector<RVec> &x_solute,
+                                       real                     rshell)
  {
-     const real                          maxRadius1
-         = *std::max_element(r->begin(), r->end());
-     const real                          maxRadius2
-         = *std::max_element(r_solute.begin(), r_solute.end());
      gmx::AtomsRemover                   remover(*atoms);
-     // If rshell is >0, the neighborhood search looks at all pairs
-     // within rshell, and unmarks those that are within the cutoff.
-     // This line marks everything, so that solvent outside rshell remains
-     // marked after the loop.
-     // Without rshell, the neighborhood search only marks the overlapping
-     // solvent atoms, and all others are left alone.
-     if (rshell > 0.0)
+     gmx::AnalysisNeighborhood           nb;
+     nb.setCutoff(rshell);
+     gmx::AnalysisNeighborhoodPositions  posSolute(x_solute);
+     gmx::AnalysisNeighborhoodSearch     search     = nb.initSearch(&pbc, posSolute);
+     gmx::AnalysisNeighborhoodPositions  pos(*x_solvent);
+     gmx::AnalysisNeighborhoodPairSearch pairSearch = search.startPairSearch(pos);
+     gmx::AnalysisNeighborhoodPair       pair;
+     // Remove everything
+     remover.markAll();
+     // Now put back those within the shell without checking for overlap
+     while (pairSearch.findNextPair(&pair))
      {
-         remover.markAll();
+         remover.markResidue(*atoms, pair.testIndex(), false);
+         pairSearch.skipRemainingPairsForTestPosition();
      }
+     remover.removeMarkedElements(x_solvent);
+     if (!v_solvent->empty())
+     {
+         remover.removeMarkedElements(v_solvent);
+     }
+     remover.removeMarkedElements(r);
+     const int originalAtomCount = atoms->nr;
+     remover.removeMarkedAtoms(atoms);
+     fprintf(stderr, "Removed %d solvent atoms more than %f nm from solute.\n",
+             originalAtomCount - atoms->nr, rshell);
+ }
+ /*! \brief
+  * Removes solvent molecules that overlap with the solute.
+  *
+  * \param[in,out] atoms    Solvent atoms.
+  * \param[in,out] x        Solvent positions.
+  * \param[in,out] v        Solvent velocities (can be empty).
+  * \param[in,out] r        Solvent exclusion radii.
+  * \param[in]     pbc      PBC information.
+  * \param[in]     x_solute Solute positions.
+  * \param[in]     r_solute Solute exclusion radii.
+  */
+ static void removeSolventOverlappingWithSolute(t_atoms                 *atoms,
+                                                std::vector<RVec>       *x,
+                                                std::vector<RVec>       *v,
+                                                std::vector<real>       *r,
+                                                const t_pbc             &pbc,
+                                                const std::vector<RVec> &x_solute,
+                                                const std::vector<real> &r_solute)
+ {
+     gmx::AtomsRemover             remover(*atoms);
+     const real                    maxRadius1
+         = *std::max_element(r->begin(), r->end());
+     const real                    maxRadius2
+         = *std::max_element(r_solute.begin(), r_solute.end());
  
+     // Now check for overlap.
      gmx::AnalysisNeighborhood           nb;
-     nb.setCutoff(std::max(maxRadius1 + maxRadius2, rshell));
+     gmx::AnalysisNeighborhoodPair       pair;
+     nb.setCutoff(maxRadius1 + maxRadius2);
      gmx::AnalysisNeighborhoodPositions  posSolute(x_solute);
      gmx::AnalysisNeighborhoodSearch     search     = nb.initSearch(&pbc, posSolute);
      gmx::AnalysisNeighborhoodPositions  pos(*x);
      gmx::AnalysisNeighborhoodPairSearch pairSearch = search.startPairSearch(pos);
-     gmx::AnalysisNeighborhoodPair       pair;
      while (pairSearch.findNextPair(&pair))
      {
          if (remover.isMarked(pair.testIndex()))
@@@ -615,7 -651,7 +651,7 @@@ static void add_solv(const char *fn, t_
  
      char             *filename = gmxlibfn(fn);
      snew(top_solvt, 1);
 -    readConformation(filename, top_solvt, &x_solvt, !v->empty() ? &v_solvt : NULL,
 +    readConformation(filename, top_solvt, &x_solvt, !v->empty() ? &v_solvt : nullptr,
                       &ePBC_solvt, box_solvt, "solvent");
      t_atoms *atoms_solvt = &top_solvt->atoms;
      if (0 == atoms_solvt->nr)
      }
      if (top->atoms.nr > 0)
      {
-         removeSoluteOverlap(atoms_solvt, &x_solvt, &v_solvt, &exclusionDistances_solvt, pbc,
-                             *x, exclusionDistances, rshell);
+         if (rshell > 0.0)
+         {
+             removeSolventOutsideShell(atoms_solvt, &x_solvt,  &v_solvt,
+                                       &exclusionDistances_solvt, pbc, *x, rshell);
+         }
+         removeSolventOverlappingWithSolute(atoms_solvt, &x_solvt, &v_solvt,
+                                            &exclusionDistances_solvt, pbc, *x,
+                                            exclusionDistances);
      }
  
      if (max_sol > 0 && atoms_solvt->nres > max_sol)
@@@ -732,7 -774,7 +774,7 @@@ static void update_top(t_atoms *atoms, 
              bSkip = FALSE;
              line++;
              strcpy(buf2, buf);
 -            if ((temp = strchr(buf2, '\n')) != NULL)
 +            if ((temp = strchr(buf2, '\n')) != nullptr)
              {
                  temp[0] = '\0';
              }
              if (buf2[0] == '[')
              {
                  buf2[0] = ' ';
 -                if ((temp = strchr(buf2, '\n')) != NULL)
 +                if ((temp = strchr(buf2, '\n')) != nullptr)
                  {
                      temp[0] = '\0';
                  }
              }
              if (bSkip)
              {
 -                if ((temp = strchr(buf, '\n')) != NULL)
 +                if ((temp = strchr(buf, '\n')) != nullptr)
                  {
                      temp[0] = '\0';
                  }
@@@ -878,8 -920,8 +920,8 @@@ int gmx_solvate(int argc, char *argv[]
      t_filenm    fnm[] = {
          { efSTX, "-cp", "protein", ffOPTRD },
          { efSTX, "-cs", "spc216",  ffLIBRD},
 -        { efSTO, NULL,  NULL,      ffWRITE},
 -        { efTOP, NULL,  NULL,      ffOPTRW},
 +        { efSTO, nullptr,  nullptr,      ffWRITE},
 +        { efTOP, nullptr,  nullptr,      ffOPTRW},
      };
  #define NFILE asize(fnm)
  
          /* Generate a solute configuration */
          conf_prot = opt2fn("-cp", NFILE, fnm);
          readConformation(conf_prot, top, &x,
 -                         bReadV ? &v : NULL, &ePBC, box, "solute");
 +                         bReadV ? &v : nullptr, &ePBC, box, "solute");
          if (bReadV && v.empty())
          {
              fprintf(stderr, "Note: no velocities found\n");
      fprintf(stderr, "Writing generated configuration to %s\n", confout);
      const char *outputTitle = (bProt ? *top->name : "Generated by gmx solvate");
      write_sto_conf(confout, outputTitle, &top->atoms, as_rvec_array(x.data()),
 -                   !v.empty() ? as_rvec_array(v.data()) : NULL, ePBC, box);
 +                   !v.empty() ? as_rvec_array(v.data()) : nullptr, ePBC, box);
  
      /* print size of generated configuration */
      fprintf(stderr, "\nOutput configuration contains %d atoms in %d residues\n",
index 3f32f7446340e4a59caa531d6f0bf95b2dec39a1,670eb6b0eb74ae570b7380dd501ac726a9cd36f3..d72582984028bc5aba19913336c80eee541da44c
@@@ -374,8 -374,8 +374,8 @@@ do_pairs_general(int ftype, int nbonds
              energygrp_vdw  = grppener->ener[egLJSR];
              break;
          default:
 -            energygrp_elec = NULL; /* Keep compiler happy */
 -            energygrp_vdw  = NULL; /* Keep compiler happy */
 +            energygrp_elec = nullptr; /* Keep compiler happy */
 +            energygrp_vdw  = nullptr; /* Keep compiler happy */
              gmx_fatal(FARGS, "Unknown function type %d in do_nonbonded14", ftype);
              break;
      }
@@@ -607,7 -607,7 +607,7 @@@ do_pairs_simple(int nbonds
          pbc_dx_aiuc(pbc, xi, xj, dr);
  
          T rsq   = dr[XX]*dr[XX] + dr[YY]*dr[YY] + dr[ZZ]*dr[ZZ];
-         T rinv  = invsqrt(rsq);
+         T rinv  = gmx::invsqrt(rsq);
          T rinv2 = rinv*rinv;
          T rinv6 = rinv2*rinv2*rinv2;
  
@@@ -668,13 -668,13 +668,13 @@@ do_pairs(int ftype, int nbonds
          t_pbc        pbc_no;
          const t_pbc *pbc_nonnull;
  
 -        if (pbc != NULL)
 +        if (pbc != nullptr)
          {
              pbc_nonnull   = pbc;
          }
          else
          {
 -            set_pbc(&pbc_no, epbcNONE, NULL);
 +            set_pbc(&pbc_no, epbcNONE, nullptr);
              pbc_nonnull   = &pbc_no;
          }
  
index 38dbbd4a67cd99ade440e11a71ca4665151e2389,4973896b2c3efd2d5eada84210a6b2b444b12a5a..f0bbde93e52d96c12ad3ad79cd3e124d9566ebc8
@@@ -45,8 -45,6 +45,8 @@@
  #include <stdio.h>
  #include <string.h>
  
 +#include <cstdint>
 +
  #include <array>
  
  #include "gromacs/domdec/domdec.h"
@@@ -54,6 -52,7 +54,6 @@@
  #include "gromacs/essentialdynamics/edsam.h"
  #include "gromacs/ewald/pme.h"
  #include "gromacs/gmxlib/chargegroup.h"
 -#include "gromacs/gmxlib/md_logging.h"
  #include "gromacs/gmxlib/network.h"
  #include "gromacs/gmxlib/nrnb.h"
  #include "gromacs/gmxlib/nonbonded/nb_free_energy.h"
  #include "gromacs/mdlib/nbnxn_kernels/simd_2xnn/nbnxn_kernel_simd_2xnn.h"
  #include "gromacs/mdlib/nbnxn_kernels/simd_4xn/nbnxn_kernel_simd_4xn.h"
  #include "gromacs/mdtypes/commrec.h"
 +#include "gromacs/mdtypes/iforceprovider.h"
  #include "gromacs/mdtypes/inputrec.h"
  #include "gromacs/mdtypes/md_enums.h"
 +#include "gromacs/mdtypes/state.h"
  #include "gromacs/pbcutil/ishift.h"
  #include "gromacs/pbcutil/mshift.h"
  #include "gromacs/pbcutil/pbc.h"
  #include "gromacs/timing/wallcycle.h"
  #include "gromacs/timing/wallcyclereporting.h"
  #include "gromacs/timing/walltime_accounting.h"
 +#include "gromacs/utility/arrayref.h"
  #include "gromacs/utility/basedefinitions.h"
  #include "gromacs/utility/cstringutil.h"
  #include "gromacs/utility/exceptions.h"
  #include "gromacs/utility/fatalerror.h"
  #include "gromacs/utility/gmxassert.h"
  #include "gromacs/utility/gmxmpi.h"
 +#include "gromacs/utility/logger.h"
  #include "gromacs/utility/pleasecite.h"
  #include "gromacs/utility/smalloc.h"
  #include "gromacs/utility/sysinfo.h"
@@@ -209,18 -204,89 +209,18 @@@ void print_start(FILE *fplog, t_commre
                          walltime_accounting_get_start_time_stamp(walltime_accounting));
  }
  
 -static void sum_forces(int start, int end, rvec f[], rvec flr[])
 +static void sum_forces(rvec f[], const PaddedRVecVector *forceToAdd)
  {
 -    int i;
 +    /* TODO: remove this - 1 when padding is properly implemented */
 +    int         end  = forceToAdd->size() - 1;
 +    const rvec *fAdd = as_rvec_array(forceToAdd->data());
  
 -    if (gmx_debug_at)
 -    {
 -        pr_rvecs(debug, 0, "fsr", f+start, end-start);
 -        pr_rvecs(debug, 0, "flr", flr+start, end-start);
 -    }
      // cppcheck-suppress unreadVariable
      int gmx_unused nt = gmx_omp_nthreads_get(emntDefault);
  #pragma omp parallel for num_threads(nt) schedule(static)
 -    for (i = start; i < end; i++)
 +    for (int i = 0; i < end; i++)
      {
 -        rvec_inc(f[i], flr[i]);
 -    }
 -}
 -
 -/*
 - * calc_f_el calculates forces due to an electric field.
 - *
 - * force is kJ mol^-1 nm^-1 = e * kJ mol^-1 nm^-1 / e
 - *
 - * Et[] contains the parameters for the time dependent
 - * part of the field.
 - * Ex[] contains the parameters for
 - * the spatial dependent part of the field.
 - * The function should return the energy due to the electric field
 - * (if any) but for now returns 0.
 - *
 - * WARNING:
 - * There can be problems with the virial.
 - * Since the field is not self-consistent this is unavoidable.
 - * For neutral molecules the virial is correct within this approximation.
 - * For neutral systems with many charged molecules the error is small.
 - * But for systems with a net charge or a few charged molecules
 - * the error can be significant when the field is high.
 - * Solution: implement a self-consistent electric field into PME.
 - */
 -static void calc_f_el(FILE *fp, int  start, int homenr,
 -                      real charge[], rvec f[],
 -                      t_cosines Ex[], t_cosines Et[], double t)
 -{
 -    rvec Ext;
 -    real t0;
 -    int  i, m;
 -
 -    for (m = 0; (m < DIM); m++)
 -    {
 -        if (Et[m].n > 0)
 -        {
 -            if (Et[m].n == 3)
 -            {
 -                t0     = Et[m].a[1];
 -                Ext[m] = std::cos(Et[m].a[0]*(t-t0))*std::exp(-gmx::square(t-t0)/(2.0*gmx::square(Et[m].a[2])));
 -            }
 -            else
 -            {
 -                Ext[m] = std::cos(Et[m].a[0]*t);
 -            }
 -        }
 -        else
 -        {
 -            Ext[m] = 1.0;
 -        }
 -        if (Ex[m].n > 0)
 -        {
 -            /* Convert the field strength from V/nm to MD-units */
 -            Ext[m] *= Ex[m].a[0]*FIELDFAC;
 -            for (i = start; (i < start+homenr); i++)
 -            {
 -                f[i][m] += charge[i]*Ext[m];
 -            }
 -        }
 -        else
 -        {
 -            Ext[m] = 0;
 -        }
 -    }
 -    if (fp != NULL)
 -    {
 -        fprintf(fp, "%10g  %10g  %10g  %10g #FIELD\n", t,
 -                Ext[XX]/FIELDFAC, Ext[YY]/FIELDFAC, Ext[ZZ]/FIELDFAC);
 +        rvec_inc(f[i], fAdd[i]);
      }
  }
  
@@@ -305,7 -371,7 +305,7 @@@ static void pme_receive_force_ener(t_co
      wallcycle_start(wcycle, ewcPP_PMEWAITRECVF);
      dvdl_q  = 0;
      dvdl_lj = 0;
 -    gmx_pme_receive_f(cr, fr->f_novirsum, fr->vir_el_recip, &e_q,
 +    gmx_pme_receive_f(cr, as_rvec_array(fr->f_novirsum->data()), fr->vir_el_recip, &e_q,
                        fr->vir_lj_recip, &e_lj, &dvdl_q, &dvdl_lj,
                        &cycles_seppme);
      enerd->term[F_COUL_RECIP] += e_q;
  }
  
  static void print_large_forces(FILE *fp, t_mdatoms *md, t_commrec *cr,
 -                               gmx_int64_t step, real pforce, rvec *x, rvec *f)
 +                               gmx_int64_t step, real forceTolerance,
 +                               const rvec *x, const rvec *f)
  {
 -    int  i;
 -    real pf2, fn2;
 -    char buf[STEPSTRSIZE];
 -
 -    pf2 = gmx::square(pforce);
 -    for (i = 0; i < md->homenr; i++)
 +    real           force2Tolerance = gmx::square(forceTolerance);
 +    std::uintmax_t numNonFinite    = 0;
 +    for (int i = 0; i < md->homenr; i++)
      {
 -        fn2 = norm2(f[i]);
 -        /* We also catch NAN, if the compiler does not optimize this away. */
 -        if (fn2 >= pf2 || fn2 != fn2)
 +        real force2    = norm2(f[i]);
 +        bool nonFinite = !std::isfinite(force2);
 +        if (force2 >= force2Tolerance || nonFinite)
 +        {
 +            fprintf(fp, "step %" GMX_PRId64 " atom %6d  x %8.3f %8.3f %8.3f  force %12.5e\n",
 +                    step,
 +                    ddglatnr(cr->dd, i), x[i][XX], x[i][YY], x[i][ZZ], std::sqrt(force2));
 +        }
 +        if (nonFinite)
          {
 -            fprintf(fp, "step %s  atom %6d  x %8.3f %8.3f %8.3f  force %12.5e\n",
 -                    gmx_step_str(step, buf),
 -                    ddglatnr(cr->dd, i), x[i][XX], x[i][YY], x[i][ZZ], std::sqrt(fn2));
 +            numNonFinite++;
          }
      }
 +    if (numNonFinite > 0)
 +    {
 +        /* Note that with MPI this fatal call on one rank might interrupt
 +         * the printing on other ranks. But we can only avoid that with
 +         * an expensive MPI barrier that we would need at each step.
 +         */
 +        gmx_fatal(FARGS, "At step %" GMX_PRId64 " detected non-finite forces on %ju atoms", step, numNonFinite);
 +    }
  }
  
  static void post_process_forces(t_commrec *cr,
               * if the constructing atoms aren't local.
               */
              wallcycle_start(wcycle, ewcVSITESPREAD);
 -            spread_vsite_f(vsite, x, fr->f_novirsum, NULL,
 +            spread_vsite_f(vsite, x, as_rvec_array(fr->f_novirsum->data()), nullptr,
                             (flags & GMX_FORCE_VIRIAL), fr->vir_el_recip,
                             nrnb,
                             &top->idef, fr->ePBC, fr->bMolPBC, graph, box, cr);
          if (flags & GMX_FORCE_VIRIAL)
          {
              /* Now add the forces, this is local */
 -            if (fr->bDomDec)
 -            {
 -                sum_forces(0, fr->f_novirsum_n, f, fr->f_novirsum);
 -            }
 -            else
 -            {
 -                sum_forces(0, mdatoms->homenr,
 -                           f, fr->f_novirsum);
 -            }
 +            sum_forces(f, fr->f_novirsum);
 +
              if (EEL_FULL(fr->eeltype))
              {
                  /* Add the mesh contribution to the virial */
@@@ -670,7 -733,7 +670,7 @@@ static void do_nb_verlet_fep(nbnxn_pair
  
  gmx_bool use_GPU(const nonbonded_verlet_t *nbv)
  {
 -    return nbv != NULL && nbv->bUseGPU;
 +    return nbv != nullptr && nbv->bUseGPU;
  }
  
  static gmx_inline void clear_rvecs_omp(int n, rvec v[])
@@@ -724,18 -787,19 +724,18 @@@ void do_force_cutsVERLET(FILE *fplog, t
                           gmx_localtop_t *top,
                           gmx_groups_t gmx_unused *groups,
                           matrix box, rvec x[], history_t *hist,
 -                         rvec f[],
 +                         PaddedRVecVector *force,
                           tensor vir_force,
                           t_mdatoms *mdatoms,
                           gmx_enerdata_t *enerd, t_fcdata *fcd,
                           real *lambda, t_graph *graph,
                           t_forcerec *fr, interaction_const_t *ic,
                           gmx_vsite_t *vsite, rvec mu_tot,
 -                         double t, FILE *field, gmx_edsam_t ed,
 +                         double t, gmx_edsam_t ed,
                           gmx_bool bBornRadii,
                           int flags)
  {
      int                 cg1, i, j;
 -    int                 start, homenr;
      double              mu[2*DIM];
      gmx_bool            bStateChanged, bNS, bFillGrid, bCalcCGCM;
      gmx_bool            bDoForces, bUseGPU, bUseOrEmulGPU;
      cycles_wait_gpu = 0;
      nbv             = fr->nbv;
  
 -    start  = 0;
 -    homenr = mdatoms->homenr;
 +    const int start  = 0;
 +    const int homenr = mdatoms->homenr;
  
      clear_mat(vir_force);
  
           * we do not need to worry about shifting.
           */
  
 -        int pme_flags = 0;
 -
          wallcycle_start(wcycle, ewcPP_PMESENDX);
  
          bBS = (inputrec->nwall == 2);
              svmul(inputrec->wall_ewald_zfac, boxs[ZZ], boxs[ZZ]);
          }
  
 -        if (EEL_PME(fr->eeltype))
 -        {
 -            pme_flags |= GMX_PME_DO_COULOMB;
 -        }
 -
 -        if (EVDW_PME(fr->vdwtype))
 -        {
 -            pme_flags |= GMX_PME_DO_LJ;
 -        }
 -
          gmx_pme_send_coordinates(cr, bBS ? boxs : box, x,
 -                                 mdatoms->nChargePerturbed, mdatoms->nTypePerturbed, lambda[efptCOUL], lambda[efptVDW],
 +                                 lambda[efptCOUL], lambda[efptVDW],
                                   (flags & (GMX_FORCE_VIRIAL | GMX_FORCE_ENERGY)),
 -                                 pme_flags, step);
 +                                 step);
  
          wallcycle_stop(wcycle, ewcPP_PMESENDX);
      }
              nbnxn_put_on_grid(nbv->nbs, fr->ePBC, box,
                                0, vzero, box_diag,
                                0, mdatoms->homenr, -1, fr->cginfo, x,
 -                              0, NULL,
 +                              0, nullptr,
                                nbv->grp[eintLocal].kernel_type,
                                nbv->grp[eintLocal].nbat);
              wallcycle_sub_stop(wcycle, ewcsNBS_GRID_LOCAL);
          wallcycle_stop(wcycle, ewcROT);
      }
  
 +    /* Temporary solution until all routines take PaddedRVecVector */
 +    rvec *f = as_rvec_array(force->data());
 +
      /* Start the force cycle counter.
       * This counter is stopped after do_force_lowlevel.
       * No parallel communication should occur while this counter is running,
          {
              if (flags & GMX_FORCE_VIRIAL)
              {
 -                fr->f_novirsum = fr->f_novirsum_alloc;
 +                fr->f_novirsum = fr->forceBufferNoVirialSummation;
              }
              else
              {
                   * a separate array for forces that do not contribute
                   * to the pressure.
                   */
 -                fr->f_novirsum = f;
 +                fr->f_novirsum = force;
              }
          }
  
          {
              if (flags & GMX_FORCE_VIRIAL)
              {
 -                if (fr->bDomDec)
 -                {
 -                    clear_rvecs_omp(fr->f_novirsum_n, fr->f_novirsum);
 -                }
 -                else
 -                {
 -                    clear_rvecs_omp(homenr, fr->f_novirsum+start);
 -                }
 +                /* TODO: remove this - 1 when padding is properly implemented */
 +                clear_rvecs_omp(fr->f_novirsum->size() - 1,
 +                                as_rvec_array(fr->f_novirsum->data()));
              }
          }
          /* Clear the short- and long-range forces */
  
      if (bDoForces)
      {
 -        if (inputrecElecField(inputrec))
 +        /* Compute forces due to electric field */
 +        if (fr->efield != nullptr)
          {
 -            /* Compute forces due to electric field */
 -            calc_f_el(MASTER(cr) ? field : NULL,
 -                      start, homenr, mdatoms->chargeA, fr->f_novirsum,
 -                      inputrec->ex, inputrec->et, t);
 +            fr->efield->calculateForces(cr, mdatoms, fr->f_novirsum, t);
          }
  
          /* If we have NoVirSum forces, but we do not calculate the virial,
          if (vsite && !(fr->bF_NoVirSum && !(flags & GMX_FORCE_VIRIAL)))
          {
              wallcycle_start(wcycle, ewcVSITESPREAD);
 -            spread_vsite_f(vsite, x, f, fr->fshift, FALSE, NULL, nrnb,
 +            spread_vsite_f(vsite, x, f, fr->fshift, FALSE, nullptr, nrnb,
                             &top->idef, fr->ePBC, fr->bMolPBC, graph, box, cr);
              wallcycle_stop(wcycle, ewcVSITESPREAD);
          }
@@@ -1466,24 -1546,25 +1466,24 @@@ void do_force_cutsGROUP(FILE *fplog, t_
                          gmx_localtop_t *top,
                          gmx_groups_t *groups,
                          matrix box, rvec x[], history_t *hist,
 -                        rvec f[],
 +                        PaddedRVecVector *force,
                          tensor vir_force,
                          t_mdatoms *mdatoms,
                          gmx_enerdata_t *enerd, t_fcdata *fcd,
                          real *lambda, t_graph *graph,
                          t_forcerec *fr, gmx_vsite_t *vsite, rvec mu_tot,
 -                        double t, FILE *field, gmx_edsam_t ed,
 +                        double t, gmx_edsam_t ed,
                          gmx_bool bBornRadii,
                          int flags)
  {
      int        cg0, cg1, i, j;
 -    int        start, homenr;
      double     mu[2*DIM];
      gmx_bool   bStateChanged, bNS, bFillGrid, bCalcCGCM;
      gmx_bool   bDoForces;
      float      cycles_pme, cycles_force;
  
 -    start  = 0;
 -    homenr = mdatoms->homenr;
 +    const int  start  = 0;
 +    const int  homenr = mdatoms->homenr;
  
      clear_mat(vir_force);
  
           * we do not need to worry about shifting.
           */
  
 -        int pme_flags = 0;
 -
          wallcycle_start(wcycle, ewcPP_PMESENDX);
  
          bBS = (inputrec->nwall == 2);
              svmul(inputrec->wall_ewald_zfac, boxs[ZZ], boxs[ZZ]);
          }
  
 -        if (EEL_PME(fr->eeltype))
 -        {
 -            pme_flags |= GMX_PME_DO_COULOMB;
 -        }
 -
 -        if (EVDW_PME(fr->vdwtype))
 -        {
 -            pme_flags |= GMX_PME_DO_LJ;
 -        }
 -
          gmx_pme_send_coordinates(cr, bBS ? boxs : box, x,
 -                                 mdatoms->nChargePerturbed, mdatoms->nTypePerturbed, lambda[efptCOUL], lambda[efptVDW],
 +                                 lambda[efptCOUL], lambda[efptVDW],
                                   (flags & (GMX_FORCE_VIRIAL | GMX_FORCE_ENERGY)),
 -                                 pme_flags, step);
 +                                 step);
  
          wallcycle_stop(wcycle, ewcPP_PMESENDX);
      }
          wallcycle_stop(wcycle, ewcROT);
      }
  
 +    /* Temporary solution until all routines take PaddedRVecVector */
 +    rvec *f = as_rvec_array(force->data());
 +
      /* Start the force cycle counter.
       * This counter is stopped after do_force_lowlevel.
       * No parallel communication should occur while this counter is running,
          {
              if (flags & GMX_FORCE_VIRIAL)
              {
 -                fr->f_novirsum = fr->f_novirsum_alloc;
 -                if (fr->bDomDec)
 -                {
 -                    clear_rvecs(fr->f_novirsum_n, fr->f_novirsum);
 -                }
 -                else
 -                {
 -                    clear_rvecs(homenr, fr->f_novirsum+start);
 -                }
 +                fr->f_novirsum = fr->forceBufferNoVirialSummation;
 +                /* TODO: remove this - 1 when padding is properly implemented */
 +                clear_rvecs(fr->f_novirsum->size() - 1,
 +                            as_rvec_array(fr->f_novirsum->data()));
              }
              else
              {
                   * a separate array for forces that do not contribute
                   * to the pressure.
                   */
 -                fr->f_novirsum = f;
 +                fr->f_novirsum = force;
              }
          }
  
  
      if (bDoForces)
      {
 -        if (inputrecElecField(inputrec))
 +        /* Compute forces due to electric field */
 +        if (fr->efield != nullptr)
          {
 -            /* Compute forces due to electric field */
 -            calc_f_el(MASTER(cr) ? field : NULL,
 -                      start, homenr, mdatoms->chargeA, fr->f_novirsum,
 -                      inputrec->ex, inputrec->et, t);
 +            fr->efield->calculateForces(cr, mdatoms, fr->f_novirsum, t);
          }
  
          /* Communicate the forces */
              if (EEL_FULL(fr->eeltype) && cr->dd->n_intercg_excl &&
                  (flags & GMX_FORCE_VIRIAL))
              {
 -                dd_move_f(cr->dd, fr->f_novirsum, NULL);
 +                dd_move_f(cr->dd, as_rvec_array(fr->f_novirsum->data()), nullptr);
              }
              wallcycle_stop(wcycle, ewcMOVEF);
          }
          if (vsite && !(fr->bF_NoVirSum && !(flags & GMX_FORCE_VIRIAL)))
          {
              wallcycle_start(wcycle, ewcVSITESPREAD);
 -            spread_vsite_f(vsite, x, f, fr->fshift, FALSE, NULL, nrnb,
 +            spread_vsite_f(vsite, x, f, fr->fshift, FALSE, nullptr, nrnb,
                             &top->idef, fr->ePBC, fr->bMolPBC, graph, box, cr);
              wallcycle_stop(wcycle, ewcVSITESPREAD);
          }
@@@ -1841,15 -1938,15 +1841,15 @@@ void do_force(FILE *fplog, t_commrec *c
                gmx_int64_t step, t_nrnb *nrnb, gmx_wallcycle_t wcycle,
                gmx_localtop_t *top,
                gmx_groups_t *groups,
 -              matrix box, rvec x[], history_t *hist,
 -              rvec f[],
 +              matrix box, PaddedRVecVector *coordinates, history_t *hist,
 +              PaddedRVecVector *force,
                tensor vir_force,
                t_mdatoms *mdatoms,
                gmx_enerdata_t *enerd, t_fcdata *fcd,
 -              real *lambda, t_graph *graph,
 +              gmx::ArrayRef<real> lambda, t_graph *graph,
                t_forcerec *fr,
                gmx_vsite_t *vsite, rvec mu_tot,
 -              double t, FILE *field, gmx_edsam_t ed,
 +              double t, gmx_edsam_t ed,
                gmx_bool bBornRadii,
                int flags)
  {
          flags &= ~GMX_FORCE_NONBONDED;
      }
  
 +    GMX_ASSERT(coordinates->size() >= static_cast<unsigned int>(fr->natoms_force + 1) ||
 +               fr->natoms_force == 0, "We might need 1 element extra for SIMD");
 +    GMX_ASSERT(force->size() >= static_cast<unsigned int>(fr->natoms_force + 1) ||
 +               fr->natoms_force == 0, "We might need 1 element extra for SIMD");
 +
 +    rvec *x = as_rvec_array(coordinates->data());
 +
      switch (inputrec->cutoff_scheme)
      {
          case ecutsVERLET:
                                  top,
                                  groups,
                                  box, x, hist,
 -                                f, vir_force,
 +                                force, vir_force,
                                  mdatoms,
                                  enerd, fcd,
 -                                lambda, graph,
 +                                lambda.data(), graph,
                                  fr, fr->ic,
                                  vsite, mu_tot,
 -                                t, field, ed,
 +                                t, ed,
                                  bBornRadii,
                                  flags);
              break;
                                 top,
                                 groups,
                                 box, x, hist,
 -                               f, vir_force,
 +                               force, vir_force,
                                 mdatoms,
                                 enerd, fcd,
 -                               lambda, graph,
 +                               lambda.data(), graph,
                                 fr, vsite, mu_tot,
 -                               t, field, ed,
 +                               t, ed,
                                 bBornRadii,
                                 flags);
              break;
@@@ -1940,22 -2030,22 +1940,22 @@@ void do_constrain_first(FILE *fplog, gm
      dvdl_dum = 0;
  
      /* constrain the current position */
 -    constrain(NULL, TRUE, FALSE, constr, &(top->idef),
 +    constrain(nullptr, TRUE, FALSE, constr, &(top->idef),
                ir, cr, step, 0, 1.0, md,
 -              state->x, state->x, NULL,
 +              as_rvec_array(state->x.data()), as_rvec_array(state->x.data()), nullptr,
                fr->bMolPBC, state->box,
                state->lambda[efptBONDED], &dvdl_dum,
 -              NULL, NULL, nrnb, econqCoord);
 +              nullptr, nullptr, nrnb, econqCoord);
      if (EI_VV(ir->eI))
      {
          /* constrain the inital velocity, and save it */
          /* also may be useful if we need the ekin from the halfstep for velocity verlet */
 -        constrain(NULL, TRUE, FALSE, constr, &(top->idef),
 +        constrain(nullptr, TRUE, FALSE, constr, &(top->idef),
                    ir, cr, step, 0, 1.0, md,
 -                  state->x, state->v, state->v,
 +                  as_rvec_array(state->x.data()), as_rvec_array(state->v.data()), as_rvec_array(state->v.data()),
                    fr->bMolPBC, state->box,
                    state->lambda[efptBONDED], &dvdl_dum,
 -                  NULL, NULL, nrnb, econqVeloc);
 +                  nullptr, nullptr, nrnb, econqVeloc);
      }
      /* constrain the inital velocities at t-dt/2 */
      if (EI_STATE_VELOCITY(ir->eI) && ir->eI != eiVV)
                      gmx_step_str(step, buf));
          }
          dvdl_dum = 0;
 -        constrain(NULL, TRUE, FALSE, constr, &(top->idef),
 +        constrain(nullptr, TRUE, FALSE, constr, &(top->idef),
                    ir, cr, step, -1, 1.0, md,
 -                  state->x, savex, NULL,
 +                  as_rvec_array(state->x.data()), savex, nullptr,
                    fr->bMolPBC, state->box,
                    state->lambda[efptBONDED], &dvdl_dum,
 -                  state->v, NULL, nrnb, econqCoord);
 +                  as_rvec_array(state->v.data()), nullptr, nrnb, econqCoord);
  
          for (i = start; i < end; i++)
          {
@@@ -2416,7 -2506,7 +2416,7 @@@ static void low_do_pbc_mtop(FILE *fplog
          else
          {
              /* Pass NULL iso fplog to avoid graph prints for each molecule type */
 -            mk_graph_ilist(NULL, mtop->moltype[molb->type].ilist,
 +            mk_graph_ilist(nullptr, mtop->moltype[molb->type].ilist,
                             0, molb->natoms_mol, FALSE, FALSE, graph);
  
              for (mol = 0; mol < molb->nmol; mol++)
@@@ -2470,25 -2560,36 +2470,35 @@@ void put_atoms_in_box_omp(int ePBC, con
  }
  
  // TODO This can be cleaned up a lot, and move back to runner.cpp
 -void finish_run(FILE *fplog, t_commrec *cr,
 +void finish_run(FILE *fplog, const gmx::MDLogger &mdlog, t_commrec *cr,
                  t_inputrec *inputrec,
                  t_nrnb nrnb[], gmx_wallcycle_t wcycle,
                  gmx_walltime_accounting_t walltime_accounting,
                  nonbonded_verlet_t *nbv,
                  gmx_bool bWriteStat)
  {
 -    t_nrnb *nrnb_tot = NULL;
 +    t_nrnb *nrnb_tot = nullptr;
      double  delta_t  = 0;
      double  nbfs     = 0, mflop = 0;
      double  elapsed_time,
              elapsed_time_over_all_ranks,
              elapsed_time_over_all_threads,
              elapsed_time_over_all_threads_over_all_ranks;
+     /* Control whether it is valid to print a report. Only the
+        simulation master may print, but it should not do so if the run
+        terminated e.g. before a scheduled reset step. This is
+        complicated by the fact that PME ranks are unaware of the
+        reason why they were sent a pmerecvqxFINISH. To avoid
+        communication deadlocks, we always do the communication for the
+        report, even if we've decided not to write the report, because
+        how long it takes to finish the run is not important when we've
+        decided not to report on the simulation performance. */
+     bool    printReport = SIMMASTER(cr);
  
      if (!walltime_accounting_get_valid_finish(walltime_accounting))
      {
 -        md_print_warn(cr, fplog,
 -                      "Simulation ended prematurely, no performance report will be written.");
 +        GMX_LOG(mdlog.warning).asParagraph().appendText("Simulation ended prematurely, no performance report will be written.");
-         return;
+         printReport = false;
      }
  
      if (cr->nnodes > 1)
      }
  #endif
  
-     if (SIMMASTER(cr))
+     if (printReport)
      {
          print_flop(fplog, nrnb_tot, &nbfs, &mflop);
      }
      wallcycle_scale_by_num_threads(wcycle, cr->duty == DUTY_PME, nthreads_pp, nthreads_pme);
      auto cycle_sum(wallcycle_sum(cr, wcycle));
  
-     if (SIMMASTER(cr))
+     if (printReport)
      {
 -        struct gmx_wallclock_gpu_t* gputimes = use_GPU(nbv) ? nbnxn_gpu_get_timings(nbv->gpu_nbv) : NULL;
 +        struct gmx_wallclock_gpu_t* gputimes = use_GPU(nbv) ? nbnxn_gpu_get_timings(nbv->gpu_nbv) : nullptr;
  
 -        wallcycle_print(fplog, cr->nnodes, cr->npmenodes, nthreads_pp, nthreads_pme,
 +        wallcycle_print(fplog, mdlog, cr->nnodes, cr->npmenodes, nthreads_pp, nthreads_pme,
                          elapsed_time_over_all_ranks,
                          wcycle, cycle_sum, gputimes);
  
      }
  }
  
 -extern void initialize_lambdas(FILE *fplog, t_inputrec *ir, int *fep_state, real *lambda, double *lam0)
 +extern void initialize_lambdas(FILE *fplog, t_inputrec *ir, int *fep_state, gmx::ArrayRef<real> lambda, double *lam0)
  {
      /* this function works, but could probably use a logic rewrite to keep all the different
         types of efep straight. */
  
 -    int       i;
 +    if ((ir->efep == efepNO) && (ir->bSimTemp == FALSE))
 +    {
 +        return;
 +    }
 +
      t_lambda *fep = ir->fepvals;
 +    *fep_state    = fep->init_fep_state; /* this might overwrite the checkpoint
 +                                            if checkpoint is set -- a kludge is in for now
 +                                            to prevent this.*/
  
 -    if ((ir->efep == efepNO) && (ir->bSimTemp == FALSE))
 +    for (int i = 0; i < efptNR; i++)
      {
 -        for (i = 0; i < efptNR; i++)
 +        /* overwrite lambda state with init_lambda for now for backwards compatibility */
 +        if (fep->init_lambda >= 0) /* if it's -1, it was never initializd */
          {
 -            lambda[i] = 0.0;
 +            lambda[i] = fep->init_lambda;
              if (lam0)
              {
 -                lam0[i] = 0.0;
 +                lam0[i] = lambda[i];
              }
          }
 -        return;
 -    }
 -    else
 -    {
 -        *fep_state = fep->init_fep_state; /* this might overwrite the checkpoint
 -                                             if checkpoint is set -- a kludge is in for now
 -                                             to prevent this.*/
 -        for (i = 0; i < efptNR; i++)
 +        else
          {
 -            /* overwrite lambda state with init_lambda for now for backwards compatibility */
 -            if (fep->init_lambda >= 0) /* if it's -1, it was never initializd */
 -            {
 -                lambda[i] = fep->init_lambda;
 -                if (lam0)
 -                {
 -                    lam0[i] = lambda[i];
 -                }
 -            }
 -            else
 +            lambda[i] = fep->all_lambda[i][*fep_state];
 +            if (lam0)
              {
 -                lambda[i] = fep->all_lambda[i][*fep_state];
 -                if (lam0)
 -                {
 -                    lam0[i] = lambda[i];
 -                }
 +                lam0[i] = lambda[i];
              }
          }
 -        if (ir->bSimTemp)
 +    }
 +    if (ir->bSimTemp)
 +    {
 +        /* need to rescale control temperatures to match current state */
 +        for (int i = 0; i < ir->opts.ngtc; i++)
          {
 -            /* need to rescale control temperatures to match current state */
 -            for (i = 0; i < ir->opts.ngtc; i++)
 +            if (ir->opts.ref_t[i] > 0)
              {
 -                if (ir->opts.ref_t[i] > 0)
 -                {
 -                    ir->opts.ref_t[i] = ir->simtempvals->temperatures[*fep_state];
 -                }
 +                ir->opts.ref_t[i] = ir->simtempvals->temperatures[*fep_state];
              }
          }
      }
  
      /* Send to the log the information on the current lambdas */
 -    if (fplog != NULL)
 +    if (fplog != nullptr)
      {
          fprintf(fplog, "Initial vector of lambda components:[ ");
 -        for (i = 0; i < efptNR; i++)
 +        for (const auto &l : lambda)
          {
 -            fprintf(fplog, "%10.4f ", lambda[i]);
 +            fprintf(fplog, "%10.4f ", l);
          }
          fprintf(fplog, "]\n");
      }
  
  
  void init_md(FILE *fplog,
 -             t_commrec *cr, t_inputrec *ir, const gmx_output_env_t *oenv,
 +             t_commrec *cr, gmx::IMDOutputProvider *outputProvider,
 +             t_inputrec *ir, const gmx_output_env_t *oenv,
               double *t, double *t0,
 -             real *lambda, int *fep_state, double *lam0,
 +             gmx::ArrayRef<real> lambda, int *fep_state, double *lam0,
               t_nrnb *nrnb, gmx_mtop_t *mtop,
               gmx_update_t **upd,
               int nfile, const t_filenm fnm[],
      }
      if (*bSimAnn)
      {
 -        update_annealing_target_temp(ir, ir->init_t, upd ? *upd : NULL);
 +        update_annealing_target_temp(ir, ir->init_t, upd ? *upd : nullptr);
      }
  
 -    if (vcm != NULL)
 +    if (vcm != nullptr)
      {
          *vcm = init_vcm(fplog, &mtop->groups, ir);
      }
              please_cite(fplog, "Goga2012");
          }
      }
 -    if ((ir->et[XX].n > 0) || (ir->et[YY].n > 0) || (ir->et[ZZ].n > 0))
 -    {
 -        please_cite(fplog, "Caleman2008a");
 -    }
      init_nrnb(nrnb);
  
      if (nfile != -1)
      {
 -        *outf = init_mdoutf(fplog, nfile, fnm, Flags, cr, ir, mtop, oenv, wcycle);
 +        *outf = init_mdoutf(fplog, nfile, fnm, Flags, cr, outputProvider, ir, mtop, oenv, wcycle);
  
 -        *mdebin = init_mdebin((Flags & MD_APPENDFILES) ? NULL : mdoutf_get_fp_ene(*outf),
 +        *mdebin = init_mdebin((Flags & MD_APPENDFILES) ? nullptr : mdoutf_get_fp_ene(*outf),
                                mtop, ir, mdoutf_get_fp_dhdl(*outf));
      }
  
index 6e9e627f12b92db476c2a59e0d3dffffa9c87bb8,9feb963c50ba81e5c961bdc5103896880830e85b..f805b1f003ee681b9dfef83d2088888755056fb2
@@@ -68,7 -68,6 +68,7 @@@
  #include <cstring>
  
  #include <algorithm>
 +#include <string>
  
  /* This file is completely threadsafe - keep it that way! */
  
  #include "gromacs/utility/path.h"
  #include "gromacs/utility/programcontext.h"
  #include "gromacs/utility/stringutil.h"
 +#include "gromacs/utility/textwriter.h"
  
  #include "cuda_version_information.h"
  
  namespace
  {
  
 +using gmx::formatString;
 +
  //! \cond Doxygen does not need to care about most of this stuff, and the macro usage is painful to document
  
  int centeringOffset(int width, int length)
      return std::max(width - length, 0) / 2;
  }
  
 -void printCentered(FILE *fp, int width, const char *text)
 +std::string formatCentered(int width, const char *text)
  {
      const int offset = centeringOffset(width, std::strlen(text));
 -    fprintf(fp, "%*s%s", offset, "", text);
 +    return formatString("%*s%s", offset, "", text);
  }
  
 -void printCopyright(FILE *fp)
 +void printCopyright(gmx::TextWriter *writer)
  {
      static const char * const Contributors[] = {
          "Emile Apol",
      };
      static const char * const CopyrightText[] = {
          "Copyright (c) 1991-2000, University of Groningen, The Netherlands.",
-         "Copyright (c) 2001-2015, The GROMACS development team at",
+         "Copyright (c) 2001-2017, The GROMACS development team at",
          "Uppsala University, Stockholm University and",
          "the Royal Institute of Technology, Sweden.",
          "check out http://www.gromacs.org for more information."
  #define NLICENSE (int)asize(LicenseText)
  #endif
  
 -    printCentered(fp, 78, "GROMACS is written by:");
 -    fprintf(fp, "\n");
 +    // TODO a centering behaviour of TextWriter could be useful here
 +    writer->writeLine(formatCentered(78, "GROMACS is written by:"));
      for (int i = 0; i < NCONTRIBUTORS; )
      {
          for (int j = 0; j < 4 && i < NCONTRIBUTORS; ++j, ++i)
                                 "Formatting buffer is not long enough");
              std::fill(buf, buf+width, ' ');
              std::strcpy(buf+offset, Contributors[i]);
 -            fprintf(fp, " %-*s", width, buf);
 +            writer->writeString(formatString(" %-*s", width, buf));
          }
 -        fprintf(fp, "\n");
 +        writer->ensureLineBreak();
      }
 -    printCentered(fp, 78, "and the project leaders:");
 -    fprintf(fp, "\n");
 -    printCentered(fp, 78, "Mark Abraham, Berk Hess, Erik Lindahl, and David van der Spoel");
 -    fprintf(fp, "\n\n");
 +    writer->writeLine(formatCentered(78, "and the project leaders:"));
 +    writer->writeLine(formatCentered(78, "Mark Abraham, Berk Hess, Erik Lindahl, and David van der Spoel"));
 +    writer->ensureEmptyLine();
      for (int i = 0; i < NCR; ++i)
      {
 -        fprintf(fp, "%s\n", CopyrightText[i]);
 +        writer->writeLine(CopyrightText[i]);
      }
 -    fprintf(fp, "\n");
 +    writer->ensureEmptyLine();
      for (int i = 0; i < NLICENSE; ++i)
      {
 -        fprintf(fp, "%s\n", LicenseText[i]);
 +        writer->writeLine(LicenseText[i]);
      }
  }
  
@@@ -217,101 -214,96 +217,101 @@@ const char *getFftDescriptionString(
  #endif
  };
  
 -void gmx_print_version_info(FILE *fp)
 +void gmx_print_version_info(gmx::TextWriter *writer)
  {
 -    fprintf(fp, "GROMACS version:    %s\n", gmx_version());
 +    writer->writeLine(formatString("GROMACS version:    %s", gmx_version()));
      const char *const git_hash = gmx_version_git_full_hash();
      if (git_hash[0] != '\0')
      {
 -        fprintf(fp, "GIT SHA1 hash:      %s\n", git_hash);
 +        writer->writeLine(formatString("GIT SHA1 hash:      %s", git_hash));
      }
      const char *const base_hash = gmx_version_git_central_base_hash();
      if (base_hash[0] != '\0')
      {
 -        fprintf(fp, "Branched from:      %s\n", base_hash);
 +        writer->writeLine(formatString("Branched from:      %s", base_hash));
      }
  
  #if GMX_DOUBLE
 -    fprintf(fp, "Precision:          double\n");
 +    writer->writeLine("Precision:          double");
  #else
 -    fprintf(fp, "Precision:          single\n");
 +    writer->writeLine("Precision:          single");
  #endif
 -    fprintf(fp, "Memory model:       %u bit\n", (unsigned)(8*sizeof(void *)));
 +    writer->writeLine(formatString("Memory model:       %u bit", (unsigned)(8*sizeof(void *))));
  
  #if GMX_THREAD_MPI
 -    fprintf(fp, "MPI library:        thread_mpi\n");
 +    writer->writeLine("MPI library:        thread_mpi");
  #elif GMX_MPI
 -    fprintf(fp, "MPI library:        MPI\n");
 +    writer->writeLine("MPI library:        MPI");
  #else
 -    fprintf(fp, "MPI library:        none\n");
 +    writer->writeLine("MPI library:        none");
  #endif
  #if GMX_OPENMP
 -    fprintf(fp, "OpenMP support:     enabled (GMX_OPENMP_MAX_THREADS = %d)\n", GMX_OPENMP_MAX_THREADS);
 +    writer->writeLine(formatString("OpenMP support:     enabled (GMX_OPENMP_MAX_THREADS = %d)", GMX_OPENMP_MAX_THREADS));
  #else
 -    fprintf(fp, "OpenMP support:     disabled\n");
 +    writer->writeLine("OpenMP support:     disabled");
  #endif
 -    fprintf(fp, "GPU support:        %s\n", getGpuImplementationString());
 -    fprintf(fp, "SIMD instructions:  %s\n", GMX_SIMD_STRING);
 -    fprintf(fp, "FFT library:        %s\n", getFftDescriptionString());
 +    writer->writeLine(formatString("GPU support:        %s", getGpuImplementationString()));
 +    writer->writeLine(formatString("SIMD instructions:  %s", GMX_SIMD_STRING));
 +    writer->writeLine(formatString("FFT library:        %s", getFftDescriptionString()));
  #ifdef HAVE_RDTSCP
 -    fprintf(fp, "RDTSCP usage:       enabled\n");
 +    writer->writeLine("RDTSCP usage:       enabled");
  #else
 -    fprintf(fp, "RDTSCP usage:       disabled\n");
 +    writer->writeLine("RDTSCP usage:       disabled");
  #endif
  #ifdef GMX_USE_TNG
 -    fprintf(fp, "TNG support:        enabled\n");
 +    writer->writeLine("TNG support:        enabled");
  #else
 -    fprintf(fp, "TNG support:        disabled\n");
 +    writer->writeLine("TNG support:        disabled");
  #endif
  #if GMX_HWLOC
 -    fprintf(fp, "Hwloc support:      hwloc-%d.%d.%d\n",
 -            HWLOC_API_VERSION>>16,
 -            (HWLOC_API_VERSION>>8) & 0xFF,
 -            HWLOC_API_VERSION & 0xFF);
 +    writer->writeLine(formatString("Hwloc support:      hwloc-%d.%d.%d",
 +                                   HWLOC_API_VERSION>>16,
 +                                   (HWLOC_API_VERSION>>8) & 0xFF,
 +                                   HWLOC_API_VERSION & 0xFF));
  #else
 -    fprintf(fp, "Hwloc support:      disabled\n");
 +    writer->writeLine("Hwloc support:      disabled");
  #endif
  #if HAVE_EXTRAE
      unsigned major, minor, revision;
      Extrae_get_version(&major, &minor, &revision);
 -    fprintf(fp, "Tracing support:    enabled. Using Extrae-%d.%d.%d\n", major, minor, revision);
 +    writer->writeLine(formatString("Tracing support:    enabled. Using Extrae-%d.%d.%d", major, minor, revision));
  #else
 -    fprintf(fp, "Tracing support:    disabled\n");
 +    writer->writeLine("Tracing support:    disabled");
  #endif
  
  
 -    fprintf(fp, "Built on:           %s\n", BUILD_TIME);
 -    fprintf(fp, "Built by:           %s\n", BUILD_USER);
 -    fprintf(fp, "Build OS/arch:      %s\n", BUILD_HOST);
 -    fprintf(fp, "Build CPU vendor:   %s\n", BUILD_CPU_VENDOR);
 -    fprintf(fp, "Build CPU brand:    %s\n", BUILD_CPU_BRAND);
 -    fprintf(fp, "Build CPU family:   %d   Model: %d   Stepping: %d\n",
 -            BUILD_CPU_FAMILY, BUILD_CPU_MODEL, BUILD_CPU_STEPPING);
 +    writer->writeLine(formatString("Built on:           %s", BUILD_TIME));
 +    writer->writeLine(formatString("Built by:           %s", BUILD_USER));
 +    writer->writeLine(formatString("Build OS/arch:      %s", BUILD_HOST));
 +    writer->writeLine(formatString("Build CPU vendor:   %s", BUILD_CPU_VENDOR));
 +    writer->writeLine(formatString("Build CPU brand:    %s", BUILD_CPU_BRAND));
 +    writer->writeLine(formatString("Build CPU family:   %d   Model: %d   Stepping: %d",
 +                                   BUILD_CPU_FAMILY, BUILD_CPU_MODEL, BUILD_CPU_STEPPING));
      /* TODO: The below strings can be quite long, so it would be nice to wrap
       * them. Can wait for later, as the master branch has ready code to do all
       * that. */
 -    fprintf(fp, "Build CPU features: %s\n", BUILD_CPU_FEATURES);
 -    fprintf(fp, "C compiler:         %s\n", BUILD_C_COMPILER);
 -    fprintf(fp, "C compiler flags:   %s\n", BUILD_CFLAGS);
 -    fprintf(fp, "C++ compiler:       %s\n", BUILD_CXX_COMPILER);
 -    fprintf(fp, "C++ compiler flags: %s\n", BUILD_CXXFLAGS);
 +    writer->writeLine(formatString("Build CPU features: %s", BUILD_CPU_FEATURES));
 +    writer->writeLine(formatString("C compiler:         %s", BUILD_C_COMPILER));
 +    writer->writeLine(formatString("C compiler flags:   %s", BUILD_CFLAGS));
 +    writer->writeLine(formatString("C++ compiler:       %s", BUILD_CXX_COMPILER));
 +    writer->writeLine(formatString("C++ compiler flags: %s", BUILD_CXXFLAGS));
  #ifdef HAVE_LIBMKL
      /* MKL might be used for LAPACK/BLAS even if FFTs use FFTW, so keep it separate */
 -    fprintf(fp, "Linked with Intel MKL version %d.%d.%d.\n",
 -            __INTEL_MKL__, __INTEL_MKL_MINOR__, __INTEL_MKL_UPDATE__);
 +    writer->writeLine(formatString("Linked with Intel MKL version %d.%d.%d.",
 +                                   __INTEL_MKL__, __INTEL_MKL_MINOR__, __INTEL_MKL_UPDATE__));
  #endif
  #if GMX_GPU == GMX_GPU_OPENCL
 -    fprintf(fp, "OpenCL include dir: %s\n", OPENCL_INCLUDE_DIR);
 -    fprintf(fp, "OpenCL library:     %s\n", OPENCL_LIBRARY);
 -    fprintf(fp, "OpenCL version:     %s\n", OPENCL_VERSION_STRING);
 +    writer->writeLine(formatString("OpenCL include dir: %s", OPENCL_INCLUDE_DIR));
 +    writer->writeLine(formatString("OpenCL library:     %s", OPENCL_LIBRARY));
 +    writer->writeLine(formatString("OpenCL version:     %s", OPENCL_VERSION_STRING));
  #endif
  #if GMX_GPU == GMX_GPU_CUDA
 -    gmx_print_version_info_cuda_gpu(fp);
 +    writer->writeLine(formatString("CUDA compiler:      %s\n", CUDA_NVCC_COMPILER_INFO));
 +    writer->writeLine(formatString("CUDA compiler flags:%s\n", CUDA_NVCC_COMPILER_FLAGS));
 +    auto driverVersion = gmx::getCudaDriverVersion();
 +    writer->writeLine(formatString("CUDA driver:        %d.%d\n", driverVersion.first, driverVersion.second));
 +    auto runtimeVersion = gmx::getCudaRuntimeVersion();
 +    writer->writeLine(formatString("CUDA runtime:       %d.%d\n", runtimeVersion.first, runtimeVersion.second));
  #endif
  }
  
@@@ -331,28 -323,13 +331,28 @@@ BinaryInformationSettings::BinaryInform
  void printBinaryInformation(FILE                  *fp,
                              const IProgramContext &programContext)
  {
 -    printBinaryInformation(fp, programContext, BinaryInformationSettings());
 +    TextWriter writer(fp);
 +    printBinaryInformation(&writer, programContext, BinaryInformationSettings());
  }
  
  void printBinaryInformation(FILE                            *fp,
                              const IProgramContext           &programContext,
                              const BinaryInformationSettings &settings)
  {
 +    try
 +    {
 +        TextWriter writer(fp);
 +        printBinaryInformation(&writer, programContext, settings);
 +    }
 +    GMX_CATCH_ALL_AND_EXIT_WITH_FATAL_ERROR;
 +}
 +
 +void printBinaryInformation(TextWriter                      *writer,
 +                            const IProgramContext           &programContext,
 +                            const BinaryInformationSettings &settings)
 +{
 +    // TODO Perhaps the writer could be configured with the prefix and
 +    // suffix strings from the settings?
      const char *prefix          = settings.prefix_;
      const char *suffix          = settings.suffix_;
      const char *precisionString = "";
      const char *const name = programContext.displayName();
      if (settings.bGeneratedByHeader_)
      {
 -        fprintf(fp, "%sCreated by:%s\n", prefix, suffix);
 +        writer->writeLine(formatString("%sCreated by:%s", prefix, suffix));
      }
      // TODO: It would be nice to know here whether we are really running a
      // Gromacs binary or some other binary that is calling Gromacs; we
          = formatString(":-) GROMACS - %s, %s%s (-:", name, gmx_version(), precisionString);
      const int   indent
          = centeringOffset(78 - std::strlen(prefix) - std::strlen(suffix), title.length()) + 1;
 -    fprintf(fp, "%s%*c%s%s\n", prefix, indent, ' ', title.c_str(), suffix);
 -    fprintf(fp, "%s%s\n", prefix, suffix);
 +    writer->writeLine(formatString("%s%*c%s%s", prefix, indent, ' ', title.c_str(), suffix));
 +    writer->writeLine(formatString("%s%s", prefix, suffix));
      if (settings.bCopyright_)
      {
          GMX_RELEASE_ASSERT(prefix[0] == '\0' && suffix[0] == '\0',
                             "Prefix/suffix not supported with copyright");
 -        printCopyright(fp);
 -        fprintf(fp, "\n");
 +        printCopyright(writer);
 +        writer->ensureEmptyLine();
          // This line is printed again after the copyright notice to make it
          // appear together with all the other information, so that it is not
          // necessary to read stuff above the copyright notice.
          // The line above the copyright notice puts the copyright notice is
          // context, though.
 -        fprintf(fp, "%sGROMACS:      %s, version %s%s%s\n", prefix, name,
 -                gmx_version(), precisionString, suffix);
 +        writer->writeLine(formatString("%sGROMACS:      %s, version %s%s%s", prefix, name,
 +                                       gmx_version(), precisionString, suffix));
      }
      const char *const binaryPath = programContext.fullBinaryPath();
      if (!gmx::isNullOrEmpty(binaryPath))
      {
 -        fprintf(fp, "%sExecutable:   %s%s\n", prefix, binaryPath, suffix);
 +        writer->writeLine(formatString("%sExecutable:   %s%s", prefix, binaryPath, suffix));
      }
      const gmx::InstallationPrefixInfo installPrefix = programContext.installationPrefix();
      if (!gmx::isNullOrEmpty(installPrefix.path))
      {
 -        fprintf(fp, "%sData prefix:  %s%s%s\n", prefix, installPrefix.path,
 -                installPrefix.bSourceLayout ? " (source tree)" : "", suffix);
 +        writer->writeLine(formatString("%sData prefix:  %s%s%s", prefix, installPrefix.path,
 +                                       installPrefix.bSourceLayout ? " (source tree)" : "", suffix));
      }
      const std::string workingDir = Path::getWorkingDirectory();
      if (!workingDir.empty())
      {
 -        fprintf(fp, "%sWorking dir:  %s%s\n", prefix, workingDir.c_str(), suffix);
 +        writer->writeLine(formatString("%sWorking dir:  %s%s", prefix, workingDir.c_str(), suffix));
      }
      const char *const commandLine = programContext.commandLine();
      if (!gmx::isNullOrEmpty(commandLine))
      {
 -        fprintf(fp, "%sCommand line:%s\n%s  %s%s\n",
 -                prefix, suffix, prefix, commandLine, suffix);
 +        writer->writeLine(formatString("%sCommand line:%s\n%s  %s%s",
 +                                       prefix, suffix, prefix, commandLine, suffix));
      }
      if (settings.bExtendedInfo_)
      {
          GMX_RELEASE_ASSERT(prefix[0] == '\0' && suffix[0] == '\0',
                             "Prefix/suffix not supported with extended info");
 -        fprintf(fp, "\n");
 -        gmx_print_version_info(fp);
 +        writer->ensureEmptyLine();
 +        gmx_print_version_info(writer);
      }
  }
  
index 2e58d0a7660190baa6d93f4e3a8281c31a03290f,8cf13863611a9927aea73f5f5e7d6303901a20da..0a52bbe3017c56414e301ec61a1d36d5d6ee39b8
@@@ -45,7 -45,6 +45,7 @@@
  #include <stdlib.h>
  
  #include <algorithm>
 +#include <memory>
  
  #include "thread_mpi/threads.h"
  
@@@ -56,6 -55,7 +56,6 @@@
  #include "gromacs/ewald/pme.h"
  #include "gromacs/ewald/pme-load-balancing.h"
  #include "gromacs/fileio/trxio.h"
 -#include "gromacs/gmxlib/md_logging.h"
  #include "gromacs/gmxlib/network.h"
  #include "gromacs/gmxlib/nrnb.h"
  #include "gromacs/gpu_utils/gpu_utils.h"
@@@ -76,7 -76,6 +76,7 @@@
  #include "gromacs/mdlib/mdebin.h"
  #include "gromacs/mdlib/mdoutf.h"
  #include "gromacs/mdlib/mdrun.h"
 +#include "gromacs/mdlib/mdsetup.h"
  #include "gromacs/mdlib/nb_verlet.h"
  #include "gromacs/mdlib/nbnxn_gpu_data_mgmt.h"
  #include "gromacs/mdlib/ns.h"
  #include "gromacs/utility/basedefinitions.h"
  #include "gromacs/utility/cstringutil.h"
  #include "gromacs/utility/fatalerror.h"
 +#include "gromacs/utility/logger.h"
  #include "gromacs/utility/real.h"
  #include "gromacs/utility/smalloc.h"
  
@@@ -155,7 -153,7 +155,7 @@@ static void checkNumberOfBondedInteract
      }
  }
  
 -static void reset_all_counters(FILE *fplog, t_commrec *cr,
 +static void reset_all_counters(FILE *fplog, const gmx::MDLogger &mdlog, t_commrec *cr,
                                 gmx_int64_t step,
                                 gmx_int64_t *step_rel, t_inputrec *ir,
                                 gmx_wallcycle_t wcycle, t_nrnb *nrnb,
      char sbuf[STEPSTRSIZE];
  
      /* Reset all the counters related to performance over the run */
 -    md_print_warn(cr, fplog, "step %s: resetting all time and cycle counters\n",
 -                  gmx_step_str(step, sbuf));
 +    GMX_LOG(mdlog.warning).asParagraph().appendTextFormatted(
 +            "step %s: resetting all time and cycle counters",
 +            gmx_step_str(step, sbuf));
  
      if (use_GPU(nbv))
      {
  }
  
  /*! \libinternal
 -    \copydoc integrator_t (FILE *fplog, t_commrec *cr,
 +    \copydoc integrator_t (FILE *fplog, t_commrec *cr, const gmx::MDLogger &mdlog,
                             int nfile, const t_filenm fnm[],
                             const gmx_output_env_t *oenv, gmx_bool bVerbose,
                             int nstglobalcomm,
                             gmx_vsite_t *vsite, gmx_constr_t constr,
                             int stepout,
 +                           gmx::IMDOutputProvider *outputProvider,
                             t_inputrec *inputrec,
                             gmx_mtop_t *top_global, t_fcdata *fcd,
                             t_state *state_global,
                             unsigned long Flags,
                             gmx_walltime_accounting_t walltime_accounting)
   */
 -double gmx::do_md(FILE *fplog, t_commrec *cr, int nfile, const t_filenm fnm[],
 +double gmx::do_md(FILE *fplog, t_commrec *cr, const gmx::MDLogger &mdlog,
 +                  int nfile, const t_filenm fnm[],
                    const gmx_output_env_t *oenv, gmx_bool bVerbose,
                    int nstglobalcomm,
                    gmx_vsite_t *vsite, gmx_constr_t constr,
 -                  int stepout, t_inputrec *ir,
 +                  int stepout, gmx::IMDOutputProvider *outputProvider,
 +                  t_inputrec *ir,
                    gmx_mtop_t *top_global,
                    t_fcdata *fcd,
                    t_state *state_global,
 +                  energyhistory_t *energyHistory,
                    t_mdatoms *mdatoms,
                    t_nrnb *nrnb, gmx_wallcycle_t wcycle,
                    gmx_edsam_t ed, t_forcerec *fr,
                    unsigned long Flags,
                    gmx_walltime_accounting_t walltime_accounting)
  {
 -    gmx_mdoutf_t    outf = NULL;
 +    gmx_mdoutf_t    outf = nullptr;
      gmx_int64_t     step, step_rel;
      double          elapsed_time;
      double          t, t0, lam0[efptNR];
      t_trxstatus      *status;
      rvec              mu_tot;
      t_vcm            *vcm;
 -    matrix            pcoupl_mu, M;
 +    matrix            parrinellorahmanMu, M;
      t_trxframe        rerun_fr;
 -    gmx_repl_ex_t     repl_ex = NULL;
 +    gmx_repl_ex_t     repl_ex = nullptr;
      int               nchkpt  = 1;
      gmx_localtop_t   *top;
 -    t_mdebin         *mdebin   = NULL;
 -    t_state          *state    = NULL;
 +    t_mdebin         *mdebin   = nullptr;
      gmx_enerdata_t   *enerd;
 -    rvec             *f = NULL;
 +    PaddedRVecVector  f {};
      gmx_global_stat_t gstat;
 -    gmx_update_t     *upd   = NULL;
 -    t_graph          *graph = NULL;
 +    gmx_update_t     *upd   = nullptr;
 +    t_graph          *graph = nullptr;
      gmx_groups_t     *groups;
      gmx_ekindata_t   *ekind;
      gmx_shellfc_t    *shellfc;
      gmx_bool          bResetCountersHalfMaxH = FALSE;
      gmx_bool          bTemp, bPres, bTrotter;
      real              dvdl_constr;
 -    rvec             *cbuf        = NULL;
 +    rvec             *cbuf        = nullptr;
      int               cbuf_nalloc = 0;
      matrix            lastbox;
      int               lamnew  = 0;
  
  
      /* PME load balancing data for GPU kernels */
 -    pme_load_balancing_t *pme_loadbal      = NULL;
 +    pme_load_balancing_t *pme_loadbal      = nullptr;
      gmx_bool              bPMETune         = FALSE;
      gmx_bool              bPMETunePrinting = FALSE;
  
          nstglobalcomm     = 1;
      }
  
 -    nstglobalcomm   = check_nstglobalcomm(fplog, cr, nstglobalcomm, ir);
 +    nstglobalcomm   = check_nstglobalcomm(mdlog, nstglobalcomm, ir);
      bGStatEveryStep = (nstglobalcomm == 1);
  
      if (bRerunMD)
      {
          /* Initialize ion swapping code */
          init_swapcoords(fplog, bVerbose, ir, opt2fn_master("-swap", nfile, fnm, cr),
 -                        top_global, state_global->x, state_global->box, &state_global->swapstate, cr, oenv, Flags);
 +                        top_global, as_rvec_array(state_global->x.data()), state_global->box, &state_global->swapstate, cr, oenv, Flags);
      }
  
      /* Initial values */
 -    init_md(fplog, cr, ir, oenv, &t, &t0, state_global->lambda,
 +    init_md(fplog, cr, outputProvider, ir, oenv, &t, &t0, state_global->lambda,
              &(state_global->fep_state), lam0,
              nrnb, top_global, &upd,
              nfile, fnm, &outf, &mdebin,
      snew(enerd, 1);
      init_enerdata(top_global->groups.grps[egcENER].nr, ir->fepvals->n_lambda,
                    enerd);
 -    if (DOMAINDECOMP(cr))
 -    {
 -        f = NULL;
 -    }
 -    else
 -    {
 -        snew(f, top_global->natoms);
 -    }
  
      /* Kinetic energy data */
      snew(ekind, 1);
          }
      }
  
 +    std::unique_ptr<t_state> stateInstance;
 +    t_state *                state;
 +
      if (DOMAINDECOMP(cr))
      {
          top = dd_init_local_top(top_global);
  
 -        snew(state, 1);
 +        stateInstance = std::unique_ptr<t_state>(new t_state);
 +        state         = stateInstance.get();
          dd_init_local_state(cr->dd, state_global, state);
      }
      else
      {
 -        top = gmx_mtop_generate_local_top(top_global, ir->efep != efepNO);
 -
 -        state    = serial_init_local_state(state_global);
 -
 -        atoms2md(top_global, ir, 0, NULL, top_global->natoms, mdatoms);
 -
 -        if (vsite)
 -        {
 -            set_vsite_top(vsite, top, mdatoms, cr);
 -        }
 -
 -        if (ir->ePBC != epbcNONE && !fr->bMolPBC)
 -        {
 -            graph = mk_graph(fplog, &(top->idef), 0, top_global->natoms, FALSE, FALSE);
 -        }
 -
 -        if (shellfc)
 -        {
 -            make_local_shells(cr, mdatoms, shellfc);
 -        }
 +        state_change_natoms(state_global, state_global->natoms);
 +        /* We need to allocate one element extra, since we might use
 +         * (unaligned) 4-wide SIMD loads to access rvec entries.
 +         */
 +        f.resize(state_global->natoms + 1);
 +        /* Copy the pointer to the global state */
 +        state = state_global;
  
 -        setup_bonded_threading(fr, &top->idef);
 +        snew(top, 1);
 +        mdAlgorithmsSetupAtomData(cr, ir, top_global, top, fr,
 +                                  &graph, mdatoms, vsite, shellfc);
  
 -        update_realloc(upd, state->nalloc);
 +        update_realloc(upd, state->natoms);
      }
  
      /* Set up interactive MD (IMD) */
 -    init_IMD(ir, cr, top_global, fplog, ir->nstcalcenergy, state_global->x,
 +    init_IMD(ir, cr, top_global, fplog, ir->nstcalcenergy, as_rvec_array(state_global->x.data()),
               nfile, fnm, oenv, imdport, Flags);
  
      if (DOMAINDECOMP(cr))
                              state_global, top_global, ir,
                              state, &f, mdatoms, top, fr,
                              vsite, constr,
 -                            nrnb, NULL, FALSE);
 +                            nrnb, nullptr, FALSE);
          shouldCheckNumberOfBondedInteractions = true;
 -        update_realloc(upd, state->nalloc);
 +        update_realloc(upd, state->natoms);
      }
  
      update_mdatoms(mdatoms, state->lambda[efptMASS]);
  
      if (ir->bExpanded)
      {
 -        init_expanded_ensemble(startingFromCheckpoint, ir, &state->dfhist);
 +        init_expanded_ensemble(startingFromCheckpoint, ir, state->dfhist);
      }
  
      if (MASTER(cr))
              /* Update mdebin with energy history if appending to output files */
              if (Flags & MD_APPENDFILES)
              {
 -                restore_energyhistory_from_state(mdebin, state_global->enerhist);
 +                restore_energyhistory_from_state(mdebin, energyHistory);
              }
              else
              {
                  /* We might have read an energy history from checkpoint,
                   * free the allocated memory and reset the counts.
                   */
 -                done_energyhistory(state_global->enerhist);
 -                init_energyhistory(state_global->enerhist);
 +                *energyHistory = {};
              }
          }
          /* Set the initial energy history in state by updating once */
 -        update_energyhistory(state_global->enerhist, mdebin);
 +        update_energyhistory(energyHistory, mdebin);
      }
  
      /* Initialize constraints */
                  !(Flags & MD_REPRODUCIBLE));
      if (bPMETune)
      {
 -        pme_loadbal_init(&pme_loadbal, cr, fplog, ir, state->box,
 +        pme_loadbal_init(&pme_loadbal, cr, mdlog, ir, state->box,
                           fr->ic, fr->pmedata, use_GPU(fr->nbv),
                           &bPMETunePrinting);
      }
  
      if (!ir->bContinuation && !bRerunMD)
      {
 -        if (mdatoms->cFREEZE && (state->flags & (1<<estV)))
 +        if (state->flags & (1 << estV))
          {
 -            /* Set the velocities of frozen particles to zero */
 +            /* Set the velocities of vsites, shells and frozen atoms to zero */
              for (i = 0; i < mdatoms->homenr; i++)
              {
 -                for (m = 0; m < DIM; m++)
 +                if (mdatoms->ptype[i] == eptVSite ||
 +                    mdatoms->ptype[i] == eptShell)
                  {
 -                    if (ir->opts.nFreeze[mdatoms->cFREEZE[i]][m])
 +                    clear_rvec(state->v[i]);
 +                }
 +                else if (mdatoms->cFREEZE)
 +                {
 +                    for (m = 0; m < DIM; m++)
                      {
 -                        state->v[i][m] = 0;
 +                        if (ir->opts.nFreeze[mdatoms->cFREEZE[i]][m])
 +                        {
 +                            state->v[i][m] = 0;
 +                        }
                      }
                  }
              }
          if (vsite)
          {
              /* Construct the virtual sites for the initial configuration */
 -            construct_vsites(vsite, state->x, ir->delta_t, NULL,
 +            construct_vsites(vsite, as_rvec_array(state->x.data()), ir->delta_t, nullptr,
                               top->idef.iparams, top->idef.il,
                               fr->ePBC, fr->bMolPBC, cr, state->box);
          }
  
      bSumEkinhOld = FALSE;
      compute_globals(fplog, gstat, cr, ir, fr, ekind, state, mdatoms, nrnb, vcm,
 -                    NULL, enerd, force_vir, shake_vir, total_vir, pres, mu_tot,
 +                    nullptr, enerd, force_vir, shake_vir, total_vir, pres, mu_tot,
                      constr, &nullSignaller, state->box,
                      &totalNumberOfBondedInteractions, &bSumEkinhOld, cglo_flags
                      | (shouldCheckNumberOfBondedInteractions ? CGLO_CHECK_NUMBER_OF_BONDED_INTERACTIONS : 0));
             perhaps loses some logic?*/
  
          compute_globals(fplog, gstat, cr, ir, fr, ekind, state, mdatoms, nrnb, vcm,
 -                        NULL, enerd, force_vir, shake_vir, total_vir, pres, mu_tot,
 +                        nullptr, enerd, force_vir, shake_vir, total_vir, pres, mu_tot,
                          constr, &nullSignaller, state->box,
 -                        NULL, &bSumEkinhOld,
 +                        nullptr, &bSumEkinhOld,
                          cglo_flags &~(CGLO_STOPCM | CGLO_PRESSURE));
      }
  
      {
          if (!multisim_int_all_are_equal(cr->ms, ir->nsteps))
          {
 -            md_print_info(cr, fplog,
 -                          "Note: The number of steps is not consistent across multi simulations,\n"
 -                          "but we are proceeding anyway!\n");
 +            GMX_LOG(mdlog.warning).appendText(
 +                    "Note: The number of steps is not consistent across multi simulations,\n"
 +                    "but we are proceeding anyway!");
          }
          if (!multisim_int_all_are_equal(cr->ms, ir->init_step))
          {
 -            md_print_info(cr, fplog,
 -                          "Note: The initial step is not consistent across multi simulations,\n"
 -                          "but we are proceeding anyway!\n");
 +            GMX_LOG(mdlog.warning).appendText(
 +                    "Note: The initial step is not consistent across multi simulations,\n"
 +                    "but we are proceeding anyway!");
          }
      }
  
          {
              /* PME grid + cut-off optimization with GPUs or PME nodes */
              pme_loadbal_do(pme_loadbal, cr,
 -                           (bVerbose && MASTER(cr)) ? stderr : NULL,
 -                           fplog,
 +                           (bVerbose && MASTER(cr)) ? stderr : nullptr,
 +                           fplog, mdlog,
                             ir, fr, state,
                             wcycle,
                             step, step_rel,
                      /* Following is necessary because the graph may get out of sync
                       * with the coordinates if we only have every N'th coordinate set
                       */
 -                    mk_mshift(fplog, graph, fr->ePBC, state->box, state->x);
 -                    shift_self(graph, state->box, state->x);
 +                    mk_mshift(fplog, graph, fr->ePBC, state->box, as_rvec_array(state->x.data()));
 +                    shift_self(graph, state->box, as_rvec_array(state->x.data()));
                  }
 -                construct_vsites(vsite, state->x, ir->delta_t, state->v,
 +                construct_vsites(vsite, as_rvec_array(state->x.data()), ir->delta_t, as_rvec_array(state->v.data()),
                                   top->idef.iparams, top->idef.il,
                                   fr->ePBC, fr->bMolPBC, cr, state->box);
                  if (graph)
                  {
 -                    unshift_self(graph, state->box, state->x);
 +                    unshift_self(graph, state->box, as_rvec_array(state->x.data()));
                  }
              }
          }
                                      nrnb, wcycle,
                                      do_verbose && !bPMETunePrinting);
                  shouldCheckNumberOfBondedInteractions = true;
 -                update_realloc(upd, state->nalloc);
 +                update_realloc(upd, state->natoms);
              }
          }
  
               * the full step kinetic energy and possibly for T-coupling.*/
              /* This may not be quite working correctly yet . . . . */
              compute_globals(fplog, gstat, cr, ir, fr, ekind, state, mdatoms, nrnb, vcm,
 -                            wcycle, enerd, NULL, NULL, NULL, NULL, mu_tot,
 +                            wcycle, enerd, nullptr, nullptr, nullptr, nullptr, mu_tot,
                              constr, &nullSignaller, state->box,
                              &totalNumberOfBondedInteractions, &bSumEkinhOld,
                              CGLO_GSTAT | CGLO_TEMPERATURE | CGLO_CHECK_NUMBER_OF_BONDED_INTERACTIONS);
              relax_shell_flexcon(fplog, cr, bVerbose, step,
                                  ir, bNS, force_flags, top,
                                  constr, enerd, fcd,
 -                                state, f, force_vir, mdatoms,
 +                                state, &f, force_vir, mdatoms,
                                  nrnb, wcycle, graph, groups,
                                  shellfc, fr, bBornRadii, t, mu_tot,
 -                                vsite, mdoutf_get_fp_field(outf));
 +                                vsite);
          }
          else
          {
               * Check comments in sim_util.c
               */
              do_force(fplog, cr, ir, step, nrnb, wcycle, top, groups,
 -                     state->box, state->x, &state->hist,
 -                     f, force_vir, mdatoms, enerd, fcd,
 +                     state->box, &state->x, &state->hist,
 +                     &f, force_vir, mdatoms, enerd, fcd,
                       state->lambda, graph,
 -                     fr, vsite, mu_tot, t, mdoutf_get_fp_field(outf), ed, bBornRadii,
 +                     fr, vsite, mu_tot, t, ed, bBornRadii,
                       (bNS ? GMX_FORCE_NS : 0) | force_flags);
          }
  
          if (EI_VV(ir->eI) && !startingFromCheckpoint && !bRerunMD)
          /*  ############### START FIRST UPDATE HALF-STEP FOR VV METHODS############### */
          {
 -            rvec *vbuf = NULL;
 +            rvec *vbuf = nullptr;
  
              wallcycle_start(wcycle, ewcUPDATE);
              if (ir->eI == eiVV && bInitStep)
                   * so that the input is actually the initial step.
                   */
                  snew(vbuf, state->natoms);
 -                copy_rvecn(state->v, vbuf, 0, state->natoms); /* should make this better for parallelizing? */
 +                copy_rvecn(as_rvec_array(state->v.data()), vbuf, 0, state->natoms); /* should make this better for parallelizing? */
              }
              else
              {
                  trotter_update(ir, step, ekind, enerd, state, total_vir, mdatoms, &MassQ, trotter_seq, ettTSEQ1);
              }
  
 -            update_coords(fplog, step, ir, mdatoms, state, f, fcd,
 +            update_coords(fplog, step, ir, mdatoms, state, &f, fcd,
                            ekind, M, upd, etrtVELOCITY1,
                            cr, constr);
  
              if (!bRerunMD || rerun_fr.bV || bForceUpdate)         /* Why is rerun_fr.bV here?  Unclear. */
              {
                  wallcycle_stop(wcycle, ewcUPDATE);
 -                update_constraints(fplog, step, NULL, ir, mdatoms,
 -                                   state, fr->bMolPBC, graph, f,
 +                update_constraints(fplog, step, nullptr, ir, mdatoms,
 +                                   state, fr->bMolPBC, graph, &f,
                                     &top->idef, shake_vir,
                                     cr, nrnb, wcycle, upd, constr,
                                     TRUE, bCalcVir);
              {
                  /* Need to unshift here if a do_force has been
                     called in the previous step */
 -                unshift_self(graph, state->box, state->x);
 +                unshift_self(graph, state->box, as_rvec_array(state->x.data()));
              }
              /* if VV, compute the pressure and constraints */
              /* For VV2, we strictly only need this if using pressure
                      m_add(force_vir, shake_vir, total_vir);     /* we need the un-dispersion corrected total vir here */
                      trotter_update(ir, step, ekind, enerd, state, total_vir, mdatoms, &MassQ, trotter_seq, ettTSEQ2);
  
-                     copy_mat(shake_vir, state->svir_prev);
-                     copy_mat(force_vir, state->fvir_prev);
+                     /* TODO This is only needed when we're about to write
+                      * a checkpoint, because we use it after the restart
+                      * (in a kludge?). But what should we be doing if
+                      * startingFromCheckpoint or bInitStep are true? */
+                     if (inputrecNptTrotter(ir) || inputrecNphTrotter(ir))
+                     {
+                         copy_mat(shake_vir, state->svir_prev);
+                         copy_mat(force_vir, state->fvir_prev);
+                     }
                      if (inputrecNvtTrotter(ir) && ir->eI == eiVV)
                      {
                          /* update temperature and kinetic energy now that step is over - this is the v(t+dt) point */
 -                        enerd->term[F_TEMP] = sum_ekin(&(ir->opts), ekind, NULL, (ir->eI == eiVV), FALSE);
 +                        enerd->term[F_TEMP] = sum_ekin(&(ir->opts), ekind, nullptr, (ir->eI == eiVV), FALSE);
                          enerd->term[F_EKIN] = trace(ekind->ekin);
                      }
                  }
                       * the full step kinetic energy and possibly for T-coupling.*/
                      /* This may not be quite working correctly yet . . . . */
                      compute_globals(fplog, gstat, cr, ir, fr, ekind, state, mdatoms, nrnb, vcm,
 -                                    wcycle, enerd, NULL, NULL, NULL, NULL, mu_tot,
 +                                    wcycle, enerd, nullptr, nullptr, nullptr, nullptr, mu_tot,
                                      constr, &nullSignaller, state->box,
 -                                    NULL, &bSumEkinhOld,
 +                                    nullptr, &bSumEkinhOld,
                                      CGLO_GSTAT | CGLO_TEMPERATURE);
                      wallcycle_start(wcycle, ewcUPDATE);
                  }
              /* if it's the initial step, we performed this first step just to get the constraint virial */
              if (ir->eI == eiVV && bInitStep)
              {
 -                copy_rvecn(vbuf, state->v, 0, state->natoms);
 +                copy_rvecn(vbuf, as_rvec_array(state->v.data()), 0, state->natoms);
                  sfree(vbuf);
              }
              wallcycle_stop(wcycle, ewcUPDATE);
          /* compute the conserved quantity */
          if (EI_VV(ir->eI))
          {
 -            saved_conserved_quantity = compute_conserved_from_auxiliary(ir, state, &MassQ);
 +            saved_conserved_quantity = NPT_energy(ir, state, &MassQ);
              if (ir->eI == eiVV)
              {
                  last_ekin = enerd->term[F_EKIN];
                 statistics, but if performing simulated tempering, we
                 do update the velocities and the tau_t. */
  
 -            lamnew = ExpandedEnsembleDynamics(fplog, ir, enerd, state, &MassQ, state->fep_state, &state->dfhist, step, state->v, mdatoms);
 +            lamnew = ExpandedEnsembleDynamics(fplog, ir, enerd, state, &MassQ, state->fep_state, state->dfhist, step, as_rvec_array(state->v.data()), mdatoms);
              /* history is maintained in state->dfhist, but state_global is what is sent to trajectory and log output */
 -            copy_df_history(&state_global->dfhist, &state->dfhist);
 +            copy_df_history(state_global->dfhist, state->dfhist);
          }
  
          /* Now we have the energies and forces corresponding to the
           * the update.
           */
          do_md_trajectory_writing(fplog, cr, nfile, fnm, step, step_rel, t,
 -                                 ir, state, state_global, top_global, fr,
 -                                 outf, mdebin, ekind, f,
 +                                 ir, state, state_global, energyHistory,
 +                                 top_global, fr,
 +                                 outf, mdebin, ekind, &f,
                                   &nchkpt,
                                   bCPT, bRerunMD, bLastStep, (Flags & MD_CONFOUT),
                                   bSumEkinhOld);
          /* Check if IMD step and do IMD communication, if bIMD is TRUE. */
 -        bIMDstep = do_IMD(ir->bIMD, step, cr, bNS, state->box, state->x, ir, t, wcycle);
 +        bIMDstep = do_IMD(ir->bIMD, step, cr, bNS, state->box, as_rvec_array(state->x.data()), ir, t, wcycle);
  
          /* kludge -- virial is lost with restart for MTTK NPT control. Must reload (saved earlier). */
-         if (startingFromCheckpoint && bTrotter)
+         if (startingFromCheckpoint && (inputrecNptTrotter(ir) || inputrecNphTrotter(ir)))
          {
              copy_mat(state->svir_prev, shake_vir);
              copy_mat(state->fvir_prev, force_vir);
              /* if we have constraints, we have to remove the kinetic energy parallel to the bonds */
              if (constr && bIfRandomize)
              {
 -                update_constraints(fplog, step, NULL, ir, mdatoms,
 -                                   state, fr->bMolPBC, graph, f,
 +                update_constraints(fplog, step, nullptr, ir, mdatoms,
 +                                   state, fr->bMolPBC, graph, &f,
                                     &top->idef, tmp_vir,
                                     cr, nrnb, wcycle, upd, constr,
                                     TRUE, bCalcVir);
              else
              {
                  update_tcouple(step, ir, state, ekind, &MassQ, mdatoms);
 -                update_pcouple(fplog, step, ir, state, pcoupl_mu, M, bInitStep);
 +                update_pcouple_before_coordinates(fplog, step, ir, state,
 +                                                  parrinellorahmanMu, M,
 +                                                  bInitStep);
              }
  
              if (EI_VV(ir->eI))
              {
                  /* velocity half-step update */
 -                update_coords(fplog, step, ir, mdatoms, state, f, fcd,
 +                update_coords(fplog, step, ir, mdatoms, state, &f, fcd,
                                ekind, M, upd, etrtVELOCITY2,
                                cr, constr);
              }
                      cbuf_nalloc = state->natoms;
                      srenew(cbuf, cbuf_nalloc);
                  }
 -                copy_rvecn(state->x, cbuf, 0, state->natoms);
 +                copy_rvecn(as_rvec_array(state->x.data()), cbuf, 0, state->natoms);
              }
  
 -            update_coords(fplog, step, ir, mdatoms, state, f, fcd,
 +            update_coords(fplog, step, ir, mdatoms, state, &f, fcd,
                            ekind, M, upd, etrtPOSITION, cr, constr);
              wallcycle_stop(wcycle, ewcUPDATE);
  
              update_constraints(fplog, step, &dvdl_constr, ir, mdatoms, state,
 -                               fr->bMolPBC, graph, f,
 +                               fr->bMolPBC, graph, &f,
                                 &top->idef, shake_vir,
                                 cr, nrnb, wcycle, upd, constr,
                                 FALSE, bCalcVir);
                  compute_globals(fplog, gstat, cr, ir, fr, ekind, state, mdatoms, nrnb, vcm,
                                  wcycle, enerd, force_vir, shake_vir, total_vir, pres, mu_tot,
                                  constr, &nullSignaller, lastbox,
 -                                NULL, &bSumEkinhOld,
 +                                nullptr, &bSumEkinhOld,
                                  (bGStat ? CGLO_GSTAT : 0) | CGLO_TEMPERATURE
                                  );
                  wallcycle_start(wcycle, ewcUPDATE);
                  trotter_update(ir, step, ekind, enerd, state, total_vir, mdatoms, &MassQ, trotter_seq, ettTSEQ4);
                  /* now we know the scaling, we can compute the positions again again */
 -                copy_rvecn(cbuf, state->x, 0, state->natoms);
 +                copy_rvecn(cbuf, as_rvec_array(state->x.data()), 0, state->natoms);
  
 -                update_coords(fplog, step, ir, mdatoms, state, f, fcd,
 +                update_coords(fplog, step, ir, mdatoms, state, &f, fcd,
                                ekind, M, upd, etrtPOSITION, cr, constr);
                  wallcycle_stop(wcycle, ewcUPDATE);
  
 -                /* do we need an extra constraint here? just need to copy out of state->v to upd->xp? */
 +                /* do we need an extra constraint here? just need to copy out of as_rvec_array(state->v.data()) to upd->xp? */
                  /* are the small terms in the shake_vir here due
                   * to numerical errors, or are they important
                   * physically? I'm thinking they are just errors, but not completely sure.
                   * For now, will call without actually constraining, constr=NULL*/
 -                update_constraints(fplog, step, NULL, ir, mdatoms,
 -                                   state, fr->bMolPBC, graph, f,
 +                update_constraints(fplog, step, nullptr, ir, mdatoms,
 +                                   state, fr->bMolPBC, graph, &f,
                                     &top->idef, tmp_vir,
 -                                   cr, nrnb, wcycle, upd, NULL,
 +                                   cr, nrnb, wcycle, upd, nullptr,
                                     FALSE, bCalcVir);
              }
              if (EI_VV(ir->eI))
          else if (graph)
          {
              /* Need to unshift here */
 -            unshift_self(graph, state->box, state->x);
 +            unshift_self(graph, state->box, as_rvec_array(state->x.data()));
          }
  
 -        if (vsite != NULL)
 +        if (vsite != nullptr)
          {
              wallcycle_start(wcycle, ewcVSITECONSTR);
 -            if (graph != NULL)
 +            if (graph != nullptr)
              {
 -                shift_self(graph, state->box, state->x);
 +                shift_self(graph, state->box, as_rvec_array(state->x.data()));
              }
 -            construct_vsites(vsite, state->x, ir->delta_t, state->v,
 +            construct_vsites(vsite, as_rvec_array(state->x.data()), ir->delta_t, as_rvec_array(state->v.data()),
                               top->idef.iparams, top->idef.il,
                               fr->ePBC, fr->bMolPBC, cr, state->box);
  
 -            if (graph != NULL)
 +            if (graph != nullptr)
              {
 -                unshift_self(graph, state->box, state->x);
 +                unshift_self(graph, state->box, as_rvec_array(state->x.data()));
              }
              wallcycle_stop(wcycle, ewcVSITECONSTR);
          }
                 Currently done every step so that dhdl is correct in the .edr */
              sum_dhdl(enerd, state->lambda, ir->fepvals);
          }
 -        update_box(fplog, step, ir, mdatoms, state, f,
 -                   pcoupl_mu, nrnb, upd);
 +
 +        update_pcouple_after_coordinates(fplog, step, ir, mdatoms,
 +                                         pres, parrinellorahmanMu,
 +                                         state, nrnb, upd);
  
          /* ################# END UPDATE STEP 2 ################# */
          /* #### We now have r(t+dt) and v(t+dt/2)  ############# */
          }
          else
          {
 -            enerd->term[F_ECONSERVED] = enerd->term[F_ETOT] + compute_conserved_from_auxiliary(ir, state, &MassQ);
 +            enerd->term[F_ECONSERVED] = enerd->term[F_ETOT] + NPT_energy(ir, state, &MassQ);
          }
          /* #########  END PREPARING EDR OUTPUT  ###########  */
  
              if (fplog && do_log && bDoExpanded)
              {
                  /* only needed if doing expanded ensemble */
 -                PrintFreeEnergyInfoToFile(fplog, ir->fepvals, ir->expandedvals, ir->bSimTemp ? ir->simtempvals : NULL,
 -                                          &state_global->dfhist, state->fep_state, ir->nstlog, step);
 +                PrintFreeEnergyInfoToFile(fplog, ir->fepvals, ir->expandedvals, ir->bSimTemp ? ir->simtempvals : nullptr,
 +                                          state_global->dfhist, state->fep_state, ir->nstlog, step);
              }
              if (bCalcEner)
              {
              gmx_bool do_dr  = do_per_step(step, ir->nstdisreout);
              gmx_bool do_or  = do_per_step(step, ir->nstorireout);
  
 -            print_ebin(mdoutf_get_fp_ene(outf), do_ene, do_dr, do_or, do_log ? fplog : NULL,
 +            print_ebin(mdoutf_get_fp_ene(outf), do_ene, do_dr, do_or, do_log ? fplog : nullptr,
                         step, t,
                         eprNORMAL, mdebin, fcd, groups, &(ir->opts));
  
              do_per_step(step, ir->swap->nstswap))
          {
              bNeedRepartition = do_swapcoords(cr, step, t, ir, wcycle,
 -                                             bRerunMD ? rerun_fr.x   : state->x,
 +                                             bRerunMD ? rerun_fr.x   : as_rvec_array(state->x.data()),
                                               bRerunMD ? rerun_fr.box : state->box,
 -                                             top_global, MASTER(cr) && bVerbose, bRerunMD);
 +                                             MASTER(cr) && bVerbose, bRerunMD);
  
              if (bNeedRepartition && DOMAINDECOMP(cr))
              {
                                  vsite, constr,
                                  nrnb, wcycle, FALSE);
              shouldCheckNumberOfBondedInteractions = true;
 -            update_realloc(upd, state->nalloc);
 +            update_realloc(upd, state->natoms);
          }
  
          bFirstStep             = FALSE;
  
          /* #######  END SET VARIABLES FOR NEXT ITERATION ###### */
  
 -        if ( (membed != NULL) && (!bLastStep) )
 +        if ( (membed != nullptr) && (!bLastStep) )
          {
 -            rescale_membed(step_rel, membed, state_global->x);
 +            rescale_membed(step_rel, membed, as_rvec_array(state_global->x.data()));
          }
  
          if (bRerunMD)
                            "resetting counters later in the run, e.g. with gmx "
                            "mdrun -resetstep.", step);
              }
 -            reset_all_counters(fplog, cr, step, &step_rel, ir, wcycle, nrnb, walltime_accounting,
 -                               use_GPU(fr->nbv) ? fr->nbv : NULL);
 +            reset_all_counters(fplog, mdlog, cr, step, &step_rel, ir, wcycle, nrnb, walltime_accounting,
 +                               use_GPU(fr->nbv) ? fr->nbv : nullptr);
              wcycle_set_reset_counters(wcycle, -1);
              if (!(cr->duty & DUTY_PME))
              {
  
      if (bPMETune)
      {
 -        pme_loadbal_done(pme_loadbal, cr, fplog, use_GPU(fr->nbv));
 +        pme_loadbal_done(pme_loadbal, fplog, mdlog, use_GPU(fr->nbv));
      }
  
      done_shellfc(fplog, shellfc, step_rel);
index 0b3dc645eed307cd51dd40afde6f009baa64fe1c,da1fc6a6603530da8b06e35fd1c8540abc21220d..0f1afabcebae4ec4b9d9835c97eca874ef83e15a
  # To help us fund GROMACS development, we humbly ask that you cite
  # the research papers on the package. Check out http://www.gromacs.org.
  
 +if (NOT GMX_BUILD_UNITTESTS)
 +    gmx_add_missing_tests_notice("Unit tests have not been run. You need to set GMX_BUILD_UNITTESTS=ON if you want to build and run them.")
 +    return()
 +endif()
 +
  include_directories(BEFORE SYSTEM ${GMOCK_INCLUDE_DIRS})
  set(TESTUTILS_SOURCES
      cmdlinetest.cpp
+     conftest.cpp
      integrationtests.cpp
      interactivetest.cpp
 +    loggertest.cpp
      mpi-printer.cpp
 +    mpitest.cpp
      refdata.cpp
      refdata-xml.cpp
      stringtest.cpp
@@@ -65,7 -59,7 +66,7 @@@ add_library(testutils STATIC ${UNITTEST
  set(TESTUTILS_LIBS testutils)
  set_property(TARGET testutils APPEND PROPERTY COMPILE_DEFINITIONS "${GMOCK_COMPILE_DEFINITIONS}")
  set_property(TARGET testutils APPEND PROPERTY COMPILE_FLAGS "${GMOCK_COMPILE_FLAGS}")
 -target_link_libraries(testutils libgromacs ${GMOCK_LIBRARIES})
 +target_link_libraries(testutils libgromacs ${GMX_COMMON_LIBRARIES} ${GMOCK_LIBRARIES})
  
  if(HAVE_TINYXML2)
      include_directories(SYSTEM ${TinyXML2_INCLUDE_DIR})
@@@ -74,35 -68,8 +75,35 @@@ else(
      include_directories(BEFORE SYSTEM "../external/tinyxml2")
  endif()
  
 +# TODO Use gmx_add_missing_tests_notice() instead of the messages below.
 +set(GMX_CAN_RUN_MPI_TESTS 1)
 +if (GMX_MPI)
 +    set(_an_mpi_variable_had_content 0)
 +    foreach(VARNAME MPIEXEC MPIEXEC_NUMPROC_FLAG MPIEXEC_PREFLAGS MPIEXEC_POSTFLAGS)
 +        # These variables need a valid value for the test to run
 +        # and pass, but conceivably any of them might be valid
 +        # with arbitrary (including empty) content. They can't be
 +        # valid if they've been populated with the CMake
 +        # find_package magic suffix/value "NOTFOUND", though.
 +        if (${VARNAME} MATCHES ".*NOTFOUND")
 +            message(STATUS "CMake variable ${VARNAME} was not detected to be a valid value. To test GROMACS correctly, check the advice in the install guide.")
 +            set(GMX_CAN_RUN_MPI_TESTS 0)
 +        endif()
 +        if (NOT VARNAME STREQUAL MPIEXEC AND ${VARNAME})
 +            set(_an_mpi_variable_had_content 1)
 +        endif()
 +    endforeach()
 +    if(_an_mpi_variable_had_content AND NOT MPIEXEC)
 +        message(STATUS "CMake variable MPIEXEC must have a valid value if one of the other related MPIEXEC variables does. To test GROMACS correctly, check the advice in the install guide.")
 +        set(GMX_CAN_RUN_MPI_TESTS 0)
 +    endif()
 +elseif (NOT GMX_THREAD_MPI)
 +    set(GMX_CAN_RUN_MPI_TESTS 0)
 +endif()
 +
  set(TESTUTILS_DIR ${CMAKE_CURRENT_SOURCE_DIR})
  set(TESTUTILS_DIR ${TESTUTILS_DIR} PARENT_SCOPE)
  set(TESTUTILS_LIBS ${TESTUTILS_LIBS} PARENT_SCOPE)
 +set(GMX_CAN_RUN_MPI_TESTS ${GMX_CAN_RUN_MPI_TESTS} PARENT_SCOPE)
  
  add_subdirectory(tests)