src/gromacs/gmxana/gmx_make_ndx.c

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