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