From: Paul Bauer <paul.bauer.q@gmail.com>
Date: Mon, 28 Oct 2019 15:37:37 +0000 (+0100)
Subject: Merge branch 'release-2019' into master
X-Git-Url: http://biod.pnpi.spb.ru/gitweb/?a=commitdiff_plain;h=506a55bfac650f7a6d9982a2d86252c4d79e9710;p=alexxy%2Fgromacs.git

Merge branch 'release-2019' into master

Updated reference data that needed updating.

Resolved Conflicts:
	cmake/gmxVersionInfo.cmake
	docs/CMakeLists.txt
	src/gromacs/ewald/pme-gpu-types-host.h
	src/gromacs/ewald/pme-only.cpp
	src/gromacs/gmxana/gmx_make_ndx.cpp
	src/gromacs/mdlib/nbnxn_gpu_common.h
	src/gromacs/mdrun/md.cpp
	src/programs/mdrun/tests/minimize.cpp

Change-Id: Ic5aeae6bf932d190b95366373b64ee3449f4a630
---

506a55bfac650f7a6d9982a2d86252c4d79e9710
diff --cc docs/CMakeLists.txt
index 46276b3401,341e14c125..b715982600
--- a/docs/CMakeLists.txt
+++ b/docs/CMakeLists.txt
@@@ -361,15 -367,7 +361,16 @@@ if (SPHINX_FOUND
          how-to/visualize.rst
          install-guide/index.rst
          release-notes/index.rst
 +        release-notes/2020/major/highlights.rst
 +        release-notes/2020/major/features.rst
 +        release-notes/2020/major/performance.rst
 +        release-notes/2020/major/tools.rst
 +        release-notes/2020/major/bugs-fixed.rst
 +        release-notes/2020/major/removed-functionality.rst
 +        release-notes/2020/major/deprecated-functionality.rst
 +        release-notes/2020/major/portability.rst
 +        release-notes/2020/major/miscellaneous.rst
+         release-notes/2019/2019.5.rst
          release-notes/2019/2019.4.rst
          release-notes/2019/2019.3.rst
          release-notes/2019/2019.2.rst
diff --cc src/gromacs/tools/make_ndx.cpp
index a3758b259a,0000000000..a97fb81ece
mode 100644,000000..100644
--- a/src/gromacs/tools/make_ndx.cpp
+++ b/src/gromacs/tools/make_ndx.cpp
@@@ -1,1649 -1,0 +1,1650 @@@
 +/*
 + * This file is part of the GROMACS molecular simulation package.
 + *
 + * Copyright (c) 1991-2000, University of Groningen, The Netherlands.
 + * Copyright (c) 2001-2004, The GROMACS development team.
 + * Copyright (c) 2012,2013,2014,2015,2016,2017,2018,2019, by the GROMACS development team, led by
 + * Mark Abraham, David van der Spoel, Berk Hess, and Erik Lindahl,
 + * and including many others, as listed in the AUTHORS file in the
 + * top-level source directory and at http://www.gromacs.org.
 + *
 + * GROMACS is free software; you can redistribute it and/or
 + * modify it under the terms of the GNU Lesser General Public License
 + * as published by the Free Software Foundation; either version 2.1
 + * of the License, or (at your option) any later version.
 + *
 + * GROMACS is distributed in the hope that it will be useful,
 + * but WITHOUT ANY WARRANTY; without even the implied warranty of
 + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
 + * Lesser General Public License for more details.
 + *
 + * You should have received a copy of the GNU Lesser General Public
 + * License along with GROMACS; if not, see
 + * http://www.gnu.org/licenses, or write to the Free Software Foundation,
 + * Inc., 51 Franklin Street, Fifth Floor, Boston, MA  02110-1301  USA.
 + *
 + * If you want to redistribute modifications to GROMACS, please
 + * consider that scientific software is very special. Version
 + * control is crucial - bugs must be traceable. We will be happy to
 + * consider code for inclusion in the official distribution, but
 + * derived work must not be called official GROMACS. Details are found
 + * in the README & COPYING files - if they are missing, get the
 + * official version at http://www.gromacs.org.
 + *
 + * To help us fund GROMACS development, we humbly ask that you cite
 + * the research papers on the package. Check out http://www.gromacs.org.
 + */
 +#include "gmxpre.h"
 +
 +#include "make_ndx.h"
 +
 +#include <cctype>
 +#include <cstring>
 +
 +#include <algorithm>
 +
 +#include "gromacs/commandline/pargs.h"
 +#include "gromacs/fileio/confio.h"
 +#include "gromacs/math/vec.h"
 +#include "gromacs/topology/block.h"
 +#include "gromacs/topology/index.h"
++#include "gromacs/topology/mtop_util.h"
 +#include "gromacs/topology/topology.h"
 +#include "gromacs/utility/arraysize.h"
 +#include "gromacs/utility/cstringutil.h"
 +#include "gromacs/utility/fatalerror.h"
 +#include "gromacs/utility/futil.h"
 +#include "gromacs/utility/smalloc.h"
 +
 +/* It's not nice to have size limits, but we should not spend more time
 + * on this ancient tool, but instead use the new selection library.
 + */
 +#define MAXNAMES 1024
 +#define NAME_LEN 1024
 +static const int NOTSET = -92637;
 +
 +static gmx_bool  bCase = FALSE;
 +
 +static int or_groups(int nr1, const int *at1, int nr2, const int *at2,
 +                     int *nr, int *at)
 +{
 +    int      i1, i2, max = 0;
 +    gmx_bool bNotIncr;
 +
 +    *nr = 0;
 +
 +    bNotIncr = FALSE;
 +    for (i1 = 0; i1 < nr1; i1++)
 +    {
 +        if ((i1 > 0) && (at1[i1] <= max))
 +        {
 +            bNotIncr = TRUE;
 +        }
 +        max = at1[i1];
 +    }
 +    for (i1 = 0; i1 < nr2; i1++)
 +    {
 +        if ((i1 > 0) && (at2[i1] <= max))
 +        {
 +            bNotIncr = TRUE;
 +        }
 +        max = at2[i1];
 +    }
 +
 +    if (bNotIncr)
 +    {
 +        printf("One of your groups is not ascending\n");
 +    }
 +    else
 +    {
 +        i1  = 0;
 +        i2  = 0;
 +        *nr = 0;
 +        while ((i1 < nr1) || (i2 < nr2))
 +        {
 +            if ((i2 == nr2) || ((i1 < nr1) && (at1[i1] < at2[i2])))
 +            {
 +                at[*nr] = at1[i1];
 +                (*nr)++;
 +                i1++;
 +            }
 +            else
 +            {
 +                if ((i2 < nr2) && ((i1 == nr1) || (at1[i1] > at2[i2])))
 +                {
 +                    at[*nr] = at2[i2];
 +                    (*nr)++;
 +                }
 +                i2++;
 +            }
 +        }
 +
 +        printf("Merged two groups with OR: %d %d -> %d\n", nr1, nr2, *nr);
 +    }
 +
 +    return *nr;
 +}
 +
 +static int and_groups(int nr1, const int *at1, int nr2, const int *at2,
 +                      int *nr, int *at)
 +{
 +    int i1, i2;
 +
 +    *nr = 0;
 +    for (i1 = 0; i1 < nr1; i1++)
 +    {
 +        for (i2 = 0; i2 < nr2; i2++)
 +        {
 +            if (at1[i1] == at2[i2])
 +            {
 +                at[*nr] = at1[i1];
 +                (*nr)++;
 +            }
 +        }
 +    }
 +
 +    printf("Merged two groups with AND: %d %d -> %d\n", nr1, nr2, *nr);
 +
 +    return *nr;
 +}
 +
 +static gmx_bool is_name_char(char c)
 +{
 +    /* This string should contain all characters that can not be
 +     * the first letter of a name due to the make_ndx syntax.
 +     */
 +    const char *spec = " !&|";
 +
 +    return (c != '\0' && std::strchr(spec, c) == nullptr);
 +}
 +
 +static int parse_names(char **string, int *n_names, char **names)
 +{
 +    int i;
 +
 +    *n_names = 0;
 +    while (is_name_char((*string)[0]) || (*string)[0] == ' ')
 +    {
 +        if (is_name_char((*string)[0]))
 +        {
 +            if (*n_names >= MAXNAMES)
 +            {
 +                gmx_fatal(FARGS, "To many names: %d\n", *n_names+1);
 +            }
 +            i = 0;
 +            while (is_name_char((*string)[i]))
 +            {
 +                names[*n_names][i] = (*string)[i];
 +                i++;
 +                if (i > NAME_LEN)
 +                {
 +                    printf("Name is too long, the maximum is %d characters\n", NAME_LEN);
 +                    return 0;
 +                }
 +            }
 +            names[*n_names][i] = '\0';
 +            if (!bCase)
 +            {
 +                upstring(names[*n_names]);
 +            }
 +            *string += i;
 +            (*n_names)++;
 +        }
 +        else
 +        {
 +            (*string)++;
 +        }
 +    }
 +
 +    return *n_names;
 +}
 +
 +static gmx_bool parse_int_char(char **string, int *nr, char *c)
 +{
 +    char    *orig;
 +    gmx_bool bRet;
 +
 +    orig = *string;
 +
 +    while ((*string)[0] == ' ')
 +    {
 +        (*string)++;
 +    }
 +
 +    bRet = FALSE;
 +
 +    *c = ' ';
 +
 +    if (std::isdigit((*string)[0]))
 +    {
 +        *nr = (*string)[0]-'0';
 +        (*string)++;
 +        while (std::isdigit((*string)[0]))
 +        {
 +            *nr = (*nr)*10+(*string)[0]-'0';
 +            (*string)++;
 +        }
 +        if (std::isalpha((*string)[0]))
 +        {
 +            *c = (*string)[0];
 +            (*string)++;
 +        }
 +        /* Check if there is at most one non-digit character */
 +        if (!std::isalnum((*string)[0]))
 +        {
 +            bRet = TRUE;
 +        }
 +        else
 +        {
 +            *string = orig;
 +        }
 +    }
 +    else
 +    {
 +        *nr = NOTSET;
 +    }
 +
 +    return bRet;
 +}
 +
 +static gmx_bool parse_int(char **string, int *nr)
 +{
 +    char    *orig, c;
 +    gmx_bool bRet;
 +
 +    orig = *string;
 +    bRet = parse_int_char(string, nr, &c);
 +    if (bRet && c != ' ')
 +    {
 +        *string = orig;
 +        bRet    = FALSE;
 +    }
 +
 +    return bRet;
 +}
 +
 +static gmx_bool isquote(char c)
 +{
 +    return (c == '\"');
 +}
 +
 +static gmx_bool parse_string(char **string, int *nr, int ngrps, char **grpname)
 +{
 +    char *s, *sp;
 +    char  c;
 +
 +    while ((*string)[0] == ' ')
 +    {
 +        (*string)++;
 +    }
 +
 +    (*nr) = NOTSET;
 +    if (isquote((*string)[0]))
 +    {
 +        c = (*string)[0];
 +        (*string)++;
 +        s  = gmx_strdup((*string));
 +        sp = std::strchr(s, c);
 +        if (sp != nullptr)
 +        {
 +            (*string) += sp-s + 1;
 +            sp[0]      = '\0';
 +            (*nr)      = find_group(s, ngrps, grpname);
 +        }
 +    }
 +
 +    return (*nr) != NOTSET;
 +}
 +
 +static int select_atomnumbers(char **string, const t_atoms *atoms, int n1,
 +                              int *nr, int *index, char *gname)
 +{
 +    char    buf[STRLEN];
 +    int     i, up;
 +
 +    *nr = 0;
 +    while ((*string)[0] == ' ')
 +    {
 +        (*string)++;
 +    }
 +    if ((*string)[0] == '-')
 +    {
 +        (*string)++;
 +        parse_int(string, &up);
 +        if ((n1 < 1) || (n1 > atoms->nr) || (up < 1) || (up > atoms->nr))
 +        {
 +            printf("Invalid atom range\n");
 +        }
 +        else
 +        {
 +            for (i = n1-1; i <= up-1; i++)
 +            {
 +                index[*nr] = i;
 +                (*nr)++;
 +            }
 +            printf("Found %d atom%s in range %d-%d\n", *nr, (*nr == 1) ? "" : "s", n1, up);
 +            if (n1 == up)
 +            {
 +                sprintf(buf, "a_%d", n1);
 +            }
 +            else
 +            {
 +                sprintf(buf, "a_%d-%d", n1, up);
 +            }
 +            std::strcpy(gname, buf);
 +        }
 +    }
 +    else
 +    {
 +        i = n1;
 +        sprintf(gname, "a");
 +        do
 +        {
 +            if ((i-1 >= 0) && (i-1 < atoms->nr))
 +            {
 +                index[*nr] = i-1;
 +                (*nr)++;
 +                sprintf(buf, "_%d", i);
 +                std::strcat(gname, buf);
 +            }
 +            else
 +            {
 +                printf("Invalid atom number %d\n", i);
 +                *nr = 0;
 +            }
 +        }
 +        while ((*nr != 0) && (parse_int(string, &i)));
 +    }
 +
 +    return *nr;
 +}
 +
 +static int select_residuenumbers(char **string, const t_atoms *atoms,
 +                                 int n1, char c,
 +                                 int *nr, int *index, char *gname)
 +{
 +    char       buf[STRLEN];
 +    int        i, j, up;
 +    t_resinfo *ri;
 +
 +    *nr = 0;
 +    while ((*string)[0] == ' ')
 +    {
 +        (*string)++;
 +    }
 +    if ((*string)[0] == '-')
 +    {
 +        /* Residue number range selection */
 +        if (c != ' ')
 +        {
 +            printf("Error: residue insertion codes can not be used with residue range selection\n");
 +            return 0;
 +        }
 +        (*string)++;
 +        parse_int(string, &up);
 +
 +        for (i = 0; i < atoms->nr; i++)
 +        {
 +            ri = &atoms->resinfo[atoms->atom[i].resind];
 +            for (j = n1; (j <= up); j++)
 +            {
 +                if (ri->nr == j && (c == ' ' || ri->ic == c))
 +                {
 +                    index[*nr] = i;
 +                    (*nr)++;
 +                }
 +            }
 +        }
 +        printf("Found %d atom%s with res.nr. in range %d-%d\n",
 +               *nr, (*nr == 1) ? "" : "s", n1, up);
 +        if (n1 == up)
 +        {
 +            sprintf(buf, "r_%d", n1);
 +        }
 +        else
 +        {
 +            sprintf(buf, "r_%d-%d", n1, up);
 +        }
 +        std::strcpy(gname, buf);
 +    }
 +    else
 +    {
 +        /* Individual residue number/insertion code selection */
 +        j = n1;
 +        sprintf(gname, "r");
 +        do
 +        {
 +            for (i = 0; i < atoms->nr; i++)
 +            {
 +                ri = &atoms->resinfo[atoms->atom[i].resind];
 +                if (ri->nr == j && ri->ic == c)
 +                {
 +                    index[*nr] = i;
 +                    (*nr)++;
 +                }
 +            }
 +            sprintf(buf, "_%d", j);
 +            std::strcat(gname, buf);
 +        }
 +        while (parse_int_char(string, &j, &c));
 +    }
 +
 +    return *nr;
 +}
 +
 +static int select_residueindices(char **string, const t_atoms *atoms,
 +                                 int n1, char c,
 +                                 int *nr, int *index, char *gname)
 +{
 +    /*this should be similar to select_residuenumbers except select by index (sequential numbering in file)*/
 +    /*resind+1 for 1-indexing*/
 +    char       buf[STRLEN];
 +    int        i, j, up;
 +    t_resinfo *ri;
 +
 +    *nr = 0;
 +    while ((*string)[0] == ' ')
 +    {
 +        (*string)++;
 +    }
 +    if ((*string)[0] == '-')
 +    {
 +        /* Residue number range selection */
 +        if (c != ' ')
 +        {
 +            printf("Error: residue insertion codes can not be used with residue range selection\n");
 +            return 0;
 +        }
 +        (*string)++;
 +        parse_int(string, &up);
 +
 +        for (i = 0; i < atoms->nr; i++)
 +        {
 +            ri = &atoms->resinfo[atoms->atom[i].resind];
 +            for (j = n1; (j <= up); j++)
 +            {
 +                if (atoms->atom[i].resind+1 == j && (c == ' ' || ri->ic == c))
 +                {
 +                    index[*nr] = i;
 +                    (*nr)++;
 +                }
 +            }
 +        }
 +        printf("Found %d atom%s with resind.+1 in range %d-%d\n",
 +               *nr, (*nr == 1) ? "" : "s", n1, up);
 +        if (n1 == up)
 +        {
 +            sprintf(buf, "r_%d", n1);
 +        }
 +        else
 +        {
 +            sprintf(buf, "r_%d-%d", n1, up);
 +        }
 +        std::strcpy(gname, buf);
 +    }
 +    else
 +    {
 +        /* Individual residue number/insertion code selection */
 +        j = n1;
 +        sprintf(gname, "r");
 +        do
 +        {
 +            for (i = 0; i < atoms->nr; i++)
 +            {
 +                ri = &atoms->resinfo[atoms->atom[i].resind];
 +                if (atoms->atom[i].resind+1 == j && ri->ic == c)
 +                {
 +                    index[*nr] = i;
 +                    (*nr)++;
 +                }
 +            }
 +            sprintf(buf, "_%d", j);
 +            std::strcat(gname, buf);
 +        }
 +        while (parse_int_char(string, &j, &c));
 +    }
 +
 +    return *nr;
 +}
 +
 +
 +static gmx_bool atoms_from_residuenumbers(const t_atoms *atoms, int group, t_blocka *block,
 +                                          int *nr, int *index, char *gname)
 +{
 +    int i, j, j0, j1, resnr, nres;
 +
 +    j0   = block->index[group];
 +    j1   = block->index[group+1];
 +    nres = atoms->nres;
 +    for (j = j0; j < j1; j++)
 +    {
 +        if (block->a[j] >= nres)
 +        {
 +            printf("Index %s contains number>nres (%d>%d)\n",
 +                   gname, block->a[j]+1, nres);
 +            return FALSE;
 +        }
 +    }
 +    for (i = 0; i < atoms->nr; i++)
 +    {
 +        resnr = atoms->resinfo[atoms->atom[i].resind].nr;
 +        for (j = j0; j < j1; j++)
 +        {
 +            if (block->a[j]+1 == resnr)
 +            {
 +                index[*nr] = i;
 +                (*nr)++;
 +                break;
 +            }
 +        }
 +    }
 +    printf("Found %d atom%s in %d residues from group %s\n",
 +           *nr, (*nr == 1) ? "" : "s", j1-j0, gname);
 +    return *nr != 0;
 +}
 +
 +static gmx_bool comp_name(const char *name, const char *search)
 +{
 +    gmx_bool matches = TRUE;
 +
 +    // Loop while name and search are not end-of-string and matches is true
 +    for (; *name && *search && matches; name++, search++)
 +    {
 +        if (*search == '?')
 +        {
 +            // still matching, continue to next character
 +            continue;
 +        }
 +        else if (*search == '*')
 +        {
 +            if (*(search+1))
 +            {
 +                printf("WARNING: Currently '*' is only supported at the end of an expression\n");
 +            }
 +            // if * is the last char in search string, we have a match,
 +            // otherwise we just failed. Return in either case, we're done.
 +            return (*(search+1) == '\0');
 +        }
 +
 +        // Compare one character
 +        if (bCase)
 +        {
 +            matches = (*name == *search);
 +        }
 +        else
 +        {
 +            matches = (std::toupper(*name) == std::toupper(*search));
 +        }
 +    }
 +
 +    matches = matches && (*name == '\0' && (*search == '\0' || *search == '*'));
 +
 +    return matches;
 +}
 +
 +static int select_chainnames(const t_atoms *atoms, int n_names, char **names,
 +                             int *nr, int *index)
 +{
 +    char    name[2];
 +    int     j;
 +    int     i;
 +
 +    name[1] = 0;
 +    *nr     = 0;
 +    for (i = 0; i < atoms->nr; i++)
 +    {
 +        name[0] = atoms->resinfo[atoms->atom[i].resind].chainid;
 +        j       = 0;
 +        while (j < n_names && !comp_name(name, names[j]))
 +        {
 +            j++;
 +        }
 +        if (j < n_names)
 +        {
 +            index[*nr] = i;
 +            (*nr)++;
 +        }
 +    }
 +    printf("Found %d atom%s with chain identifier%s",
 +           *nr, (*nr == 1) ? "" : "s", (n_names == 1) ? "" : "s");
 +    for (j = 0; (j < n_names); j++)
 +    {
 +        printf(" %s", names[j]);
 +    }
 +    printf("\n");
 +
 +    return *nr;
 +}
 +
 +static int select_atomnames(const t_atoms *atoms, int n_names, char **names,
 +                            int *nr, int *index, gmx_bool bType)
 +{
 +    char   *name;
 +    int     j;
 +    int     i;
 +
 +    *nr = 0;
 +    for (i = 0; i < atoms->nr; i++)
 +    {
 +        if (bType)
 +        {
 +            name = *(atoms->atomtype[i]);
 +        }
 +        else
 +        {
 +            name = *(atoms->atomname[i]);
 +        }
 +        j = 0;
 +        while (j < n_names && !comp_name(name, names[j]))
 +        {
 +            j++;
 +        }
 +        if (j < n_names)
 +        {
 +            index[*nr] = i;
 +            (*nr)++;
 +        }
 +    }
 +    printf("Found %d atoms with %s%s",
 +           *nr, bType ? "type" : "name", (n_names == 1) ? "" : "s");
 +    for (j = 0; (j < n_names); j++)
 +    {
 +        printf(" %s", names[j]);
 +    }
 +    printf("\n");
 +
 +    return *nr;
 +}
 +
 +static int select_residuenames(const t_atoms *atoms, int n_names, char **names,
 +                               int *nr, int *index)
 +{
 +    char   *name;
 +    int     j;
 +    int     i;
 +
 +    *nr = 0;
 +    for (i = 0; i < atoms->nr; i++)
 +    {
 +        name = *(atoms->resinfo[atoms->atom[i].resind].name);
 +        j    = 0;
 +        while (j < n_names && !comp_name(name, names[j]))
 +        {
 +            j++;
 +        }
 +        if (j < n_names)
 +        {
 +            index[*nr] = i;
 +            (*nr)++;
 +        }
 +    }
 +    printf("Found %d atoms with residue name%s", *nr, (n_names == 1) ? "" : "s");
 +    for (j = 0; (j < n_names); j++)
 +    {
 +        printf(" %s", names[j]);
 +    }
 +    printf("\n");
 +
 +    return *nr;
 +}
 +
 +static void copy2block(int n, const int *index, t_blocka *block)
 +{
 +    int i, n0;
 +
 +    block->nr++;
 +    n0         = block->nra;
 +    block->nra = n0+n;
 +    srenew(block->index, block->nr+1);
 +    block->index[block->nr] = n0+n;
 +    srenew(block->a, n0+n);
 +    for (i = 0; (i < n); i++)
 +    {
 +        block->a[n0+i] = index[i];
 +    }
 +}
 +
 +static void make_gname(int n, char **names, char *gname)
 +{
 +    int i;
 +
 +    std::strcpy(gname, names[0]);
 +    for (i = 1; i < n; i++)
 +    {
 +        std::strcat(gname, "_");
 +        std::strcat(gname, names[i]);
 +    }
 +}
 +
 +static void copy_group(int g, t_blocka *block, int *nr, int *index)
 +{
 +    int i, i0;
 +
 +    i0  = block->index[g];
 +    *nr = block->index[g+1]-i0;
 +    for (i = 0; i < *nr; i++)
 +    {
 +        index[i] = block->a[i0+i];
 +    }
 +}
 +
 +static void remove_group(int nr, int nr2, t_blocka *block, char ***gn)
 +{
 +    int   i, j, shift;
 +    char *name;
 +
 +    if (nr2 == NOTSET)
 +    {
 +        nr2 = nr;
 +    }
 +
 +    for (j = 0; j <= nr2-nr; j++)
 +    {
 +        if ((nr < 0) || (nr >= block->nr))
 +        {
 +            printf("Group %d does not exist\n", nr+j);
 +        }
 +        else
 +        {
 +            shift = block->index[nr+1]-block->index[nr];
 +            for (i = block->index[nr+1]; i < block->nra; i++)
 +            {
 +                block->a[i-shift] = block->a[i];
 +            }
 +
 +            for (i = nr; i < block->nr; i++)
 +            {
 +                block->index[i] = block->index[i+1]-shift;
 +            }
 +            name = gmx_strdup((*gn)[nr]);
 +            sfree((*gn)[nr]);
 +            for (i = nr; i < block->nr-1; i++)
 +            {
 +                (*gn)[i] = (*gn)[i+1];
 +            }
 +            block->nr--;
 +            block->nra = block->index[block->nr];
 +            printf("Removed group %d '%s'\n", nr+j, name);
 +            sfree(name);
 +        }
 +    }
 +}
 +
 +static void split_group(const t_atoms *atoms, int sel_nr, t_blocka *block, char ***gn,
 +                        gmx_bool bAtom)
 +{
 +    char    buf[STRLEN], *name;
 +    int     i, resind;
 +    int     a, n0, n1;
 +
 +    printf("Splitting group %d '%s' into %s\n", sel_nr, (*gn)[sel_nr],
 +           bAtom ? "atoms" : "residues");
 +
 +    n0 = block->index[sel_nr];
 +    n1 = block->index[sel_nr+1];
 +    srenew(block->a, block->nra+n1-n0);
 +    for (i = n0; i < n1; i++)
 +    {
 +        a      = block->a[i];
 +        resind = atoms->atom[a].resind;
 +        name   = *(atoms->resinfo[resind].name);
 +        if (bAtom || (i == n0) || (atoms->atom[block->a[i-1]].resind != resind))
 +        {
 +            if (i > n0)
 +            {
 +                block->index[block->nr] = block->nra;
 +            }
 +            block->nr++;
 +            srenew(block->index, block->nr+1);
 +            srenew(*gn, block->nr);
 +            if (bAtom)
 +            {
 +                sprintf(buf, "%s_%s_%d", (*gn)[sel_nr], *atoms->atomname[a], a+1);
 +            }
 +            else
 +            {
 +                sprintf(buf, "%s_%s_%d", (*gn)[sel_nr], name, atoms->resinfo[resind].nr);
 +            }
 +            (*gn)[block->nr-1] = gmx_strdup(buf);
 +        }
 +        block->a[block->nra] = a;
 +        block->nra++;
 +    }
 +    block->index[block->nr] = block->nra;
 +}
 +
 +static int split_chain(const t_atoms *atoms, const rvec *x,
 +                       int sel_nr, t_blocka *block, char ***gn)
 +{
 +    char    buf[STRLEN];
 +    int     j, nchain;
 +    int     i, a, natoms, *start = nullptr, *end = nullptr, ca_start, ca_end;
 +    rvec    vec;
 +
 +    natoms   = atoms->nr;
 +    nchain   = 0;
 +    ca_start = 0;
 +
 +    while (ca_start < natoms)
 +    {
 +        while ((ca_start < natoms) && std::strcmp(*atoms->atomname[ca_start], "CA") != 0)
 +        {
 +            ca_start++;
 +        }
 +        if (ca_start < natoms)
 +        {
 +            srenew(start, nchain+1);
 +            srenew(end, nchain+1);
 +            start[nchain] = ca_start;
 +            while ((start[nchain] > 0) &&
 +                   (atoms->atom[start[nchain]-1].resind ==
 +                    atoms->atom[ca_start].resind))
 +            {
 +                start[nchain]--;
 +            }
 +
 +            i = ca_start;
 +            do
 +            {
 +                ca_end = i;
 +                do
 +                {
 +                    i++;
 +                }
 +                while ((i < natoms) && std::strcmp(*atoms->atomname[i], "CA") != 0);
 +                if (i < natoms)
 +                {
 +                    rvec_sub(x[ca_end], x[i], vec);
 +                }
 +                else
 +                {
 +                    break;
 +                }
 +            }
 +            while (norm(vec) < 0.45);
 +
 +            end[nchain] = ca_end;
 +            while ((end[nchain]+1 < natoms) &&
 +                   (atoms->atom[end[nchain]+1].resind == atoms->atom[ca_end].resind))
 +            {
 +                end[nchain]++;
 +            }
 +            ca_start = end[nchain]+1;
 +            nchain++;
 +        }
 +    }
 +    if (nchain == 1)
 +    {
 +        printf("Found 1 chain, will not split\n");
 +    }
 +    else
 +    {
 +        printf("Found %d chains\n", nchain);
 +    }
 +    for (j = 0; j < nchain; j++)
 +    {
 +        printf("%d:%6d atoms (%d to %d)\n",
 +               j+1, end[j]-start[j]+1, start[j]+1, end[j]+1);
 +    }
 +
 +    if (nchain > 1)
 +    {
 +        srenew(block->a, block->nra+block->index[sel_nr+1]-block->index[sel_nr]);
 +        for (j = 0; j < nchain; j++)
 +        {
 +            block->nr++;
 +            srenew(block->index, block->nr+1);
 +            srenew(*gn, block->nr);
 +            sprintf(buf, "%s_chain%d", (*gn)[sel_nr], j+1);
 +            (*gn)[block->nr-1] = gmx_strdup(buf);
 +            for (i = block->index[sel_nr]; i < block->index[sel_nr+1]; i++)
 +            {
 +                a = block->a[i];
 +                if ((a >= start[j]) && (a <= end[j]))
 +                {
 +                    block->a[block->nra] = a;
 +                    block->nra++;
 +                }
 +            }
 +            block->index[block->nr] = block->nra;
 +            if (block->index[block->nr-1] == block->index[block->nr])
 +            {
 +                remove_group(block->nr-1, NOTSET, block, gn);
 +            }
 +        }
 +    }
 +    sfree(start);
 +    sfree(end);
 +
 +    return nchain;
 +}
 +
 +static gmx_bool check_have_atoms(const t_atoms *atoms, char *string)
 +{
 +    if (atoms == nullptr)
 +    {
 +        printf("Can not process '%s' without atom info, use option -f\n", string);
 +        return FALSE;
 +    }
 +    else
 +    {
 +        return TRUE;
 +    }
 +}
 +
 +static gmx_bool parse_entry(char **string, int natoms, const t_atoms *atoms,
 +                            t_blocka *block, char ***gn,
 +                            int *nr, int *index, char *gname)
 +{
 +    static char   **names, *ostring;
 +    static gmx_bool bFirst = TRUE;
 +    int             j, n_names, sel_nr1;
 +    int             i, nr1, *index1;
 +    char            c;
 +    gmx_bool        bRet, bCompl;
 +
 +    if (bFirst)
 +    {
 +        bFirst = FALSE;
 +        snew(names, MAXNAMES);
 +        for (i = 0; i < MAXNAMES; i++)
 +        {
 +            snew(names[i], NAME_LEN+1);
 +        }
 +    }
 +
 +    bRet    = FALSE;
 +    sel_nr1 = NOTSET;
 +
 +    while (*string[0] == ' ')
 +    {
 +        (*string)++;
 +    }
 +
 +    if ((*string)[0] == '!')
 +    {
 +        bCompl = TRUE;
 +        (*string)++;
 +        while (*string[0] == ' ')
 +        {
 +            (*string)++;
 +        }
 +    }
 +    else
 +    {
 +        bCompl = FALSE;
 +    }
 +
 +    ostring = *string;
 +
 +    if (parse_int(string, &sel_nr1) ||
 +        parse_string(string, &sel_nr1, block->nr, *gn))
 +    {
 +        if ((sel_nr1 >= 0) && (sel_nr1 < block->nr))
 +        {
 +            copy_group(sel_nr1, block, nr, index);
 +            std::strcpy(gname, (*gn)[sel_nr1]);
 +            printf("Copied index group %d '%s'\n", sel_nr1, (*gn)[sel_nr1]);
 +            bRet = TRUE;
 +        }
 +        else
 +        {
 +            printf("Group %d does not exist\n", sel_nr1);
 +        }
 +    }
 +    else if ((*string)[0] == 'a')
 +    {
 +        (*string)++;
 +        if (check_have_atoms(atoms, ostring))
 +        {
 +            if (parse_int(string, &sel_nr1))
 +            {
 +                bRet = (select_atomnumbers(string, atoms, sel_nr1, nr, index, gname) != 0);
 +            }
 +            else if (parse_names(string, &n_names, names))
 +            {
 +                bRet = (select_atomnames(atoms, n_names, names, nr, index, FALSE) != 0);
 +                make_gname(n_names, names, gname);
 +            }
 +        }
 +    }
 +    else if ((*string)[0] == 't')
 +    {
 +        (*string)++;
 +        if (check_have_atoms(atoms, ostring) &&
 +            parse_names(string, &n_names, names))
 +        {
 +            if (!(atoms->haveType))
 +            {
 +                printf("Need a run input file to select atom types\n");
 +            }
 +            else
 +            {
 +                bRet = (select_atomnames(atoms, n_names, names, nr, index, TRUE) != 0);
 +                make_gname(n_names, names, gname);
 +            }
 +        }
 +    }
 +    else if (std::strncmp(*string, "res", 3) == 0)
 +    {
 +        (*string) += 3;
 +        if (check_have_atoms(atoms, ostring) &&
 +            parse_int(string, &sel_nr1) &&
 +            (sel_nr1 >= 0) && (sel_nr1 < block->nr) )
 +        {
 +            bRet = atoms_from_residuenumbers(atoms,
 +                                             sel_nr1, block, nr, index, (*gn)[sel_nr1]);
 +            sprintf(gname, "atom_%s", (*gn)[sel_nr1]);
 +        }
 +    }
 +    else if (std::strncmp(*string, "ri", 2) == 0)
 +    {
 +        (*string) += 2;
 +        if (check_have_atoms(atoms, ostring) &&
 +            parse_int_char(string, &sel_nr1, &c))
 +        {
 +            bRet = (select_residueindices(string, atoms, sel_nr1, c, nr, index, gname) != 0);
 +        }
 +    }
 +    else if ((*string)[0] == 'r')
 +    {
 +        (*string)++;
 +        if (check_have_atoms(atoms, ostring))
 +        {
 +            if (parse_int_char(string, &sel_nr1, &c))
 +            {
 +                bRet = (select_residuenumbers(string, atoms, sel_nr1, c, nr, index, gname) != 0);
 +            }
 +            else if (parse_names(string, &n_names, names))
 +            {
 +                bRet = (select_residuenames(atoms, n_names, names, nr, index) != 0);
 +                make_gname(n_names, names, gname);
 +            }
 +        }
 +    }
 +    else if (std::strncmp(*string, "chain", 5) == 0)
 +    {
 +        (*string) += 5;
 +        if (check_have_atoms(atoms, ostring) &&
 +            parse_names(string, &n_names, names))
 +        {
 +            bRet = (select_chainnames(atoms, n_names, names, nr, index) != 0);
 +            sprintf(gname, "ch%s", names[0]);
 +            for (i = 1; i < n_names; i++)
 +            {
 +                std::strcat(gname, names[i]);
 +            }
 +        }
 +    }
 +    if (bRet && bCompl)
 +    {
 +        snew(index1, natoms-*nr);
 +        nr1 = 0;
 +        for (i = 0; i < natoms; i++)
 +        {
 +            j = 0;
 +            while ((j < *nr) && (index[j] != i))
 +            {
 +                j++;
 +            }
 +            if (j == *nr)
 +            {
 +                if (nr1 >= natoms-*nr)
 +                {
 +                    printf("There are double atoms in your index group\n");
 +                    break;
 +                }
 +                index1[nr1] = i;
 +                nr1++;
 +            }
 +        }
 +        *nr = nr1;
 +        for (i = 0; i < nr1; i++)
 +        {
 +            index[i] = index1[i];
 +        }
 +        sfree(index1);
 +
 +        for (i = std::strlen(gname)+1; i > 0; i--)
 +        {
 +            gname[i] = gname[i-1];
 +        }
 +        gname[0] = '!';
 +        printf("Complemented group: %d atoms\n", *nr);
 +    }
 +
 +    return bRet;
 +}
 +
 +static void list_residues(const t_atoms *atoms)
 +{
 +    int      i, j, start, end, prev_resind, resind;
 +    gmx_bool bDiff;
 +
 +    /* Print all the residues, assuming continuous resnr count */
 +    start       = atoms->atom[0].resind;
 +    prev_resind = start;
 +    for (i = 0; i < atoms->nr; i++)
 +    {
 +        resind = atoms->atom[i].resind;
 +        if ((resind != prev_resind) || (i == atoms->nr-1))
 +        {
 +            if ((bDiff = (std::strcmp(*atoms->resinfo[resind].name,
 +                                      *atoms->resinfo[start].name) != 0)) ||
 +                (i == atoms->nr-1))
 +            {
 +                if (bDiff)
 +                {
 +                    end = prev_resind;
 +                }
 +                else
 +                {
 +                    end = resind;
 +                }
 +                if (end < start+3)
 +                {
 +                    for (j = start; j <= end; j++)
 +                    {
 +                        printf("%4d %-5s",
 +                               j+1, *(atoms->resinfo[j].name));
 +                    }
 +                }
 +                else
 +                {
 +                    printf(" %4d - %4d %-5s  ",
 +                           start+1, end+1, *(atoms->resinfo[start].name));
 +                }
 +                start = resind;
 +            }
 +        }
 +        prev_resind = resind;
 +    }
 +    printf("\n");
 +}
 +
 +static void edit_index(int natoms, const t_atoms *atoms, const rvec *x, t_blocka *block, char ***gn, gmx_bool bVerbose)
 +{
 +    static char   **atnames, *ostring;
 +    static gmx_bool bFirst = TRUE;
 +    char            inp_string[STRLEN], *string;
 +    char            gname[STRLEN*3], gname1[STRLEN], gname2[STRLEN];
 +    int             i, i0, i1, sel_nr, sel_nr2, newgroup;
 +    int             nr, nr1, nr2, *index, *index1, *index2;
 +    gmx_bool        bAnd, bOr, bPrintOnce;
 +
 +    if (bFirst)
 +    {
 +        bFirst = FALSE;
 +        snew(atnames, MAXNAMES);
 +        for (i = 0; i < MAXNAMES; i++)
 +        {
 +            snew(atnames[i], NAME_LEN+1);
 +        }
 +    }
 +
 +    string = nullptr;
 +
 +    snew(index, natoms);
 +    snew(index1, natoms);
 +    snew(index2, natoms);
 +
 +    newgroup   = NOTSET;
 +    bPrintOnce = TRUE;
 +    do
 +    {
 +        gname1[0] = '\0';
 +        if (bVerbose || bPrintOnce || newgroup != NOTSET)
 +        {
 +            printf("\n");
 +            if (bVerbose || bPrintOnce || newgroup == NOTSET)
 +            {
 +                i0 = 0;
 +                i1 = block->nr;
 +            }
 +            else
 +            {
 +                i0 = newgroup;
 +                i1 = newgroup+1;
 +            }
 +            for (i = i0; i < i1; i++)
 +            {
 +                printf("%3d %-20s: %5d atoms\n", i, (*gn)[i],
 +                       block->index[i+1]-block->index[i]);
 +            }
 +            newgroup = NOTSET;
 +        }
 +        if (bVerbose || bPrintOnce)
 +        {
 +            printf("\n");
 +            printf(" nr : group      '!': not  'name' nr name   'splitch' nr    Enter: list groups\n");
 +            printf(" 'a': atom       '&': and  'del' nr         'splitres' nr   'l': list residues\n");
 +            printf(" 't': atom type  '|': or   'keep' nr        'splitat' nr    'h': help\n");
 +            printf(" 'r': residue              'res' nr         'chain' char\n");
 +            printf(" \"name\": group             'case': case %s         'q': save and quit\n",
 +                   bCase ? "insensitive" : "sensitive  ");
 +            printf(" 'ri': residue index\n");
 +            bPrintOnce = FALSE;
 +        }
 +        printf("\n");
 +        printf("> ");
 +        if (nullptr == fgets(inp_string, STRLEN, stdin))
 +        {
 +            gmx_fatal(FARGS, "Error reading user input");
 +        }
 +        inp_string[std::strlen(inp_string)-1] = 0;
 +        printf("\n");
 +        string = inp_string;
 +        while (string[0] == ' ')
 +        {
 +            string++;
 +        }
 +
 +        ostring = string;
 +        nr      = 0;
 +        if (string[0] == 'h')
 +        {
 +            printf(" nr                : selects an index group by number or quoted string.\n");
 +            printf("                     The string is first matched against the whole group name,\n");
 +            printf("                     then against the beginning and finally against an\n");
 +            printf("                     arbitrary substring. A multiple match is an error.\n");
 +
 +            printf(" 'a' nr1 [nr2 ...] : selects atoms, atom numbering starts at 1.\n");
 +            printf(" 'a' nr1 - nr2     : selects atoms in the range from nr1 to nr2.\n");
 +            printf(" 'a' name1[*] [name2[*] ...] : selects atoms by name(s), '?' matches any char,\n");
 +            printf("                               wildcard '*' allowed at the end of a name.\n");
 +            printf(" 't' type1[*] [type2[*] ...] : as 'a', but for type, run input file required.\n");
 +            printf(" 'r' nr1[ic1] [nr2[ic2] ...] : selects residues by number and insertion code.\n");
 +            printf(" 'r' nr1 - nr2               : selects residues in the range from nr1 to nr2.\n");
 +            printf(" 'r' name1[*] [name2[*] ...] : as 'a', but for residue names.\n");
 +            printf(" 'ri' nr1 - nr2              : selects residue indices, 1-indexed, (as opposed to numbers) in the range from nr1 to nr2.\n");
 +            printf(" 'chain' ch1 [ch2 ...]       : selects atoms by chain identifier(s),\n");
 +            printf("                               not available with a .gro file as input.\n");
 +            printf(" !                 : takes the complement of a group with respect to all\n");
 +            printf("                     the atoms in the input file.\n");
 +            printf(" & |               : AND and OR, can be placed between any of the options\n");
 +            printf("                     above, the input is processed from left to right.\n");
 +            printf(" 'name' nr name    : rename group nr to name.\n");
 +            printf(" 'del' nr1 [- nr2] : deletes one group or groups in the range from nr1 to nr2.\n");
 +            printf(" 'keep' nr         : deletes all groups except nr.\n");
 +            printf(" 'case'            : make all name compares case (in)sensitive.\n");
 +            printf(" 'splitch' nr      : split group into chains using CA distances.\n");
 +            printf(" 'splitres' nr     : split group into residues.\n");
 +            printf(" 'splitat' nr      : split group into atoms.\n");
 +            printf(" 'res' nr          : interpret numbers in group as residue numbers\n");
 +            printf(" Enter             : list the currently defined groups and commands\n");
 +            printf(" 'l'               : list the residues.\n");
 +            printf(" 'h'               : show this help.\n");
 +            printf(" 'q'               : save and quit.\n");
 +            printf("\n");
 +            printf(" Examples:\n");
 +            printf(" > 2 | 4 & r 3-5\n");
 +            printf(" selects all atoms from group 2 and 4 that have residue numbers 3, 4 or 5\n");
 +            printf(" > a C* & !a C CA\n");
 +            printf(" selects all atoms starting with 'C' but not the atoms 'C' and 'CA'\n");
 +            printf(" > \"protein\" & ! \"backb\"\n");
 +            printf(" selects all atoms that are in group 'protein' and not in group 'backbone'\n");
 +            if (bVerbose)
 +            {
 +                printf("\npress Enter ");
 +                getchar();
 +            }
 +        }
 +        else if (std::strncmp(string, "del", 3) == 0)
 +        {
 +            string += 3;
 +            if (parse_int(&string, &sel_nr))
 +            {
 +                while (string[0] == ' ')
 +                {
 +                    string++;
 +                }
 +                if (string[0] == '-')
 +                {
 +                    string++;
 +                    parse_int(&string, &sel_nr2);
 +                }
 +                else
 +                {
 +                    sel_nr2 = NOTSET;
 +                }
 +                while (string[0] == ' ')
 +                {
 +                    string++;
 +                }
 +                if (string[0] == '\0')
 +                {
 +                    remove_group(sel_nr, sel_nr2, block, gn);
 +                }
 +                else
 +                {
 +                    printf("\nSyntax error: \"%s\"\n", string);
 +                }
 +            }
 +        }
 +        else if (std::strncmp(string, "keep", 4) == 0)
 +        {
 +            string += 4;
 +            if (parse_int(&string, &sel_nr))
 +            {
 +                remove_group(sel_nr+1, block->nr-1, block, gn);
 +                remove_group(0, sel_nr-1, block, gn);
 +            }
 +        }
 +        else if (std::strncmp(string, "name", 4) == 0)
 +        {
 +            string += 4;
 +            if (parse_int(&string, &sel_nr))
 +            {
 +                if ((sel_nr >= 0) && (sel_nr < block->nr))
 +                {
 +                    sscanf(string, "%s", gname);
 +                    sfree((*gn)[sel_nr]);
 +                    (*gn)[sel_nr] = gmx_strdup(gname);
 +                }
 +            }
 +        }
 +        else if (std::strncmp(string, "case", 4) == 0)
 +        {
 +            bCase = !bCase;
 +            printf("Switched to case %s\n", bCase ? "sensitive" : "insensitive");
 +        }
 +        else if (string[0] == 'v')
 +        {
 +            bVerbose = !bVerbose;
 +            printf("Turned verbose %s\n", bVerbose ? "on" : "off");
 +        }
 +        else if (string[0] == 'l')
 +        {
 +            if (check_have_atoms(atoms, ostring) )
 +            {
 +                list_residues(atoms);
 +            }
 +        }
 +        else if (std::strncmp(string, "splitch", 7) == 0)
 +        {
 +            string += 7;
 +            if (check_have_atoms(atoms, ostring) &&
 +                parse_int(&string, &sel_nr) &&
 +                (sel_nr >= 0) && (sel_nr < block->nr))
 +            {
 +                split_chain(atoms, x, sel_nr, block, gn);
 +            }
 +        }
 +        else if (std::strncmp(string, "splitres", 8) == 0)
 +        {
 +            string += 8;
 +            if (check_have_atoms(atoms, ostring) &&
 +                parse_int(&string, &sel_nr) &&
 +                (sel_nr >= 0) && (sel_nr < block->nr))
 +            {
 +                split_group(atoms, sel_nr, block, gn, FALSE);
 +            }
 +        }
 +        else if (std::strncmp(string, "splitat", 7) == 0)
 +        {
 +            string += 7;
 +            if (check_have_atoms(atoms, ostring) &&
 +                parse_int(&string, &sel_nr) &&
 +                (sel_nr >= 0) && (sel_nr < block->nr))
 +            {
 +                split_group(atoms, sel_nr, block, gn, TRUE);
 +            }
 +        }
 +        else if (string[0] == '\0')
 +        {
 +            bPrintOnce = TRUE;
 +        }
 +        else if (string[0] != 'q')
 +        {
 +            nr2 = -1;
 +            if (parse_entry(&string, natoms, atoms, block, gn, &nr, index, gname))
 +            {
 +                do
 +                {
 +                    while (string[0] == ' ')
 +                    {
 +                        string++;
 +                    }
 +
 +                    bAnd = FALSE;
 +                    bOr  = FALSE;
 +                    if (string[0] == '&')
 +                    {
 +                        bAnd = TRUE;
 +                    }
 +                    else if (string[0] == '|')
 +                    {
 +                        bOr = TRUE;
 +                    }
 +
 +                    if (bAnd || bOr)
 +                    {
 +                        string++;
 +                        nr1 = nr;
 +                        for (i = 0; i < nr; i++)
 +                        {
 +                            index1[i] = index[i];
 +                        }
 +                        std::strcpy(gname1, gname);
 +                        if (parse_entry(&string, natoms, atoms, block, gn, &nr2, index2, gname2))
 +                        {
 +                            if (bOr)
 +                            {
 +                                or_groups(nr1, index1, nr2, index2, &nr, index);
 +                                sprintf(gname, "%s_%s", gname1, gname2);
 +                            }
 +                            else
 +                            {
 +                                and_groups(nr1, index1, nr2, index2, &nr, index);
 +                                sprintf(gname, "%s_&_%s", gname1, gname2);
 +                            }
 +                        }
 +                    }
 +                }
 +                while (bAnd || bOr);
 +            }
 +            while (string[0] == ' ')
 +            {
 +                string++;
 +            }
 +            if (string[0])
 +            {
 +                printf("\nSyntax error: \"%s\"\n", string);
 +            }
 +            else if (nr > 0)
 +            {
 +                copy2block(nr, index, block);
 +                srenew(*gn, block->nr);
 +                newgroup        = block->nr-1;
 +                (*gn)[newgroup] = gmx_strdup(gname);
 +            }
 +            else
 +            {
 +                printf("Group is empty\n");
 +            }
 +        }
 +    }
 +    while (string[0] != 'q');
 +
 +    sfree(index);
 +    sfree(index1);
 +    sfree(index2);
 +}
 +
 +static int block2natoms(t_blocka *block)
 +{
 +    int i, natoms;
 +
 +    natoms = 0;
 +    for (i = 0; i < block->nra; i++)
 +    {
 +        natoms = std::max(natoms, block->a[i]);
 +    }
 +    natoms++;
 +
 +    return natoms;
 +}
 +
 +static void merge_blocks(t_blocka *dest, t_blocka *source)
 +{
 +    int     i, nra0, i0;
 +
 +    /* count groups, srenew and fill */
 +    i0        = dest->nr;
 +    nra0      = dest->nra;
 +    dest->nr += source->nr;
 +    srenew(dest->index, dest->nr+1);
 +    for (i = 0; i < source->nr; i++)
 +    {
 +        dest->index[i0+i] = nra0 + source->index[i];
 +    }
 +    /* count atoms, srenew and fill */
 +    dest->nra += source->nra;
 +    srenew(dest->a, dest->nra);
 +    for (i = 0; i < source->nra; i++)
 +    {
 +        dest->a[nra0+i] = source->a[i];
 +    }
 +
 +    /* terminate list */
 +    dest->index[dest->nr] = dest->nra;
 +
 +}
 +
 +int gmx_make_ndx(int argc, char *argv[])
 +{
 +    const char     *desc[] = {
 +        "Index groups are necessary for almost every GROMACS program.",
 +        "All these programs can generate default index groups. You ONLY",
 +        "have to use [THISMODULE] when you need SPECIAL index groups.",
 +        "There is a default index group for the whole system, 9 default",
 +        "index groups for proteins, and a default index group",
 +        "is generated for every other residue name.",
 +        "",
 +        "When no index file is supplied, also [THISMODULE] will generate the",
 +        "default groups.",
 +        "With the index editor you can select on atom, residue and chain names",
 +        "and numbers.",
 +        "When a run input file is supplied you can also select on atom type.",
 +        "You can use boolean operations, you can split groups",
 +        "into chains, residues or atoms. You can delete and rename groups.",
 +        "Type 'h' in the editor for more details.",
 +        "",
 +        "The atom numbering in the editor and the index file starts at 1.",
 +        "",
 +        "The [TT]-twin[tt] switch duplicates all index groups with an offset of",
 +        "[TT]-natoms[tt], which is useful for Computational Electrophysiology",
 +        "double-layer membrane setups.",
 +        "",
 +        "See also [gmx-select] [TT]-on[tt], which provides an alternative way",
 +        "for constructing index groups.  It covers nearly all of [THISMODULE]",
 +        "functionality, and in many cases much more."
 +    };
 +
 +    static int      natoms     = 0;
 +    static gmx_bool bVerbose   = FALSE;
 +    static gmx_bool bDuplicate = FALSE;
 +    t_pargs         pa[]       = {
 +        { "-natoms",  FALSE, etINT, {&natoms},
 +          "set number of atoms (default: read from coordinate or index file)" },
 +        { "-twin",     FALSE, etBOOL, {&bDuplicate},
 +          "Duplicate all index groups with an offset of -natoms" },
 +        { "-verbose", FALSE, etBOOL, {&bVerbose},
 +          "HIDDENVerbose output" }
 +    };
 +#define NPA asize(pa)
 +
 +    gmx_output_env_t *oenv;
 +    const char       *stxfile;
 +    const char       *ndxoutfile;
 +    gmx_bool          bNatoms;
 +    int               j;
-     t_atoms          *atoms;
++    t_atoms           atoms;
 +    rvec             *x, *v;
 +    int               ePBC;
 +    matrix            box;
 +    t_blocka         *block, *block2;
 +    char            **gnames, **gnames2;
 +    t_filenm          fnm[] = {
 +        { efSTX, "-f", nullptr,     ffOPTRD  },
 +        { efNDX, "-n", nullptr,     ffOPTRDMULT },
 +        { efNDX, "-o", nullptr,     ffWRITE }
 +    };
 +#define NFILE asize(fnm)
 +
 +    if (!parse_common_args(&argc, argv, 0, NFILE, fnm, NPA, pa, asize(desc), desc,
 +                           0, nullptr, &oenv))
 +    {
 +        return 0;
 +    }
 +
 +    stxfile = ftp2fn_null(efSTX, NFILE, fnm);
 +    gmx::ArrayRef<const std::string> ndxInFiles = opt2fnsIfOptionSet("-n", NFILE, fnm);
 +    ndxoutfile = opt2fn("-o", NFILE, fnm);
 +    bNatoms    = opt2parg_bSet("-natoms", NPA, pa);
 +
 +    if (!stxfile && ndxInFiles.empty())
 +    {
 +        gmx_fatal(FARGS, "No input files (structure or index)");
 +    }
 +
++    gmx_mtop_t mtop;
 +    if (stxfile)
 +    {
-         t_topology *top;
-         snew(top, 1);
++        bool haveFullTopology = false;
 +        fprintf(stderr, "\nReading structure file\n");
-         read_tps_conf(stxfile, top, &ePBC, &x, &v, box, FALSE);
-         atoms = &top->atoms;
-         if (atoms->pdbinfo == nullptr)
++        readConfAndTopology(stxfile, &haveFullTopology, &mtop,
++                            &ePBC, &x, &v, box);
++        atoms = gmx_mtop_global_atoms(&mtop);
++        if (atoms.pdbinfo == nullptr)
 +        {
-             snew(atoms->pdbinfo, atoms->nr);
++            snew(atoms.pdbinfo, atoms.nr);
 +        }
-         natoms  = atoms->nr;
++        natoms  = atoms.nr;
 +        bNatoms = TRUE;
 +    }
 +    else
 +    {
-         atoms = nullptr;
 +        x     = nullptr;
 +    }
 +
 +    /* read input file(s) */
 +    block  = new_blocka();
 +    gnames = nullptr;
 +    printf("Going to read %td old index file(s)\n", ndxInFiles.ssize());
 +    if (!ndxInFiles.empty())
 +    {
 +        for (const std::string &ndxInFile : ndxInFiles)
 +        {
 +            block2 = init_index(ndxInFile.c_str(), &gnames2);
 +            srenew(gnames, block->nr+block2->nr);
 +            for (j = 0; j < block2->nr; j++)
 +            {
 +                gnames[block->nr+j] = gnames2[j];
 +            }
 +            sfree(gnames2);
 +            merge_blocks(block, block2);
 +            sfree(block2->a);
 +            sfree(block2->index);
 +/*       done_block(block2); */
 +            sfree(block2);
 +        }
 +    }
 +    else
 +    {
 +        snew(gnames, 1);
-         analyse(atoms, block, &gnames, FALSE, TRUE);
++        analyse(&atoms, block, &gnames, FALSE, TRUE);
 +    }
 +
 +    if (!bNatoms)
 +    {
 +        natoms = block2natoms(block);
 +        printf("Counted atom numbers up to %d in index file\n", natoms);
 +    }
 +
-     edit_index(natoms, atoms, x, block, &gnames, bVerbose);
++    edit_index(natoms, &atoms, x, block, &gnames, bVerbose);
 +
 +    write_index(ndxoutfile, block, gnames, bDuplicate, natoms);
 +
 +    return 0;
 +}
diff --cc src/programs/mdrun/tests/refdata/Angles1_SimpleMdrunTest_WithinTolerances_0.xml
index e6b945db3c,0000000000..3e723a8d3e
mode 100644,000000..100644
--- a/src/programs/mdrun/tests/refdata/Angles1_SimpleMdrunTest_WithinTolerances_0.xml
+++ b/src/programs/mdrun/tests/refdata/Angles1_SimpleMdrunTest_WithinTolerances_0.xml
@@@ -1,316 -1,0 +1,316 @@@
 +<?xml version="1.0"?>
 +<?xml-stylesheet type="text/xsl" href="referencedata.xsl"?>
 +<ReferenceData>
 +  <Simulation Name="angles1">
 +    <Mdrun Name="md">
++      <Energy Name="Pressure">
++        <Real Name="Time 0.000000 Step 0 in frame 0">0.050796915464453078</Real>
++        <Real Name="Time 0.004000 Step 4 in frame 1">0.049668265065138127</Real>
++        <Real Name="Time 0.008000 Step 8 in frame 2">0.051740106421168469</Real>
++        <Real Name="Time 0.012000 Step 12 in frame 3">0.05656710035689768</Real>
++        <Real Name="Time 0.016000 Step 16 in frame 4">0.062737560105270665</Real>
++        <Real Name="Time 0.020000 Step 20 in frame 5">0.068449719298513415</Real>
++        <Real Name="Time 0.024000 Step 24 in frame 6">0.072170660789040553</Real>
++        <Real Name="Time 0.028000 Step 28 in frame 7">0.073111436629933149</Real>
++        <Real Name="Time 0.032000 Step 32 in frame 8">0.0713787995074181</Real>
++        <Real Name="Time 0.036000 Step 36 in frame 9">0.067818279739478599</Real>
++        <Real Name="Time 0.040000 Step 40 in frame 10">0.06366651387852211</Real>
++        <Real Name="Time 0.044000 Step 44 in frame 11">0.060156895239239429</Real>
++        <Real Name="Time 0.048000 Step 48 in frame 12">0.05819249706054485</Real>
++        <Real Name="Time 0.050000 Step 50 in frame 13">0.057930280864272023</Real>
++      </Energy>
 +      <Energy Name="Kinetic En.">
 +        <Real Name="Time 0.000000 Step 0 in frame 0">8.8827049229659867</Real>
-         <Real Name="Time 0.004000 Step 4 in frame 1">8.6853411978925212</Real>
-         <Real Name="Time 0.008000 Step 8 in frame 2">9.0476379091109891</Real>
-         <Real Name="Time 0.012000 Step 12 in frame 3">9.8917199248001868</Real>
-         <Real Name="Time 0.016000 Step 16 in frame 4">10.970729795432797</Real>
-         <Real Name="Time 0.020000 Step 20 in frame 5">11.969598016517697</Real>
-         <Real Name="Time 0.024000 Step 24 in frame 6">12.620267943889669</Real>
++        <Real Name="Time 0.004000 Step 4 in frame 1">8.6853411978925283</Real>
++        <Real Name="Time 0.008000 Step 8 in frame 2">9.0476379091110388</Real>
++        <Real Name="Time 0.012000 Step 12 in frame 3">9.8917199248002774</Real>
++        <Real Name="Time 0.016000 Step 16 in frame 4">10.970729795432909</Real>
++        <Real Name="Time 0.020000 Step 20 in frame 5">11.969598016517836</Real>
++        <Real Name="Time 0.024000 Step 24 in frame 6">12.620267943889758</Real>
 +        <Real Name="Time 0.028000 Step 28 in frame 7">12.784778605942707</Real>
-         <Real Name="Time 0.032000 Step 32 in frame 8">12.481797526143797</Real>
-         <Real Name="Time 0.036000 Step 36 in frame 9">11.859180066366713</Real>
-         <Real Name="Time 0.040000 Step 40 in frame 10">11.133173167819335</Real>
-         <Real Name="Time 0.044000 Step 44 in frame 11">10.519456636414786</Real>
-         <Real Name="Time 0.048000 Step 48 in frame 12">10.175948192781643</Real>
-         <Real Name="Time 0.050000 Step 50 in frame 13">10.141536444106846</Real>
-       </Energy>
-       <Energy Name="Pressure">
-         <Real Name="Time 0.000000 Step 0 in frame 0">0.050796915464453078</Real>
-         <Real Name="Time 0.004000 Step 4 in frame 1">0.049668265065138079</Real>
-         <Real Name="Time 0.008000 Step 8 in frame 2">0.051740106421168025</Real>
-         <Real Name="Time 0.012000 Step 12 in frame 3">0.056567100356897236</Real>
-         <Real Name="Time 0.016000 Step 16 in frame 4">0.062737560105270318</Real>
-         <Real Name="Time 0.020000 Step 20 in frame 5">0.06844971929851254</Real>
-         <Real Name="Time 0.024000 Step 24 in frame 6">0.072170660789039956</Real>
-         <Real Name="Time 0.028000 Step 28 in frame 7">0.073111436629933163</Real>
-         <Real Name="Time 0.032000 Step 32 in frame 8">0.071378799507418461</Real>
-         <Real Name="Time 0.036000 Step 36 in frame 9">0.067818279739479667</Real>
-         <Real Name="Time 0.040000 Step 40 in frame 10">0.063666513878523082</Real>
-         <Real Name="Time 0.044000 Step 44 in frame 11">0.060156895239239817</Real>
-         <Real Name="Time 0.048000 Step 48 in frame 12">0.058192497060544725</Real>
-         <Real Name="Time 0.050000 Step 50 in frame 13">0.057995708953366192</Real>
++        <Real Name="Time 0.032000 Step 32 in frame 8">12.481797526143708</Real>
++        <Real Name="Time 0.036000 Step 36 in frame 9">11.85918006636658</Real>
++        <Real Name="Time 0.040000 Step 40 in frame 10">11.133173167819194</Real>
++        <Real Name="Time 0.044000 Step 44 in frame 11">10.51945663641469</Real>
++        <Real Name="Time 0.048000 Step 48 in frame 12">10.17594819278162</Real>
++        <Real Name="Time 0.050000 Step 50 in frame 13">10.130095229541357</Real>
 +      </Energy>
 +      <Energy Name="Potential">
 +        <Real Name="Time 0.000000 Step 0 in frame 0">4.1472047198540372</Real>
 +        <Real Name="Time 0.004000 Step 4 in frame 1">4.3456029525254678</Real>
-         <Real Name="Time 0.008000 Step 8 in frame 2">3.9826669068573985</Real>
-         <Real Name="Time 0.012000 Step 12 in frame 3">3.1365732905082528</Real>
-         <Real Name="Time 0.016000 Step 16 in frame 4">2.0549883062523469</Real>
-         <Real Name="Time 0.020000 Step 20 in frame 5">1.053909305545295</Real>
-         <Real Name="Time 0.024000 Step 24 in frame 6">0.40206460374619313</Real>
-         <Real Name="Time 0.028000 Step 28 in frame 7">0.23763430982967265</Real>
-         <Real Name="Time 0.032000 Step 32 in frame 8">0.54173749646261815</Real>
-         <Real Name="Time 0.036000 Step 36 in frame 9">1.1660367172497461</Real>
-         <Real Name="Time 0.040000 Step 40 in frame 10">1.8937346262503456</Real>
-         <Real Name="Time 0.044000 Step 44 in frame 11">2.5086918336925299</Real>
-         <Real Name="Time 0.048000 Step 48 in frame 12">2.852711964582892</Real>
-         <Real Name="Time 0.050000 Step 50 in frame 13">2.8985177691300135</Real>
++        <Real Name="Time 0.008000 Step 8 in frame 2">3.9826669068572804</Real>
++        <Real Name="Time 0.012000 Step 12 in frame 3">3.1365732905081485</Real>
++        <Real Name="Time 0.016000 Step 16 in frame 4">2.0549883062522021</Real>
++        <Real Name="Time 0.020000 Step 20 in frame 5">1.0539093055450899</Real>
++        <Real Name="Time 0.024000 Step 24 in frame 6">0.40206460374611186</Real>
++        <Real Name="Time 0.028000 Step 28 in frame 7">0.23763430982970496</Real>
++        <Real Name="Time 0.032000 Step 32 in frame 8">0.54173749646274727</Real>
++        <Real Name="Time 0.036000 Step 36 in frame 9">1.1660367172499504</Real>
++        <Real Name="Time 0.040000 Step 40 in frame 10">1.8937346262505983</Real>
++        <Real Name="Time 0.044000 Step 44 in frame 11">2.5086918336925828</Real>
++        <Real Name="Time 0.048000 Step 48 in frame 12">2.8527119645828178</Real>
++        <Real Name="Time 0.050000 Step 50 in frame 13">2.8985177691298398</Real>
 +      </Energy>
 +      <Vector Name="Time 0.000000 Step 0 F[0]">
 +        <Real Name="X">10.938095074317319</Real>
 +        <Real Name="Y">-9.6849636829809604</Real>
 +        <Real Name="Z">14.9477673456121</Real>
 +      </Vector>
 +      <Vector Name="Time 0.000000 Step 0 F[1]">
 +        <Real Name="X">-358.95983701890032</Real>
 +        <Real Name="Y">-72.184737465987098</Real>
 +        <Real Name="Z">5.2427273882850791</Real>
 +      </Vector>
 +      <Vector Name="Time 0.000000 Step 0 F[2]">
 +        <Real Name="X">504.44011431612421</Real>
 +        <Real Name="Y">-30.492879467234754</Real>
 +        <Real Name="Z">132.02668359243958</Real>
 +      </Vector>
 +      <Vector Name="Time 0.000000 Step 0 F[3]">
 +        <Real Name="X">-156.41837237154124</Real>
 +        <Real Name="Y">112.36258061620281</Real>
 +        <Real Name="Z">-152.21717832633675</Real>
 +      </Vector>
 +      <Vector Name="Time 0.004000 Step 4 F[0]">
 +        <Real Name="X">20.435838614386633</Real>
 +        <Real Name="Y">-17.686356823439912</Real>
 +        <Real Name="Z">28.095606243064644</Real>
 +      </Vector>
 +      <Vector Name="Time 0.004000 Step 4 F[1]">
 +        <Real Name="X">-393.58451331352791</Real>
 +        <Real Name="Y">-74.845578585694341</Real>
 +        <Real Name="Z">-13.392242553088089</Real>
 +      </Vector>
 +      <Vector Name="Time 0.004000 Step 4 F[2]">
 +        <Real Name="X">529.98619683151344</Real>
 +        <Real Name="Y">-19.588183147046106</Real>
 +        <Real Name="Z">140.73458818251191</Real>
 +      </Vector>
 +      <Vector Name="Time 0.004000 Step 4 F[3]">
 +        <Real Name="X">-156.83752213237216</Real>
 +        <Real Name="Y">112.12011855618036</Real>
 +        <Real Name="Z">-155.43795187248847</Real>
 +      </Vector>
 +      <Vector Name="Time 0.008000 Step 8 F[0]">
-         <Real Name="X">23.175091525779919</Real>
-         <Real Name="Y">-19.684699303050994</Real>
-         <Real Name="Z">32.068217603600317</Real>
++        <Real Name="X">23.175091525778825</Real>
++        <Real Name="Y">-19.684699303050085</Real>
++        <Real Name="Z">32.068217603598796</Real>
 +      </Vector>
 +      <Vector Name="Time 0.008000 Step 8 F[1]">
-         <Real Name="X">-384.67295235829749</Real>
-         <Real Name="Y">-74.107890253924694</Real>
-         <Real Name="Z">-21.330768782037424</Real>
++        <Real Name="X">-384.67295235828954</Real>
++        <Real Name="Y">-74.107890253923756</Real>
++        <Real Name="Z">-21.33076878203557</Real>
 +      </Vector>
 +      <Vector Name="Time 0.008000 Step 8 F[2]">
-         <Real Name="X">508.23023246865938</Real>
-         <Real Name="Y">-11.323306055502048</Real>
-         <Real Name="Z">138.18426885013943</Real>
++        <Real Name="X">508.23023246865029</Real>
++        <Real Name="Y">-11.323306055502478</Real>
++        <Real Name="Z">138.18426885013693</Real>
 +      </Vector>
 +      <Vector Name="Time 0.008000 Step 8 F[3]">
-         <Real Name="X">-146.73237163614181</Real>
-         <Real Name="Y">105.11589561247773</Real>
-         <Real Name="Z">-148.92171767170231</Real>
++        <Real Name="X">-146.73237163613959</Real>
++        <Real Name="Y">105.11589561247632</Real>
++        <Real Name="Z">-148.92171767170015</Real>
 +      </Vector>
 +      <Vector Name="Time 0.012000 Step 12 F[0]">
-         <Real Name="X">19.606862325228409</Real>
-         <Real Name="Y">-16.414806546801842</Real>
-         <Real Name="Z">27.313596732132655</Real>
++        <Real Name="X">19.606862325227326</Real>
++        <Real Name="Y">-16.414806546800946</Real>
++        <Real Name="Z">27.313596732131138</Real>
 +      </Vector>
 +      <Vector Name="Time 0.012000 Step 12 F[1]">
-         <Real Name="X">-335.63542363030876</Real>
-         <Real Name="Y">-68.955808848547065</Real>
-         <Real Name="Z">-18.44424868720855</Real>
++        <Real Name="X">-335.63542363030092</Real>
++        <Real Name="Y">-68.955808848546042</Real>
++        <Real Name="Z">-18.444248687206709</Real>
 +      </Vector>
 +      <Vector Name="Time 0.012000 Step 12 F[2]">
-         <Real Name="X">443.40685745535319</Real>
-         <Real Name="Y">-6.6985698054608456</Real>
-         <Real Name="Z">124.1640033656024</Real>
++        <Real Name="X">443.40685745534432</Real>
++        <Real Name="Y">-6.6985698054614105</Real>
++        <Real Name="Z">124.16400336559995</Real>
 +      </Vector>
 +      <Vector Name="Time 0.012000 Step 12 F[3]">
-         <Real Name="X">-127.37829615027287</Real>
-         <Real Name="Y">92.069185200809756</Real>
-         <Real Name="Z">-133.03335141052651</Real>
++        <Real Name="X">-127.37829615027078</Real>
++        <Real Name="Y">92.069185200808391</Real>
++        <Real Name="Z">-133.03335141052438</Real>
 +      </Vector>
 +      <Vector Name="Time 0.016000 Step 16 F[0]">
-         <Real Name="X">10.936701481888193</Real>
-         <Real Name="Y">-9.0594222620341984</Real>
-         <Real Name="Z">15.337114133364635</Real>
++        <Real Name="X">10.936701481886285</Real>
++        <Real Name="Y">-9.0594222620326246</Real>
++        <Real Name="Z">15.337114133361936</Real>
 +      </Vector>
 +      <Vector Name="Time 0.016000 Step 16 F[1]">
-         <Real Name="X">-254.29369565940516</Real>
-         <Real Name="Y">-58.536578272057184</Real>
-         <Real Name="Z">-6.7734573229946236</Real>
++        <Real Name="X">-254.29369565939146</Real>
++        <Real Name="Y">-58.536578272055259</Real>
++        <Real Name="Z">-6.7734573229913586</Real>
 +      </Vector>
 +      <Vector Name="Time 0.016000 Step 16 F[2]">
-         <Real Name="X">344.4276213065869</Real>
-         <Real Name="Y">-6.5739101403710727</Real>
-         <Real Name="Z">100.47506026189879</Real>
++        <Real Name="X">344.4276213065715</Real>
++        <Real Name="Y">-6.5739101403722024</Real>
++        <Real Name="Z">100.47506026189447</Real>
 +      </Vector>
 +      <Vector Name="Time 0.016000 Step 16 F[3]">
-         <Real Name="X">-101.0706271290699</Real>
-         <Real Name="Y">74.169910674462457</Real>
-         <Real Name="Z">-109.0387170722688</Real>
++        <Real Name="X">-101.07062712906634</Real>
++        <Real Name="Y">74.169910674460084</Real>
++        <Real Name="Z">-109.03871707226504</Real>
 +      </Vector>
 +      <Vector Name="Time 0.020000 Step 20 F[0]">
-         <Real Name="X">-1.1487885143590644</Real>
-         <Real Name="Y">0.94425702845751291</Real>
-         <Real Name="Z">-1.6210190752303417</Real>
++        <Real Name="X">-1.148788514363047</Real>
++        <Real Name="Y">0.94425702846078807</Real>
++        <Real Name="Z">-1.621019075235957</Real>
 +      </Vector>
 +      <Vector Name="Time 0.020000 Step 20 F[1]">
-         <Real Name="X">-151.33525622409871</Real>
-         <Real Name="Y">-42.597492262017155</Real>
-         <Real Name="Z">10.276839538853924</Real>
++        <Real Name="X">-151.33525622407006</Real>
++        <Real Name="Y">-42.597492262012373</Real>
++        <Real Name="Z">10.276839538860115</Real>
 +      </Vector>
 +      <Vector Name="Time 0.020000 Step 20 F[2]">
-         <Real Name="X">223.13701175223176</Real>
-         <Real Name="Y">-11.261469518677318</Real>
-         <Real Name="Z">70.282144298824107</Real>
++        <Real Name="X">223.13701175219981</Real>
++        <Real Name="Y">-11.26146951868024</Real>
++        <Real Name="Z">70.282144298815567</Real>
 +      </Vector>
 +      <Vector Name="Time 0.020000 Step 20 F[3]">
-         <Real Name="X">-70.652967013773974</Real>
-         <Real Name="Y">52.914704752236958</Real>
-         <Real Name="Z">-78.937964762447692</Real>
++        <Real Name="X">-70.652967013766684</Real>
++        <Real Name="Y">52.914704752231827</Real>
++        <Real Name="Z">-78.937964762439719</Real>
 +      </Vector>
 +      <Vector Name="Time 0.024000 Step 24 F[0]">
-         <Real Name="X">-14.792631934835677</Real>
-         <Real Name="Y">12.084501645102721</Real>
-         <Real Name="Z">-20.985856455180208</Real>
++        <Real Name="X">-14.792631934839411</Real>
++        <Real Name="Y">12.084501645105796</Real>
++        <Real Name="Z">-20.985856455185452</Real>
 +      </Vector>
 +      <Vector Name="Time 0.024000 Step 24 F[1]">
-         <Real Name="X">-38.508823233299573</Real>
-         <Real Name="Y">-21.773206178588957</Real>
-         <Real Name="Z">28.978581607086344</Real>
++        <Real Name="X">-38.508823233273233</Real>
++        <Real Name="Y">-21.773206178584019</Real>
++        <Real Name="Z">28.978581607091613</Real>
 +      </Vector>
 +      <Vector Name="Time 0.024000 Step 24 F[2]">
-         <Real Name="X">92.311360568250592</Real>
-         <Real Name="Y">-20.243825211273709</Real>
-         <Real Name="Z">37.199039838943072</Real>
++        <Real Name="X">92.311360568221545</Real>
++        <Real Name="Y">-20.243825211276956</Real>
++        <Real Name="Z">37.199039838935668</Real>
 +      </Vector>
 +      <Vector Name="Time 0.024000 Step 24 F[3]">
-         <Real Name="X">-39.009905400115343</Real>
-         <Real Name="Y">29.932529744759947</Real>
-         <Real Name="Z">-45.191764990849208</Real>
++        <Real Name="X">-39.009905400108892</Real>
++        <Real Name="Y">29.932529744755179</Real>
++        <Real Name="Z">-45.191764990841826</Real>
 +      </Vector>
 +      <Vector Name="Time 0.028000 Step 28 F[0]">
-         <Real Name="X">-28.233093664967807</Real>
-         <Real Name="Y">22.92859444472731</Real>
-         <Real Name="Z">-40.223321948721754</Real>
++        <Real Name="X">-28.233093664971499</Real>
++        <Real Name="Y">22.928594444730347</Real>
++        <Real Name="Z">-40.223321948726912</Real>
 +      </Vector>
 +      <Vector Name="Time 0.028000 Step 28 F[1]">
-         <Real Name="X">73.021460992109311</Real>
-         <Real Name="Y">2.4029740292514781</Real>
-         <Real Name="Z">46.15342995506451</Real>
++        <Real Name="X">73.021460992134578</Real>
++        <Real Name="Y">2.4029740292567254</Real>
++        <Real Name="Z">46.153429955069178</Real>
 +      </Vector>
 +      <Vector Name="Time 0.028000 Step 28 F[2]">
-         <Real Name="X">-36.130408162563846</Real>
-         <Real Name="Y">-32.154073025311916</Real>
-         <Real Name="Z">4.4691445584323635</Real>
++        <Real Name="X">-36.130408162591316</Real>
++        <Real Name="Y">-32.154073025315604</Real>
++        <Real Name="Z">4.4691445584258105</Real>
 +      </Vector>
 +      <Vector Name="Time 0.028000 Step 28 F[3]">
-         <Real Name="X">-8.6579591645776546</Real>
-         <Real Name="Y">6.8225045513331253</Real>
-         <Real Name="Z">-10.399252564775118</Real>
++        <Real Name="X">-8.6579591645717695</Real>
++        <Real Name="Y">6.8225045513285316</Real>
++        <Real Name="Z">-10.399252564768073</Real>
 +      </Vector>
 +      <Vector Name="Time 0.032000 Step 32 F[0]">
-         <Real Name="X">-39.97789984097686</Real>
-         <Real Name="Y">32.241671576985929</Real>
-         <Real Name="Z">-57.119576941001135</Real>
++        <Real Name="X">-39.977899840980314</Real>
++        <Real Name="Y">32.24167157698875</Real>
++        <Real Name="Z">-57.119576941005896</Real>
 +      </Vector>
 +      <Vector Name="Time 0.032000 Step 32 F[1]">
-         <Real Name="X">173.85168027656647</Real>
-         <Real Name="Y">27.749078001711162</Real>
-         <Real Name="Z">59.637350464835151</Real>
++        <Real Name="X">173.85168027659</Real>
++        <Real Name="Y">27.749078001716573</Real>
++        <Real Name="Z">59.637350464838917</Real>
 +      </Vector>
 +      <Vector Name="Time 0.032000 Step 32 F[2]">
-         <Real Name="X">-152.36496859639158</Real>
-         <Real Name="Y">-45.012821768681185</Real>
-         <Real Name="Z">-25.518087185941891</Real>
++        <Real Name="X">-152.36496859641699</Real>
++        <Real Name="Y">-45.012821768685001</Real>
++        <Real Name="Z">-25.5180871859476</Real>
 +      </Vector>
 +      <Vector Name="Time 0.032000 Step 32 F[3]">
-         <Real Name="X">18.49118816080194</Real>
-         <Real Name="Y">-14.977927810015906</Real>
-         <Real Name="Z">23.000313662107878</Real>
++        <Real Name="X">18.491188160807297</Real>
++        <Real Name="Y">-14.977927810020326</Real>
++        <Real Name="Z">23.000313662114582</Real>
 +      </Vector>
 +      <Vector Name="Time 0.036000 Step 36 F[0]">
-         <Real Name="X">-48.881612139197181</Real>
-         <Real Name="Y">39.065295762634229</Real>
-         <Real Name="Z">-69.938279010926436</Real>
++        <Real Name="X">-48.881612139200257</Real>
++        <Real Name="Y">39.065295762636694</Real>
++        <Real Name="Z">-69.938279010930557</Real>
 +      </Vector>
 +      <Vector Name="Time 0.036000 Step 36 F[1]">
-         <Real Name="X">256.91883657637214</Real>
-         <Real Name="Y">51.827136837357685</Real>
-         <Real Name="Z">68.316578380756738</Real>
++        <Real Name="X">256.91883657639391</Real>
++        <Real Name="Y">51.827136837363227</Real>
++        <Real Name="Z">68.316578380759381</Real>
 +      </Vector>
 +      <Vector Name="Time 0.036000 Step 36 F[2]">
-         <Real Name="X">-249.22935104064919</Real>
-         <Real Name="Y">-56.609707958967739</Real>
-         <Real Name="Z">-51.324820635516417</Real>
++        <Real Name="X">-249.22935104067298</Real>
++        <Real Name="Y">-56.609707958971349</Real>
++        <Real Name="Z">-51.324820635521533</Real>
 +      </Vector>
 +      <Vector Name="Time 0.036000 Step 36 F[3]">
-         <Real Name="X">41.192126603474236</Real>
-         <Real Name="Y">-34.282724641024174</Real>
-         <Real Name="Z">52.946521265686115</Real>
++        <Real Name="X">41.192126603479288</Real>
++        <Real Name="Y">-34.282724641028572</Real>
++        <Real Name="Z">52.946521265692709</Real>
 +      </Vector>
 +      <Vector Name="Time 0.040000 Step 40 F[0]">
-         <Real Name="X">-54.149503868648566</Real>
-         <Real Name="Y">42.755819150349467</Real>
-         <Real Name="Z">-77.469738353419245</Real>
++        <Real Name="X">-54.149503868651053</Real>
++        <Real Name="Y">42.755819150351471</Real>
++        <Real Name="Z">-77.469738353422414</Real>
 +      </Vector>
 +      <Vector Name="Time 0.040000 Step 40 F[1]">
-         <Real Name="X">317.55427307025371</Real>
-         <Real Name="Y">72.337791383415436</Real>
-         <Real Name="Z">71.871971514297229</Real>
++        <Real Name="X">317.55427307027304</Real>
++        <Real Name="Y">72.337791383420608</Real>
++        <Real Name="Z">71.871971514298565</Real>
 +      </Vector>
 +      <Vector Name="Time 0.040000 Step 40 F[2]">
-         <Real Name="X">-322.20372287496099</Real>
-         <Real Name="Y">-64.891170578042846</Real>
-         <Real Name="Z">-72.27849537849221</Real>
++        <Real Name="X">-322.2037228749827</Real>
++        <Real Name="Y">-64.891170578045518</Real>
++        <Real Name="Z">-72.278495378496842</Real>
 +      </Vector>
 +      <Vector Name="Time 0.040000 Step 40 F[3]">
-         <Real Name="X">58.798953673355882</Real>
-         <Real Name="Y">-50.202439955722056</Real>
-         <Real Name="Z">77.876262217614226</Real>
++        <Real Name="X">58.798953673360693</Real>
++        <Real Name="Y">-50.202439955726554</Real>
++        <Real Name="Z">77.876262217620692</Real>
 +      </Vector>
 +      <Vector Name="Time 0.044000 Step 44 F[0]">
-         <Real Name="X">-55.30214765388515</Real>
-         <Real Name="Y">42.989286333517626</Real>
-         <Real Name="Z">-79.007731097948096</Real>
++        <Real Name="X">-55.302147653884383</Real>
++        <Real Name="Y">42.989286333516887</Real>
++        <Real Name="Z">-79.00773109794666</Real>
 +      </Vector>
 +      <Vector Name="Time 0.044000 Step 44 F[1]">
-         <Real Name="X">353.20119046090213</Real>
-         <Real Name="Y">87.420405578975277</Real>
-         <Real Name="Z">70.415575063717398</Real>
++        <Real Name="X">353.2011904609036</Real>
++        <Real Name="Y">87.420405578976244</Real>
++        <Real Name="Z">70.41557506371511</Real>
 +      </Vector>
 +      <Vector Name="Time 0.044000 Step 44 F[2]">
-         <Real Name="X">-369.03341302494488</Real>
-         <Real Name="Y">-68.257569667105031</Real>
-         <Real Name="Z">-88.161363599526965</Real>
++        <Real Name="X">-369.03341302494869</Real>
++        <Real Name="Y">-68.257569667103468</Real>
++        <Real Name="Z">-88.161363599528485</Real>
 +      </Vector>
 +      <Vector Name="Time 0.044000 Step 44 F[3]">
-         <Real Name="X">71.134370217927881</Real>
-         <Real Name="Y">-62.152122245387858</Real>
-         <Real Name="Z">96.753519633757662</Real>
++        <Real Name="X">71.134370217929529</Real>
++        <Real Name="Y">-62.152122245389663</Real>
++        <Real Name="Z">96.753519633760035</Real>
 +      </Vector>
 +      <Vector Name="Time 0.048000 Step 48 F[0]">
-         <Real Name="X">-52.13141853615042</Real>
-         <Real Name="Y">39.743487788303739</Real>
-         <Real Name="Z">-74.292830164081394</Real>
++        <Real Name="X">-52.131418536148871</Real>
++        <Real Name="Y">39.743487788302431</Real>
++        <Real Name="Z">-74.292830164078936</Real>
 +      </Vector>
 +      <Vector Name="Time 0.048000 Step 48 F[1]">
-         <Real Name="X">363.0139691466951</Real>
-         <Real Name="Y">95.826941958577009</Real>
-         <Real Name="Z">64.166722501371126</Real>
++        <Real Name="X">363.01396914668874</Real>
++        <Real Name="Y">95.826941958575318</Real>
++        <Real Name="Z">64.166722501368355</Real>
 +      </Vector>
 +      <Vector Name="Time 0.048000 Step 48 F[2]">
-         <Real Name="X">-389.2319923561339</Real>
-         <Real Name="Y">-65.735297034854213</Real>
-         <Real Name="Z">-98.914481786706801</Real>
++        <Real Name="X">-389.23199235612861</Real>
++        <Real Name="Y">-65.73529703485147</Real>
++        <Real Name="Z">-98.914481786705664</Real>
 +      </Vector>
 +      <Vector Name="Time 0.048000 Step 48 F[3]">
-         <Real Name="X">78.349441745589218</Real>
-         <Real Name="Y">-69.835132712026535</Real>
-         <Real Name="Z">109.04058944941707</Real>
++        <Real Name="X">78.349441745588734</Real>
++        <Real Name="Y">-69.835132712026279</Real>
++        <Real Name="Z">109.04058944941626</Real>
 +      </Vector>
 +    </Mdrun>
 +  </Simulation>
 +</ReferenceData>