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 #include "gmxpre.h"
  38
  39 #include <ctype.h>
  40 #include <string.h>
  41
  42 #include "gromacs/utility/futil.h"
  43 #include "gromacs/legacyheaders/macros.h"
  44 #include "gromacs/fileio/confio.h"
  45 #include "gromacs/legacyheaders/typedefs.h"
  46 #include "gromacs/topology/index.h"
  47 #include "gromacs/math/vec.h"
  48
  49 #include "gromacs/commandline/pargs.h"
  50 #include "gromacs/topology/block.h"
  51 #include "gromacs/utility/cstringutil.h"
  52 #include "gromacs/utility/fatalerror.h"
  53 #include "gromacs/utility/smalloc.h"
  54
  55 /* It's not nice to have size limits, but we should not spend more time
  56  * on this ancient tool, but instead use the new selection library.
  57  */
  58 #define MAXNAMES 1024
  59 #define NAME_LEN 1024
  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  = gmx_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 = gmx_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] = gmx_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                 else
 851                 {
 852                     break;
 853                 }
 854             }
 855             while (norm(vec) < 0.45);
 856
 857             end[nchain] = ca_end;
 858             while ((end[nchain]+1 < natoms) &&
 859                    (atoms->atom[end[nchain]+1].resind == atoms->atom[ca_end].resind))
 860             {
 861                 end[nchain]++;
 862             }
 863             ca_start = end[nchain]+1;
 864             nchain++;
 865         }
 866     }
 867     if (nchain == 1)
 868     {
 869         printf("Found 1 chain, will not split\n");
 870     }
 871     else
 872     {
 873         printf("Found %d chains\n", nchain);
 874     }
 875     for (j = 0; j < nchain; j++)
 876     {
 877         printf("%d:%6d atoms (%d to %d)\n",
 878                j+1, end[j]-start[j]+1, start[j]+1, end[j]+1);
 879     }
 880
 881     if (nchain > 1)
 882     {
 883         srenew(block->a, block->nra+block->index[sel_nr+1]-block->index[sel_nr]);
 884         for (j = 0; j < nchain; j++)
 885         {
 886             block->nr++;
 887             srenew(block->index, block->nr+1);
 888             srenew(*gn, block->nr);
 889             sprintf(buf, "%s_chain%d", (*gn)[sel_nr], j+1);
 890             (*gn)[block->nr-1] = gmx_strdup(buf);
 891             for (i = block->index[sel_nr]; i < block->index[sel_nr+1]; i++)
 892             {
 893                 a = block->a[i];
 894                 if ((a >= start[j]) && (a <= end[j]))
 895                 {
 896                     block->a[block->nra] = a;
 897                     block->nra++;
 898                 }
 899             }
 900             block->index[block->nr] = block->nra;
 901             if (block->index[block->nr-1] == block->index[block->nr])
 902             {
 903                 remove_group(block->nr-1, NOTSET, block, gn);
 904             }
 905         }
 906     }
 907     sfree(start);
 908     sfree(end);
 909
 910     return nchain;
 911 }
 912
 913 static gmx_bool check_have_atoms(t_atoms *atoms, char *string)
 914 {
 915     if (atoms == NULL)
 916     {
 917         printf("Can not process '%s' without atom info, use option -f\n", string);
 918         return FALSE;
 919     }
 920     else
 921     {
 922         return TRUE;
 923     }
 924 }
 925
 926 static gmx_bool parse_entry(char **string, int natoms, t_atoms *atoms,
 927                             t_blocka *block, char ***gn,
 928                             atom_id *nr, atom_id *index, char *gname)
 929 {
 930     static char   **names, *ostring;
 931     static gmx_bool bFirst = TRUE;
 932     int             j, n_names, sel_nr1;
 933     atom_id         i, nr1, *index1;
 934     char            c;
 935     gmx_bool        bRet, bCompl;
 936
 937     if (bFirst)
 938     {
 939         bFirst = FALSE;
 940         snew(names, MAXNAMES);
 941         for (i = 0; i < MAXNAMES; i++)
 942         {
 943             snew(names[i], NAME_LEN+1);
 944         }
 945     }
 946
 947     bRet    = FALSE;
 948     sel_nr1 = NOTSET;
 949
 950     while (*string[0] == ' ')
 951     {
 952         (*string)++;
 953     }
 954
 955     if ((*string)[0] == '!')
 956     {
 957         bCompl = TRUE;
 958         (*string)++;
 959         while (*string[0] == ' ')
 960         {
 961             (*string)++;
 962         }
 963     }
 964     else
 965     {
 966         bCompl = FALSE;
 967     }
 968
 969     ostring = *string;
 970
 971     if (parse_int(string, &sel_nr1) ||
 972         parse_string(string, &sel_nr1, block->nr, *gn))
 973     {
 974         if ((sel_nr1 >= 0) && (sel_nr1 < block->nr))
 975         {
 976             copy_group(sel_nr1, block, nr, index);
 977             strcpy(gname, (*gn)[sel_nr1]);
 978             printf("Copied index group %d '%s'\n", sel_nr1, (*gn)[sel_nr1]);
 979             bRet = TRUE;
 980         }
 981         else
 982         {
 983             printf("Group %d does not exist\n", sel_nr1);
 984         }
 985     }
 986     else if ((*string)[0] == 'a')
 987     {
 988         (*string)++;
 989         if (check_have_atoms(atoms, ostring))
 990         {
 991             if (parse_int(string, &sel_nr1))
 992             {
 993                 bRet = select_atomnumbers(string, atoms, sel_nr1, nr, index, gname);
 994             }
 995             else if (parse_names(string, &n_names, names))
 996             {
 997                 bRet = select_atomnames(atoms, n_names, names, nr, index, FALSE);
 998                 make_gname(n_names, names, gname);
 999             }
1000         }
1001     }
1002     else if ((*string)[0] == 't')
1003     {
1004         (*string)++;
1005         if (check_have_atoms(atoms, ostring) &&
1006             parse_names(string, &n_names, names))
1007         {
1008             if (atoms->atomtype == NULL)
1009             {
1010                 printf("Need a run input file to select atom types\n");
1011             }
1012             else
1013             {
1014                 bRet = select_atomnames(atoms, n_names, names, nr, index, TRUE);
1015                 make_gname(n_names, names, gname);
1016             }
1017         }
1018     }
1019     else if (strncmp(*string, "res", 3) == 0)
1020     {
1021         (*string) += 3;
1022         if (check_have_atoms(atoms, ostring) &&
1023             parse_int(string, &sel_nr1) &&
1024             (sel_nr1 >= 0) && (sel_nr1 < block->nr) )
1025         {
1026             bRet = atoms_from_residuenumbers(atoms,
1027                                              sel_nr1, block, nr, index, (*gn)[sel_nr1]);
1028             sprintf(gname, "atom_%s", (*gn)[sel_nr1]);
1029         }
1030     }
1031     else if (strncmp(*string, "ri", 2) == 0)
1032     {
1033         (*string) += 2;
1034         if (check_have_atoms(atoms, ostring) &&
1035             parse_int_char(string, &sel_nr1, &c))
1036         {
1037             bRet = select_residueindices(string, atoms, sel_nr1, c, nr, index, gname);
1038         }
1039     }
1040     else if ((*string)[0] == 'r')
1041     {
1042         (*string)++;
1043         if (check_have_atoms(atoms, ostring))
1044         {
1045             if (parse_int_char(string, &sel_nr1, &c))
1046             {
1047                 bRet = select_residuenumbers(string, atoms, sel_nr1, c, nr, index, gname);
1048             }
1049             else if (parse_names(string, &n_names, names))
1050             {
1051                 bRet = select_residuenames(atoms, n_names, names, nr, index);
1052                 make_gname(n_names, names, gname);
1053             }
1054         }
1055     }
1056     else if (strncmp(*string, "chain", 5) == 0)
1057     {
1058         (*string) += 5;
1059         if (check_have_atoms(atoms, ostring) &&
1060             parse_names(string, &n_names, names))
1061         {
1062             bRet = select_chainnames(atoms, n_names, names, nr, index);
1063             sprintf(gname, "ch%s", names[0]);
1064             for (i = 1; i < n_names; i++)
1065             {
1066                 strcat(gname, names[i]);
1067             }
1068         }
1069     }
1070     if (bRet && bCompl)
1071     {
1072         snew(index1, natoms-*nr);
1073         nr1 = 0;
1074         for (i = 0; i < natoms; i++)
1075         {
1076             j = 0;
1077             while ((j < *nr) && (index[j] != i))
1078             {
1079                 j++;
1080             }
1081             if (j == *nr)
1082             {
1083                 if (nr1 >= natoms-*nr)
1084                 {
1085                     printf("There are double atoms in your index group\n");
1086                     break;
1087                 }
1088                 index1[nr1] = i;
1089                 nr1++;
1090             }
1091         }
1092         *nr = nr1;
1093         for (i = 0; i < nr1; i++)
1094         {
1095             index[i] = index1[i];
1096         }
1097         sfree(index1);
1098
1099         for (i = strlen(gname)+1; i > 0; i--)
1100         {
1101             gname[i] = gname[i-1];
1102         }
1103         gname[0] = '!';
1104         printf("Complemented group: %d atoms\n", *nr);
1105     }
1106
1107     return bRet;
1108 }
1109
1110 static void list_residues(t_atoms *atoms)
1111 {
1112     int      i, j, start, end, prev_resind, resind;
1113     gmx_bool bDiff;
1114
1115     /* Print all the residues, assuming continuous resnr count */
1116     start       = atoms->atom[0].resind;
1117     prev_resind = start;
1118     for (i = 0; i < atoms->nr; i++)
1119     {
1120         resind = atoms->atom[i].resind;
1121         if ((resind != prev_resind) || (i == atoms->nr-1))
1122         {
1123             if ((bDiff = strcmp(*atoms->resinfo[resind].name,
1124                                 *atoms->resinfo[start].name)) ||
1125                 (i == atoms->nr-1))
1126             {
1127                 if (bDiff)
1128                 {
1129                     end = prev_resind;
1130                 }
1131                 else
1132                 {
1133                     end = resind;
1134                 }
1135                 if (end < start+3)
1136                 {
1137                     for (j = start; j <= end; j++)
1138                     {
1139                         printf("%4d %-5s",
1140                                j+1, *(atoms->resinfo[j].name));
1141                     }
1142                 }
1143                 else
1144                 {
1145                     printf(" %4d - %4d %-5s  ",
1146                            start+1, end+1, *(atoms->resinfo[start].name));
1147                 }
1148                 start = resind;
1149             }
1150         }
1151         prev_resind = resind;
1152     }
1153     printf("\n");
1154 }
1155
1156 static void edit_index(int natoms, t_atoms *atoms, rvec *x, t_blocka *block, char ***gn, gmx_bool bVerbose)
1157 {
1158     static char   **atnames, *ostring;
1159     static gmx_bool bFirst = TRUE;
1160     char            inp_string[STRLEN], *string;
1161     char            gname[STRLEN], gname1[STRLEN], gname2[STRLEN];
1162     int             i, i0, i1, sel_nr, sel_nr2, newgroup;
1163     atom_id         nr, nr1, nr2, *index, *index1, *index2;
1164     gmx_bool        bAnd, bOr, bPrintOnce;
1165
1166     if (bFirst)
1167     {
1168         bFirst = FALSE;
1169         snew(atnames, MAXNAMES);
1170         for (i = 0; i < MAXNAMES; i++)
1171         {
1172             snew(atnames[i], NAME_LEN+1);
1173         }
1174     }
1175
1176     string = NULL;
1177
1178     snew(index, natoms);
1179     snew(index1, natoms);
1180     snew(index2, natoms);
1181
1182     newgroup   = NOTSET;
1183     bPrintOnce = TRUE;
1184     do
1185     {
1186         gname1[0] = '\0';
1187         if (bVerbose || bPrintOnce || newgroup != NOTSET)
1188         {
1189             printf("\n");
1190             if (bVerbose || bPrintOnce || newgroup == NOTSET)
1191             {
1192                 i0 = 0;
1193                 i1 = block->nr;
1194             }
1195             else
1196             {
1197                 i0 = newgroup;
1198                 i1 = newgroup+1;
1199             }
1200             for (i = i0; i < i1; i++)
1201             {
1202                 printf("%3d %-20s: %5d atoms\n", i, (*gn)[i],
1203                        block->index[i+1]-block->index[i]);
1204             }
1205             newgroup = NOTSET;
1206         }
1207         if (bVerbose || bPrintOnce)
1208         {
1209             printf("\n");
1210             printf(" nr : group       !   'name' nr name   'splitch' nr    Enter: list groups\n");
1211             printf(" 'a': atom        &   'del' nr         'splitres' nr   'l': list residues\n");
1212             printf(" 't': atom type   |   'keep' nr        'splitat' nr    'h': help\n");
1213             printf(" 'r': residue         'res' nr         'chain' char\n");
1214             printf(" \"name\": group        'case': case %s         'q': save and quit\n",
1215                    bCase ? "insensitive" : "sensitive  ");
1216             printf(" 'ri': residue index\n");
1217             bPrintOnce = FALSE;
1218         }
1219         printf("\n");
1220         printf("> ");
1221         if (NULL == fgets(inp_string, STRLEN, stdin))
1222         {
1223             gmx_fatal(FARGS, "Error reading user input");
1224         }
1225         inp_string[strlen(inp_string)-1] = 0;
1226         printf("\n");
1227         string = inp_string;
1228         while (string[0] == ' ')
1229         {
1230             string++;
1231         }
1232
1233         ostring = string;
1234         nr      = 0;
1235         if (string[0] == 'h')
1236         {
1237             printf(" nr                : selects an index group by number or quoted string.\n");
1238             printf("                     The string is first matched against the whole group name,\n");
1239             printf("                     then against the beginning and finally against an\n");
1240             printf("                     arbitrary substring. A multiple match is an error.\n");
1241
1242             printf(" 'a' nr1 [nr2 ...] : selects atoms, atom numbering starts at 1.\n");
1243             printf(" 'a' nr1 - nr2     : selects atoms in the range from nr1 to nr2.\n");
1244             printf(" 'a' name1[*] [name2[*] ...] : selects atoms by name(s), '?' matches any char,\n");
1245             printf("                               wildcard '*' allowed at the end of a name.\n");
1246             printf(" 't' type1[*] [type2[*] ...] : as 'a', but for type, run input file required.\n");
1247             printf(" 'r' nr1[ic1] [nr2[ic2] ...] : selects residues by number and insertion code.\n");
1248             printf(" 'r' nr1 - nr2               : selects residues in the range from nr1 to nr2.\n");
1249             printf(" 'r' name1[*] [name2[*] ...] : as 'a', but for residue names.\n");
1250             printf(" 'ri' nr1 - nr2              : selects residue indices, 1-indexed, (as opposed to numbers) in the range from nr1 to nr2.\n");
1251             printf(" 'chain' ch1 [ch2 ...]       : selects atoms by chain identifier(s),\n");
1252             printf("                               not available with a .gro file as input.\n");
1253             printf(" !                 : takes the complement of a group with respect to all\n");
1254             printf("                     the atoms in the input file.\n");
1255             printf(" & |               : AND and OR, can be placed between any of the options\n");
1256             printf("                     above, the input is processed from left to right.\n");
1257             printf(" 'name' nr name    : rename group nr to name.\n");
1258             printf(" 'del' nr1 [- nr2] : deletes one group or groups in the range from nr1 to nr2.\n");
1259             printf(" 'keep' nr         : deletes all groups except nr.\n");
1260             printf(" 'case'            : make all name compares case (in)sensitive.\n");
1261             printf(" 'splitch' nr      : split group into chains using CA distances.\n");
1262             printf(" 'splitres' nr     : split group into residues.\n");
1263             printf(" 'splitat' nr      : split group into atoms.\n");
1264             printf(" 'res' nr          : interpret numbers in group as residue numbers\n");
1265             printf(" Enter             : list the currently defined groups and commands\n");
1266             printf(" 'l'               : list the residues.\n");
1267             printf(" 'h'               : show this help.\n");
1268             printf(" 'q'               : save and quit.\n");
1269             printf("\n");
1270             printf(" Examples:\n");
1271             printf(" > 2 | 4 & r 3-5\n");
1272             printf(" selects all atoms from group 2 and 4 that have residue numbers 3, 4 or 5\n");
1273             printf(" > a C* & !a C CA\n");
1274             printf(" selects all atoms starting with 'C' but not the atoms 'C' and 'CA'\n");
1275             printf(" > \"protein\" & ! \"backb\"\n");
1276             printf(" selects all atoms that are in group 'protein' and not in group 'backbone'\n");
1277             if (bVerbose)
1278             {
1279                 printf("\npress Enter ");
1280                 getchar();
1281             }
1282         }
1283         else if (strncmp(string, "del", 3) == 0)
1284         {
1285             string += 3;
1286             if (parse_int(&string, &sel_nr))
1287             {
1288                 while (string[0] == ' ')
1289                 {
1290                     string++;
1291                 }
1292                 if (string[0] == '-')
1293                 {
1294                     string++;
1295                     parse_int(&string, &sel_nr2);
1296                 }
1297                 else
1298                 {
1299                     sel_nr2 = NOTSET;
1300                 }
1301                 while (string[0] == ' ')
1302                 {
1303                     string++;
1304                 }
1305                 if (string[0] == '\0')
1306                 {
1307                     remove_group(sel_nr, sel_nr2, block, gn);
1308                 }
1309                 else
1310                 {
1311                     printf("\nSyntax error: \"%s\"\n", string);
1312                 }
1313             }
1314         }
1315         else if (strncmp(string, "keep", 4) == 0)
1316         {
1317             string += 4;
1318             if (parse_int(&string, &sel_nr))
1319             {
1320                 remove_group(sel_nr+1, block->nr-1, block, gn);
1321                 remove_group(0, sel_nr-1, block, gn);
1322             }
1323         }
1324         else if (strncmp(string, "name", 4) == 0)
1325         {
1326             string += 4;
1327             if (parse_int(&string, &sel_nr))
1328             {
1329                 if ((sel_nr >= 0) && (sel_nr < block->nr))
1330                 {
1331                     sscanf(string, "%s", gname);
1332                     sfree((*gn)[sel_nr]);
1333                     (*gn)[sel_nr] = gmx_strdup(gname);
1334                 }
1335             }
1336         }
1337         else if (strncmp(string, "case", 4) == 0)
1338         {
1339             bCase = !bCase;
1340             printf("Switched to case %s\n", bCase ? "sensitive" : "insensitive");
1341         }
1342         else if (string[0] == 'v')
1343         {
1344             bVerbose = !bVerbose;
1345             printf("Turned verbose %s\n", bVerbose ? "on" : "off");
1346         }
1347         else if (string[0] == 'l')
1348         {
1349             if (check_have_atoms(atoms, ostring) )
1350             {
1351                 list_residues(atoms);
1352             }
1353         }
1354         else if (strncmp(string, "splitch", 7) == 0)
1355         {
1356             string += 7;
1357             if (check_have_atoms(atoms, ostring) &&
1358                 parse_int(&string, &sel_nr) &&
1359                 (sel_nr >= 0) && (sel_nr < block->nr))
1360             {
1361                 split_chain(atoms, x, sel_nr, block, gn);
1362             }
1363         }
1364         else if (strncmp(string, "splitres", 8) == 0)
1365         {
1366             string += 8;
1367             if (check_have_atoms(atoms, ostring) &&
1368                 parse_int(&string, &sel_nr) &&
1369                 (sel_nr >= 0) && (sel_nr < block->nr))
1370             {
1371                 split_group(atoms, sel_nr, block, gn, FALSE);
1372             }
1373         }
1374         else if (strncmp(string, "splitat", 7) == 0)
1375         {
1376             string += 7;
1377             if (check_have_atoms(atoms, ostring) &&
1378                 parse_int(&string, &sel_nr) &&
1379                 (sel_nr >= 0) && (sel_nr < block->nr))
1380             {
1381                 split_group(atoms, sel_nr, block, gn, TRUE);
1382             }
1383         }
1384         else if (string[0] == '\0')
1385         {
1386             bPrintOnce = TRUE;
1387         }
1388         else if (string[0] != 'q')
1389         {
1390             nr1 = -1;
1391             nr2 = -1;
1392             if (parse_entry(&string, natoms, atoms, block, gn, &nr, index, gname))
1393             {
1394                 do
1395                 {
1396                     while (string[0] == ' ')
1397                     {
1398                         string++;
1399                     }
1400
1401                     bAnd = FALSE;
1402                     bOr  = FALSE;
1403                     if (string[0] == '&')
1404                     {
1405                         bAnd = TRUE;
1406                     }
1407                     else if (string[0] == '|')
1408                     {
1409                         bOr = TRUE;
1410                     }
1411
1412                     if (bAnd || bOr)
1413                     {
1414                         string++;
1415                         nr1 = nr;
1416                         for (i = 0; i < nr; i++)
1417                         {
1418                             index1[i] = index[i];
1419                         }
1420                         strcpy(gname1, gname);
1421                         if (parse_entry(&string, natoms, atoms, block, gn, &nr2, index2, gname2))
1422                         {
1423                             if (bOr)
1424                             {
1425                                 or_groups(nr1, index1, nr2, index2, &nr, index);
1426                                 sprintf(gname, "%s_%s", gname1, gname2);
1427                             }
1428                             else
1429                             {
1430                                 and_groups(nr1, index1, nr2, index2, &nr, index);
1431                                 sprintf(gname, "%s_&_%s", gname1, gname2);
1432                             }
1433                         }
1434                     }
1435                 }
1436                 while (bAnd || bOr);
1437             }
1438             while (string[0] == ' ')
1439             {
1440                 string++;
1441             }
1442             if (string[0])
1443             {
1444                 printf("\nSyntax error: \"%s\"\n", string);
1445             }
1446             else if (nr > 0)
1447             {
1448                 copy2block(nr, index, block);
1449                 srenew(*gn, block->nr);
1450                 newgroup        = block->nr-1;
1451                 (*gn)[newgroup] = gmx_strdup(gname);
1452             }
1453             else
1454             {
1455                 printf("Group is empty\n");
1456             }
1457         }
1458     }
1459     while (string[0] != 'q');
1460
1461     sfree(index);
1462     sfree(index1);
1463     sfree(index2);
1464 }
1465
1466 static int block2natoms(t_blocka *block)
1467 {
1468     int i, natoms;
1469
1470     natoms = 0;
1471     for (i = 0; i < block->nra; i++)
1472     {
1473         natoms = max(natoms, block->a[i]);
1474     }
1475     natoms++;
1476
1477     return natoms;
1478 }
1479
1480 void merge_blocks(t_blocka *dest, t_blocka *source)
1481 {
1482     int     i, nra0, i0;
1483
1484     /* count groups, srenew and fill */
1485     i0        = dest->nr;
1486     nra0      = dest->nra;
1487     dest->nr += source->nr;
1488     srenew(dest->index, dest->nr+1);
1489     for (i = 0; i < source->nr; i++)
1490     {
1491         dest->index[i0+i] = nra0 + source->index[i];
1492     }
1493     /* count atoms, srenew and fill */
1494     dest->nra += source->nra;
1495     srenew(dest->a, dest->nra);
1496     for (i = 0; i < source->nra; i++)
1497     {
1498         dest->a[nra0+i] = source->a[i];
1499     }
1500
1501     /* terminate list */
1502     dest->index[dest->nr] = dest->nra;
1503
1504 }
1505
1506 int gmx_make_ndx(int argc, char *argv[])
1507 {
1508     const char     *desc[] = {
1509         "Index groups are necessary for almost every GROMACS program.",
1510         "All these programs can generate default index groups. You ONLY",
1511         "have to use [THISMODULE] when you need SPECIAL index groups.",
1512         "There is a default index group for the whole system, 9 default",
1513         "index groups for proteins, and a default index group",
1514         "is generated for every other residue name.[PAR]",
1515         "When no index file is supplied, also [THISMODULE] will generate the",
1516         "default groups.",
1517         "With the index editor you can select on atom, residue and chain names",
1518         "and numbers.",
1519         "When a run input file is supplied you can also select on atom type.",
1520         "You can use NOT, AND and OR, you can split groups",
1521         "into chains, residues or atoms. You can delete and rename groups.[PAR]",
1522         "The atom numbering in the editor and the index file starts at 1.[PAR]",
1523         "The [TT]-twin[tt] switch duplicates all index groups with an offset of",
1524         "[TT]-natoms[tt], which is useful for Computational Electrophysiology",
1525         "double-layer membrane setups."
1526     };
1527
1528     static int      natoms     = 0;
1529     static gmx_bool bVerbose   = FALSE;
1530     static gmx_bool bDuplicate = FALSE;
1531     t_pargs         pa[]       = {
1532         { "-natoms",  FALSE, etINT, {&natoms},
1533           "set number of atoms (default: read from coordinate or index file)" },
1534         { "-twin",     FALSE, etBOOL, {&bDuplicate},
1535           "Duplicate all index groups with an offset of -natoms" },
1536         { "-verbose", FALSE, etBOOL, {&bVerbose},
1537           "HIDDENVerbose output" }
1538     };
1539 #define NPA asize(pa)
1540
1541     output_env_t oenv;
1542     char         title[STRLEN];
1543     int          nndxin;
1544     const char  *stxfile;
1545     char       **ndxinfiles;
1546     const char  *ndxoutfile;
1547     gmx_bool     bNatoms;
1548     int          i, j;
1549     t_atoms     *atoms;
1550     rvec        *x, *v;
1551     int          ePBC;
1552     matrix       box;
1553     t_blocka    *block, *block2;
1554     char       **gnames, **gnames2;
1555     t_filenm     fnm[] = {
1556         { efSTX, "-f", NULL,     ffOPTRD  },
1557         { efNDX, "-n", NULL,     ffOPTRDMULT },
1558         { efNDX, "-o", NULL,     ffWRITE }
1559     };
1560 #define NFILE asize(fnm)
1561
1562     if (!parse_common_args(&argc, argv, 0, NFILE, fnm, NPA, pa, asize(desc), desc,
1563                            0, NULL, &oenv))
1564     {
1565         return 0;
1566     }
1567
1568     stxfile = ftp2fn_null(efSTX, NFILE, fnm);
1569     if (opt2bSet("-n", NFILE, fnm))
1570     {
1571         nndxin = opt2fns(&ndxinfiles, "-n", NFILE, fnm);
1572     }
1573     else
1574     {
1575         nndxin = 0;
1576     }
1577     ndxoutfile = opt2fn("-o", NFILE, fnm);
1578     bNatoms    = opt2parg_bSet("-natoms", NPA, pa);
1579
1580     if (!stxfile && !nndxin)
1581     {
1582         gmx_fatal(FARGS, "No input files (structure or index)");
1583     }
1584
1585     if (stxfile)
1586     {
1587         snew(atoms, 1);
1588         get_stx_coordnum(stxfile, &(atoms->nr));
1589         init_t_atoms(atoms, atoms->nr, TRUE);
1590         snew(x, atoms->nr);
1591         snew(v, atoms->nr);
1592         fprintf(stderr, "\nReading structure file\n");
1593         read_stx_conf(stxfile, title, atoms, x, v, &ePBC, box);
1594         natoms  = atoms->nr;
1595         bNatoms = TRUE;
1596     }
1597     else
1598     {
1599         atoms = NULL;
1600         x     = NULL;
1601     }
1602
1603     /* read input file(s) */
1604     block  = new_blocka();
1605     gnames = NULL;
1606     printf("Going to read %d old index file(s)\n", nndxin);
1607     if (nndxin)
1608     {
1609         for (i = 0; i < nndxin; i++)
1610         {
1611             block2 = init_index(ndxinfiles[i], &gnames2);
1612             srenew(gnames, block->nr+block2->nr);
1613             for (j = 0; j < block2->nr; j++)
1614             {
1615                 gnames[block->nr+j] = gnames2[j];
1616             }
1617             sfree(gnames2);
1618             merge_blocks(block, block2);
1619             sfree(block2->a);
1620             sfree(block2->index);
1621 /*       done_block(block2); */
1622             sfree(block2);
1623         }
1624     }
1625     else
1626     {
1627         snew(gnames, 1);
1628         analyse(atoms, block, &gnames, FALSE, TRUE);
1629     }
1630
1631     if (!bNatoms)
1632     {
1633         natoms = block2natoms(block);
1634         printf("Counted atom numbers up to %d in index file\n", natoms);
1635     }
1636
1637     edit_index(natoms, atoms, x, block, &gnames, bVerbose);
1638
1639     write_index(ndxoutfile, block, gnames, bDuplicate, natoms);
1640
1641     return 0;
1642 }