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