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