src/gromacs/gmxana/gmx_make_ndx.c

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