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