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