src/gromacs/tools/make_ndx.cpp

   1 /*
   2  * This file is part of the GROMACS molecular simulation package.
   3  *
   4  * Copyright (c) 1991-2000, University of Groningen, The Netherlands.
   5  * Copyright (c) 2001-2004, The GROMACS development team.
   6  * Copyright (c) 2012,2013,2014,2015,2016 by the GROMACS development team.
   7  * Copyright (c) 2017,2018,2019,2020,2021, by the GROMACS development team, led by
   8  * Mark Abraham, David van der Spoel, Berk Hess, and Erik Lindahl,
   9  * and including many others, as listed in the AUTHORS file in the
  10  * top-level source directory and at http://www.gromacs.org.
  11  *
  12  * GROMACS is free software; you can redistribute it and/or
  13  * modify it under the terms of the GNU Lesser General Public License
  14  * as published by the Free Software Foundation; either version 2.1
  15  * of the License, or (at your option) any later version.
  16  *
  17  * GROMACS is distributed in the hope that it will be useful,
  18  * but WITHOUT ANY WARRANTY; without even the implied warranty of
  19  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
  20  * Lesser General Public License for more details.
  21  *
  22  * You should have received a copy of the GNU Lesser General Public
  23  * License along with GROMACS; if not, see
  24  * http://www.gnu.org/licenses, or write to the Free Software Foundation,
  25  * Inc., 51 Franklin Street, Fifth Floor, Boston, MA  02110-1301  USA.
  26  *
  27  * If you want to redistribute modifications to GROMACS, please
  28  * consider that scientific software is very special. Version
  29  * control is crucial - bugs must be traceable. We will be happy to
  30  * consider code for inclusion in the official distribution, but
  31  * derived work must not be called official GROMACS. Details are found
  32  * in the README & COPYING files - if they are missing, get the
  33  * official version at http://www.gromacs.org.
  34  *
  35  * To help us fund GROMACS development, we humbly ask that you cite
  36  * the research papers on the package. Check out http://www.gromacs.org.
  37  */
  38 #include "gmxpre.h"
  39
  40 #include "gromacs/utility/arrayref.h"
  41 #include "make_ndx.h"
  42
  43 #include <array>
  44 #include <cctype>
  45 #include <cstring>
  46
  47 #include <algorithm>
  48 #include <array>
  49 #include <vector>
  50
  51 #include "gromacs/commandline/pargs.h"
  52 #include "gromacs/fileio/confio.h"
  53 #include "gromacs/math/vec.h"
  54 #include "gromacs/topology/block.h"
  55 #include "gromacs/topology/index.h"
  56 #include "gromacs/topology/mtop_util.h"
  57 #include "gromacs/topology/topology.h"
  58 #include "gromacs/utility/arrayref.h"
  59 #include "gromacs/utility/arraysize.h"
  60 #include "gromacs/utility/cstringutil.h"
  61 #include "gromacs/utility/fatalerror.h"
  62 #include "gromacs/utility/futil.h"
  63 #include "gromacs/utility/smalloc.h"
  64
  65 /* It's not nice to have size limits, but we should not spend more time
  66  * on this ancient tool, but instead use the new selection library.
  67  */
  68 #define MAXNAMES 1024
  69 #define NAME_LEN 1024
  70 static const int NOTSET = -92637;
  71
  72 // NOLINTNEXTLINE(cppcoreguidelines-avoid-non-const-global-variables)
  73 static gmx_bool bCase = FALSE;
  74
  75 static int or_groups(int nr1, const int* at1, int nr2, const int* at2, int* nr, int* at)
  76 {
  77     int      i1, i2, max = 0;
  78     gmx_bool bNotIncr;
  79
  80     *nr = 0;
  81
  82     bNotIncr = FALSE;
  83     for (i1 = 0; i1 < nr1; i1++)
  84     {
  85         if ((i1 > 0) && (at1[i1] <= max))
  86         {
  87             bNotIncr = TRUE;
  88         }
  89         max = at1[i1];
  90     }
  91     for (i1 = 0; i1 < nr2; i1++)
  92     {
  93         if ((i1 > 0) && (at2[i1] <= max))
  94         {
  95             bNotIncr = TRUE;
  96         }
  97         max = at2[i1];
  98     }
  99
 100     if (bNotIncr)
 101     {
 102         printf("One of your groups is not ascending\n");
 103     }
 104     else
 105     {
 106         i1  = 0;
 107         i2  = 0;
 108         *nr = 0;
 109         while ((i1 < nr1) || (i2 < nr2))
 110         {
 111             if ((i2 == nr2) || ((i1 < nr1) && (at1[i1] < at2[i2])))
 112             {
 113                 at[*nr] = at1[i1];
 114                 (*nr)++;
 115                 i1++;
 116             }
 117             else
 118             {
 119                 if ((i2 < nr2) && ((i1 == nr1) || (at1[i1] > at2[i2])))
 120                 {
 121                     at[*nr] = at2[i2];
 122                     (*nr)++;
 123                 }
 124                 i2++;
 125             }
 126         }
 127
 128         printf("Merged two groups with OR: %d %d -> %d\n", nr1, nr2, *nr);
 129     }
 130
 131     return *nr;
 132 }
 133
 134 static int and_groups(int nr1, const int* at1, int nr2, const int* at2, int* nr, int* at)
 135 {
 136     int i1, i2;
 137
 138     *nr = 0;
 139     for (i1 = 0; i1 < nr1; i1++)
 140     {
 141         for (i2 = 0; i2 < nr2; i2++)
 142         {
 143             if (at1[i1] == at2[i2])
 144             {
 145                 at[*nr] = at1[i1];
 146                 (*nr)++;
 147             }
 148         }
 149     }
 150
 151     printf("Merged two groups with AND: %d %d -> %d\n", nr1, nr2, *nr);
 152
 153     return *nr;
 154 }
 155
 156 static gmx_bool is_name_char(char c)
 157 {
 158     /* This string should contain all characters that can not be
 159      * the first letter of a name due to the make_ndx syntax.
 160      */
 161     const char* spec = " !&|";
 162
 163     return (c != '\0' && std::strchr(spec, c) == nullptr);
 164 }
 165
 166 static int parse_names(char** string, int* n_names, gmx::ArrayRef<char*> names)
 167 {
 168     int i;
 169
 170     *n_names = 0;
 171     while (is_name_char((*string)[0]) || (*string)[0] == ' ')
 172     {
 173         if (is_name_char((*string)[0]))
 174         {
 175             if (*n_names >= MAXNAMES)
 176             {
 177                 gmx_fatal(FARGS, "To many names: %d\n", *n_names + 1);
 178             }
 179             i = 0;
 180             while (is_name_char((*string)[i]))
 181             {
 182                 names[*n_names][i] = (*string)[i];
 183                 i++;
 184                 if (i > NAME_LEN)
 185                 {
 186                     printf("Name is too long, the maximum is %d characters\n", NAME_LEN);
 187                     return 0;
 188                 }
 189             }
 190             names[*n_names][i] = '\0';
 191             if (!bCase)
 192             {
 193                 upstring(names[*n_names]);
 194             }
 195             *string += i;
 196             (*n_names)++;
 197         }
 198         else
 199         {
 200             (*string)++;
 201         }
 202     }
 203
 204     return *n_names;
 205 }
 206
 207 static gmx_bool parse_int_char(char** string, int* nr, unsigned char* c)
 208 {
 209     char*    orig;
 210     gmx_bool bRet;
 211
 212     orig = *string;
 213
 214     while ((*string)[0] == ' ')
 215     {
 216         (*string)++;
 217     }
 218
 219     bRet = FALSE;
 220
 221     *c = ' ';
 222
 223     if (std::isdigit((*string)[0]))
 224     {
 225         *nr = (*string)[0] - '0';
 226         (*string)++;
 227         while (std::isdigit((*string)[0]))
 228         {
 229             *nr = (*nr) * 10 + (*string)[0] - '0';
 230             (*string)++;
 231         }
 232         if (std::isalpha((*string)[0]))
 233         {
 234             *c = (*string)[0];
 235             (*string)++;
 236         }
 237         /* Check if there is at most one non-digit character */
 238         if (!std::isalnum((*string)[0]))
 239         {
 240             bRet = TRUE;
 241         }
 242         else
 243         {
 244             *string = orig;
 245         }
 246     }
 247     else
 248     {
 249         *nr = NOTSET;
 250     }
 251
 252     return bRet;
 253 }
 254
 255 static gmx_bool parse_int(char** string, int* nr)
 256 {
 257     char*         orig;
 258     unsigned char c;
 259     gmx_bool      bRet;
 260
 261     orig = *string;
 262     bRet = parse_int_char(string, nr, &c);
 263     if (bRet && c != ' ')
 264     {
 265         *string = orig;
 266         bRet    = FALSE;
 267     }
 268
 269     return bRet;
 270 }
 271
 272 static gmx_bool isquote(char c)
 273 {
 274     return (c == '\"');
 275 }
 276
 277 static gmx_bool parse_string(char** string, int* nr, int ngrps, char** grpname)
 278 {
 279     char *s, *sp;
 280     char  c;
 281
 282     while ((*string)[0] == ' ')
 283     {
 284         (*string)++;
 285     }
 286
 287     (*nr) = NOTSET;
 288     if (isquote((*string)[0]))
 289     {
 290         c = (*string)[0];
 291         (*string)++;
 292         s  = gmx_strdup((*string));
 293         sp = std::strchr(s, c);
 294         if (sp != nullptr)
 295         {
 296             (*string) += sp - s + 1;
 297             sp[0] = '\0';
 298             (*nr) = find_group(s, ngrps, grpname);
 299         }
 300     }
 301
 302     return (*nr) != NOTSET;
 303 }
 304
 305 static int select_atomnumbers(char** string, const t_atoms* atoms, int n1, int* nr, int* index, char* gname)
 306 {
 307     char buf[STRLEN];
 308     int  i, up;
 309
 310     *nr = 0;
 311     while ((*string)[0] == ' ')
 312     {
 313         (*string)++;
 314     }
 315     if ((*string)[0] == '-')
 316     {
 317         (*string)++;
 318         parse_int(string, &up);
 319         if ((n1 < 1) || (n1 > atoms->nr) || (up < 1) || (up > atoms->nr))
 320         {
 321             printf("Invalid atom range\n");
 322         }
 323         else
 324         {
 325             for (i = n1 - 1; i <= up - 1; i++)
 326             {
 327                 index[*nr] = i;
 328                 (*nr)++;
 329             }
 330             printf("Found %d atom%s in range %d-%d\n", *nr, (*nr == 1) ? "" : "s", n1, up);
 331             if (n1 == up)
 332             {
 333                 sprintf(buf, "a_%d", n1);
 334             }
 335             else
 336             {
 337                 sprintf(buf, "a_%d-%d", n1, up);
 338             }
 339             std::strcpy(gname, buf);
 340         }
 341     }
 342     else
 343     {
 344         i = n1;
 345         sprintf(gname, "a");
 346         do
 347         {
 348             if ((i - 1 >= 0) && (i - 1 < atoms->nr))
 349             {
 350                 index[*nr] = i - 1;
 351                 (*nr)++;
 352                 sprintf(buf, "_%d", i);
 353                 std::strcat(gname, buf);
 354             }
 355             else
 356             {
 357                 printf("Invalid atom number %d\n", i);
 358                 *nr = 0;
 359             }
 360         } while ((*nr != 0) && (parse_int(string, &i)));
 361     }
 362
 363     return *nr;
 364 }
 365
 366 static int
 367 select_residuenumbers(char** string, const t_atoms* atoms, int n1, unsigned char c, int* nr, int* index, char* gname)
 368 {
 369     char       buf[STRLEN];
 370     int        i, j, up;
 371     t_resinfo* ri;
 372
 373     *nr = 0;
 374     while ((*string)[0] == ' ')
 375     {
 376         (*string)++;
 377     }
 378     if ((*string)[0] == '-')
 379     {
 380         /* Residue number range selection */
 381         if (c != ' ')
 382         {
 383             printf("Error: residue insertion codes can not be used with residue range selection\n");
 384             return 0;
 385         }
 386         (*string)++;
 387         parse_int(string, &up);
 388
 389         for (i = 0; i < atoms->nr; i++)
 390         {
 391             ri = &atoms->resinfo[atoms->atom[i].resind];
 392             for (j = n1; (j <= up); j++)
 393             {
 394                 if (ri->nr == j && (c == ' ' || ri->ic == c))
 395                 {
 396                     index[*nr] = i;
 397                     (*nr)++;
 398                 }
 399             }
 400         }
 401         printf("Found %d atom%s with res.nr. in range %d-%d\n", *nr, (*nr == 1) ? "" : "s", n1, up);
 402         if (n1 == up)
 403         {
 404             sprintf(buf, "r_%d", n1);
 405         }
 406         else
 407         {
 408             sprintf(buf, "r_%d-%d", n1, up);
 409         }
 410         std::strcpy(gname, buf);
 411     }
 412     else
 413     {
 414         /* Individual residue number/insertion code selection */
 415         j = n1;
 416         sprintf(gname, "r");
 417         do
 418         {
 419             for (i = 0; i < atoms->nr; i++)
 420             {
 421                 ri = &atoms->resinfo[atoms->atom[i].resind];
 422                 if (ri->nr == j && ri->ic == c)
 423                 {
 424                     index[*nr] = i;
 425                     (*nr)++;
 426                 }
 427             }
 428             sprintf(buf, "_%d", j);
 429             std::strcat(gname, buf);
 430         } while (parse_int_char(string, &j, &c));
 431     }
 432
 433     return *nr;
 434 }
 435
 436 static int
 437 select_residueindices(char** string, const t_atoms* atoms, int n1, unsigned char c, int* nr, int* index, char* gname)
 438 {
 439     /*this should be similar to select_residuenumbers except select by index (sequential numbering in file)*/
 440     /*resind+1 for 1-indexing*/
 441     char       buf[STRLEN];
 442     int        i, j, up;
 443     t_resinfo* ri;
 444
 445     *nr = 0;
 446     while ((*string)[0] == ' ')
 447     {
 448         (*string)++;
 449     }
 450     if ((*string)[0] == '-')
 451     {
 452         /* Residue number range selection */
 453         if (c != ' ')
 454         {
 455             printf("Error: residue insertion codes can not be used with residue range selection\n");
 456             return 0;
 457         }
 458         (*string)++;
 459         parse_int(string, &up);
 460
 461         for (i = 0; i < atoms->nr; i++)
 462         {
 463             ri = &atoms->resinfo[atoms->atom[i].resind];
 464             for (j = n1; (j <= up); j++)
 465             {
 466                 if (atoms->atom[i].resind + 1 == j && (c == ' ' || ri->ic == c))
 467                 {
 468                     index[*nr] = i;
 469                     (*nr)++;
 470                 }
 471             }
 472         }
 473         printf("Found %d atom%s with resind.+1 in range %d-%d\n", *nr, (*nr == 1) ? "" : "s", n1, up);
 474         if (n1 == up)
 475         {
 476             sprintf(buf, "r_%d", n1);
 477         }
 478         else
 479         {
 480             sprintf(buf, "r_%d-%d", n1, up);
 481         }
 482         std::strcpy(gname, buf);
 483     }
 484     else
 485     {
 486         /* Individual residue number/insertion code selection */
 487         j = n1;
 488         sprintf(gname, "r");
 489         do
 490         {
 491             for (i = 0; i < atoms->nr; i++)
 492             {
 493                 ri = &atoms->resinfo[atoms->atom[i].resind];
 494                 if (atoms->atom[i].resind + 1 == j && ri->ic == c)
 495                 {
 496                     index[*nr] = i;
 497                     (*nr)++;
 498                 }
 499             }
 500             sprintf(buf, "_%d", j);
 501             std::strcat(gname, buf);
 502         } while (parse_int_char(string, &j, &c));
 503     }
 504
 505     return *nr;
 506 }
 507
 508
 509 static gmx_bool
 510 atoms_from_residuenumbers(const t_atoms* atoms, int group, t_blocka* block, int* nr, int* index, char* gname)
 511 {
 512     int i, j, j0, j1, resnr, nres;
 513
 514     j0   = block->index[group];
 515     j1   = block->index[group + 1];
 516     nres = atoms->nres;
 517     for (j = j0; j < j1; j++)
 518     {
 519         if (block->a[j] >= nres)
 520         {
 521             printf("Index %s contains number>nres (%d>%d)\n", gname, block->a[j] + 1, nres);
 522             return FALSE;
 523         }
 524     }
 525     for (i = 0; i < atoms->nr; i++)
 526     {
 527         resnr = atoms->resinfo[atoms->atom[i].resind].nr;
 528         for (j = j0; j < j1; j++)
 529         {
 530             if (block->a[j] + 1 == resnr)
 531             {
 532                 index[*nr] = i;
 533                 (*nr)++;
 534                 break;
 535             }
 536         }
 537     }
 538     printf("Found %d atom%s in %d residues from group %s\n", *nr, (*nr == 1) ? "" : "s", j1 - j0, gname);
 539     return *nr != 0;
 540 }
 541
 542 static gmx_bool comp_name(const char* name, const char* search)
 543 {
 544     gmx_bool matches = TRUE;
 545
 546     // Loop while name and search are not end-of-string and matches is true
 547     for (; *name && *search && matches; name++, search++)
 548     {
 549         if (*search == '?')
 550         {
 551             // still matching, continue to next character
 552             continue;
 553         }
 554         else if (*search == '*')
 555         {
 556             if (*(search + 1))
 557             {
 558                 printf("WARNING: Currently '*' is only supported at the end of an expression\n");
 559             }
 560             // if * is the last char in search string, we have a match,
 561             // otherwise we just failed. Return in either case, we're done.
 562             return (*(search + 1) == '\0');
 563         }
 564
 565         // Compare one character
 566         if (bCase)
 567         {
 568             matches = (*name == *search);
 569         }
 570         else
 571         {
 572             matches = (std::toupper(*name) == std::toupper(*search));
 573         }
 574     }
 575
 576     matches = matches && (*name == '\0' && (*search == '\0' || *search == '*'));
 577
 578     return matches;
 579 }
 580
 581 static int select_chainnames(const t_atoms* atoms, int n_names, gmx::ArrayRef<char*> names, int* nr, int* index)
 582 {
 583     char name[2];
 584     int  j;
 585     int  i;
 586
 587     name[1] = 0;
 588     *nr     = 0;
 589     for (i = 0; i < atoms->nr; i++)
 590     {
 591         name[0] = atoms->resinfo[atoms->atom[i].resind].chainid;
 592         j       = 0;
 593         while (j < n_names && !comp_name(name, names[j]))
 594         {
 595             j++;
 596         }
 597         if (j < n_names)
 598         {
 599             index[*nr] = i;
 600             (*nr)++;
 601         }
 602     }
 603     printf("Found %d atom%s with chain identifier%s", *nr, (*nr == 1) ? "" : "s", (n_names == 1) ? "" : "s");
 604     for (j = 0; (j < n_names); j++)
 605     {
 606         printf(" %s", names[j]);
 607     }
 608     printf("\n");
 609
 610     return *nr;
 611 }
 612
 613 static int select_atomnames(const t_atoms* atoms, int n_names, gmx::ArrayRef<char*> names, int* nr, int* index, gmx_bool bType)
 614 {
 615     char* name;
 616     int   j;
 617     int   i;
 618
 619     *nr = 0;
 620     for (i = 0; i < atoms->nr; i++)
 621     {
 622         if (bType)
 623         {
 624             name = *(atoms->atomtype[i]);
 625         }
 626         else
 627         {
 628             name = *(atoms->atomname[i]);
 629         }
 630         j = 0;
 631         while (j < n_names && !comp_name(name, names[j]))
 632         {
 633             j++;
 634         }
 635         if (j < n_names)
 636         {
 637             index[*nr] = i;
 638             (*nr)++;
 639         }
 640     }
 641     printf("Found %d atoms with %s%s", *nr, bType ? "type" : "name", (n_names == 1) ? "" : "s");
 642     for (j = 0; (j < n_names); j++)
 643     {
 644         printf(" %s", names[j]);
 645     }
 646     printf("\n");
 647
 648     return *nr;
 649 }
 650
 651 static int select_residuenames(const t_atoms* atoms, int n_names, gmx::ArrayRef<char*> names, int* nr, int* index)
 652 {
 653     char* name;
 654     int   j;
 655     int   i;
 656
 657     *nr = 0;
 658     for (i = 0; i < atoms->nr; i++)
 659     {
 660         name = *(atoms->resinfo[atoms->atom[i].resind].name);
 661         j    = 0;
 662         while (j < n_names && !comp_name(name, names[j]))
 663         {
 664             j++;
 665         }
 666         if (j < n_names)
 667         {
 668             index[*nr] = i;
 669             (*nr)++;
 670         }
 671     }
 672     printf("Found %d atoms with residue name%s", *nr, (n_names == 1) ? "" : "s");
 673     for (j = 0; (j < n_names); j++)
 674     {
 675         printf(" %s", names[j]);
 676     }
 677     printf("\n");
 678
 679     return *nr;
 680 }
 681
 682 static void copy2block(int n, const int* index, t_blocka* block)
 683 {
 684     int i, n0;
 685
 686     block->nr++;
 687     n0         = block->nra;
 688     block->nra = n0 + n;
 689     srenew(block->index, block->nr + 1);
 690     block->index[block->nr] = n0 + n;
 691     srenew(block->a, n0 + n);
 692     for (i = 0; (i < n); i++)
 693     {
 694         block->a[n0 + i] = index[i];
 695     }
 696 }
 697
 698 static void make_gname(int n, gmx::ArrayRef<char*> names, char* gname)
 699 {
 700     int i;
 701
 702     std::strcpy(gname, names[0]);
 703     for (i = 1; i < n; i++)
 704     {
 705         std::strcat(gname, "_");
 706         std::strcat(gname, names[i]);
 707     }
 708 }
 709
 710 static void copy_group(int g, t_blocka* block, int* nr, int* index)
 711 {
 712     int i, i0;
 713
 714     i0  = block->index[g];
 715     *nr = block->index[g + 1] - i0;
 716     for (i = 0; i < *nr; i++)
 717     {
 718         index[i] = block->a[i0 + i];
 719     }
 720 }
 721
 722 static void remove_group(int nr, int nr2, t_blocka* block, char*** gn)
 723 {
 724     int   i, j, shift;
 725     char* name;
 726
 727     if (nr2 == NOTSET)
 728     {
 729         nr2 = nr;
 730     }
 731
 732     for (j = 0; j <= nr2 - nr; j++)
 733     {
 734         if ((nr < 0) || (nr >= block->nr))
 735         {
 736             printf("Group %d does not exist\n", nr + j);
 737         }
 738         else
 739         {
 740             shift = block->index[nr + 1] - block->index[nr];
 741             for (i = block->index[nr + 1]; i < block->nra; i++)
 742             {
 743                 block->a[i - shift] = block->a[i];
 744             }
 745
 746             for (i = nr; i < block->nr; i++)
 747             {
 748                 block->index[i] = block->index[i + 1] - shift;
 749             }
 750             name = gmx_strdup((*gn)[nr]);
 751             sfree((*gn)[nr]);
 752             for (i = nr; i < block->nr - 1; i++)
 753             {
 754                 (*gn)[i] = (*gn)[i + 1];
 755             }
 756             block->nr--;
 757             block->nra = block->index[block->nr];
 758             printf("Removed group %d '%s'\n", nr + j, name);
 759             sfree(name);
 760         }
 761     }
 762 }
 763
 764 static void split_group(const t_atoms* atoms, int sel_nr, t_blocka* block, char*** gn, gmx_bool bAtom)
 765 {
 766     char buf[STRLEN], *name;
 767     int  i, resind;
 768     int  a, n0, n1;
 769
 770     printf("Splitting group %d '%s' into %s\n", sel_nr, (*gn)[sel_nr], bAtom ? "atoms" : "residues");
 771
 772     n0 = block->index[sel_nr];
 773     n1 = block->index[sel_nr + 1];
 774     srenew(block->a, block->nra + n1 - n0);
 775     for (i = n0; i < n1; i++)
 776     {
 777         a      = block->a[i];
 778         resind = atoms->atom[a].resind;
 779         name   = *(atoms->resinfo[resind].name);
 780         if (bAtom || (i == n0) || (atoms->atom[block->a[i - 1]].resind != resind))
 781         {
 782             if (i > n0)
 783             {
 784                 block->index[block->nr] = block->nra;
 785             }
 786             block->nr++;
 787             srenew(block->index, block->nr + 1);
 788             srenew(*gn, block->nr);
 789             if (bAtom)
 790             {
 791                 sprintf(buf, "%s_%s_%d", (*gn)[sel_nr], *atoms->atomname[a], a + 1);
 792             }
 793             else
 794             {
 795                 sprintf(buf, "%s_%s_%d", (*gn)[sel_nr], name, atoms->resinfo[resind].nr);
 796             }
 797             (*gn)[block->nr - 1] = gmx_strdup(buf);
 798         }
 799         block->a[block->nra] = a;
 800         block->nra++;
 801     }
 802     block->index[block->nr] = block->nra;
 803 }
 804
 805 static int split_chain(const t_atoms* atoms, const rvec* x, int sel_nr, t_blocka* block, char*** gn)
 806 {
 807     char buf[STRLEN];
 808     int  j, nchain;
 809     int  i, a, natoms, *start = nullptr, *end = nullptr, ca_start, ca_end;
 810     rvec vec;
 811
 812     natoms   = atoms->nr;
 813     nchain   = 0;
 814     ca_start = 0;
 815
 816     while (ca_start < natoms)
 817     {
 818         while ((ca_start < natoms) && std::strcmp(*atoms->atomname[ca_start], "CA") != 0)
 819         {
 820             ca_start++;
 821         }
 822         if (ca_start < natoms)
 823         {
 824             srenew(start, nchain + 1);
 825             srenew(end, nchain + 1);
 826             start[nchain] = ca_start;
 827             while ((start[nchain] > 0)
 828                    && (atoms->atom[start[nchain] - 1].resind == atoms->atom[ca_start].resind))
 829             {
 830                 start[nchain]--;
 831             }
 832
 833             i = ca_start;
 834             do
 835             {
 836                 ca_end = i;
 837                 do
 838                 {
 839                     i++;
 840                 } while ((i < natoms) && std::strcmp(*atoms->atomname[i], "CA") != 0);
 841                 if (i < natoms)
 842                 {
 843                     rvec_sub(x[ca_end], x[i], vec);
 844                 }
 845                 else
 846                 {
 847                     break;
 848                 }
 849             } while (norm(vec) < 0.45);
 850
 851             end[nchain] = ca_end;
 852             while ((end[nchain] + 1 < natoms)
 853                    && (atoms->atom[end[nchain] + 1].resind == atoms->atom[ca_end].resind))
 854             {
 855                 end[nchain]++;
 856             }
 857             ca_start = end[nchain] + 1;
 858             nchain++;
 859         }
 860     }
 861     if (nchain == 1)
 862     {
 863         printf("Found 1 chain, will not split\n");
 864     }
 865     else
 866     {
 867         printf("Found %d chains\n", nchain);
 868     }
 869     for (j = 0; j < nchain; j++)
 870     {
 871         printf("%d:%6d atoms (%d to %d)\n", j + 1, end[j] - start[j] + 1, start[j] + 1, end[j] + 1);
 872     }
 873
 874     if (nchain > 1)
 875     {
 876         srenew(block->a, block->nra + block->index[sel_nr + 1] - block->index[sel_nr]);
 877         for (j = 0; j < nchain; j++)
 878         {
 879             block->nr++;
 880             srenew(block->index, block->nr + 1);
 881             srenew(*gn, block->nr);
 882             sprintf(buf, "%s_chain%d", (*gn)[sel_nr], j + 1);
 883             (*gn)[block->nr - 1] = gmx_strdup(buf);
 884             for (i = block->index[sel_nr]; i < block->index[sel_nr + 1]; i++)
 885             {
 886                 a = block->a[i];
 887                 if ((a >= start[j]) && (a <= end[j]))
 888                 {
 889                     block->a[block->nra] = a;
 890                     block->nra++;
 891                 }
 892             }
 893             block->index[block->nr] = block->nra;
 894             if (block->index[block->nr - 1] == block->index[block->nr])
 895             {
 896                 remove_group(block->nr - 1, NOTSET, block, gn);
 897             }
 898         }
 899     }
 900     sfree(start);
 901     sfree(end);
 902
 903     return nchain;
 904 }
 905
 906 static gmx_bool check_have_atoms(const t_atoms* atoms, char* string)
 907 {
 908     if (atoms == nullptr)
 909     {
 910         printf("Can not process '%s' without atom info, use option -f\n", string);
 911         return FALSE;
 912     }
 913     else
 914     {
 915         return TRUE;
 916     }
 917 }
 918
 919 static gmx_bool parse_entry(char**               string,
 920                             int                  natoms,
 921                             const t_atoms*       atoms,
 922                             t_blocka*            block,
 923                             char***              gn,
 924                             int*                 nr,
 925                             int*                 index,
 926                             char*                gname,
 927                             gmx::ArrayRef<char*> entryNames)
 928 {
 929     char*         ostring;
 930     int           j, n_names, sel_nr1;
 931     int           i, nr1, *index1;
 932     unsigned char c;
 933     gmx_bool      bRet, bCompl;
 934
 935     bRet    = FALSE;
 936     sel_nr1 = NOTSET;
 937
 938     while (*string[0] == ' ')
 939     {
 940         (*string)++;
 941     }
 942
 943     if ((*string)[0] == '!')
 944     {
 945         bCompl = TRUE;
 946         (*string)++;
 947         while (*string[0] == ' ')
 948         {
 949             (*string)++;
 950         }
 951     }
 952     else
 953     {
 954         bCompl = FALSE;
 955     }
 956
 957     ostring = *string;
 958
 959     if (parse_int(string, &sel_nr1) || parse_string(string, &sel_nr1, block->nr, *gn))
 960     {
 961         if ((sel_nr1 >= 0) && (sel_nr1 < block->nr))
 962         {
 963             copy_group(sel_nr1, block, nr, index);
 964             std::strcpy(gname, (*gn)[sel_nr1]);
 965             printf("Copied index group %d '%s'\n", sel_nr1, (*gn)[sel_nr1]);
 966             bRet = TRUE;
 967         }
 968         else
 969         {
 970             printf("Group %d does not exist\n", sel_nr1);
 971         }
 972     }
 973     else if ((*string)[0] == 'a')
 974     {
 975         (*string)++;
 976         if (check_have_atoms(atoms, ostring))
 977         {
 978             if (parse_int(string, &sel_nr1))
 979             {
 980                 bRet = (select_atomnumbers(string, atoms, sel_nr1, nr, index, gname) != 0);
 981             }
 982             else if (parse_names(string, &n_names, entryNames))
 983             {
 984                 bRet = (select_atomnames(atoms, n_names, entryNames, nr, index, FALSE) != 0);
 985                 make_gname(n_names, entryNames, gname);
 986             }
 987         }
 988     }
 989     else if ((*string)[0] == 't')
 990     {
 991         (*string)++;
 992         if (check_have_atoms(atoms, ostring) && parse_names(string, &n_names, entryNames))
 993         {
 994             if (!(atoms->haveType))
 995             {
 996                 printf("Need a run input file to select atom types\n");
 997             }
 998             else
 999             {
1000                 bRet = (select_atomnames(atoms, n_names, entryNames, nr, index, TRUE) != 0);
1001                 make_gname(n_names, entryNames, gname);
1002             }
1003         }
1004     }
1005     else if (std::strncmp(*string, "res", 3) == 0)
1006     {
1007         (*string) += 3;
1008         if (check_have_atoms(atoms, ostring) && parse_int(string, &sel_nr1) && (sel_nr1 >= 0)
1009             && (sel_nr1 < block->nr))
1010         {
1011             bRet = atoms_from_residuenumbers(atoms, sel_nr1, block, nr, index, (*gn)[sel_nr1]);
1012             sprintf(gname, "atom_%s", (*gn)[sel_nr1]);
1013         }
1014     }
1015     else if (std::strncmp(*string, "ri", 2) == 0)
1016     {
1017         (*string) += 2;
1018         if (check_have_atoms(atoms, ostring) && parse_int_char(string, &sel_nr1, &c))
1019         {
1020             bRet = (select_residueindices(string, atoms, sel_nr1, c, nr, index, gname) != 0);
1021         }
1022     }
1023     else if ((*string)[0] == 'r')
1024     {
1025         (*string)++;
1026         if (check_have_atoms(atoms, ostring))
1027         {
1028             if (parse_int_char(string, &sel_nr1, &c))
1029             {
1030                 bRet = (select_residuenumbers(string, atoms, sel_nr1, c, nr, index, gname) != 0);
1031             }
1032             else if (parse_names(string, &n_names, entryNames))
1033             {
1034                 bRet = (select_residuenames(atoms, n_names, entryNames, nr, index) != 0);
1035                 make_gname(n_names, entryNames, gname);
1036             }
1037         }
1038     }
1039     else if (std::strncmp(*string, "chain", 5) == 0)
1040     {
1041         (*string) += 5;
1042         if (check_have_atoms(atoms, ostring) && parse_names(string, &n_names, entryNames))
1043         {
1044             bRet = (select_chainnames(atoms, n_names, entryNames, nr, index) != 0);
1045             sprintf(gname, "ch%s", entryNames[0]);
1046             for (i = 1; i < n_names; i++)
1047             {
1048                 std::strcat(gname, entryNames[i]);
1049             }
1050         }
1051     }
1052     if (bRet && bCompl)
1053     {
1054         snew(index1, natoms - *nr);
1055         nr1 = 0;
1056         for (i = 0; i < natoms; i++)
1057         {
1058             j = 0;
1059             while ((j < *nr) && (index[j] != i))
1060             {
1061                 j++;
1062             }
1063             if (j == *nr)
1064             {
1065                 if (nr1 >= natoms - *nr)
1066                 {
1067                     printf("There are double atoms in your index group\n");
1068                     break;
1069                 }
1070                 index1[nr1] = i;
1071                 nr1++;
1072             }
1073         }
1074         *nr = nr1;
1075         for (i = 0; i < nr1; i++)
1076         {
1077             index[i] = index1[i];
1078         }
1079         sfree(index1);
1080
1081         for (i = std::strlen(gname) + 1; i > 0; i--)
1082         {
1083             gname[i] = gname[i - 1];
1084         }
1085         gname[0] = '!';
1086         printf("Complemented group: %d atoms\n", *nr);
1087     }
1088
1089     return bRet;
1090 }
1091
1092 static void list_residues(const t_atoms* atoms)
1093 {
1094     int      i, j, start, end, prev_resind, resind;
1095     gmx_bool bDiff;
1096
1097     /* Print all the residues, assuming continuous resnr count */
1098     start       = atoms->atom[0].resind;
1099     prev_resind = start;
1100     for (i = 0; i < atoms->nr; i++)
1101     {
1102         resind = atoms->atom[i].resind;
1103         if ((resind != prev_resind) || (i == atoms->nr - 1))
1104         {
1105             if ((bDiff = (std::strcmp(*atoms->resinfo[resind].name, *atoms->resinfo[start].name) != 0))
1106                 || (i == atoms->nr - 1))
1107             {
1108                 if (bDiff)
1109                 {
1110                     end = prev_resind;
1111                 }
1112                 else
1113                 {
1114                     end = resind;
1115                 }
1116                 if (end < start + 3)
1117                 {
1118                     for (j = start; j <= end; j++)
1119                     {
1120                         printf("%4d %-5s", j + 1, *(atoms->resinfo[j].name));
1121                     }
1122                 }
1123                 else
1124                 {
1125                     printf(" %4d - %4d %-5s  ", start + 1, end + 1, *(atoms->resinfo[start].name));
1126                 }
1127                 start = resind;
1128             }
1129         }
1130         prev_resind = resind;
1131     }
1132     printf("\n");
1133 }
1134
1135 static void edit_index(int natoms, const t_atoms* atoms, const rvec* x, t_blocka* block, char*** gn, gmx_bool bVerbose)
1136 {
1137     char*    ostring;
1138     char     inp_string[STRLEN], *string;
1139     char     gname[STRLEN * 3], gname1[STRLEN], gname2[STRLEN];
1140     int      i, i0, i1, sel_nr, sel_nr2, newgroup;
1141     int      nr, nr1, nr2, *index, *index1, *index2;
1142     gmx_bool bAnd, bOr, bPrintOnce;
1143
1144     string = nullptr;
1145
1146     snew(index, natoms);
1147     snew(index1, natoms);
1148     snew(index2, natoms);
1149
1150     newgroup   = NOTSET;
1151     bPrintOnce = TRUE;
1152     std::array<char*, MAXNAMES> entryNames;
1153     for (auto& name : entryNames)
1154     {
1155         snew(name, NAME_LEN + 1);
1156     }
1157
1158     do
1159     {
1160         gname1[0] = '\0';
1161         if (bVerbose || bPrintOnce || newgroup != NOTSET)
1162         {
1163             printf("\n");
1164             if (bVerbose || bPrintOnce || newgroup == NOTSET)
1165             {
1166                 i0 = 0;
1167                 i1 = block->nr;
1168             }
1169             else
1170             {
1171                 i0 = newgroup;
1172                 i1 = newgroup + 1;
1173             }
1174             for (i = i0; i < i1; i++)
1175             {
1176                 printf("%3d %-20s: %5d atoms\n", i, (*gn)[i], block->index[i + 1] - block->index[i]);
1177             }
1178             newgroup = NOTSET;
1179         }
1180         if (bVerbose || bPrintOnce)
1181         {
1182             printf("\n");
1183             printf(" nr : group      '!': not  'name' nr name   'splitch' nr    Enter: list "
1184                    "groups\n");
1185             printf(" 'a': atom       '&': and  'del' nr         'splitres' nr   'l': list "
1186                    "residues\n");
1187             printf(" 't': atom type  '|': or   'keep' nr        'splitat' nr    'h': help\n");
1188             printf(" 'r': residue              'res' nr         'chain' char\n");
1189             printf(" \"name\": group             'case': case %s         'q': save and quit\n",
1190                    bCase ? "insensitive" : "sensitive  ");
1191             printf(" 'ri': residue index\n");
1192             bPrintOnce = FALSE;
1193         }
1194         printf("\n");
1195         printf("> ");
1196         if (nullptr == fgets(inp_string, STRLEN, stdin))
1197         {
1198             gmx_fatal(FARGS, "Error reading user input");
1199         }
1200         inp_string[std::strlen(inp_string) - 1] = 0;
1201         printf("\n");
1202         string = inp_string;
1203         while (string[0] == ' ')
1204         {
1205             string++;
1206         }
1207
1208         ostring = string;
1209         nr      = 0;
1210         if (string[0] == 'h')
1211         {
1212             printf(" nr                : selects an index group by number or quoted string.\n");
1213             printf("                     The string is first matched against the whole group "
1214                    "name,\n");
1215             printf("                     then against the beginning and finally against an\n");
1216             printf("                     arbitrary substring. A multiple match is an error.\n");
1217
1218             printf(" 'a' nr1 [nr2 ...] : selects atoms, atom numbering starts at 1.\n");
1219             printf(" 'a' nr1 - nr2     : selects atoms in the range from nr1 to nr2.\n");
1220             printf(" 'a' name1[*] [name2[*] ...] : selects atoms by name(s), '?' matches any "
1221                    "char,\n");
1222             printf("                               wildcard '*' allowed at the end of a name.\n");
1223             printf(" 't' type1[*] [type2[*] ...] : as 'a', but for type, run input file "
1224                    "required.\n");
1225             printf(" 'r' nr1[ic1] [nr2[ic2] ...] : selects residues by number and insertion "
1226                    "code.\n");
1227             printf(" 'r' nr1 - nr2               : selects residues in the range from nr1 to "
1228                    "nr2.\n");
1229             printf(" 'r' name1[*] [name2[*] ...] : as 'a', but for residue names.\n");
1230             printf(" 'ri' nr1 - nr2              : selects residue indices, 1-indexed, (as opposed "
1231                    "to numbers) in the range from nr1 to nr2.\n");
1232             printf(" 'chain' ch1 [ch2 ...]       : selects atoms by chain identifier(s),\n");
1233             printf("                               not available with a .gro file as input.\n");
1234             printf(" !                 : takes the complement of a group with respect to all\n");
1235             printf("                     the atoms in the input file.\n");
1236             printf(" & |               : AND and OR, can be placed between any of the options\n");
1237             printf("                     above, the input is processed from left to right.\n");
1238             printf(" 'name' nr name    : rename group nr to name.\n");
1239             printf(" 'del' nr1 [- nr2] : deletes one group or groups in the range from nr1 to "
1240                    "nr2.\n");
1241             printf(" 'keep' nr         : deletes all groups except nr.\n");
1242             printf(" 'case'            : make all name compares case (in)sensitive.\n");
1243             printf(" 'splitch' nr      : split group into chains using CA distances.\n");
1244             printf(" 'splitres' nr     : split group into residues.\n");
1245             printf(" 'splitat' nr      : split group into atoms.\n");
1246             printf(" 'res' nr          : interpret numbers in group as residue numbers\n");
1247             printf(" Enter             : list the currently defined groups and commands\n");
1248             printf(" 'l'               : list the residues.\n");
1249             printf(" 'h'               : show this help.\n");
1250             printf(" 'q'               : save and quit.\n");
1251             printf("\n");
1252             printf(" Examples:\n");
1253             printf(" > 2 | 4 & r 3-5\n");
1254             printf(" selects all atoms from group 2 and 4 that have residue numbers 3, 4 or 5\n");
1255             printf(" > a C* & !a C CA\n");
1256             printf(" selects all atoms starting with 'C' but not the atoms 'C' and 'CA'\n");
1257             printf(" > \"protein\" & ! \"backb\"\n");
1258             printf(" selects all atoms that are in group 'protein' and not in group 'backbone'\n");
1259             if (bVerbose)
1260             {
1261                 printf("\npress Enter ");
1262                 getchar();
1263             }
1264         }
1265         else if (std::strncmp(string, "del", 3) == 0)
1266         {
1267             string += 3;
1268             if (parse_int(&string, &sel_nr))
1269             {
1270                 while (string[0] == ' ')
1271                 {
1272                     string++;
1273                 }
1274                 if (string[0] == '-')
1275                 {
1276                     string++;
1277                     parse_int(&string, &sel_nr2);
1278                 }
1279                 else
1280                 {
1281                     sel_nr2 = NOTSET;
1282                 }
1283                 while (string[0] == ' ')
1284                 {
1285                     string++;
1286                 }
1287                 if (string[0] == '\0')
1288                 {
1289                     remove_group(sel_nr, sel_nr2, block, gn);
1290                 }
1291                 else
1292                 {
1293                     printf("\nSyntax error: \"%s\"\n", string);
1294                 }
1295             }
1296         }
1297         else if (std::strncmp(string, "keep", 4) == 0)
1298         {
1299             string += 4;
1300             if (parse_int(&string, &sel_nr))
1301             {
1302                 remove_group(sel_nr + 1, block->nr - 1, block, gn);
1303                 remove_group(0, sel_nr - 1, block, gn);
1304             }
1305         }
1306         else if (std::strncmp(string, "name", 4) == 0)
1307         {
1308             string += 4;
1309             if (parse_int(&string, &sel_nr))
1310             {
1311                 if ((sel_nr >= 0) && (sel_nr < block->nr))
1312                 {
1313                     sscanf(string, "%s", gname);
1314                     sfree((*gn)[sel_nr]);
1315                     (*gn)[sel_nr] = gmx_strdup(gname);
1316                 }
1317             }
1318         }
1319         else if (std::strncmp(string, "case", 4) == 0)
1320         {
1321             bCase = !bCase;
1322             printf("Switched to case %s\n", bCase ? "sensitive" : "insensitive");
1323         }
1324         else if (string[0] == 'v')
1325         {
1326             bVerbose = !bVerbose;
1327             printf("Turned verbose %s\n", bVerbose ? "on" : "off");
1328         }
1329         else if (string[0] == 'l')
1330         {
1331             if (check_have_atoms(atoms, ostring))
1332             {
1333                 list_residues(atoms);
1334             }
1335         }
1336         else if (std::strncmp(string, "splitch", 7) == 0)
1337         {
1338             string += 7;
1339             if (check_have_atoms(atoms, ostring) && parse_int(&string, &sel_nr) && (sel_nr >= 0)
1340                 && (sel_nr < block->nr))
1341             {
1342                 split_chain(atoms, x, sel_nr, block, gn);
1343             }
1344         }
1345         else if (std::strncmp(string, "splitres", 8) == 0)
1346         {
1347             string += 8;
1348             if (check_have_atoms(atoms, ostring) && parse_int(&string, &sel_nr) && (sel_nr >= 0)
1349                 && (sel_nr < block->nr))
1350             {
1351                 split_group(atoms, sel_nr, block, gn, FALSE);
1352             }
1353         }
1354         else if (std::strncmp(string, "splitat", 7) == 0)
1355         {
1356             string += 7;
1357             if (check_have_atoms(atoms, ostring) && parse_int(&string, &sel_nr) && (sel_nr >= 0)
1358                 && (sel_nr < block->nr))
1359             {
1360                 split_group(atoms, sel_nr, block, gn, TRUE);
1361             }
1362         }
1363         else if (string[0] == '\0')
1364         {
1365             bPrintOnce = TRUE;
1366         }
1367         else if (string[0] != 'q')
1368         {
1369             nr2 = -1;
1370             if (parse_entry(&string, natoms, atoms, block, gn, &nr, index, gname, entryNames))
1371             {
1372                 do
1373                 {
1374                     while (string[0] == ' ')
1375                     {
1376                         string++;
1377                     }
1378
1379                     bAnd = FALSE;
1380                     bOr  = FALSE;
1381                     if (string[0] == '&')
1382                     {
1383                         bAnd = TRUE;
1384                     }
1385                     else if (string[0] == '|')
1386                     {
1387                         bOr = TRUE;
1388                     }
1389
1390                     if (bAnd || bOr)
1391                     {
1392                         string++;
1393                         nr1 = nr;
1394                         for (i = 0; i < nr; i++)
1395                         {
1396                             index1[i] = index[i];
1397                         }
1398                         std::strcpy(gname1, gname);
1399                         if (parse_entry(&string, natoms, atoms, block, gn, &nr2, index2, gname2, entryNames))
1400                         {
1401                             if (bOr)
1402                             {
1403                                 or_groups(nr1, index1, nr2, index2, &nr, index);
1404                                 sprintf(gname, "%s_%s", gname1, gname2);
1405                             }
1406                             else
1407                             {
1408                                 and_groups(nr1, index1, nr2, index2, &nr, index);
1409                                 sprintf(gname, "%s_&_%s", gname1, gname2);
1410                             }
1411                         }
1412                     }
1413                 } while (bAnd || bOr);
1414             }
1415             while (string[0] == ' ')
1416             {
1417                 string++;
1418             }
1419             if (string[0])
1420             {
1421                 printf("\nSyntax error: \"%s\"\n", string);
1422             }
1423             else if (nr > 0)
1424             {
1425                 copy2block(nr, index, block);
1426                 srenew(*gn, block->nr);
1427                 newgroup        = block->nr - 1;
1428                 (*gn)[newgroup] = gmx_strdup(gname);
1429             }
1430             else
1431             {
1432                 printf("Group is empty\n");
1433             }
1434         }
1435     } while (string[0] != 'q');
1436
1437     for (auto& name : entryNames)
1438     {
1439         sfree(name);
1440     }
1441     sfree(index);
1442     sfree(index1);
1443     sfree(index2);
1444 }
1445
1446 static int block2natoms(t_blocka* block)
1447 {
1448     int i, natoms;
1449
1450     natoms = 0;
1451     for (i = 0; i < block->nra; i++)
1452     {
1453         natoms = std::max(natoms, block->a[i]);
1454     }
1455     natoms++;
1456
1457     return natoms;
1458 }
1459
1460 static void merge_blocks(t_blocka* dest, t_blocka* source)
1461 {
1462     int i, nra0, i0;
1463
1464     /* count groups, srenew and fill */
1465     i0   = dest->nr;
1466     nra0 = dest->nra;
1467     dest->nr += source->nr;
1468     srenew(dest->index, dest->nr + 1);
1469     for (i = 0; i < source->nr; i++)
1470     {
1471         dest->index[i0 + i] = nra0 + source->index[i];
1472     }
1473     /* count atoms, srenew and fill */
1474     dest->nra += source->nra;
1475     srenew(dest->a, dest->nra);
1476     for (i = 0; i < source->nra; i++)
1477     {
1478         dest->a[nra0 + i] = source->a[i];
1479     }
1480
1481     /* terminate list */
1482     dest->index[dest->nr] = dest->nra;
1483 }
1484
1485 int gmx_make_ndx(int argc, char* argv[])
1486 {
1487     const char* desc[] = { "Index groups are necessary for almost every GROMACS program.",
1488                            "All these programs can generate default index groups. You ONLY",
1489                            "have to use [THISMODULE] when you need SPECIAL index groups.",
1490                            "There is a default index group for the whole system, 9 default",
1491                            "index groups for proteins, and a default index group",
1492                            "is generated for every other residue name.",
1493                            "",
1494                            "When no index file is supplied, also [THISMODULE] will generate the",
1495                            "default groups.",
1496                            "With the index editor you can select on atom, residue and chain names",
1497                            "and numbers.",
1498                            "When a run input file is supplied you can also select on atom type.",
1499                            "You can use boolean operations, you can split groups",
1500                            "into chains, residues or atoms. You can delete and rename groups.",
1501                            "Type 'h' in the editor for more details.",
1502                            "",
1503                            "The atom numbering in the editor and the index file starts at 1.",
1504                            "",
1505                            "The [TT]-twin[tt] switch duplicates all index groups with an offset of",
1506                            "[TT]-natoms[tt], which is useful for Computational Electrophysiology",
1507                            "double-layer membrane setups.",
1508                            "",
1509                            "See also [gmx-select] [TT]-on[tt], which provides an alternative way",
1510                            "for constructing index groups.  It covers nearly all of [THISMODULE]",
1511                            "functionality, and in many cases much more." };
1512
1513     static int      natoms     = 0;
1514     static gmx_bool bVerbose   = FALSE;
1515     static gmx_bool bDuplicate = FALSE;
1516     t_pargs         pa[]       = { { "-natoms",
1517                        FALSE,
1518                        etINT,
1519                        { &natoms },
1520                        "set number of atoms (default: read from coordinate or index file)" },
1521                      { "-twin",
1522                        FALSE,
1523                        etBOOL,
1524                        { &bDuplicate },
1525                        "Duplicate all index groups with an offset of -natoms" },
1526                      { "-verbose", FALSE, etBOOL, { &bVerbose }, "HIDDENVerbose output" } };
1527 #define NPA asize(pa)
1528
1529     gmx_output_env_t* oenv;
1530     const char*       stxfile;
1531     const char*       ndxoutfile;
1532     gmx_bool          bNatoms;
1533     int               j;
1534     t_atoms           atoms;
1535     rvec *            x, *v;
1536     PbcType           pbcType;
1537     matrix            box;
1538     t_blocka *        block, *block2;
1539     char **           gnames, **gnames2;
1540     t_filenm          fnm[] = { { efSTX, "-f", nullptr, ffOPTRD },
1541                        { efNDX, "-n", nullptr, ffOPTRDMULT },
1542                        { efNDX, "-o", nullptr, ffWRITE } };
1543 #define NFILE asize(fnm)
1544
1545     if (!parse_common_args(&argc, argv, 0, NFILE, fnm, NPA, pa, asize(desc), desc, 0, nullptr, &oenv))
1546     {
1547         return 0;
1548     }
1549
1550     stxfile                                     = ftp2fn_null(efSTX, NFILE, fnm);
1551     gmx::ArrayRef<const std::string> ndxInFiles = opt2fnsIfOptionSet("-n", NFILE, fnm);
1552     ndxoutfile                                  = opt2fn("-o", NFILE, fnm);
1553     bNatoms                                     = opt2parg_bSet("-natoms", NPA, pa);
1554
1555     if (!stxfile && ndxInFiles.empty())
1556     {
1557         gmx_fatal(FARGS, "No input files (structure or index)");
1558     }
1559
1560     gmx_mtop_t mtop;
1561     if (stxfile)
1562     {
1563         bool haveFullTopology = false;
1564         fprintf(stderr, "\nReading structure file\n");
1565         readConfAndTopology(stxfile, &haveFullTopology, &mtop, &pbcType, &x, &v, box);
1566         atoms = gmx_mtop_global_atoms(mtop);
1567         if (atoms.pdbinfo == nullptr)
1568         {
1569             snew(atoms.pdbinfo, atoms.nr);
1570         }
1571         natoms  = atoms.nr;
1572         bNatoms = TRUE;
1573     }
1574     else
1575     {
1576         x = nullptr;
1577     }
1578
1579     /* read input file(s) */
1580     block  = new_blocka();
1581     gnames = nullptr;
1582     printf("Going to read %td old index file(s)\n", ndxInFiles.ssize());
1583     if (!ndxInFiles.empty())
1584     {
1585         for (const std::string& ndxInFile : ndxInFiles)
1586         {
1587             block2 = init_index(ndxInFile.c_str(), &gnames2);
1588             srenew(gnames, block->nr + block2->nr);
1589             for (j = 0; j < block2->nr; j++)
1590             {
1591                 gnames[block->nr + j] = gnames2[j];
1592             }
1593             sfree(gnames2);
1594             merge_blocks(block, block2);
1595             sfree(block2->a);
1596             sfree(block2->index);
1597             /*       done_block(block2); */
1598             sfree(block2);
1599         }
1600     }
1601     else
1602     {
1603         snew(gnames, 1);
1604         analyse(&atoms, block, &gnames, FALSE, TRUE);
1605     }
1606
1607     if (!bNatoms)
1608     {
1609         natoms = block2natoms(block);
1610         printf("Counted atom numbers up to %d in index file\n", natoms);
1611     }
1612     edit_index(natoms, &atoms, x, block, &gnames, bVerbose);
1613
1614     write_index(ndxoutfile, block, gnames, bDuplicate, natoms);
1615
1616     return 0;
1617 }