src/gromacs/gmxana/gmx_make_ndx.c

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