src/gromacs/gmxana/gmx_make_ndx.c

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