src/gromacs/tools/make_ndx.cpp

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