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